energy.cc

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

//
// energy.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 <stdlib.h>
#include <math.h>
#include <stdexcept>

#include <util/misc/formio.h>
#include <util/state/stateio.h>
#include <math/scmat/local.h>
#include <util/keyval/keyval.h>
#include <chemistry/molecule/energy.h>

using namespace std;
using namespace sc;

/////////////////////////////////////////////////////////////////
// MolecularEnergy

static ClassDesc MolecularEnergy_cd(
  typeid(MolecularEnergy),"MolecularEnergy",6,"public Function",
  0, 0, 0);

MolecularEnergy::MolecularEnergy(const MolecularEnergy& mole):
  Function(mole)
{
  print_molecule_when_changed_ = mole.print_molecule_when_changed_;
  mc_ = mole.mc_;
  moldim_ = mole.moldim_;
  mol_ = mole.mol_;
  initial_pg_ = new PointGroup(mol_->point_group());
  ckpt_ = mole.ckpt_;
  ckpt_file_ = strdup(mole.ckpt_file_);
  ckpt_freq_ = mole.ckpt_freq_;
}

MolecularEnergy::MolecularEnergy(const Ref<KeyVal>&keyval):
  Function(keyval,1.0e-6,1.0e-6,1.0e-4)
{
  // The following code is a solaris workshop 5.0 hack, since it doesn't
  // seem to pass the right arguments in the Function CTOR.  This code can
  // be deleted with other C++ compilers.
  KeyValValuedouble funcaccval(1.0e-6);
  value_.set_desired_accuracy(keyval->doublevalue("value_accuracy",
                                                  funcaccval));
  if (value_.desired_accuracy() < DBL_EPSILON)
    value_.set_desired_accuracy(DBL_EPSILON);
  KeyValValuedouble gradaccval(1.0e-6);
  gradient_.set_desired_accuracy(keyval->doublevalue("gradient_accuracy",
                                                     gradaccval));
  if (gradient_.desired_accuracy() < DBL_EPSILON)
    gradient_.set_desired_accuracy(DBL_EPSILON);
  KeyValValuedouble hessaccval(1.0e-4);
  hessian_.set_desired_accuracy(keyval->doublevalue("hessian_accuracy",
                                                    hessaccval));
  if (hessian_.desired_accuracy() < DBL_EPSILON)
    hessian_.set_desired_accuracy(DBL_EPSILON);
  // end of solaris workshop 5.0 hack

  print_molecule_when_changed_
      = keyval->booleanvalue("print_molecule_when_changed");
  if (keyval->error() != KeyVal::OK) print_molecule_when_changed_ = 1;

  mol_ << keyval->describedclassvalue("molecule");
  if (mol_.null()) {
      ExEnv::errn() << indent << "MolecularEnergy(Keyval): no molecule found"
           << endl;
      abort();
    }

  initial_pg_ = new PointGroup(mol_->point_group());

  hess_ << keyval->describedclassvalue("hessian");

  guesshess_ << keyval->describedclassvalue("guess_hessian");

  moldim_ = new SCDimension(3 * mol_->natom(), "3Natom");

  // the molecule coordinate object needs moldim_
  // so constract a keyval that has it
  Ref<AssignedKeyVal> assignedkeyval = new AssignedKeyVal;
  Ref<DescribedClass> dc = moldim_.pointer();
  assignedkeyval->assign("natom3", dc);
  dc = matrixkit();
  assignedkeyval->assign("matrixkit", dc);
  Ref<KeyVal> asskeyval(assignedkeyval.pointer());
  Ref<KeyVal> aggkeyval = new AggregateKeyVal(asskeyval, keyval);

  // Don't bother with internal coordinates if there is only 1 atom
  if (mol_->natom() > 1) {
      mc_ << aggkeyval->describedclassvalue("coor");
    }

  RefSCDimension dim;
  if (mc_.null()) {
      dim = moldim_;
    }
  else {
      dim = mc_->dim();
    }
  set_dimension(dim);

  ckpt_ = keyval->booleanvalue("checkpoint");
  if (keyval->error() != KeyVal::OK) ckpt_ = false;
  ckpt_file_ = keyval->pcharvalue("checkpoint_file");
  if (keyval->error() != KeyVal::OK) {
    char* filename = SCFormIO::fileext_to_filename(".wfn.ckpt");
    ckpt_file_ = strdup(filename);
    delete[] filename;
  }
  ckpt_freq_ = keyval->intvalue("checkpoint_freq");
  if (keyval->error() != KeyVal::OK) {
    ckpt_freq_ = 1;
  }
  
  do_value(1);
  do_gradient(0);
  do_hessian(0);

  molecule_to_x();
}

MolecularEnergy::~MolecularEnergy()
{
  if (ckpt_file_) free(ckpt_file_);
  ckpt_file_ = 0;
}

MolecularEnergy::MolecularEnergy(StateIn&s):
  SavableState(s),
  Function(s)
{
  mc_ << SavableState::restore_state(s);
  moldim_ << SavableState::restore_state(s);
  mol_ << SavableState::restore_state(s);
  if (s.version(::class_desc<MolecularEnergy>()) >= 2)
      s.get(print_molecule_when_changed_);
  else print_molecule_when_changed_ = 1;
  if (s.version(::class_desc<MolecularEnergy>()) >= 3) {
      hess_ << SavableState::restore_state(s);
      guesshess_ << SavableState::restore_state(s);
    }
  if (s.version(::class_desc<MolecularEnergy>()) >= 4)
      initial_pg_ << SavableState::restore_state(s);
  else initial_pg_ = new PointGroup(mol_->point_group());
  if (s.version(::class_desc<MolecularEnergy>()) >= 5) {
    int ckpt; s.get(ckpt); ckpt_ = (bool)ckpt;
    s.getstring(ckpt_file_);
  }
  else {
    ckpt_ = false;
    char* filename = SCFormIO::fileext_to_filename(".wfn.ckpt");
    ckpt_file_ = strdup(filename);
    delete[] filename;
  }
  if (s.version(::class_desc<MolecularEnergy>()) >= 6)
    s.get(ckpt_freq_);
  else
    ckpt_freq_ = 1;
}

MolecularEnergy&
MolecularEnergy::operator=(const MolecularEnergy& mole)
{
  Function::operator=(mole);
  mc_ = mole.mc_;
  moldim_ = mole.moldim_;
  mol_ = mole.mol_;
  print_molecule_when_changed_ = mole.print_molecule_when_changed_;
  initial_pg_ = new PointGroup(mole.initial_pg_);
  return *this;
}

void
MolecularEnergy::save_data_state(StateOut&s)
{
  Function::save_data_state(s);
  SavableState::save_state(mc_.pointer(),s);
  SavableState::save_state(moldim_.pointer(),s);
  SavableState::save_state(mol_.pointer(),s);
  s.put(print_molecule_when_changed_);
  SavableState::save_state(hess_.pointer(),s);
  SavableState::save_state(guesshess_.pointer(),s);
  SavableState::save_state(initial_pg_.pointer(),s);
  s.put((int)ckpt_);
  s.putstring(ckpt_file_);
  s.put(ckpt_freq_);
}

void
MolecularEnergy::set_checkpoint()
{
  ckpt_ = true;
}

void
MolecularEnergy::set_checkpoint_file(const char *path)
{
  if (ckpt_file_) free(ckpt_file_);
  if (path) {
    ckpt_file_ = strdup(path);
  } else
    ckpt_file_ = 0;
}

void
MolecularEnergy::set_checkpoint_freq(int freq)
{
  if (freq >= 1)
    ckpt_freq_ = freq;
  else
    throw std::runtime_error("MolecularEnergy::set_checkpoint_freq() -- invalid checkpointing frequency");
}

bool
MolecularEnergy::if_to_checkpoint() const
{
  return ckpt_;
}

const char*
MolecularEnergy::checkpoint_file() const
{
  return strdup(ckpt_file_);
}

int
MolecularEnergy::checkpoint_freq() const
{
  return ckpt_freq_;
}

void
MolecularEnergy::failure(const char * msg)
{
  ExEnv::err0() << indent << "MolecularEnergy::failure: " << msg << endl;
  abort();
}

void
MolecularEnergy::set_energy(double e)
{
  set_value(e);
}

double
MolecularEnergy::energy()
{
  return value();
}

void
MolecularEnergy::set_gradient(RefSCVector&g)
{
  cartesian_gradient_ = g.copy();
  if (mc_.null()) {
    Function::set_gradient(g);
  } else {
    RefSCVector grad(dimension(), matrixkit());
    mc_->to_internal(grad,g);
    Function::set_gradient(grad);
  }
}

void
MolecularEnergy::set_hessian(RefSymmSCMatrix&h)
{
  cartesian_hessian_ = h.copy();
  if (mc_.null()) {
    Function::set_hessian(h);
  } else {
    RefSymmSCMatrix hess(dimension(), matrixkit());
    mc_->to_internal(hess,h);
    Function::set_hessian(hess);
  }
}

void
MolecularEnergy::x_to_molecule()
{
  RefSCVector x = get_x_no_copy();

  if (mc_.null()) {
    int c = 0;
    
    for (int i=0; i<mol_->natom(); i++) {
      mol_->r(i,0) = x(c); c++;
      mol_->r(i,1) = x(c); c++;
      mol_->r(i,2) = x(c); c++;
    }
  } else {
    mc_->to_cartesian(get_x_no_copy());
  }
}

void
MolecularEnergy::molecule_to_x()
{
  if (mc_.null()) {
    RefSCVector cartesian(moldim(),matrixkit());
    int c = 0;
    for (int i=0; i < mol_->natom(); i++) {
      cartesian(c) = mol_->r(i,0); c++;
      cartesian(c) = mol_->r(i,1); c++;
      cartesian(c) = mol_->r(i,2); c++;
    }
    Function::set_x(cartesian);
  } else {
    mc_->to_internal(get_x_reference());
  }
}

void
MolecularEnergy::set_x(const RefSCVector&v)
{
  Function::set_x(v);
  x_to_molecule();
  if (print_molecule_when_changed_) {
      ExEnv::out0() << endl << indent << class_name()
           << ": changing atomic coordinates:" << endl;
      molecule()->print();
    }
}

RefSCVector
MolecularEnergy::get_cartesian_x()
{
  RefSCVector cartesian(moldim(),matrixkit());
  int c = 0;
  for (int i=0; i < mol_->natom(); i++) {
      cartesian(c) = mol_->r(i,0); c++;
      cartesian(c) = mol_->r(i,1); c++;
      cartesian(c) = mol_->r(i,2); c++;
    }
  return cartesian;
}

RefSCVector
MolecularEnergy::get_cartesian_gradient()
{
  gradient();
  if (cartesian_gradient_.null()) {
      ExEnv::errn() << "MolecularEnergy::get_cartesian_gradient(): "
           << "cartesian gradient not available"
           << endl;
      abort();
    }
  return cartesian_gradient_;
}

RefSymmSCMatrix
MolecularEnergy::get_cartesian_hessian()
{
  hessian();
  if (cartesian_hessian_.null()) {
      ExEnv::errn() << "MolecularEnergy::get_cartesian_hessian(): "
           << "cartesian hessian not available"
           << endl;
      abort();
    }
  return cartesian_hessian_;
}

RefSCDimension
MolecularEnergy::moldim() const
{
  return moldim_;
}

Ref<Molecule>
MolecularEnergy::molecule() const
{
  return mol_;
}

void
MolecularEnergy::guess_hessian(RefSymmSCMatrix&hessian)
{
  if (guesshess_.nonnull()) {
      int nullmole = (guesshess_->energy() == 0);
      this->reference();
      if (nullmole) guesshess_->set_energy(this);
      RefSymmSCMatrix xhess = guesshess_->cartesian_hessian();
      if (nullmole) guesshess_->set_energy(0);
      this->dereference();
      if (mc_.nonnull()) {
          mc_->to_internal(hessian, xhess);
        }
      else {
          hessian.assign(xhess);
        }
    }
  else if (mc_.nonnull()) {
      mc_->guess_hessian(hessian);
    }
  else {
      Function::guess_hessian(hessian);
    }
}

RefSymmSCMatrix
MolecularEnergy::inverse_hessian(RefSymmSCMatrix&hessian)
{
  if (mc_.nonnull()) {
      return mc_->inverse_hessian(hessian);
    }
  else {
      return Function::inverse_hessian(hessian);
    }