molfreq.cc

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

//
// molfreq.cc
//
// Copyright (C) 1996 Limit Point Systems, Inc.
//
// Author: Curtis Janssen <[e-mail]>
// Maintainer: LPS
//
// This file is part of the SC Toolkit.
//
// The SC Toolkit is free software; you can redistribute it and/or modify
// it under the terms of the GNU Library General Public License as published by
// the Free Software Foundation; either version 2, or (at your option)
// any later version.
//
// The SC Toolkit is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
//
// The U.S. Government is granted a limited license as per AL 91-7.
//

#ifdef __GNUC__
#pragma implementation
#endif

#include <math.h>
#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <util/group/message.h>
#include <math/symmetry/corrtab.h>
#include <math/scmat/local.h>
#include <math/scmat/blocked.h>
#include <chemistry/molecule/molfreq.h>
#include <chemistry/molecule/molrender.h>

using namespace std;
using namespace sc;

#undef DEBUG

static ClassDesc MolecularFrequencies_cd(
  typeid(MolecularFrequencies),"MolecularFrequencies",3,"public SavableState",
  0, create<MolecularFrequencies>, create<MolecularFrequencies>);

MolecularFrequencies::MolecularFrequencies(const Ref<KeyVal>& keyval)
{
  mol_ << keyval->describedclassvalue("molecule");
  if (mol_.null()) {
      ExEnv::err0() << "MolecularFrequencies: KeyVal CTOR requires molecule"
           << endl;
      abort();
    }
  KeyValValueRefDescribedClass def_pg(mol_->point_group().pointer());
  pg_ << keyval->describedclassvalue("point_group", def_pg);
  nirrep_ = pg_->char_table().nirrep();
  debug_ = keyval->booleanvalue("debug");
  nfreq_ = 0;
  freq_ = 0;
}

MolecularFrequencies::~MolecularFrequencies()
{
  delete[] nfreq_;
  if (freq_) {
      for (int i=0; i<nirrep_; i++) {
          delete[] freq_[i];
        }
      delete[] freq_;
    }
}

MolecularFrequencies::MolecularFrequencies(StateIn& si):
  SavableState(si)
{
  int i;

  if (si.version(::class_desc<MolecularFrequencies>()) < 3) {
      ExEnv::errn() << "MolecularFrequencies: cannot restore from old version" << endl;
      abort();
    }

  mol_ << SavableState::restore_state(si);
  pg_ << SavableState::restore_state(si);

  si.get(nirrep_);
  si.get(nfreq_);
  for (i=0; i<nirrep_; i++) si.get(freq_[i]);
}

void
MolecularFrequencies::save_data_state(StateOut& so)
{
  int i;

  SavableState::save_state(mol_.pointer(),so);
  SavableState::save_state(pg_.pointer(),so);

  so.put(nirrep_);
  so.put(nfreq_,nirrep_);
  for (i=0; i<nirrep_; i++) so.put(freq_[i],nfreq_[i]);
}

void
MolecularFrequencies::compute_frequencies(const RefSymmSCMatrix &xhessian)
{
  int i, coor;

  RefSCMatrix symmbasis
      = MolecularHessian::cartesian_to_symmetry(mol_,pg_);
  BlockedSCMatrix *bsymmbasis = dynamic_cast<BlockedSCMatrix*>(symmbasis.pointer());

  kit_ = xhessian->kit();
  d3natom_ = xhessian->dim();
  symkit_ = symmbasis->kit();
  bd3natom_ = symmbasis->coldim();
  disym_ = symmbasis->rowdim();

  ExEnv::out0() << endl
               << indent << "Frequencies (cm-1; negative is imaginary):"
               << endl;

  // initialize the frequency tables
  if (nfreq_) delete[] nfreq_;
  nfreq_ = new int[nirrep_];
  if (freq_) delete[] freq_;
  freq_ = new double*[nirrep_];

  // initialize normal coordinate matrix
  normco_ = symmatrixkit()->matrix(bd3natom_, disym_);

  // find the inverse sqrt mass matrix
  RefDiagSCMatrix m(d3natom_, matrixkit());
  for (i=0,coor=0; i<mol_->natom(); i++) {
      for (int j=0; j<3; j++, coor++) {
          m(coor) = 1.0/sqrt(mol_->mass(i)*(1.0/5.48579903e-4));
        }
    }

  RefSymmSCMatrix dhessian;

  for (int irrep=0; irrep<nirrep_; irrep++) {
      RefSCMatrix dtranst = bsymmbasis->block(irrep);
      RefSCDimension ddim = dtranst.rowdim();
      nfreq_[irrep] = ddim.n();
      freq_[irrep] = new double[nfreq_[irrep]];
      if (ddim.n() == 0) continue;
      dhessian = matrixkit()->symmmatrix(ddim);
      dhessian.assign(0.0);
      dhessian.accumulate_transform(dtranst,xhessian);
      do_freq_for_irrep(irrep, m, dhessian, dtranst);
    }
}

void
MolecularFrequencies::do_freq_for_irrep(
    int irrep,
    const RefDiagSCMatrix &m,
    const RefSymmSCMatrix &dhessian,
    const RefSCMatrix &dtranst)
{
  int i;
  RefSCMatrix dtrans = dtranst.t();
  RefSCDimension ddim = dtrans.coldim();
  if (ddim.n() == 0) return;
  if (debug_) {
      dhessian.print("dhessian");
      dtrans.print("dtrans");
    }
  // find the basis for the normal coordinates
  RefSCMatrix ncbasis = m * dtrans;
  // use the SVD to orthogonalize and check this basis
  RefSCMatrix basU(d3natom_, d3natom_, matrixkit());
  RefSCMatrix basV(ddim, ddim, matrixkit());
  RefDiagSCMatrix bassigma(ddim, matrixkit());
  ncbasis.svd(basU, bassigma, basV);
  for (i=0; i<ddim.n(); i++) {
      if (bassigma(i) < 1.e0-3) {
          ExEnv::err0() << indent
               << "MolecularFrequencies: displacements don't span"
               << " normal coordinates"
               << endl;
          abort();
        }
    }
  ncbasis.assign_subblock(basU, 0, d3natom_.n()-1, 0, ddim.n()-1, 0, 0);
  // a transform from disp to x to q (mass weighted x) to disp
  RefSCMatrix dxqd = ncbasis.t() * m * dtrans;
  // transform the dhessian to the mass weighted dhessian
  RefSymmSCMatrix mdhessian = matrixkit()->symmmatrix(dxqd.rowdim());
  mdhessian.assign(0.0);
  mdhessian.accumulate_transform(dxqd, dhessian);
  if (debug_) {
      mdhessian.print("mass weighted dhessian");
    }
  // diagonalize the hessian
  RefDiagSCMatrix freqs(ddim,matrixkit());
  RefSCMatrix eigvecs(ddim,ddim,matrixkit());
  mdhessian.diagonalize(freqs,eigvecs);
  // convert the eigvals to frequencies in wavenumbers
  for (i=0; i<freqs.n(); i++) {
      if (freqs(i) >=0.0) freqs(i) = sqrt(freqs(i));
      else freqs(i) = -sqrt(-freqs(i));
      freq_[irrep][i] = freqs(i);
      freqs(i) = freqs->get_element(i) * 219474.63;
    }

  ExEnv::out0() << indent
               << pg_->char_table().gamma(irrep).symbol() << endl;
  int ifreqoff = 1;
  for (i=0; i<irrep; i++) ifreqoff += nfreq_[i];
  for (i=0; i<freqs.n(); i++) {
      double freq = freqs(freqs.n()-i-1);
      ExEnv::out0() << indent
                   << scprintf("%4d % 8.2f",i+ifreqoff,freq)
                   << endl;
    }
  ExEnv::out0() << endl;

  if (debug_) {
      eigvecs.print("eigenvectors");
      ncbasis.print("ncbasis");
      (ncbasis*eigvecs).print("ncbasis*eigvecs");
    }
  dynamic_cast<BlockedSCMatrix*>(
      normco_.pointer())->block(irrep).assign(ncbasis*eigvecs);
}

void
MolecularFrequencies::thermochemistry(int degeneracy, double T, double P)
{
  int i;
  double tmpvar;

  if (!nfreq_) return;

  // default values for temperature T and pressure P are
  // 298.15 K and 1 atm (=101325.0 Pa), respectively 

  // 1986 CODATA
  const double NA = 6.0221367e23;  // Avogadro's number
  const double k  = 1.380658e-23;  // Boltzmann's constant (J/K)
  const double h  = 6.6260755e-34; // Planck's constant (J*s)
  const double R  = 8.314510;      // gas constant (J/(mol*K))  (R=k*NA)
  const double pi = 3.14159265358979323846;
  const double hartree_to_hertz = 6.5796838e15; // (hertz/hartree)

  const double hartree_to_joule = 4.3597482e-18; // (J/hartree)
  const double hartree_to_joule_per_mol = hartree_to_joule*NA;
                                            // (J/(mol*hartree))
  const double amu_to_kg = 1.6605402e-27; // (kg/amu)
  const double angstrom_to_meter = 1.0e-10;
  const double atm_to_Pa = 101325.0; // (Pa/atm)


  ////////////////////////////////////////////////////////////////////////
  // compute the molar entropy using formulas for ideal polyatomic gasses
  // from McQuarrie, Statistical Mechanics, 1976, Ch. 8; [use (8-27) for
  // linear and (8-33) for non-linear molecules]
  // S = S_trans + S_rot + S_vib + S_el
  ////////////////////////////////////////////////////////////////////////

  // compute the mass of the molecule (in kg)
  double mass = 0.0;
  for (i=0; i<mol_->natom(); i++) {
      mass += mol_->mass(i);
      }
  mass *= amu_to_kg;

  // compute principal moments of inertia (pmi) in amu*angstrom^2
  double pmi[3];
  mol_->principal_moments_of_inertia(pmi);

  // find out if molecule is linear (if smallest pmi < 1.0e-5 amu angstrom^2)
  // (elements of pmi are sorted in order smallest to largest)
  int linear = 0;
  if (pmi[0] < 1.0e-5) linear = 1;

  // compute the symmetry number sigma;
  // for linear molecules: sigma = 2 (D_inf_h), sigma = 1 (C_inf_v)
  // for non-linear molecules: sigma = # of rot. in pt. grp, including E
  int sigma;
  CharacterTable ct = pg_->char_table();