imcoor.cc

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

  // compute the internal coordinate displacements
  RefSCVector old_internal(dvc_,matrixkit_);

  RefSCMatrix internal_to_cart_disp;
  double maxabs_cart_diff = 0.0;
  for (int step = 0; step < max_update_steps_; step++) {
      // compute the old internal coordinates
      all_to_internal(mol, old_internal);

      if (debug_) {
          ExEnv::out0()
               << indent << "Coordinates on step " << step << ":" << endl;
          variable_->print_details(0,ExEnv::out0());
        }

      // the displacements
      RefSCVector displacement = new_internal - old_internal;
      if (debug_ && step == 0) {
          displacement.print("Step 0 Internal Coordinate Displacement");
        }

      if ((update_bmat_ && maxabs_cart_diff>update_tolerance)
          || internal_to_cart_disp.null()) {
          if (debug_) {
              ExEnv::out0() << indent << "updating bmatrix" << endl;
            }

          int i;
          RefSCMatrix bmat(dvc_,dnatom3_,matrixkit_);

          // form the set of all coordinates
          Ref<SetIntCoor> variable_and_constant = new SetIntCoor();
          variable_and_constant->add(variable_);
          variable_and_constant->add(constant_);

          // form the bmatrix
          variable_and_constant->bmat(mol,bmat);

          // Compute the singular value decomposition of B
          RefSCMatrix U(dvc_,dvc_,matrixkit_);
          RefSCMatrix V(dnatom3_,dnatom3_,matrixkit_);
          RefSCDimension min;
          if (dnatom3_.n()<dvc_.n()) min = dnatom3_;
          else min = dvc_;
          int nmin = min.n();
          RefDiagSCMatrix sigma(min,matrixkit_);
          bmat.svd(U,sigma,V);

          // compute the epsilon rank of B
          int rank = 0;
          for (i=0; i<nmin; i++) {
              if (fabs(sigma(i)) > 0.0001) rank++;
            }

          RefSCDimension drank = new SCDimension(rank);
          RefDiagSCMatrix sigma_i(drank,matrixkit_);
          for (i=0; i<rank; i++) {
              sigma_i(i) = 1.0/sigma(i);
            }
          RefSCMatrix Ur(dvc_, drank, matrixkit_);
          RefSCMatrix Vr(dnatom3_, drank, matrixkit_);
          Ur.assign_subblock(U,0, dvc_.n()-1, 0, drank.n()-1, 0, 0);
          Vr.assign_subblock(V,0, dnatom3_.n()-1, 0, drank.n()-1, 0, 0);
          internal_to_cart_disp = Vr * sigma_i * Ur.t();

        }

      // compute the cartesian displacements
      RefSCVector cartesian_displacement = internal_to_cart_disp*displacement;
      if (debug_ && step == 0) {
          internal_to_cart_disp.print("Internal to Cartesian Transform");
          cartesian_displacement.print("Step 0 Cartesian Displacment");
        }
      // update the geometry
      for (int i=0; i < dnatom3_.n(); i++) {
#if OLD_BMAT
          molecule.r(i/3,i%3) += cartesian_displacement(i) * 1.88972666;
#else        
          molecule.r(i/3,i%3) += cartesian_displacement(i);
#endif          
        }

      // fix symmetry breaking due to numerical round-off
      molecule.cleanup_molecule();

      // check for convergence
      Ref<SCElementMaxAbs> maxabs = new SCElementMaxAbs();
      Ref<SCElementOp> op = maxabs.pointer();
      cartesian_displacement.element_op(op);
      maxabs_cart_diff = maxabs->result();
      if (maxabs_cart_diff < cartesian_tolerance_) {

          constant_->update_values(mol);
          variable_->update_values(mol);

          return 0;
        }
    }

  ExEnv::err0() << indent
       << "WARNING: IntMolecularCoor::all_to_cartesian(RefSCVector&):"
       << " too many iterations in geometry update" << endl;

  new_internal.print("desired internal coordinates");
  (new_internal
   - old_internal).print("difference of desired and actual coordinates");

  return -1;
}

int
IntMolecularCoor::to_cartesian(const Ref<Molecule> &mol,
                               const RefSCVector&new_internal)
{
  if (new_internal.dim().n() != dim_.n()
      || dvc_.n() != variable_->n() + constant_->n()
      || new_internal.dim().n() != variable_->n()) {
      ExEnv::err0() << indent << "to_cartesian: internal error in dim\n";
      abort();
    }

  RefSCVector all_internal(dvc_,matrixkit_);

  int i,j;

  for (i=0; i<variable_->n(); i++) all_internal(i) = new_internal(i);
  for (j=0; j<constant_->n(); i++,j++) {
      all_internal(i) = constant_->coor(j)->value();
    }

  int ret = all_to_cartesian(mol, all_internal);

  if (watched_.nonnull()) {
      ExEnv::out0() << endl
           << indent << "Watched coordinate(s):\n" << incindent;
      watched_->update_values(mol);
      watched_->print_details(mol,ExEnv::out0());
      ExEnv::out0() << decindent;
    }
  
  return ret;
}

int
IntMolecularCoor::all_to_internal(const Ref<Molecule> &mol,RefSCVector&internal)
{
  if (internal.dim().n() != dvc_.n()
      || dim_.n() != variable_->n()
      || dvc_.n() != variable_->n() + constant_->n()) {
      ExEnv::err0() << indent << "all_to_internal: internal error in dim\n";
      abort();
    }

  variable_->update_values(mol);
  constant_->update_values(mol);
   
  int n = dim_.n();
  int i;
  for (i=0; i<n; i++) {
      internal(i) = variable_->coor(i)->value();
    }
  n = dvc_.n();
  for (int j=0; i<n; i++,j++) {
      internal(i) = constant_->coor(j)->value();
    }

  return 0;
}

int
IntMolecularCoor::to_internal(RefSCVector&internal)
{
  if (internal.dim().n() != dim_.n()
      || dim_.n() != variable_->n()) {
      ExEnv::err0() << indent << "to_internal: internal error in dim\n";
      abort();
    }

  variable_->update_values(molecule_);
   
  int n = dim_.n();
  for (int i=0; i<n; i++) {
      internal(i) = variable_->coor(i)->value();
    }

  return 0;
}

int
IntMolecularCoor::to_cartesian(RefSCVector&gradient,RefSCVector&internal)
{
  RefSCMatrix bmat(dim_,gradient.dim(),matrixkit_);
  variable_->bmat(molecule_,bmat);

  gradient = bmat.t() * internal;
  
  return 0;
}

// converts the gradient in cartesian coordinates to internal coordinates
int
IntMolecularCoor::to_internal(RefSCVector&internal,RefSCVector&gradient)
{
  RefSCMatrix bmat(dvc_,gradient.dim(),matrixkit_);
  RefSymmSCMatrix bmbt(dvc_,matrixkit_);

  Ref<SetIntCoor> variable_and_constant = new SetIntCoor();
  variable_and_constant->add(variable_);
  variable_and_constant->add(constant_);

  // form the bmatrix
  variable_and_constant->bmat(molecule_,bmat);
  
  // form the inverse of bmatrix * bmatrix_t
  bmbt.assign(0.0);
  bmbt.accumulate_symmetric_product(bmat);
  bmbt = bmbt.gi();

#if OLD_BMAT  
  RefSCVector all_internal = bmbt*bmat*(gradient*8.2388575);
#else
  RefSCVector all_internal = bmbt*bmat*gradient;
#endif  

  // put the variable coordinate gradients into internal
  int n = variable_->n();
  for (int i=0; i<n; i++) {
      internal.set_element(i,all_internal.get_element(i));
    }

  return 0;
}

int
IntMolecularCoor::to_cartesian(RefSymmSCMatrix&cart,RefSymmSCMatrix&internal)
{
  cart.assign(0.0);
  RefSCMatrix bmat(dim_,cart.dim(),matrixkit_);
  variable_->bmat(molecule_,bmat);
  cart.accumulate_transform(bmat.t(),internal);
  return 0;
}

int
IntMolecularCoor::to_internal(RefSymmSCMatrix&internal,RefSymmSCMatrix&cart)
{
  // form bmat
  RefSCMatrix bmat(dim_,cart.dim(),matrixkit_);
  variable_->bmat(molecule_,bmat);
  // and (B*B+)^-1
  RefSymmSCMatrix bmbt(dim_,matrixkit_);
  bmbt.assign(0.0);
  bmbt.accumulate_symmetric_product(bmat);
  bmbt = bmbt.gi();

  internal.assign(0.0);
  internal.accumulate_transform(bmbt*bmat,cart);
  return 0;
}

int
IntMolecularCoor::nconstrained()
{
  return fixed_->n();
}

void
IntMolecularCoor::print(ostream& os) const
{
  all_->update_values(molecule_);

  os << indent << "IntMolecularCoor Parameters:\n" << incindent
     << indent << "update_bmat = " << (update_bmat_?"yes":"no") << endl
     << indent << "scale_bonds = " << scale_bonds_ << endl
     << indent << "scale_bends = " << scale_bends_ << endl
     << indent << "scale_tors = " << scale_tors_ << endl
     << indent << "scale_outs = " << scale_outs_ << endl
     << indent << scprintf("symmetry_tolerance = %e\n",symmetry_tolerance_)
     << indent << scprintf("simple_tolerance = %e\n",simple_tolerance_)
     << indent << scprintf("coordinate_tolerance = %e\n",coordinate_tolerance_)
     << indent << "have_fixed_values = " << given_fixed_values_ << endl
     << indent << "max_update_steps = " << max_update_steps_ << endl
     << indent << scprintf("max_update_disp = %f\n",max_update_disp_)
     << indent << "have_fixed_values = " << given_fixed_values_ << endl
     << decindent << endl;

  molecule_->print(os);
  os << endl;

  print_simples(os);
  os << endl;

  if (form_print_variable_) {
      print_variable(os);
      os << endl;
    }

  if (form_print_constant_) {
      print_constant(os);
      os << endl;
    }
}

void
IntMolecularCoor::print_simples(ostream& os) const
{
  if (matrixkit()->messagegrp()->me()==0) {
    if (bonds_->n()) {
      os << indent << "Bonds:\n" << incindent;
      bonds_->print_details(molecule_,os);
      os << decindent;
    }
    if (bends_->n()) {
      os << indent << "Bends:\n" << incindent;
      bends_->print_details(molecule_,os);
      os << decindent;
    }
    if (tors_->n()) {
      os << indent << "Torsions:\n" << incindent;
      tors_->print_details(molecule_,os);
      os << decindent;
    }
    if (outs_->n()) {
      os << indent << "Out of Plane:\n" << incindent;
      outs_->print_details(molecule_,os);
      os << decindent;
    }
    if (extras_->n()) {
      os << indent << "Extras:\n" << incindent;
      extras_->print_details(molecule_,os);
      os << decindent;
    }
    if (fixed_->n()) {
      os << indent << "Fixed:\n" << incindent;
      fixed_->print_details(molecule_,os);
      os << decindent;
    }
    if (followed_.nonnull()) {
      os << indent << "Followed:\n" << incindent;
      followed_->print_details(molecule_,os);
      os << decindent;
    }
    if (watched_.nonnull()) {
      os << indent << "Watched:\n" << incindent;
      watched_->print_details(molecule_,os);
      os << decindent;
    }
  }
}

void
IntMolecularCoor::print_variable(ostream& os) const
{
  if (variable_->n() == 0) return;
  os << indent
     << "Variable Coordinates:" << endl;
  os << incindent;
  variable_->print_details(molecule_,os);
  os << decindent;
}

void
IntMolecularCoor::print_constant(ostream& os) const
{
  if (constant_->n() == 0) return;
  os << indent
     << "Constant Coordinates:" << endl;
  os << incindent;
  constant_->print_details(molecule_,os);
  os << decindent;
}

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

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