symmcoor.cc

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

{
  IntMolecularCoor::save_data_state(s);
  s.put(change_coordinates_);
  s.put(transform_hessian_);
  s.put(max_kappa2_);
}

void
SymmMolecularCoor::init()
{
  IntMolecularCoor::init();
  change_coordinates_ = 0;
  max_kappa2_ = 10.0;
  transform_hessian_ = 1;
}

void
SymmMolecularCoor::form_coordinates(int keep_variable)
{
  int i;
  int nbonds = bonds_->n();
  int nbends = bends_->n();
  int ntors = tors_->n();
  int nouts = outs_->n();
  int nextras = extras_->n();

  Ref<SetIntCoor> saved_fixed_ = fixed_;
  fixed_ = new SetIntCoor;
  fixed_->add(saved_fixed_);
  // if we're following coordinates, add them to the fixed list
  if (followed_.nonnull())
    fixed_->add(followed_);
  
  int nredundant = nbonds + nbends + ntors + nouts + nextras;
  int nfixed = fixed_->n();

  // see how many coords we expect
  int n3 = molecule_->natom()*3;
  int nunique = n3 - 6; // need to detect linear

  if (nredundant < nunique) {
      ExEnv::err0() << indent
           << "SymmMolecularCoor::form_coordinates: "
           << "found too few redundant coordinates\n"
           << indent << scprintf("nredundant = %d, 3n-6 = %d\n",
                                 nredundant, nunique)
           << indent << "  (the geometry is probably bad)\n";
      molecule_->print(ExEnv::err0());
      abort();
    }

  RefSCDimension dredundant = new SCDimension(nredundant, "Nredund");
  RefSCDimension dfixed = new SCDimension(nfixed, "Nfixed");
  RefSCMatrix K; // nredundant x nnonzero
  int* is_totally_symmetric; // nnonzero; if 1 coor has tot. symm. component

  form_K_matrix(dredundant, dfixed, K, is_totally_symmetric);

  RefSCDimension dnonzero = K.coldim();
  int nnonzero = dnonzero.n();

  if (!keep_variable) variable_->clear();
  constant_->clear();

  // now remove followed coords from the fixed list, and add to the
  // variable list
  if (followed_.nonnull()) {
    fixed_->pop();
    variable_->add(followed_);
  }
  
  // put the fixed coordinates into the constant list
  nfixed = fixed_->n();
  for (i=0; i<nfixed; i++) {
      constant_->add(fixed_->coor(i));
    }

  // ok, now we have the K matrix, the columns of which give us the
  // contribution from each red. coord to the ith non-red. coord.
  // this gets a little hairy since the red coords can themselves be
  // linear combinations of simple coords
  for(i=0; i < nnonzero; i++) {
      // construct the new linear combination coordinate
      char label[80];
      if (is_totally_symmetric[i]) {
          sprintf(label,"symm_coord_%03d",i+1);
        }
      else {
          sprintf(label,"asymm_coord_%03d",i+1);
        }
      SumIntCoor* coordinate = new SumIntCoor(label);

      int j;
      for(j=0; j < nredundant; j++) {
          if(pow(K(j,i),2.0) > simple_tolerance_) {
              Ref<IntCoor> c = all_->coor(j);
              coordinate->add(c,K(j,i));
              if (debug_) {
                  ExEnv::out0() << "added redund coor "
                       << j << " to coor " << i << ":" << endl;
                  c->print();
                }
            }
        }

      // normalize the coordinate
      coordinate->normalize();

      if (only_totally_symmetric_ && !is_totally_symmetric[i]) {
          // Don't put nonsymmetric coordinates into the
          // constant_ coordinate set.  This causes problems
          // when coordinates with small coefficients are eliminated
          // since they can then acquire symmetric components.
          // constant_->add(coordinate);
          delete coordinate;
        }
      else {
          if (!keep_variable) variable_->add(coordinate);
        }
    }

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

  ExEnv::out0() << incindent << indent
       << "SymmMolecularCoor::form_variable_coordinates()\n" << incindent
       << indent << "expected " << nunique << " coordinates\n"
       << indent << "found " << variable_->n() << " variable coordinates\n"
       << indent << "found " << constant_->n() << " constant coordinates\n"
       << decindent << decindent << flush;

  delete[] is_totally_symmetric;
  fixed_ = saved_fixed_;

  if (form_print_molecule_) molecule_->print();
  if (form_print_simples_) print_simples(ExEnv::out0());
  if (form_print_variable_) print_variable(ExEnv::out0());
  if (form_print_constant_) print_constant(ExEnv::out0());
}

void
SymmMolecularCoor::guess_hessian(RefSymmSCMatrix&hessian)
{
  // first form diagonal hessian in redundant internal coordinates
  RefSCDimension rdim = new SCDimension(all_->n(), "Nall");
  RefSymmSCMatrix rhessian(rdim,matrixkit_);
  rhessian.assign(0.0);
  all_->guess_hessian(molecule_,rhessian);

  // create redundant coordinate bmat
  RefSCDimension dn3 = dnatom3_;
  RefSCMatrix bmatr(rdim,dn3,matrixkit_);
  all_->bmat(molecule_,bmatr);

  // then form the variable coordinate bmat
  RefSCDimension dredundant = new SCDimension(variable_->n(), "Nvar");
  RefSCMatrix bmat(dredundant,dn3,matrixkit_);
  variable_->bmat(molecule_,bmat);

  // and (B*B+)^-1
  RefSymmSCMatrix bmbt(dredundant,matrixkit_);
  bmbt.assign(0.0);
  bmbt.accumulate_symmetric_product(bmat);
  bmbt = bmbt.gi();

  // now transform redundant hessian to internal coordinates
  // Hc = Br+ * Hr * Br
  // Hi = (B*B+)^-1 * B * Hc * B+ * (B*B+)^-1+
  //    = bmbt_inv*B*Br+ * Hr * Br*B+*bmbt_inv+
  //    = b * Hr * b+  (b = (B*B+)^-1 * B * Br+)
  RefSCMatrix b = bmbt * bmat * bmatr.t();
  
  hessian.assign(0.0);
  hessian.accumulate_transform(b,rhessian);
}

RefSymmSCMatrix
SymmMolecularCoor::inverse_hessian(RefSymmSCMatrix& hessian)
{
  return hessian.gi();
}

// Possibly change to a new coordinate system
Ref<NonlinearTransform>
SymmMolecularCoor::change_coordinates()
{
  if (dim_.n() == 0 || !change_coordinates_) return 0;

  const double epsilon = 0.001;

  // compute the condition number of the old coordinate system at the
  // current point
  RefSCMatrix B(dim_, dnatom3_, matrixkit_);
  variable_->bmat(molecule_, B);

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

  // Compute the epsilon rank of B
  int i, rank = 0;
  for (i=0; i<nmin; i++) {
      if (sigma(i) > epsilon) rank++;
    }
  // the rank could get bigger if there is a fixed coordinate
  if (rank < dim_.n() || ((fixed_.null()
                           || fixed_->n() == 0) && rank != dim_.n())) {
      ExEnv::err0() << indent
           << "SymmMolecularCoor::change_coordinates: "
           << "disallowed rank change\n";
      abort();
    }
  if (rank != dim_.n()) {
      ExEnv::out0() << indent
           << "SymmMolecularCoor::change_coordinates: rank changed\n";
    }

  double kappa2 = sigma(0)/sigma(dim_.n()-1);

  ExEnv::out0() << indent
       << scprintf(
           "SymmMolecularCoor: condition number = %14.8f (max = %14.8f)\n",
           kappa2, max_kappa2_);

  if (kappa2 > max_kappa2_) {
      Ref<SetIntCoor> oldvariable = new SetIntCoor;
      oldvariable->add(variable_);

      // form the new variable coordinates
      form_coordinates();

      SymmCoorTransform *trans = new SymmCoorTransform(molecule_,
                                                       dnatom3_,
                                                       matrixkit_,
                                                       oldvariable,
                                                       variable_,
                                                       transform_hessian_);
      return trans;
    }

  return 0;
}

void
SymmMolecularCoor::print(ostream& os) const
{
  IntMolecularCoor::print(os);
  
  os << indent << "SymmMolecularCoor Parameters:\n" << incindent
     << indent << "change_coordinates = "
     << (change_coordinates_?"yes":"no") << endl
     << indent << "transform_hessian = "
     << (transform_hessian_?"yes":"no") << endl
     << indent << scprintf("max_kappa2 = %f",max_kappa2_) << endl
     << decindent << endl;
}

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

}

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