fdhess.cc

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

//
// fdhess.cc
//
// Copyright (C) 1997 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 <stdlib.h>
#include <sys/stat.h>

#include <util/misc/formio.h>
#include <util/misc/timer.h>
#include <util/state/stateio.h>
#include <util/state/state_bin.h>
#include <util/group/mstate.h>
#include <util/keyval/keyval.h>
#include <math/scmat/blocked.h>
#include <math/symmetry/corrtab.h>
#include <chemistry/molecule/fdhess.h>

using namespace std;
using namespace sc;

#define DEFAULT_CHECKPOINT  1
#define DEFAULT_RESTART     1

/////////////////////////////////////////////////////////////////
// FinDispMolecularHessian

static ClassDesc FinDispMolecularHessian_cd(
  typeid(FinDispMolecularHessian),"FinDispMolecularHessian",1,"public MolecularHessian",
  0, create<FinDispMolecularHessian>, create<FinDispMolecularHessian>);

FinDispMolecularHessian::FinDispMolecularHessian(const Ref<MolecularEnergy> &e):
  mole_(e)
{
  only_totally_symmetric_ = 0;
  eliminate_cubic_terms_ = 1;
  do_null_displacement_ = 1;
  disp_ = 1.0e-2;
  ndisp_ = 0;
  debug_ = 0;
  gradients_ = 0;
  accuracy_ = disp_/1000;
  restart_ = DEFAULT_RESTART;
  checkpoint_ = DEFAULT_CHECKPOINT;
  checkpoint_file_ = 0;
  restart_file_ = 0;
  checkpoint_file_ = SCFormIO::fileext_to_filename(".ckpt.hess");
  restart_file_ = SCFormIO::fileext_to_filename(".ckpt.hess");
}

FinDispMolecularHessian::FinDispMolecularHessian(const Ref<KeyVal>&keyval):
  MolecularHessian(keyval)
{
  mole_ << keyval->describedclassvalue("energy");

  debug_ = keyval->booleanvalue("debug");

  displacement_point_group_ << keyval->describedclassvalue("point_group");

  disp_ = keyval->doublevalue("displacement",KeyValValuedouble(1.0e-2));

  KeyValValueboolean def_restart(DEFAULT_RESTART);
  KeyValValueboolean def_checkpoint(DEFAULT_CHECKPOINT);
  KeyValValueboolean truevalue(1);
  KeyValValueboolean falsevalue(0);

  restart_ = keyval->booleanvalue("restart", def_restart);
  char *def_file = SCFormIO::fileext_to_filename(".ckpt.hess");
  KeyValValueString def_restart_file(def_file,KeyValValueString::Steal);
  restart_file_ = keyval->pcharvalue("restart_file", def_restart_file);
  checkpoint_ = keyval->booleanvalue("checkpoint", def_checkpoint);
  checkpoint_file_ = keyval->pcharvalue("checkpoint_file", def_restart_file);
  only_totally_symmetric_ = keyval->booleanvalue("only_totally_symmetric",
                                                 falsevalue);
  eliminate_cubic_terms_ = keyval->booleanvalue("eliminate_cubic_terms",
                                                truevalue);

  do_null_displacement_ = keyval->booleanvalue("do_null_displacement",
                                               truevalue);

  accuracy_ = keyval->doublevalue("gradient_accuracy",
                                  KeyValValuedouble(disp_/1000));

  gradients_ = 0;
  ndisp_ = 0;
}

FinDispMolecularHessian::FinDispMolecularHessian(StateIn&s):
  SavableState(s),
  MolecularHessian(s)
{
  mole_ << SavableState::restore_state(s);
  s.get(checkpoint_);
  s.get(debug_);
  s.get(accuracy_);
  s.getstring(checkpoint_file_);
  s.getstring(restart_file_);

  gradients_ = 0;
  restore_displacements(s);
}

FinDispMolecularHessian::~FinDispMolecularHessian()
{
  delete[] gradients_;
  delete[] checkpoint_file_;
  delete[] restart_file_;
}

void
FinDispMolecularHessian::save_data_state(StateOut&s)
{
  MolecularHessian::save_data_state(s);

  SavableState::save_state(mole_.pointer(),s);
  s.put(checkpoint_);
  s.put(debug_);
  s.put(accuracy_);
  s.putstring(checkpoint_file_);
  s.putstring(restart_file_);

  checkpoint_displacements(s);
}

void
FinDispMolecularHessian::init()
{
  if (mole_.null()) return;

  mol_ = mole_->molecule();

  if (displacement_point_group_.null()) {
    displacement_point_group_
      = new PointGroup(*mol_->point_group().pointer());
    }

  nirrep_ = displacement_point_group_->char_table().nirrep();

  original_point_group_ = mol_->point_group();
  original_geometry_ = matrixkit()->vector(d3natom());

  int i, coor;
  for (i=0, coor=0; i<mol_->natom(); i++) {
      for (int j=0; j<3; j++, coor++) {
          original_geometry_(coor) = mol_->r(i,j);
        }
    }

  ndisp_ = 0;
  symbasis_ = cartesian_to_symmetry(mol_,
                                      displacement_point_group_,
                                      matrixkit());
  delete[] gradients_;
  gradients_ = new RefSCVector[ndisplace()];
}

void
FinDispMolecularHessian::set_energy(const Ref<MolecularEnergy> &energy)
{
  mole_ = energy;
}

MolecularEnergy*
FinDispMolecularHessian::energy() const
{
  return mole_.pointer();
}

void
FinDispMolecularHessian::restart()
{
  int statresult, statsize;
  Ref<MessageGrp> grp = MessageGrp::get_default_messagegrp();
  if (grp->me() == 0) {
    struct stat sb;
    statresult = stat(restart_file_,&sb);
    statsize = (statresult==0) ? sb.st_size : 0;
    }
  grp->bcast(statsize);
  if (statsize) {
    BcastStateInBin si(grp,restart_file_);
    restore_displacements(si);
    mol_ = mole_->molecule();

    if (ndisplacements_done() >= ndisplace()) {
        restart_=0;
        return;
      }
    }

  if (ndisp_) {
    int irrep, index;
    double coef;
    get_disp(ndisplacements_done(), irrep, index, coef);
    if (irrep != 0 && index != 0) {
        displace(ndisplacements_done());
        mole_->symmetry_changed();
      }
    }
  else {
    init();
    }
  restart_ = 0;
}

void
FinDispMolecularHessian::restore_displacements(StateIn& s)
{
  int i;
  displacement_point_group_ << SavableState::restore_state(s);
  original_point_group_ << SavableState::restore_state(s);
  original_geometry_ = matrixkit()->vector(d3natom());
  original_geometry_.restore(s);

  s.get(disp_);
  s.get(ndisp_);
  s.get(nirrep_);
  s.get(only_totally_symmetric_);
  s.get(eliminate_cubic_terms_);
  s.get(do_null_displacement_);

  if (ndisp_) {
    RefSCDimension symrow, symcol;
    symrow << SavableState::restore_state(s);
    symcol << SavableState::restore_state(s);
    Ref<SCMatrixKit> symkit = new BlockedSCMatrixKit(matrixkit());
    symbasis_ = symkit->matrix(symrow,symcol);
    symbasis_.restore(s);

    delete[] gradients_;
    gradients_ = new RefSCVector[ndisplace()];
    for (i=0; i < ndisp_; i++) {
      int ndisp;
      s.get(ndisp);
      RefSCDimension ddisp = new SCDimension(ndisp);
      gradients_[i] = matrixkit()->vector(ddisp);
      gradients_[i].restore(s);
      }
    }
}

void
FinDispMolecularHessian::checkpoint_displacements(StateOut& s)
{
  int i;
  SavableState::save_state(displacement_point_group_.pointer(),s);
  SavableState::save_state(original_point_group_.pointer(),s);
  original_geometry_.save(s);

  s.put(disp_);
  s.put(ndisp_);
  s.put(nirrep_);
  s.put(only_totally_symmetric_);
  s.put(eliminate_cubic_terms_);
  s.put(do_null_displacement_);

  if (ndisp_) {
    SavableState::save_state(symbasis_.rowdim().pointer(),s);
    SavableState::save_state(symbasis_.coldim().pointer(),s);
    symbasis_.save(s);

    for (i=0; i < ndisp_; i++) {
      s.put(gradients_[i].n());
      gradients_[i].save(s);
      }
    }
}

RefSCMatrix
FinDispMolecularHessian::displacements(int irrep) const
{
  BlockedSCMatrix *bsymbasis = dynamic_cast<BlockedSCMatrix*>(symbasis_.pointer());
  RefSCMatrix block = bsymbasis->block(irrep);
  if (block.null() || (only_totally_symmetric_ && irrep > 0)) {
    RefSCDimension zero = new SCDimension(0);
    block = matrixkit()->matrix(zero,zero);
    return block;
    }