imcoor.cc

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

      // 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++;
        }

      if (rank != dim_.n()) {
          ExEnv::out0() << indent << "IntMolecularCoor::init: rank changed\n";
          sigma.print("sigma");
        }

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

      ExEnv::out0() << indent
           << scprintf("IntMolecularCoor::init: condition number = %14.8f\n",
                       kappa2);
    }
#endif

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

static int
count_nonzero(const RefSCVector &vec, double eps)
{
  int nz=0, i, n=vec.n();
  for (i=0; i<n; i++) {
      if (fabs(vec(i)) > eps) nz++;
    }
  return nz;
}

static RefSymmSCMatrix
form_partial_K(const Ref<SetIntCoor>& coor, Ref<Molecule>& molecule,
               const RefSCVector& geom,
               double epsilon,
               const RefSCDimension& dnatom3,
               const Ref<SCMatrixKit>& matrixkit,
               RefSCMatrix& projection,
               RefSCVector& totally_symmetric,
               RefSCMatrix& K,int debug)
{
  if (debug) {
      ExEnv::out0() << indent << "form_partial_K:" << endl;
      ExEnv::out0() << incindent;
    }

  // Compute the B matrix for the coordinates
  RefSCDimension dcoor = new SCDimension(coor->n());
  RefSCMatrix B(dcoor, dnatom3,matrixkit);
  coor->bmat(molecule, B);

  if (debug) B.print("B");

  // Project out the previously discovered internal coordinates
  if (projection.nonnull()) {
      B = B * projection;
      if (debug) B.print("Projected B");
    }

  // Compute the singular value decomposition of B
  RefSCMatrix U(dcoor,dcoor,matrixkit);
  RefSCMatrix V(dnatom3,dnatom3,matrixkit);
  RefSCDimension min;
  if (dnatom3.n()<dcoor.n()) min = dnatom3;
  else min = dcoor;
  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++;
    }

  if (debug)
      ExEnv::out0() << indent << "rank(" << epsilon << ",B) = " << rank
           << endl;

  RefSCMatrix SIGMA(dcoor, dnatom3,matrixkit);
  SIGMA.assign(0.0);
  for (i=0; i<nmin; i++) {
      SIGMA(i,i) = sigma(i);
    }

  // return if there are no new coordinates
  if (rank==0) {
      if (debug) ExEnv::out0() << decindent;
      return 0;
    }

  // Find an orthogonal matrix that spans the range of B
  RefSCMatrix Ur;
  RefSCDimension drank = new SCDimension(rank);
  if (rank) {
      Ur = matrixkit->matrix(dcoor,drank);
      Ur.assign_subblock(U,0, dcoor.n()-1, 0, drank.n()-1, 0, 0);
    }

  // Find an orthogonal matrix that spans the null space of B
  int rank_tilde = dnatom3.n() - rank;
  RefSCMatrix Vr_tilde;
  RefSCDimension drank_tilde = new SCDimension(rank_tilde);
  if (rank_tilde) {
      Vr_tilde = matrixkit->matrix(dnatom3,drank_tilde);
      Vr_tilde.assign_subblock(V,0, dnatom3.n()-1, 0, drank_tilde.n()-1,
                               0, drank.n());
    }

  // Find an orthogonal matrix that spans the null(B) perp
  RefSCMatrix Vr;
  if (rank) {
      Vr = matrixkit->matrix(dnatom3,drank);
      Vr.assign_subblock(V,0, dnatom3.n()-1, 0, drank.n()-1, 0, 0);
    }

  // compute the projection into the null space of B
  RefSymmSCMatrix proj_nullspace_B;
  if (rank_tilde) {
      proj_nullspace_B = matrixkit->symmmatrix(dnatom3);
      proj_nullspace_B.assign(0.0);
      proj_nullspace_B.accumulate_symmetric_product(Vr_tilde);
    }

  // compute the projection into the null(B) perp
  RefSymmSCMatrix proj_nullspace_B_perp;
  if (rank) {
      proj_nullspace_B_perp = matrixkit->symmmatrix(dnatom3);
      proj_nullspace_B_perp.assign(0.0);
      proj_nullspace_B_perp.accumulate_symmetric_product(Vr);
    }

  if (Ur.nonnull()) {
      // totally_symmetric will be nonzero for totally symmetric coordinates
      totally_symmetric = Ur.t() * B * geom;

      if (debug) {
          Ur.print("Ur");
          geom.print("geom");
          totally_symmetric.print("totally_symmetric = Ur.t()*B*geom");

          int ntotally_symmetric = count_nonzero(totally_symmetric,0.001);
          ExEnv::out0() << indent << "found " << ntotally_symmetric
               << " totally symmetric coordinates\n";
        }

      // compute the cumulative projection
      if (projection.null()) {
          projection = matrixkit->matrix(dnatom3,dnatom3);
          projection->unit();
        }
      projection = projection * proj_nullspace_B;
    }

  // give Ur to caller
  K = Ur;

  if (debug) ExEnv::out0() << decindent;

  return proj_nullspace_B_perp;
}

// this allocates storage for and computes K and is_totally_symmetric
void
IntMolecularCoor::form_K_matrix(RefSCDimension& dredundant,
                                RefSCDimension& dfixed,
                                RefSCMatrix& K,
                                int*& is_totally_symmetric)
{
  int i,j;

  // The cutoff for whether or not a coordinate is considered totally symmetric
  double ts_eps = 0.0001;

  // The geometry will be needed to check for totally symmetric
  // coordinates
  RefSCVector geom(dnatom3_,matrixkit_);
  for(i=0; i < geom.n()/3; i++) {
      geom(3*i  ) = molecule_->r(i,0);
      geom(3*i+1) = molecule_->r(i,1);
      geom(3*i+2) = molecule_->r(i,2);
    }

  RefSCDimension dcoor = new SCDimension(all_->n());

  // this keeps track of the total projection for the b matrices
  RefSCMatrix projection;
  if (dfixed.n()) {
      ExEnv::out0() << indent
           << "Forming fixed optimization coordinates:" << endl;
      RefSCMatrix Ktmp;
      RefSCVector totally_symmetric_fixed;
      RefSymmSCMatrix null_bfixed_perp
          = form_partial_K(fixed_, molecule_, geom, 0.001, dnatom3_,
                           matrixkit_, projection, totally_symmetric_fixed,
                           Ktmp,debug_);
      // require that the epsilon rank equal the number of fixed coordinates
      if (Ktmp.nrow() != dfixed.n()) {
          ExEnv::err0() << indent
               << scprintf("ERROR: IntMolecularCoor: nfixed = %d rank = %d\n",
                           dfixed.n(), Ktmp.ncol());
          abort();
        }
      // check that fixed coordinates be totally symmetric
      //if (Ktmp.nrow() != count_nonzero(totally_symmetric_fixed, ts_eps)) {
      //    ExEnv::err0() << indent
      //         << scprintf("WARNING: only %d of %d fixed coordinates are"
      //                     " totally symmetric\n",
      //                     count_nonzero(totally_symmetric_fixed, ts_eps),
      //                     dfixed.n());
      //  }
    }

  ExEnv::out0() << indent << "Forming optimization coordinates:" << endl;

  int n_total = 0;

  RefSCVector totally_symmetric_bond;
  RefSCMatrix Kbond;
  if (decouple_bonds_) {
      ExEnv::out0() << indent << "looking for bonds" << endl;
      form_partial_K(bonds_, molecule_, geom, 0.1, dnatom3_, matrixkit_,
                     projection, totally_symmetric_bond, Kbond, debug_);
      if (Kbond.nonnull()) n_total += Kbond.ncol();
    }

  RefSCVector totally_symmetric_bend;
  RefSCMatrix Kbend;
  if (decouple_bends_) {
      ExEnv::out0() << indent << "looking for bends" << endl;
      form_partial_K(bends_, molecule_, geom, 0.1, dnatom3_, matrixkit_,
                     projection, totally_symmetric_bend, Kbend, debug_);
      if (Kbend.nonnull()) n_total += Kbend.ncol();
    }

  if (decouple_bonds_ || decouple_bends_) {
      ExEnv::out0() << indent << "looking for remaining coordinates" << endl;
    }
  RefSCVector totally_symmetric_all;
  RefSCMatrix Kall;
  // I hope the IntCoorSet keeps the ordering
  form_partial_K(all_, molecule_, geom, 0.001, dnatom3_, matrixkit_,
                 projection, totally_symmetric_all, Kall, debug_);
  if (Kall.nonnull()) n_total += Kall.ncol();

  // This requires that all_ coordinates is made up of first bonds,
  // bends, and finally the rest of the coordinates.
  RefSCDimension dtot = new SCDimension(n_total);
  K = matrixkit_->matrix(dcoor, dtot);
  K.assign(0.0);
  int istart=0, jstart=0;
  if (Kbond.nonnull()) {
      if (debug_) Kbond.print("Kbond");
      K.assign_subblock(Kbond, 0, Kbond.nrow()-1, 0, Kbond.ncol()-1, 0, 0);
      istart += Kbond.nrow();
      jstart += Kbond.ncol();
    }
  if (Kbend.nonnull()) {
      if (debug_) Kbend.print("Kbend");
      K.assign_subblock(Kbend, istart, istart+Kbend.nrow()-1,
                        jstart, jstart+Kbend.ncol()-1, 0, 0);
      istart += Kbend.nrow();
      jstart += Kbend.ncol();
    }
  if (Kall.nonnull()) {
      if (debug_) Kall.print("Kall");
      K.assign_subblock(Kall, 0, Kall.nrow()-1,
                        jstart, jstart+Kall.ncol()-1, 0, 0);
    }
  if (debug_) K.print("K");

  is_totally_symmetric = new int[K.ncol()];
  j=0;
  if (Kbond.nonnull()) {
      for (i=0; i<Kbond.ncol(); i++,j++) {
          if (fabs(totally_symmetric_bond(i)) > ts_eps)
              is_totally_symmetric[j] = 1;
          else is_totally_symmetric[j] = 0;
        }
    }
  if (Kbend.nonnull()) {
      for (i=0; i<Kbend.ncol(); i++,j++) {
          if (fabs(totally_symmetric_bend(i)) > ts_eps)
              is_totally_symmetric[j] = 1;
          else is_totally_symmetric[j] = 0;
        }
    }
  if (Kall.nonnull()) {
      for (i=0; i<Kall.ncol(); i++,j++) {
          if (fabs(totally_symmetric_all(i)) > ts_eps)
              is_totally_symmetric[j] = 1;
          else is_totally_symmetric[j] = 0;
        }
    }
}

IntMolecularCoor::~IntMolecularCoor()
{
}

void
IntMolecularCoor::save_data_state(StateOut&s)
{
  MolecularCoor::save_data_state(s);

  SavableState::save_state(generator_.pointer(),s);

  s.put(decouple_bonds_);
  s.put(decouple_bends_);

  s.put(max_update_steps_);
  s.put(max_update_disp_);
  s.put(given_fixed_values_);

  s.put(form_print_simples_);
  s.put(form_print_variable_);
  s.put(form_print_constant_);
  s.put(form_print_molecule_);

  SavableState::save_state(dim_.pointer(),s);
  SavableState::save_state(dvc_.pointer(),s);

  SavableState::save_state(all_.pointer(),s);
  
  SavableState::save_state(variable_.pointer(),s);
  SavableState::save_state(constant_.pointer(),s);

  SavableState::save_state(fixed_.pointer(),s);
  SavableState::save_state(followed_.pointer(),s);
  SavableState::save_state(watched_.pointer(),s);

  SavableState::save_state(bonds_.pointer(),s);
  SavableState::save_state(bends_.pointer(),s);
  SavableState::save_state(tors_.pointer(),s);
  SavableState::save_state(outs_.pointer(),s);
  SavableState::save_state(extras_.pointer(),s);

  s.put(update_bmat_);
  s.put(only_totally_symmetric_);
  s.put(scale_bonds_);
  s.put(scale_bends_);
  s.put(scale_tors_);
  s.put(scale_outs_);
  s.put(simple_tolerance_);
  s.put(symmetry_tolerance_);
  s.put(coordinate_tolerance_);
  s.put(cartesian_tolerance_);
}

RefSCDimension
IntMolecularCoor::dim()
{
  return dim_;
}

int
IntMolecularCoor::all_to_cartesian(const Ref<Molecule> &mol,
                                   RefSCVector&new_internal)
{
  // get a reference to Molecule for convenience
  Molecule& molecule = *(mol.pointer());

  // don't bother updating the bmatrix when the error is less than this
  const double update_tolerance = 1.0e-6;