energy.cc

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

}

RefSymmSCMatrix
MolecularEnergy::hessian()
{
  if (hess_.null()) return hessian_.result();

  if (hessian_.computed()) return hessian_.result();

  int nullmole = (hess_->energy() == 0);
  this->reference();
  if (nullmole) hess_->set_energy(this);
  RefSymmSCMatrix xhess = hess_->cartesian_hessian();
  if (nullmole) hess_->set_energy(0);
  this->dereference();
  set_hessian(xhess);
  return hessian_.result();
}

int
MolecularEnergy::hessian_implemented() const
{
  return hess_.nonnull();
}

void
MolecularEnergy::symmetry_changed()
{
  obsolete();
}

Ref<NonlinearTransform>
MolecularEnergy::change_coordinates()
{
  if (!mc_) return 0;
  Ref<NonlinearTransform> t = mc_->change_coordinates();
  do_change_coordinates(t);
  return t;
}

void
MolecularEnergy::print_natom_3(const RefSCVector &v,
                               const char *title, ostream&o) const
{
  int precision = 10;
  int lwidth = precision + 4;
  int n = v.n()/3;
  if (title) {
    o << indent << title << endl;
    o << incindent;
  }
  for (int i=0,ii=0; i<n; i++) {
    o << indent
      << scprintf("%4d %3s",
                  i+1,AtomInfo::symbol(molecule()->Z(i)));
    for (int j=0; j<3; j++,ii++) {
      o << scprintf(" % *.*f", lwidth,precision,double(v(ii)));
    }
    o << endl;
  }
  if (title) {
    o << decindent;
  }
  o.flush();
}

void
MolecularEnergy::print_natom_3(double **vn3,
                               const char *title, ostream&o) const
{
  int precision = 10;
  int lwidth = precision + 4;
  int n = molecule()->natom();
  if (title) {
    o << indent << title << endl;
    o << incindent;
  }
  for (int i=0; i<n; i++) {
    o << indent
      << scprintf("%4d %3s",
                  i+1,AtomInfo::symbol(molecule()->Z(i)));
    for (int j=0; j<3; j++) {
      o << scprintf(" % *.*f", lwidth,precision,double(vn3[i][j]));
    }
    o << endl;
  }
  if (title) {
    o << decindent;
  }
  o.flush();
}

void
MolecularEnergy::print_natom_3(double *vn3,
                               const char *title, ostream&o) const
{
  int precision = 10;
  int lwidth = precision + 4;
  int n = molecule()->natom();
  if (title) {
    o << indent << title << endl;
    o << incindent;
  }
  for (int i=0; i<n; i++) {
    o << indent
      << scprintf("%4d %3s",
                  i+1,AtomInfo::symbol(molecule()->Z(i)));
    for (int j=0; j<3; j++) {
      o << scprintf(" % *.*f", lwidth,precision,double(vn3[3*i+j]));
    }
    o << endl;
  }
  if (title) {
    o << decindent;
  }
  o.flush();
}

void
MolecularEnergy::print(ostream&o) const
{
  Function::print(o);
  if (mc_.nonnull()) {
      o << indent << "Molecular Coordinates:\n" << incindent;
      mc_->print(o);
      o << decindent;
    }
  else {
      o << indent << "Molecule:\n" << incindent;
      mol_->print(o);
      o << decindent << endl;
    }
}

/////////////////////////////////////////////////////////////////
// SumMolecularEnergy

static ClassDesc SumMolecularEnergy_cd(
  typeid(SumMolecularEnergy),"SumMolecularEnergy",1,"public MolecularEnergy",
  0, create<SumMolecularEnergy>, create<SumMolecularEnergy>);

SumMolecularEnergy::SumMolecularEnergy(const Ref<KeyVal> &keyval):
  MolecularEnergy(keyval)
{
  n_ = keyval->count("mole");
  mole_ = new Ref<MolecularEnergy>[n_];
  coef_ = new double[n_];
  for (int i=0; i<n_; i++) {
      mole_[i] << keyval->describedclassvalue("mole",i);
      coef_[i] = keyval->intvalue("coef",i);
      if (mole_[i].null()
          || mole_[i]->molecule()->natom() != molecule()->natom()) {
          ExEnv::errn() << "SumMolecularEnergy: a mole is null or has a molecule"
               << " with the wrong number of atoms" << endl;
          abort();
        }
    }
}

SumMolecularEnergy::SumMolecularEnergy(StateIn&s):
  MolecularEnergy(s)
{
  s.get(n_);
  coef_ = new double[n_];
  mole_ = new Ref<MolecularEnergy>[n_];
  s.get_array_double(coef_,n_);
  for (int i=0; i<n_; i++) {
      mole_[i] << SavableState::restore_state(s);
    }
}

void
SumMolecularEnergy::save_data_state(StateOut&s)
{
  MolecularEnergy::save_data_state(s);
  s.put(n_);
  s.put_array_double(coef_,n_);
  for (int i=0; i<n_; i++) {
      SavableState::save_state(mole_[i].pointer(),s);
    }
}

SumMolecularEnergy::~SumMolecularEnergy()
{
  delete[] mole_;
  delete[] coef_;
}

int
SumMolecularEnergy::value_implemented() const
{
  for (int i=0; i<n_; i++) {
      if (!mole_[i]->value_implemented()) return 0;
    }
  return 1;
}

int
SumMolecularEnergy::gradient_implemented() const
{
  for (int i=0; i<n_; i++) {
      if (!mole_[i]->gradient_implemented()) return 0;
    }
  return 1;
}

int
SumMolecularEnergy::hessian_implemented() const
{
  for (int i=0; i<n_; i++) {
      if (!mole_[i]->hessian_implemented()) return 0;
    }
  return 1;
}

void
SumMolecularEnergy::set_x(const RefSCVector&v)
{
  MolecularEnergy::set_x(v);
  for (int i=0; i<n_; i++) {
      mole_[i]->set_x(v);
    }
}

void
SumMolecularEnergy::compute()
{
  int i;

  int *old_do_value = new int[n_];
  int *old_do_gradient = new int[n_];
  int *old_do_hessian = new int[n_];

  for (i=0; i<n_; i++)
      old_do_value[i] = mole_[i]->do_value(value_.compute());
  for (i=0; i<n_; i++)
      old_do_gradient[i]=mole_[i]->do_gradient(gradient_.compute());
  for (i=0; i<n_; i++)
      old_do_hessian[i] = mole_[i]->do_hessian(hessian_.compute());

  ExEnv::out0() << indent
       << "SumMolecularEnergy: compute" << endl;

  ExEnv::out0() << incindent;

  if (value_needed()) {
      double val = 0.0;
      for (i=0; i<n_; i++) {
          val += coef_[i] * mole_[i]->value();
        }
      ExEnv::out0() << endl << indent
           << "SumMolecularEnergy =" << endl;
      for (i=0; i<n_; i++) {
          ExEnv::out0() << indent
               << scprintf("  %c % 16.12f * % 16.12f",
                           (i==0?' ':'+'),
                           coef_[i], mole_[i]->value())
               << endl;
        }
      ExEnv::out0() << indent
           << scprintf("  = % 16.12f", val) << endl;
      set_energy(val);
    }
  if (gradient_needed()) {
      RefSCVector gradientvec = matrixkit()->vector(moldim());
      gradientvec->assign(0.0);
      for (i=0; i<n_; i++)
          gradientvec.accumulate(coef_[i] * mole_[i]->gradient());
      set_gradient(gradientvec);
    }
  if (hessian_needed()) {
      RefSymmSCMatrix hessianmat = matrixkit()->symmmatrix(moldim());
      hessianmat->assign(0.0);
      for (i=0; i<n_; i++)
          hessianmat.accumulate(coef_[i] * mole_[i]->hessian());
      set_hessian(hessianmat);
    }

  ExEnv::out0() << decindent;

  for (i=0; i<n_; i++) mole_[i]->do_value(old_do_value[i]);
  for (i=0; i<n_; i++) mole_[i]->do_gradient(old_do_gradient[i]);
  for (i=0; i<n_; i++) mole_[i]->do_hessian(old_do_hessian[i]);

  delete[] old_do_value;
  delete[] old_do_gradient;
  delete[] old_do_hessian;
}

/////////////////////////////////////////////////////////////////
// MolEnergyConvergence

static ClassDesc MolEnergyConvergence_cd(
  typeid(MolEnergyConvergence),"MolEnergyConvergence",3,"public Convergence",
  0, create<MolEnergyConvergence>, create<MolEnergyConvergence>);

MolEnergyConvergence::MolEnergyConvergence()
{
  set_defaults();
}

MolEnergyConvergence::MolEnergyConvergence(StateIn&s):
  SavableState(s),
  Convergence(s)
{
  if (s.version(::class_desc<MolEnergyConvergence>()) >= 2)
      s.get(cartesian_);
  if (s.version(::class_desc<MolEnergyConvergence>()) >= 3)
      mole_ << SavableState::restore_state(s);
}

MolEnergyConvergence::MolEnergyConvergence(const Ref<KeyVal>&keyval)
{
  mole_ << keyval->describedclassvalue("energy");
  if (mole_.null()) {
      ExEnv::errn() << "MolEnergyConvergence(const Ref<KeyVal>&keyval): "
           << "require an energy keyword of type MolecularEnergy"
           << endl;
      abort();
    }

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

  use_max_disp_ = keyval->exists("max_disp");
  use_max_grad_ = keyval->exists("max_grad");
  use_rms_disp_ = keyval->exists("rms_disp");
  use_rms_grad_ = keyval->exists("rms_grad");
  use_graddisp_ = keyval->exists("graddisp");
  if (use_max_disp_) max_disp_ = keyval->doublevalue("max_disp");
  if (use_max_grad_) max_grad_ = keyval->doublevalue("max_grad");
  if (use_rms_disp_) rms_disp_ = keyval->doublevalue("rms_disp");
  if (use_rms_grad_) rms_grad_ = keyval->doublevalue("rms_grad");
  if (use_graddisp_) graddisp_ = keyval->doublevalue("graddisp");

  if (!use_max_disp_ && !use_max_grad_
      && !use_rms_disp_ && !use_rms_grad_
      && !use_graddisp_) {
      set_defaults();
    }
}

MolEnergyConvergence::~MolEnergyConvergence()
{
}

void
MolEnergyConvergence::save_data_state(StateOut&s)
{
  Convergence::save_data_state(s);
  s.put(cartesian_);
  SavableState::save_state(mole_.pointer(),s);
}

void
MolEnergyConvergence::set_defaults()
{
  use_max_disp_ = 1;
  use_max_grad_ = 1;
  use_rms_disp_ = 0;
  use_rms_grad_ = 0;
  use_graddisp_ = 1;
  max_disp_ = 1.0e-4;
  max_grad_ = 1.0e-4;
  graddisp_ = 1.0e-4;
}

void
MolEnergyConvergence::get_x(const Ref<Function> &f)
{
  Ref<MolecularEnergy> m; m << f;
  if (cartesian_ && m.nonnull() && m->molecularcoor().nonnull()) {
      x_ = m->get_cartesian_x();
    }
  else {
      x_ = f->get_x();
    }
}


void
MolEnergyConvergence::set_nextx(const RefSCVector& x)
{
  if (cartesian_ && mole_.nonnull() && mole_->molecularcoor().nonnull()) {
      Ref<Molecule> mol = new Molecule(*(mole_->molecule().pointer()));
      mole_->molecularcoor()->to_cartesian(mol, x);
      nextx_ = mole_->matrixkit()->vector(mole_->moldim());
      int c = 0;
      for (int i=0; i < mol->natom(); i++) {
          nextx_(c) = mol->r(i,0); c++;
          nextx_(c) = mol->r(i,1); c++;
          nextx_(c) = mol->r(i,2); c++;
        }
    }
  else if (mole_.null()) {
      // this only happens after restoring state from old versions
      // of MolEnergyConvergence
      nextx_ = 0;
    }
  else {
      nextx_ = x.copy();
    }
}

void
MolEnergyConvergence::get_grad(const Ref<Function> &f)
{
  Ref<MolecularEnergy> m; m << f;
  if (cartesian_ && m.nonnull() && m->molecularcoor().nonnull()) {
      RefSCVector cartesian_grad = m->get_cartesian_gradient()->copy();
      if (m->molecularcoor()->nconstrained()) {
          // convert the gradient to internal coordinates and back
          // this will project out the fixed coordinates
          RefSCVector internal_grad(m->dimension(), m->matrixkit());
          m->molecularcoor()->to_internal(internal_grad,cartesian_grad);
          m->molecularcoor()->to_cartesian(cartesian_grad,internal_grad);
        }
      grad_ = cartesian_grad;
    }
  else {
      grad_ = f->gradient();
    }
}

int
MolEnergyConvergence::converged()
{
  return Convergence::converged();
}

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

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