coor.cc

来自「大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CH」· CC 代码 · 共 1,155 行 · 第 1/2 页

//
// coor.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 <set>

#include <math.h>

#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <math/scmat/matrix.h>
#include <chemistry/molecule/molecule.h>
#include <chemistry/molecule/coor.h>
#include <chemistry/molecule/simple.h>
#include <chemistry/molecule/localdef.h>

#include <util/container/bitarray.h>

using namespace std;
using namespace sc;

///////////////////////////////////////////////////////////////////////////
// members of IntCoor

double IntCoor::bohr_conv = 0.52917706;
double IntCoor::radian_conv = 180.0/3.14159265358979323846;

static ClassDesc IntCoor_cd(
  typeid(IntCoor),"IntCoor",1,"public SavableState",
  0, 0, 0);

IntCoor::IntCoor(const char *re):
  label_(0), value_(0.0)
{
  if (!re) re = "noname";
  label_=new char[strlen(re)+1]; strcpy(label_,re);
}

IntCoor::IntCoor(const IntCoor& c):
  label_(0)
{
  value_ = c.value_;
  if (c.label_) label_ = strcpy(new char[strlen(c.label_)+1],c.label_);
}

IntCoor::IntCoor(const Ref<KeyVal>&keyval)
{
  label_ = keyval->pcharvalue("label");
  value_ = keyval->doublevalue("value");

  char* unit = keyval->pcharvalue("unit");
  if (unit) {
      if (!strcmp(unit, "bohr")) {
        }
      else if (!strcmp(unit, "angstrom")) {
          value_ /= bohr_conv;
        }
      else if (!strcmp(unit, "radian")) {
        }
      else if (!strcmp(unit, "degree")) {
          value_ *= M_PI/180.0;
        }
      else {
          ExEnv::err0() << indent
               << "IntCoor::IntCoor(KeyVal): unknown unit = \""
               << unit << "\"\n";
          abort();
        }
      delete[] unit;
    }
}

IntCoor::IntCoor(StateIn& si):
  SavableState(si)
{
  si.get(value_);
  si.getstring(label_);
}

IntCoor::~IntCoor()
{
  if (label_) delete[] label_;
}

void
IntCoor::save_data_state(StateOut& so)
{
  so.put(value_);
  so.putstring(label_);
}

const char*
IntCoor::label() const
{
  return label_;
}

double
IntCoor::value() const
{
  return value_;
}

void
IntCoor::set_value(double v)
{
  value_ = v;
}

void
IntCoor::print(ostream &o) const
{
  print_details(0,o);
}

void
IntCoor::print_details(const Ref<Molecule> &mol, ostream& os) const
{
  os.setf(ios::fixed,ios::floatfield);
  os.precision(10);
  os.setf(ios::left,ios::adjustfield);
  os.width(10);

  os << indent
     << scprintf("%-5s \"%10s\" %15.10f\n",ctype(),label(),preferred_value());
}

double
IntCoor::preferred_value() const
{
  return value_;
}

///////////////////////////////////////////////////////////////////////////
// members of SetIntCoor

static ClassDesc SetIntCoor_cd(
  typeid(SetIntCoor),"SetIntCoor",1,"public SavableState",
  create<SetIntCoor>, create<SetIntCoor>, create<SetIntCoor>);

SetIntCoor::SetIntCoor()
{
}

SetIntCoor::SetIntCoor(const Ref<KeyVal>& keyval)
{
  int n = keyval->count();

  Ref<IntCoorGen> gen; gen << keyval->describedclassvalue("generator");

  if (gen.null() && !n) {
      ExEnv::err0() << indent << "SetIntCoor::SetIntCoor: bad input\n";
      abort();
    }

  if (gen.nonnull()) {
      // Make sure that gen doesn't delete me before my reference
      // count gets incremented.
      this->reference();
      gen->generate(this);
      // Now it is safe to decrement my reference count back down to zero.
      this->dereference();
    }

  for (int i=0; i<n; i++) {
      Ref<IntCoor> coori; coori << keyval->describedclassvalue(i);
      coor_.push_back(coori);
    }
}

SetIntCoor::SetIntCoor(StateIn& s):
  SavableState(s)
{
  int n;
  s.get(n);

  Ref<IntCoor> tmp;
  for (int i=0; i<n; i++) {
      tmp << SavableState::restore_state(s);
      coor_.push_back(tmp);
    }
}

SetIntCoor::~SetIntCoor()
{
}

void
SetIntCoor::save_data_state(StateOut& s)
{
  int n = coor_.size();
  s.put(n);

  for (int i=0; i<n; i++) {
      SavableState::save_state(coor_[i].pointer(),s);
    }
}

void
SetIntCoor::add(const Ref<IntCoor>& coor)
{
  coor_.push_back(coor);
}

void
SetIntCoor::add(const Ref<SetIntCoor>& coor)
{
  for (int i=0; i<coor->n(); i++) {
      coor_.push_back(coor->coor(i));
    }
}

void
SetIntCoor::pop()
{
  coor_.pop_back();
}

int
SetIntCoor::n() const
{
  return coor_.size();
}

Ref<IntCoor>
SetIntCoor::coor(int i) const
{
  return coor_[i];
}

// compute the bmatrix by finite displacements
void
SetIntCoor::fd_bmat(const Ref<Molecule>& mol,RefSCMatrix& fd_bmatrix)
{
  Ref<SCMatrixKit> kit = fd_bmatrix.kit();

  fd_bmatrix.assign(0.0);
  
  int i;
  Molecule& m = * mol.pointer();

  const double cart_disp = 0.01;

  RefSCDimension dn3(fd_bmatrix.coldim());
  RefSCDimension dnc(fd_bmatrix.rowdim());
  int n3 = dn3.n();
  int nc = dnc.n();
  RefSCVector internal(dnc,kit);
  RefSCVector internal_p(dnc,kit);
  RefSCVector internal_m(dnc,kit);

  // the internal coordinates
  update_values(mol);
  for (i=0; i<nc; i++) {
      internal(i) = coor_[i]->value();
    }

  // the finite displacement bmat
  for (i=0; i<n3; i++) {
      // the plus displacement
      m.r(i/3,i%3) += cart_disp;
      update_values(mol);
      int j;
      for (j=0; j<nc; j++) {
          internal_p(j) = coor_[j]->value();
        }
      // the minus displacement
      m.r(i/3,i%3) -= 2.0*cart_disp;
      update_values(mol);
      for (j=0; j<nc; j++) {
          internal_m(j) = coor_[j]->value();
        }
      // reset the cartesian coordinate to its original value
      m.r(i/3,i%3) += cart_disp;

      // construct the entries in the finite displacement bmat
      for (j=0; j<nc; j++) {
          fd_bmatrix(j,i) = (internal_p(j)-internal_m(j))/(2.0*cart_disp);
        }
    }
}

void
SetIntCoor::bmat(const Ref<Molecule>& mol, RefSCMatrix& bmat)
{
  bmat.assign(0.0);

  int i, ncoor = n();

  RefSCVector bmatrow(bmat.coldim(),bmat.kit());
  // send the rows of the b matrix to each of the coordinates
  for (i=0; i<ncoor; i++) {
      bmatrow.assign(0.0);
      coor_[i]->bmat(mol,bmatrow);
      bmat.assign_row(bmatrow,i);
    }
}

void
SetIntCoor::guess_hessian(Ref<Molecule>& mol,RefSymmSCMatrix& hessian)
{
  int ncoor = hessian.n();

  hessian.assign(0.0);
  for (int i=0; i<ncoor; i++) {
      hessian(i,i) = coor_[i]->force_constant(mol);
    }
}

void
SetIntCoor::print_details(const Ref<Molecule> &mol, ostream& os) const
{
  int i;

  for(i=0; i<coor_.size(); i++) {
      coor_[i]->print_details(mol,os);
    }
}

void
SetIntCoor::update_values(const Ref<Molecule>&mol)
{
  for (int i=0; i<coor_.size(); i++) {
      coor_[i]->update_value(mol);
    }
}

void
SetIntCoor::values_to_vector(const RefSCVector&v)
{
  for (int i=0; i<coor_.size(); i++) {
      v(i) = coor_[i]->value();
    }
}  

void
SetIntCoor::clear()
{
  coor_.clear();
}

///////////////////////////////////////////////////////////////////////////
// members of SumIntCoor

static ClassDesc SumIntCoor_cd(
  typeid(SumIntCoor),"SumIntCoor",1,"public IntCoor",
  0, create<SumIntCoor>, create<SumIntCoor>);

SumIntCoor::SumIntCoor(const char* label):
  IntCoor(label)
{
}

SumIntCoor::SumIntCoor(const Ref<KeyVal>&keyval):
  IntCoor(keyval)
{
  static const char* coor = "coor";
  static const char* coef = "coef";
  int n = keyval->count(coor);
  int ncoef = keyval->count(coef);
  if (n != ncoef || !n) {
      ExEnv::err0() << indent << "SumIntCoor::SumIntCoor: bad input\n";
      abort();
    }

  for (int i=0; i<n; i++) {
      double coe = keyval->doublevalue(coef,i);
      Ref<IntCoor> coo; coo << keyval->describedclassvalue(coor,i);
      add(coo,coe);
    }
}

SumIntCoor::SumIntCoor(StateIn&s):
  IntCoor(s)
{
  int n;
  s.get(n);

  coef_.resize(n);
  coor_.resize(n);
  for (int i=0; i<n; i++) {
      s.get(coef_[i]);
      coor_[i] << SavableState::restore_state(s);
    }
}

SumIntCoor::~SumIntCoor()
{
}

void
SumIntCoor::save_data_state(StateOut&s)
{
  int n = coef_.size();
  IntCoor::save_data_state(s);
  s.put(int(coef_.size()));

  for (int i=0; i<n; i++) {
      s.put(coef_[i]);
      SavableState::save_state(coor_[i].pointer(),s);
    }
}

int
SumIntCoor::n()
{
  return coef_.size();
}

void
SumIntCoor::add(Ref<IntCoor>&coor,double coef)
{
  // if a sum is added to a sum, unfold the nested sum
  SumIntCoor* scoor = dynamic_cast<SumIntCoor*>(coor.pointer());
  if (scoor) {
      int l = scoor->coor_.size();
      for (int i=0; i<l; i++) {
          add(scoor->coor_[i],coef * scoor->coef_[i]);
        }
    }
  else {
      int l = coef_.size();
      for (int i=0; i<l; i++) {
          if (coor_[i]->equivalent(coor)) {
              coef_[i] += coef;
              return;
            }
        }
      coef_.resize(l+1);
      coor_.resize(l+1);
      coef_[l] = coef;
      coor_[l] = coor;
    }
}

int
SumIntCoor::equivalent(Ref<IntCoor>&c)
{
  return 0;
}

// this normalizes and makes the biggest coordinate positive
void
SumIntCoor::normalize()
{
  int i;
  int n = coef_.size();
  double norm = 0.0;

  double biggest = 0.0;
  for (i=0; i<n; i++) {
      norm += coef_[i] * coef_[i];
      if (fabs(biggest) < fabs(coef_[i])) biggest = coef_[i];
    }
  norm = (biggest < 0.0? -1.0:1.0)/sqrt(norm);

  for (i=0; i<n; i++) {
      coef_[i] = coef_[i]*norm;
    }
}

double
SumIntCoor::preferred_value() const
{
  return value_;
}

const char*
SumIntCoor::ctype() const
{
  return "SUM";
}

void
SumIntCoor::print_details(const Ref<Molecule> &mol, ostream& os) const
{
  int initial_indent = SCFormIO::getindent(os);
  int i;

  os << indent
     << scprintf("%-5s %10s %14.10f\n",ctype(),
                 (label()?label():""), preferred_value());

  for(i=0; i<coor_.size(); i++) {
      os << incindent
         << indent << scprintf("%14.10f ",coef_[i]);

      SCFormIO::setindent(os, SCFormIO::getindent(os) + 15);
      os << skipnextindent;
      coor_[i]->print_details(mol,os);
      SCFormIO::setindent(os, initial_indent);
    }
}

// the SumIntCoor should be normalized before this is called.
double
SumIntCoor::force_constant(Ref<Molecule>&molecule)
{
  double fc = 0.0;
  
  for (int i=0; i<n(); i++) {
      fc += coef_[i] * coef_[i] * coor_[i]->force_constant(molecule);
    }

  return fc;
}

void
SumIntCoor::update_value(const Ref<Molecule>&molecule)
{
  int i, l = n();

  value_ = 0.0;
  for (i=0; i<l; i++) {
      coor_[i]->update_value(molecule);
#if OLD_BMAT
      if (dynamic_cast<StreSimpleCo*>(coor_[i]))
        value_ += coef_[i] * dynamic_cast<StreSimpleCo*>(coor_[i])->angstrom();
      else
#endif        
      value_ += coef_[i] * coor_[i]->value();
    }
}

void
SumIntCoor::bmat(const Ref<Molecule>&molecule,RefSCVector&bmat,double coef)
{
  int i, l = n();
  
  for (i=0; i<l; i++) {
      coor_[i]->bmat(molecule,bmat,coef*coef_[i]);
    }
}

///////////////////////////////////////////////////////////////////////////
// members of MolecularCoor

static ClassDesc MolecularCoor_cd(
  typeid(MolecularCoor),"MolecularCoor",1,"public SavableState",
  0, 0, 0);

MolecularCoor::MolecularCoor(Ref<Molecule>&mol):
  molecule_(mol)
{
  debug_ = 0;
  matrixkit_ = SCMatrixKit::default_matrixkit();
  dnatom3_ = new SCDimension(3*molecule_->natom());
}

MolecularCoor::MolecularCoor(const Ref<KeyVal>&keyval)
{