coor.cc

来自「大型并行量子化学软件；支持密度泛函（DFT）。可以进行各种量子化学计算。支持CH」· CC 代码 · 共 1,155 行 · 第 1/2 页

  molecule_ << keyval->describedclassvalue("molecule");

  if (molecule_.null()) {
      ExEnv::err0() << indent
           << "MolecularCoor(const Ref<KeyVal>&keyval): molecule not found\n";
      abort();
    }

  debug_ = keyval->intvalue("debug");

  matrixkit_ << keyval->describedclassvalue("matrixkit");
  dnatom3_ << keyval->describedclassvalue("natom3");

  if (matrixkit_.null()) matrixkit_ = SCMatrixKit::default_matrixkit();

  if (dnatom3_.null()) dnatom3_ = new SCDimension(3*molecule_->natom());
  else if (dnatom3_->n() != 3 * molecule_->natom()) {
      ExEnv::err0() << indent << "MolecularCoor(KeyVal): bad dnatom3 value\n";
      abort();
    }
}

MolecularCoor::MolecularCoor(StateIn&s):
  SavableState(s)
{
  debug_ = 0;
  matrixkit_ = SCMatrixKit::default_matrixkit();
  molecule_ << SavableState::restore_state(s);
  dnatom3_ << SavableState::restore_state(s);
}

MolecularCoor::~MolecularCoor()
{
}

void
MolecularCoor::save_data_state(StateOut&s)
{
  SavableState::save_state(molecule_.pointer(),s);
  SavableState::save_state(dnatom3_.pointer(),s);
}

int
MolecularCoor::nconstrained()
{
  return 0;
}

// The default action is to never change the coordinates.
Ref<NonlinearTransform>
MolecularCoor::change_coordinates()
{
  return 0;
}

int
MolecularCoor::to_cartesian(const RefSCVector&internal)
{
  return to_cartesian(molecule_, internal);
}

///////////////////////////////////////////////////////////////////////////
// members of IntCoorGen

static ClassDesc IntCoorGen_cd(
  typeid(IntCoorGen),"IntCoorGen",2,"public SavableState",
  0, create<IntCoorGen>, create<IntCoorGen>);

IntCoorGen::IntCoorGen(const Ref<Molecule>& mol,
                       int nextra_bonds, int *extra_bonds)
{
  init_constants();

  molecule_ = mol;
  nextra_bonds_ = nextra_bonds;
  extra_bonds_ = extra_bonds;
}

IntCoorGen::IntCoorGen(const Ref<KeyVal>& keyval)
{
  init_constants();

  molecule_ << keyval->describedclassvalue("molecule");

  radius_scale_factor_
    = keyval->doublevalue("radius_scale_factor",
                          KeyValValuedouble(radius_scale_factor_));

  // degrees
  linear_bend_thres_
    = keyval->doublevalue("linear_bend_threshold",
                          KeyValValuedouble(linear_bend_thres_));

  // entered in degrees; stored as cos(theta)
  linear_tors_thres_
    = keyval->doublevalue("linear_tors_threshold",
                          KeyValValuedouble(linear_tors_thres_));

  linear_bends_
    = keyval->booleanvalue("linear_bend",
                           KeyValValueboolean(linear_bends_));

  linear_lbends_
    = keyval->booleanvalue("linear_lbend",
                           KeyValValueboolean(linear_lbends_));

  linear_tors_
    = keyval->booleanvalue("linear_tors",
                           KeyValValueboolean(linear_tors_));

  linear_stors_
    = keyval->booleanvalue("linear_stors",
                           KeyValValueboolean(linear_stors_));

  // the extra_bonds list is given as a vector of atom numbers
  // (atom numbering starts at 1)
  nextra_bonds_ = keyval->count("extra_bonds");
  nextra_bonds_ /= 2;
  if (nextra_bonds_) {
      extra_bonds_ = new int[nextra_bonds_*2];
      for (int i=0; i<nextra_bonds_*2; i++) {
          extra_bonds_[i] = keyval->intvalue("extra_bonds",i);
          if (keyval->error() != KeyVal::OK) {
              ExEnv::err0() << indent
                   << "IntCoorGen:: keyval CTOR: problem reading "
                   << "\"extra_bonds:" << i << "\"\n";
              abort();
            }
        }
    }
  else {
      extra_bonds_ = 0;
    }
}

IntCoorGen::IntCoorGen(StateIn& s):
  SavableState(s)
{
  molecule_ << SavableState::restore_state(s);
  s.get(linear_bends_);
  if (s.version(::class_desc<IntCoorGen>()) >= 2) {
      s.get(linear_lbends_);
    }
  s.get(linear_tors_);
  s.get(linear_stors_);
  s.get(linear_bend_thres_);
  s.get(linear_tors_thres_);
  s.get(nextra_bonds_);
  s.get(extra_bonds_);
  s.get(radius_scale_factor_);
}

void
IntCoorGen::init_constants()
{
  nextra_bonds_ = 0;
  extra_bonds_ = 0;
  radius_scale_factor_ = 1.1;
  linear_bend_thres_ = 1.0;
  linear_tors_thres_ = 1.0;
  linear_bends_ = 0;
  linear_lbends_ = 1;
  linear_tors_ = 0;
  linear_stors_ = 1;
}

IntCoorGen::~IntCoorGen()
{
  if (extra_bonds_) delete[] extra_bonds_;
}

void
IntCoorGen::save_data_state(StateOut& s)
{
  SavableState::save_state(molecule_.pointer(),s);
  s.put(linear_bends_);
  s.put(linear_lbends_);
  s.put(linear_tors_);
  s.put(linear_stors_);
  s.put(linear_bend_thres_);
  s.put(linear_tors_thres_);
  s.put(nextra_bonds_);
  s.put(extra_bonds_,2*nextra_bonds_);
  s.put(radius_scale_factor_);
}

void
IntCoorGen::print(ostream& out) const
{
  out << indent << "IntCoorGen:" << endl << incindent
      << indent << "linear_bends = " << linear_bends_ << endl
      << indent << "linear_lbends = " << linear_lbends_ << endl
      << indent << "linear_tors = " << linear_tors_ << endl
      << indent << "linear_stors = " << linear_stors_ << endl
      << indent << scprintf("linear_bend_threshold = %f\n",linear_bend_thres_)
      << indent << scprintf("linear_tors_threshold = %f\n",linear_tors_thres_)
      << indent << scprintf("radius_scale_factor = %f\n",radius_scale_factor_)
      << indent << "nextra_bonds = " << nextra_bonds_ << endl
      << decindent;
}

static void
find_bonds(Molecule &m, BitArrayLTri &bonds,
           double radius_scale_factor_)
{
  int i, j;
  for(i=0; i < m.natom(); i++) {
    double at_rad_i = m.atominfo()->atomic_radius(m.Z(i));
    SCVector3 ri(m.r(i));

    for(j=0; j < i; j++) {
      double at_rad_j = m.atominfo()->atomic_radius(m.Z(j));
      SCVector3 rj(m.r(j));

      if (ri.dist(rj)
          < radius_scale_factor_*(at_rad_i+at_rad_j))
        bonds.set(i,j);
      }
    }

  // check for groups of atoms bound to nothing
  std::set<int> boundatoms;
  std::set<int> newatoms, nextnewatoms;
  // start out with atom 0
  newatoms.insert(0);
  std::set<int>::iterator iatom;
  int warning_printed = 0;
  while (newatoms.size() > 0) {
    while (newatoms.size() > 0) {
      // newatoms gets merged into boundatoms
      for (iatom=newatoms.begin(); iatom!=newatoms.end(); iatom++) {
        boundatoms.insert(*iatom);
        }
      // set nextnewatoms to atoms bound to boundatoms that are not already
      // in boundatoms
      nextnewatoms.clear();
      for (iatom=newatoms.begin(); iatom!=newatoms.end(); iatom++) {
        int atom = *iatom;
        for (i=0; i<m.natom(); i++) {
          if (bonds(i,atom) && boundatoms.find(i) == boundatoms.end()) {
            nextnewatoms.insert(i);
            }
          }
        }
      // set newatoms to nextnewatoms to start off the next iteration
      newatoms.clear();
      for (iatom=nextnewatoms.begin(); iatom!=nextnewatoms.end(); iatom++) {
        newatoms.insert(*iatom);
        }
      }

    if (boundatoms.size() != m.natom()) {
      if (!warning_printed) {
        warning_printed = 1;
        ExEnv::out0()
             << indent << "WARNING: two unbound groups of atoms" << endl
             << indent << "         consider using extra_bonds input" << endl
             << endl;
        }
      // find an unbound group
      double nearest_dist;
      int nearest_bound = -1, nearest_unbound = -1;
      for(i=0; i < m.natom(); i++) {
        if (boundatoms.find(i) == boundatoms.end()) {
          SCVector3 ri(m.r(i));
          for (iatom=boundatoms.begin(); iatom!=boundatoms.end(); iatom++) {
            SCVector3 rj(m.r(*iatom));
            double d = ri.dist(rj);
            if (nearest_unbound == -1 || d < nearest_dist) {
              nearest_dist = d;
              nearest_bound = *iatom;
              nearest_unbound = i;
              }
            }
          }
        }
      if (nearest_bound == -1) {
        ExEnv::out0() << indent
             << "ERROR: impossible error generating coordinates"
             << endl;
        abort();
        }
      // add all bound atoms within a certain distance of nearest_unbound
      // --- should really do this for all atoms that nearest_unbound
      // --- is already bound to, too
      SCVector3 rnearest_unbound(m.r(nearest_unbound));
      for (iatom=boundatoms.begin(); iatom!=boundatoms.end(); iatom++) {
        SCVector3 r(m.r(*iatom));
        if (*iatom == nearest_bound
            || rnearest_unbound.dist(r) < 1.1 * nearest_dist) {
          ExEnv::out0() << indent
               << "         adding bond between "
               << *iatom+1 << " and " << nearest_unbound+1 << endl;
          bonds.set(*iatom,nearest_unbound);
          }
        }
      newatoms.insert(nearest_unbound);
      }
    }
}

void
IntCoorGen::generate(const Ref<SetIntCoor>& sic)
{
  int i;
  Molecule& m = *molecule_.pointer();

  // let's go through the geometry and find all the close contacts
  // bonds is a lower triangle matrix of 1's and 0's indicating whether
  // there is a bond between atoms i and j

  BitArrayLTri bonds(m.natom(),m.natom());

  for (i=0; i<nextra_bonds_; i++) {
    bonds.set(extra_bonds_[i*2]-1,extra_bonds_[i*2+1]-1);
  }

  find_bonds(m, bonds, radius_scale_factor_);
      
  // compute the simple internal coordinates by type
  add_bonds(sic,bonds,m);
  add_bends(sic,bonds,m);
  add_tors(sic,bonds,m);
  add_out(sic,bonds,m);

  ExEnv::out0() << endl << indent
       << "IntCoorGen: generated " << sic->n() << " coordinates." << endl;
}

///////////////////////////////////////////////////////////////////////////
// auxillary functions of IntCoorGen

/*
 * the following are translations of functions written by Gregory Humphreys
 * at the NIH
 */

/*
 * for each bonded pair, add an entry to the simple coord list
 */

void
IntCoorGen::add_bonds(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
{
  int i,j,ij;
  int labelc=0;
  char label[80];

  for(i=ij=0; i < m.natom(); i++) {
    for(j=0; j <= i; j++,ij++) {
      if(bonds[ij]) {
        labelc++;
        sprintf(label,"s%d",labelc);
        list->add(new Stre(label,j+1,i+1));
        }
      }
    }
  }

/*
 * return 1 if all three atoms are nearly on the same line.
 */

// returns fabs(cos(theta_ijk))
double
IntCoorGen::cos_ijk(Molecule& m, int i, int j, int k)
{
  SCVector3 a, b, c;
  int xyz;
  for (xyz=0; xyz<3; xyz++) {
      a[xyz] = m.r(i,xyz);
      b[xyz] = m.r(j,xyz);
      c[xyz] = m.r(k,xyz);
    }
  SCVector3 ab = a - b;
  SCVector3 cb = c - b;
  return fabs(ab.dot(cb)/(ab.norm()*cb.norm()));
}

void
IntCoorGen::add_bends(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
{
  int i,j,k;
  int labelc=0;
  char label[80];

  int n = m.natom();

  double thres = cos(linear_bend_thres_*M_PI/180.0);

  for(i=0; i < n; i++) {
    SCVector3 ri(m.r(i));
    for(j=0; j < n; j++) {
      if(bonds(i,j)) {
        SCVector3 rj(m.r(j));
        for(k=0; k < i; k++) {
          if(bonds(j,k)) {
            SCVector3 rk(m.r(k));
            int is_linear = (cos_ijk(m,i,j,k) >= thres);
            if (linear_bends_ || !is_linear) {
              labelc++;
              sprintf(label,"b%d",labelc);
              list->add(new Bend(label,k+1,j+1,i+1));
              }
            if (linear_lbends_ && is_linear) {
              // find a unit vector roughly perp to the bonds
              SCVector3 u;
              // first try to find another atom, that'll help keep one of
              // the coordinates totally symmetric in some cases
              int most_perp_atom = -1;
              double cos_most_perp = thres;
              for (int l=0; l < n; l++) {
                if (l == i || l == j || l == k) continue;
                double tmp = cos_ijk(m,i,j,l);
                if (tmp < cos_most_perp) {
                  cos_most_perp = tmp;
                  most_perp_atom = l;
                  }
                }
              if (most_perp_atom != -1) {
                SCVector3 rmpa(m.r(most_perp_atom));
                u = rj-rmpa;
                u.normalize();
                }
              else {
                SCVector3 b1, b2;
                b1 = ri-rj;
                b2 = rk-rj;
                u = b1.perp_unit(b2);
                }
              labelc++;
              sprintf(label,"b%d",labelc);
              list->add(new LinIP(label,k+1,j+1,i+1,u));
              labelc++;
              sprintf(label,"b%d",labelc);
              list->add(new LinOP(label,k+1,j+1,i+1,u));
              }
	    }
          }
	}
      }
    }
  }

/*
 * for each pair of bends which share a common bond, add a torsion
 */

/*
 * just look at the heavy-atom skeleton. return true if i is a terminal
 * atom.
 */

int
IntCoorGen::hterminal(Molecule& m, BitArrayLTri& bonds, int i)
{
  int nh=0;
  for (int j=0; j < m.natom(); j++)
    if (bonds(i,j) && m.Z(j) > 1) nh++;
  return (nh==1);
}

void
IntCoorGen::add_tors(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
{
  int i,j,k,l;
  int labelc=0;
  char label[80];

  int n = m.natom();

  double thres = cos(linear_tors_thres_*M_PI/180.0);

  for(j=0; j < n; j++) {
    for(k=0; k < j; k++) {
      if(bonds(j,k)) {
        for(i=0; i < n; i++) {
          if(k==i) continue;

         // no hydrogen torsions, ok?
	  if (m.Z(i) == 1 && !hterminal(m,bonds,j)) continue;

          if (bonds(j,i)) {
            int is_linear = 0;
	    if (cos_ijk(m,i,j,k)>=thres) is_linear = 1;

            for (l=0; l < n; l++) {
              if (l==j || l==i) continue;

             // no hydrogen torsions, ok?
	      if (m.Z(l) == 1 && !hterminal(m,bonds,k))
                continue;

              if (bonds(k,l)) {
		if(cos_ijk(m,j,k,l)>=thres) is_linear = 1;

                if (is_linear && linear_stors_) {
                    labelc++;
                    sprintf(label,"st%d",labelc);
                    list->add(new ScaledTors(label,l+1,k+1,j+1,i+1));
                  }
                if (!is_linear || linear_tors_) {
                    labelc++;
                    sprintf(label,"t%d",labelc);
                    list->add(new Tors(label,l+1,k+1,j+1,i+1));
                  }
		}
	      }
	    }
          }
	}
      }
    }
  }

void
IntCoorGen::add_out(const Ref<SetIntCoor>& list, BitArrayLTri& bonds, Molecule& m)
{
  int i,j,k,l;
  int labelc=0;
  char label[80];

  int n = m.natom();

 // first find all tri-coordinate atoms
  for(i=0; i < n; i++) {
    if(bonds.degree(i)!=3) continue;

   // then look for terminal atoms connected to i
    for(j=0; j < n; j++) {
      if(bonds(i,j) && bonds.degree(j)==1) {

        for(k=0; k < n; k++) {
          if(k!=j && bonds(i,k)) {
            for(l=0; l < k; l++) {
              if(l!=j && bonds(i,l)) {
		labelc++;
		sprintf(label,"o%d",labelc);
		list->add(new Out(label,j+1,i+1,l+1,k+1));
		}
	      }
	    }
          }
	}
      }
    }
  }

int
IntCoorGen::nearest_contact(int i, Molecule& m)
{
  double d=-1.0;
  int n=0;

  SCVector3 ri(m.r(i));
  
  for (int j=0; j < m.natom(); j++) {
    SCVector3 rj(m.r(j));
    double td = ri.dist(rj);
    if (j==i)
      continue;
    else if (d < 0 || td < d) {
      d = td;
      n = j;
    }
  }
  
  return n;
}

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

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