fdhess.cc

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

  return block.t();
}

void
FinDispMolecularHessian::get_disp(int disp, int &irrep,
                                  int &index, double &coef)
{
  int disp_offset = 0;

  if (do_null_displacement_ && disp == 0) {
    irrep = 0;
    coef = 0.0;
    index = -1;
    return;
    }
  disp_offset++;
  // check for +ve totally symmetric displacements
  if (disp < disp_offset + displacements(0).ncol()) {
    irrep = 0;
    coef = 1.0;
    index = disp - disp_offset;
    return;
    }
  disp_offset += displacements(0).ncol();
  // check for -ve totally symmetric displacements
  if (eliminate_cubic_terms_) {
    if (disp < disp_offset + displacements(0).ncol()) {
      irrep = 0;
      coef = -1.0;
      index = disp - disp_offset;
      return;
      }
    disp_offset += displacements(0).ncol();
    }
  for (int i=1; i<nirrep_; i++) {
    if (disp < disp_offset + displacements(i).ncol()) {
      irrep = i;
      coef = 1.0;
      index = disp - disp_offset;
      return;
      }
    disp_offset += displacements(i).ncol();
    }
  ExEnv::err0() << indent
       << "FinDispMolecularHessian::get_disp: bad disp number" << endl;
  abort();
}

int
FinDispMolecularHessian::ndisplace() const
{
  int ndisp = displacements(0).ncol();
  if (eliminate_cubic_terms_) {
    ndisp *= 2;
    }
  for (int i=1; i<nirrep_; i++) {
    ndisp += displacements(i).ncol();
    }
  if (do_null_displacement_) ndisp++;
  return ndisp;
}

void
FinDispMolecularHessian::displace(int disp)
{
  int irrep, index;
  double coef;
  get_disp(disp, irrep, index, coef);

  if (mole_.nonnull()) mole_->obsolete();

  for (int i=0, coor=0; i<mol_->natom(); i++) {
    for (int j=0; j<3; j++, coor++) {
      if (index >= 0) {
        mol_->r(i,j) = original_geometry_(coor)
                       + coef * disp_
                       * displacements(irrep)->get_element(coor,index);

        }
      else {
        mol_->r(i,j) = original_geometry_(coor);
        }
      }
    }

  if (irrep == 0) {
    mol_->set_point_group(original_point_group_);
    }
  else {
    Ref<PointGroup> oldpg = mol_->point_group();
    Ref<PointGroup> newpg = mol_->highest_point_group();
    CorrelationTable corrtab;
    if (corrtab.initialize_table(original_point_group_, newpg)) {
      // something went wrong so use c1 symmetry
      newpg = new PointGroup("c1");
      }
    if (!oldpg->equiv(newpg)) {
      mol_->set_point_group(newpg);
      mole_->symmetry_changed();
      }
    }

#ifdef DEBUG
  ExEnv::out0() << indent
       << "Displacement point group: " << endl
       << incindent << displacement_point_group_ << decindent;
  ExEnv::out0() << indent
       << "Displaced molecule: " << endl
       << incindent << mol_ << decindent;
#endif

  ExEnv::out0() << indent
       << "Displacement is "
       << displacement_point_group_->char_table().gamma(irrep).symbol()
       << " in " << displacement_point_group_->symbol()
       << ".  Using point group "
       << mol_->point_group()->symbol()
       << " for displaced molecule."
       << endl;
}

void
FinDispMolecularHessian::original_geometry()
{
  if (mole_.nonnull()) mole_->obsolete();

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

  if (!mol_->point_group()->equiv(original_point_group_)) {
    mol_->set_point_group(original_point_group_);
    mole_->symmetry_changed();
    }
}

void
FinDispMolecularHessian::set_gradient(int disp, const RefSCVector &grad)
{
  int irrep, index;
  double coef;
  get_disp(disp, irrep, index, coef);

  // transform the gradient into symmetrized coordinates
  gradients_[disp] = displacements(irrep).t() * grad;
  if (debug_) {
    grad.print("cartesian gradient");
    gradients_[disp].print("internal gradient");
    }

  ndisp_++;
}

RefSymmSCMatrix
FinDispMolecularHessian::compute_hessian_from_gradients()
{
  int i;

  RefSymmSCMatrix dhessian;

  RefSymmSCMatrix xhessian = matrixkit()->symmmatrix(d3natom());
  xhessian.assign(0.0);

  // start with the totally symmetric displacments
  int offset = 0;
  if (do_null_displacement_) offset++;
  RefSCMatrix dtrans = displacements(0);
  RefSCDimension ddim = dtrans.coldim();
  dhessian = matrixkit()->symmmatrix(ddim);
  for (i=0; i<ddim.n(); i++) {
    for (int j=0; j<=i; j++) {
      double hij = gradients_[i+offset](j) + gradients_[j+offset](i);
      double ncontrib = 2.0;
      if (do_null_displacement_) {
        hij -= gradients_[0](j) + gradients_[0](i);
        }
      if (eliminate_cubic_terms_) {
        hij -=   gradients_[i+ddim.n()+offset](j)
                 + gradients_[j+ddim.n()+offset](i);
        ncontrib += 2.0;
        if (do_null_displacement_) {
          hij += gradients_[0](j) + gradients_[0](i);
          }
        }
      hij /= ncontrib*disp_;
      dhessian(i,j) = hij;
      }
    }
  do_hess_for_irrep(0, dhessian, xhessian);

  offset += ddim.n();
  if (eliminate_cubic_terms_) offset += ddim.n();
  for (int irrep=1; irrep<nirrep_; irrep++) {
    dtrans = displacements(irrep);
    ddim = dtrans.coldim();
    if (ddim.n() == 0) continue;
    dhessian = matrixkit()->symmmatrix(ddim);
    for (i=0; i<ddim.n(); i++) {
      for (int j=0; j<=i; j++) {
        dhessian(i,j) = (gradients_[i+offset](j)
                         + gradients_[j+offset](i))
                        /(2.0*disp_);
        }
      }
    do_hess_for_irrep(irrep, dhessian, xhessian);
    offset += ddim.n();
    }

  if (debug_) {
    xhessian.print("xhessian");
    }

  return xhessian;
}

void
FinDispMolecularHessian::do_hess_for_irrep(int irrep,
                                        const RefSymmSCMatrix &dhessian,
                                        const RefSymmSCMatrix &xhessian)
{
  RefSCMatrix dtrans = displacements(irrep);
  RefSCDimension ddim = dtrans.coldim();
  if (ddim.n() == 0) return;
  if (debug_) {
    dhessian.print("dhessian");
    dtrans.print("dtrans");
    }
  xhessian.accumulate_transform(dtrans, dhessian);
}

RefSymmSCMatrix
FinDispMolecularHessian::cartesian_hessian()
{
  tim_enter("hessian");

  if (restart_) restart();
  else init();

  ExEnv::out0() << indent
       << "Computing molecular hessian from "
       << ndisplace() << " displacements:" << endl
       << indent << "Starting at displacement: "
       << ndisplacements_done() << endl;
  ExEnv::out0() << indent << "Hessian options: " << endl;
  ExEnv::out0() << indent << "  displacement: " << disp_
               << " bohr" << endl;
  ExEnv::out0() << indent << "  gradient_accuracy: "
               << accuracy_ << " au" << endl;
  ExEnv::out0() << indent << "  eliminate_cubic_terms: "
               << (eliminate_cubic_terms_==0?"no":"yes") << endl;
  ExEnv::out0() << indent << "  only_totally_symmetric: "
               << (only_totally_symmetric_==0?"no":"yes") << endl;

  for (int i=ndisplacements_done(); i<ndisplace(); i++) {
    // This produces side-effects in mol and may even change
    // its symmetry.
    ExEnv::out0() << endl << indent
         << "Beginning displacement " << i << ":" << endl;
    displace(i);

    mole_->obsolete();
    double original_accuracy;
    original_accuracy = mole_->desired_gradient_accuracy();
    if (accuracy_ > 0.0)
      mole_->set_desired_gradient_accuracy(accuracy_);
    else
      mole_->set_desired_gradient_accuracy(disp_/1000.0);
    RefSCVector gradv = mole_->get_cartesian_gradient();
    mole_->set_desired_gradient_accuracy(original_accuracy);
    set_gradient(i, gradv);

    if (checkpoint_) {
      const char *hessckptfile;
      if (MessageGrp::get_default_messagegrp()->me() == 0) {
        hessckptfile = checkpoint_file_;
        }
      else {
        hessckptfile = "/dev/null";
        }
      StateOutBin so(hessckptfile);
      checkpoint_displacements(so);
      }
    }
  original_geometry();
  RefSymmSCMatrix xhessian = compute_hessian_from_gradients();
  tim_exit("hessian");

  symbasis_ = 0;
  delete[] gradients_;
  gradients_ = 0;
  ndisp_ = 0;

  return xhessian;
}

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

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