molfreq.cc

大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CHARMM并行计算。非常具有应用价值。

  if (linear) {
      //if (D_inf_h) sigma = 2;
      if (ct.symbol()[0] == 'D' ||
          ct.symbol()[0] == 'd') sigma = 2;
      else if (ct.symbol()[0] == 'C' ||
               ct.symbol()[0] == 'c') sigma = 1;
      else {
          ExEnv::errn() << "MolecularFrequencies: For linear molecules"
               << " the specified point group must be Cnv or Dnh"
               << endl;
          abort();
          }
      }
  else if ((ct.symbol()[0] == 'C' ||
            ct.symbol()[0] == 'c') &&
           (ct.symbol()[1] >= '1'  &&
            ct.symbol()[1] <= '8') &&
            ct.symbol()[2] == '\0') {
      sigma = ct.order();  // group is a valid CN
      }
  else if ((ct.symbol()[0] == 'D' ||
            ct.symbol()[0] == 'd') &&
           (ct.symbol()[1] >= '2'  &&
            ct.symbol()[1] <= '6') &&
            ct.symbol()[2] == '\0') {
      sigma = ct.order();  // group is a valid DN
      }
  else if ((ct.symbol()[0] == 'T' ||
            ct.symbol()[0] == 't') &&
            ct.symbol()[1] == '\0') {
      sigma = ct.order();  // group is T
      }
  else sigma = (int)(0.5*ct.order()); // group is not pure rot. group (CN, DN, or T)

  // compute S_trans
  double S_trans;
  tmpvar = pow(2*pi*mass*k*T/(h*h),1.5);
  S_trans = R*(log(tmpvar*R*T/(P*atm_to_Pa)) + 2.5 - log(NA));

  // compute S_rot
  double S_rot;
  double theta[3]; // rotational temperatures (K)
  if (linear) {
      theta[1] = h*h/(8*pi*pi*pmi[1]*amu_to_kg*pow(angstrom_to_meter,2.0)*k);
      S_rot = log(T/(sigma*theta[1])) + 1.0;
      }
  else {
      theta[0] = h*h/(8*pi*pi*pmi[0]*amu_to_kg*pow(angstrom_to_meter,2.0)*k);
      theta[1] = h*h/(8*pi*pi*pmi[1]*amu_to_kg*pow(angstrom_to_meter,2.0)*k);
      theta[2] = h*h/(8*pi*pi*pmi[2]*amu_to_kg*pow(angstrom_to_meter,2.0)*k);
      tmpvar = theta[0]*theta[1]*theta[2];
      S_rot = log(pow(pi*T*T*T/tmpvar,0.5)/sigma) + 1.5;
      }
  S_rot *= R;

  // compute S_vib
  double S_vib = 0.0;
  for (i=0; i<nirrep_; i++) {
      for (int j=0; j<nfreq_[i]; j++) {
          if (freq_[i][j] > 0.0) {
              tmpvar = hartree_to_hertz*h*freq_[i][j]/(k*T);
              double expval = exp(-tmpvar);
              S_vib += tmpvar*expval/(1.0-expval) - log(1.0-expval);
              }
          }
      }
   S_vib *= R;

  // compute S_el
  double S_el;
  S_el = R*log(double(degeneracy));

  // compute total molar entropy S (in J/(mol*K))
  double S;
  S = S_trans + S_rot + S_vib + S_el;

  
  //////////////////////////////////////////////
  // compute the molar enthalpy (nonelectronic) 
  //////////////////////////////////////////////

  int n_zero_or_imaginary = 0;
  double E0vib = 0.0;
  for (i=0; i<nirrep_; i++) {
      for (int j=0; j<nfreq_[i]; j++) {
          if (freq_[i][j] > 0.0) E0vib += freq_[i][j] * hartree_to_joule_per_mol;
          else n_zero_or_imaginary++;
        }
    }
  E0vib *= 0.5;

  double EvibT = 0.0;
  for (i=0; i<nirrep_; i++) {
      for (int j=0; j<nfreq_[i]; j++) {
          if (freq_[i][j] > 0.0) {
              double expval = exp(-freq_[i][j]*hartree_to_joule/(k*T));
              EvibT += freq_[i][j] * hartree_to_joule_per_mol
                     * expval/(1.0-expval);
            }
        }
    }

  double EPV = NA*k*T;

  int nexternal = 6;
  if (mol_->natom() == 1) nexternal = 3;
  else if (mol_->is_linear()) nexternal = 5;

  double Erot;
  if (nexternal == 3) {
      // atom
      Erot = 0.0;
    }
  else if (nexternal == 5) {
      // linear
      Erot = EPV;
    }
  else if (nexternal == 6) {
      // nonlinear
      Erot = 1.5 * EPV;
    }
  else {
      ExEnv::errn() << "Strange number of external coordinates: " << nexternal
           << ".  Setting Erot to 0.0" << endl;
      Erot = 0.0;
    }

  double Etrans = 1.5 * EPV;

  ////////////////////////////////////////////////
  // Print out results of thermodynamic analysis
  ////////////////////////////////////////////////

  ExEnv::out0() << indent << "THERMODYNAMIC ANALYSIS:" << endl << endl
       << indent << scprintf("Contributions to the nonelectronic enthalpy at %.2lf K:\n",T)
       << indent << "                   kJ/mol       kcal/mol"<< endl
       << indent << scprintf("  E0vib        = %9.4lf    %9.4lf\n",
          E0vib/1000, E0vib/(4.184*1000))
       << indent << scprintf("  Evib(T)      = %9.4lf    %9.4lf\n",
          EvibT/1000, EvibT/(4.184*1000))
       << indent << scprintf("  Erot(T)      = %9.4lf    %9.4lf\n",
          Erot/1000, Erot/(4.184*1000))
       << indent << scprintf("  Etrans(T)    = %9.4lf    %9.4lf\n",
          Etrans/1000, Etrans/(4.184*1000))
       << indent << scprintf("  PV(T)        = %9.4lf    %9.4lf\n",
          EPV/1000, EPV/(4.184*1000))
       << indent << scprintf("  Total nonelectronic enthalpy:\n")
       << indent << scprintf("  H_nonel(T)   = %9.4lf    %9.4lf\n",
         (E0vib+EvibT+Erot+Etrans+EPV)/1000,
         (E0vib+EvibT+Erot+Etrans+EPV)/(4.184*1000))
       << endl
       << indent
       << scprintf("Contributions to the entropy at %.2lf K and %.1lf atm:\n",
                   T, P)
       << indent << "                   J/(mol*K)    cal/(mol*K)"<< endl
       << indent
       << scprintf("  S_trans(T,P) = %9.4lf    %9.4lf\n",
                   S_trans, S_trans/4.184)
       << indent
       << scprintf("  S_rot(T)     = %9.4lf    %9.4lf\n", S_rot,S_rot/4.184)
       << indent
       << scprintf("  S_vib(T)     = %9.4lf    %9.4lf\n", S_vib,S_vib/4.184)
       << indent
       << scprintf("  S_el         = %9.4lf    %9.4lf\n", S_el,S_el/4.184)
       << indent << scprintf("  Total entropy:\n")
       << indent << scprintf("  S_total(T,P) = %9.4lf    %9.4lf\n", S, S/4.184)
       << indent << endl

       << indent << "Various data used for thermodynamic analysis:" << endl
       << indent << endl;

  if (linear) ExEnv::out0() << indent << "Linear molecule" << endl;
  else ExEnv::out0() << indent << "Nonlinear molecule" << endl;

  ExEnv::out0() << indent
       << scprintf("Principal moments of inertia (amu*angstrom^2):"
          " %.5lf, %.5lf, %.5lf\n", pmi[0], pmi[1], pmi[2])
       << indent << "Point group: " << ct.symbol()
       << endl
       << indent << "Order of point group: " << ct.order() << endl
       << indent << "Rotational symmetry number: " << sigma << endl;

  if (linear) {
      ExEnv::out0() << indent
           << scprintf("Rotational temperature (K): %.4lf\n", theta[1]);
    }
  else {
      ExEnv::out0() << indent
           << scprintf("Rotational temperatures (K): %.4lf, %.4lf, %.4lf\n",
                       theta[0], theta[1], theta[2]);
    }

  ExEnv::out0() << indent << "Electronic degeneracy: " << degeneracy
       << endl << endl;
}

void
MolecularFrequencies::animate(const Ref<Render>& render,
                              const Ref<MolFreqAnimate>& anim)
{
  int i,j, symoff = 0;
  for (i=0; i<nirrep_; i++) {
      int nfreq = disym_->blocks()->size(i);
      for (j=0; j<nfreq; j++) {
          char name[128];
          sprintf(name,"%02d.%s",
                  nfreq-j+symoff, pg_->char_table().gamma(i).symbol_ns());
          anim->set_name(name);
          anim->set_mode(i,j);
          render->animate(anim.pointer());
        }
      symoff += nfreq;
    }
}

/////////////////////////////////////////////////////////////////////////////
// MolFreqAnimate

static ClassDesc MolFreqAnimate_cd(
  typeid(MolFreqAnimate),"MolFreqAnimate",1,"public AnimatedObject",
  0, create<MolFreqAnimate>, 0);

MolFreqAnimate::MolFreqAnimate(const Ref<KeyVal> &keyval):
  AnimatedObject(keyval)
{
  renmol_ << keyval->describedclassvalue("rendered");
  molfreq_ << keyval->describedclassvalue("freq");
  dependent_mole_ << keyval->describedclassvalue("dependent_mole");
  irrep_ = keyval->intvalue("irrep");
  mode_ = keyval->intvalue("mode");
  KeyValValueint default_nframe(10);
  nframe_ = keyval->intvalue("nframe",default_nframe);
  KeyValValuedouble default_disp(0.2);
  disp_ = keyval->doublevalue("displacement", default_disp);
}

MolFreqAnimate::~MolFreqAnimate()
{
}

int
MolFreqAnimate::nobject()
{
  return nframe_;
}

Ref<RenderedObject>
MolFreqAnimate::object(int iobject)
{
  BlockedSCMatrix *normco
      = dynamic_cast<BlockedSCMatrix*>(molfreq_->normal_coordinates().pointer());
  Ref<Molecule> mol = renmol_->molecule();
  Ref<Molecule> molcopy = new Molecule(*mol.pointer());

  double scale = disp_ * cos(M_PI*(iobject+0.5)/(double)nframe_);

  RefSCMatrix irrepblock = normco->block(irrep_);
  int ixyz, iatom, icoor=0;
  for (iatom=0; iatom<mol->natom(); iatom++) {
      for (ixyz=0; ixyz<3; ixyz++, icoor++) {
          mol->r(iatom,ixyz) += scale
                                   * irrepblock->get_element(icoor,mode_);
        }
    }

  if (dependent_mole_.nonnull()) dependent_mole_->obsolete();
  renmol_->init();

  char name[64];
  sprintf(name,"%02d",iobject);
  renmol_->set_name(name);

  // restore the original molecule
  mol->operator = (*molcopy.pointer());
  if (dependent_mole_.nonnull()) dependent_mole_->obsolete();

  return renmol_.pointer();
}


/////////////////////////////////////////////////////////////////////////////

// Local Variables:
// mode: c++
// c-file-style: "CLJ"
// End: