molecule.cc

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

    }
  }

  cmat_diag(inert, evals, evecs, 3, 1, 1e-14);

 // cleanup evecs
  for (i=0; i < 3; i++) {
    for (int j=0; j < 3; j++) {
      if (fabs(evecs[i][j]) < 1.0e-5) {
        evecs[i][j]=0.0;
      }
    }
  }

  for (i=0; i<natom(); i++) {
      double a[3];
      a[0] = r(i,0); a[1] = r(i,1); a[2] = r(i,2);
      for (j=0; j<3; j++) {
          double e = 0.0;
          for (k=0; k<3; k++) {
              e += a[k] * evecs[k][j];
            }
          r_[i][j] = e;
        }
    }

  if (!trans_frame) return;
  
  SymmetryOperation tso=point_group()->symm_frame();

  for (i=0; i < 3; i++) {
    for (int j=0; j < 3; j++) {
      double t=0;
      for (int k=0; k < 3; k++) t += tso[k][j]*evecs[k][i];
      pg_->symm_frame()[i][j] = t;
    }
  }
}

void
Molecule::transform_to_symmetry_frame()
{
  int i,j,k;
  double t[3][3];

  SymmetryOperation tso=point_group()->symm_frame();

  for (i=0; i<3; i++) {
      for (j=0; j<3; j++) {
          t[i][j] = tso[i][j];
        }
    }

  for (i=0; i<natom(); i++) {
      double a[3];
      a[0] = r(i,0); a[1] = r(i,1); a[2] = r(i,2);
      for (j=0; j<3; j++) {
          double e = 0.0;
          for (k=0; k<3; k++) {
              e += a[k] * t[k][j];
            }
          r_[i][j] = e;
        }
    }

  for (i=0; i<3; i++) {
      for (j=0; j<3; j++) {
          double e=0;
          for (k=0; k<3; k++) e += tso[k][j]*t[k][i];
          pg_->symm_frame()[i][j] = e;
        }
    }
}

// given a molecule, make sure that equivalent centers have coordinates
// that really map into each other

void
Molecule::cleanup_molecule(double tol)
{
  // if symmetry is c1, do nothing else
  if (!strcmp(point_group()->symbol(),"c1")) return;

  int i;
  SCVector3 up,np,ap;
  SymmetryOperation so;
  CharacterTable ct = point_group()->char_table();

  // first clean up the unique atoms by replacing each coordinate with the
  // average of coordinates obtained by applying all symmetry operations to
  // the original atom, iff the new atom ends up near the original atom
  for (i=0; i < nunique(); i++) {
      // up will store the original coordinates of unique atom i
      up = r(unique(i));
      // ap will hold the average coordinate (times the number of coordinates)
      // initialize it to the E result
      ap = up;
      int ncoor = 1;
      // loop through all sym ops except E
      for (int g=1; g < ct.order(); g++) {
          so = ct.symm_operation(g);
          for (int ii=0; ii < 3; ii++) {
              np[ii]=0;
              for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * up[jj];
            }
          if (np.dist(up) < 0.1) {
              for (int jj=0; jj < 3; jj++) ap[jj] += np[jj];
              ncoor++;
            }
        }
      // replace the unique coordinate with the average coordinate
      r_[unique(i)][0] = ap[0] / ncoor;
      r_[unique(i)][1] = ap[1] / ncoor;
      r_[unique(i)][2] = ap[2] / ncoor;
    }

  // find the atoms equivalent to each unique atom and eliminate
  // numerical errors that may be in the equivalent atom's coordinates

  // loop through unique atoms
  for (i=0; i < nunique(); i++) {
      // up will store the coordinates of unique atom i
      up = r(unique(i));

      // loop through all sym ops except E
      for (int g=1; g < ct.order(); g++) {
          so = ct.symm_operation(g);
          for (int ii=0; ii < 3; ii++) {
              np[ii]=0;
              for (int jj=0; jj < 3; jj++) np[ii] += so(ii,jj) * up[jj];
            }

          // loop through equivalent atoms
          int found = 0;
          for (int j=0; j < natom(); j++) {
              // see if j is generated from i
              if (np.dist(SCVector3(r(j))) < tol) {
                  r_[j][0] = np[0];
                  r_[j][1] = np[1];
                  r_[j][2] = np[2];
                  found = 1;
                }
            }
          if (!found) {
              ExEnv::err0()
                   << "Molecule: cleanup: couldn't find atom at " << np << endl
                   << "  transforming uniq atom " << i << " at " << up << endl
                   << "  with symmetry op " << g << ":" << endl;
              so.print(ExEnv::err0());
              abort();
            }
        }
    }

}

///////////////////////////////////////////////////////////////////
// Compute the principal axes and the principal moments of inertia
///////////////////////////////////////////////////////////////////

void
Molecule::principal_moments_of_inertia(double *evals, double **evecs) const
{

  // The principal moments of inertia are computed in amu*angstrom^2
  // evals: principal moments of inertia
  // evecs: principal axes (optional argument)

  const double au_to_angs = 0.2800283608302436; // for moments of inertia

  double *inert[3];  // inertia tensor

  int i, j;
  int delete_evecs = 0;

  // (allocate and) initialize evecs, evals, and inert
  if (!evecs) {
    evecs = new double*[3];
    for (i=0; i<3; i++) evecs[i] = new double[3];
    delete_evecs = 1;
    }
  for (i=0; i<3; i++) {
    inert[i] = new double[3];
    memset(inert[i],'\0',sizeof(double)*3);
    memset(evecs[i],'\0',sizeof(double)*3);
    }
  memset(evals,'\0',sizeof(double)*3);

  SCVector3 com = center_of_mass();

  // compute inertia tensor
  SCVector3 ac;
  for (i=0; i<natom(); i++) {
    ac = r(i);
    // compute moments of inertia wrt center of mass
    for (j=0; j<3; j++) ac(j) -= com(j);
    double m=au_to_angs*mass(i);
    inert[0][0] += m * (ac[1]*ac[1] + ac[2]*ac[2]);
    inert[1][0] -= m * ac[0]*ac[1];
    inert[1][1] += m * (ac[0]*ac[0] + ac[2]*ac[2]);
    inert[2][0] -= m * ac[0]*ac[2];
    inert[2][1] -= m * ac[1]*ac[2];
    inert[2][2] += m * (ac[0]*ac[0] + ac[1]*ac[1]);
    }
  inert[0][1] = inert[1][0];
  inert[0][2] = inert[2][0];
  inert[1][2] = inert[2][1];

  cmat_diag(inert, evals, evecs, 3, 1, 1e-14);

  if (delete_evecs) {
    for (i=0; i<3; i++) delete[] evecs[i];
    delete[] evecs;
    }
  for (i=0; i<3; i++) {
    delete[] inert[i];
    }
}

int
Molecule::n_core_electrons()
{
  int i,n=0;
  for (i=0; i<natom(); i++) {
      if (charge(i) == 0.0) continue;
      int z = Z_[i];
      if (z > 2) n += 2;
      if (z > 10) n += 8;
      if (z > 18) n += 8;
      if (z > 30) n += 10;
      if (z > 36) n += 8;
      if (z > 48) n += 10;
      if (z > 54) n += 8;
      if (z > 72) {
          ExEnv::errn() << "Molecule::n_core_electrons: atomic number too large"
               << endl;
          abort();
        }
    }
  return n;
}

int
Molecule::max_z()
{
  int i, maxz=0;
  for (i=0; i<natom(); i++) {
      int z = Z_[i];
      if (z>maxz) maxz = z;
    }
  return maxz;
}

void
Molecule::read_pdb(const char *filename)
{
  clear();
  ifstream in(filename);
  Ref<Units> units = new Units("angstrom");
  while (in.good()) {
      const int max_line = 80;
      char line[max_line];
      in.getline(line,max_line);
      char *endofline = (char*) memchr(line, 0, max_line);
      if (endofline) memset(endofline, ' ', &line[max_line-1] - endofline);
      if (!in.good()) break;
      if (strncmp(line,"ATOM  ",6) == 0
          ||strncmp(line,"HETATM",6) == 0) {

          char element[3];
          strncpy(element,&line[76],2); element[2] = '\0';
          char name[5];
          strncpy(name,&line[12],4); name[4] = '\0';
          if (element[0]==' '&&element[1]==' ') {
              // no element was given so get the element from the atom name
              if (name[0]!=' '&&name[3]!=' ') {
                  // some of the atom label may have been
                  // pushed over into the element fields
                  // so check the residue
                  char resName[4];
                  strncpy(resName,&line[17],3); resName[3] = '\0';
                  if (strncmp(line,"ATOM  ",6)==0&&(0
                      ||strcmp(resName,"ALA")==0||strcmp(resName,"A  ")==0
                      ||strcmp(resName,"ARG")==0||strcmp(resName,"R  ")==0
                      ||strcmp(resName,"ASN")==0||strcmp(resName,"N  ")==0
                      ||strcmp(resName,"ASP")==0||strcmp(resName,"D  ")==0
                      ||strcmp(resName,"ASX")==0||strcmp(resName,"B  ")==0
                      ||strcmp(resName,"CYS")==0||strcmp(resName,"C  ")==0
                      ||strcmp(resName,"GLN")==0||strcmp(resName,"Q  ")==0
                      ||strcmp(resName,"GLU")==0||strcmp(resName,"E  ")==0
                      ||strcmp(resName,"GLX")==0||strcmp(resName,"Z  ")==0
                      ||strcmp(resName,"GLY")==0||strcmp(resName,"G  ")==0
                      ||strcmp(resName,"HIS")==0||strcmp(resName,"H  ")==0
                      ||strcmp(resName,"ILE")==0||strcmp(resName,"I  ")==0
                      ||strcmp(resName,"LEU")==0||strcmp(resName,"L  ")==0
                      ||strcmp(resName,"LYS")==0||strcmp(resName,"K  ")==0
                      ||strcmp(resName,"MET")==0||strcmp(resName,"M  ")==0
                      ||strcmp(resName,"PHE")==0||strcmp(resName,"F  ")==0
                      ||strcmp(resName,"PRO")==0||strcmp(resName,"P  ")==0
                      ||strcmp(resName,"SER")==0||strcmp(resName,"S  ")==0
                      ||strcmp(resName,"THR")==0||strcmp(resName,"T  ")==0
                      ||strcmp(resName,"TRP")==0||strcmp(resName,"W  ")==0
                      ||strcmp(resName,"TYR")==0||strcmp(resName,"Y  ")==0
                      ||strcmp(resName,"VAL")==0||strcmp(resName,"V  ")==0
                      ||strcmp(resName,"A  ")==0
                      ||strcmp(resName,"+A ")==0
                      ||strcmp(resName,"C  ")==0
                      ||strcmp(resName,"+C ")==0
                      ||strcmp(resName,"G  ")==0
                      ||strcmp(resName,"+G ")==0
                      ||strcmp(resName,"I  ")==0
                      ||strcmp(resName,"+I ")==0
                      ||strcmp(resName,"T  ")==0
                      ||strcmp(resName,"+T ")==0
                      ||strcmp(resName,"U  ")==0
                      ||strcmp(resName,"+U ")==0
                      )) {
                      // there no two letter elements for these cases
                      element[0] = name[0];
                      element[1] = '\0';
                    }
                }
              else {
                  strncpy(element,name,2); element[2] = '\0';
                }
            }
          if (element[0] == ' ') {
              element[0] = element[1];
              element[1] = '\0';
            }

          int Z = AtomInfo::string_to_Z(element);

          char field[9];
          strncpy(field,&line[30],8); field[8] = '\0';
          double x = atof(field);
          strncpy(field,&line[38],8); field[8] = '\0';
          double y = atof(field);
          strncpy(field,&line[46],8); field[8] = '\0';
          double z = atof(field);
          add_atom(Z,
                   x*units->to_atomic_units(),
                   y*units->to_atomic_units(),
                   z*units->to_atomic_units());
        }
      else {
        // skip to next record
        }
    }
}

void
Molecule::print_pdb(ostream& os, char *title) const
{
  Ref<Units> u = new Units("angstrom");
  double bohr = u->from_atomic_units();

  if (title)
      os << scprintf("%-10s%-60s\n","COMPND",title);
  else
      os << scprintf("%-10s%-60s\n","COMPND","Title");

  if (title)
      os << scprintf("REMARK   %s\n", title);

  int i;
  for (i=0; i < natom(); i++) {
    char symb[4];
    sprintf(symb,"%s1",AtomInfo::symbol(Z_[i]));

    os << scprintf(
        "HETATM%5d  %-3s UNK %5d    %8.3f%8.3f%8.3f  0.00  0.00   0\n",
        i+1, symb, 0, r(i,0)*bohr, r(i,1)*bohr, r(i,2)*bohr);
  }

  for (i=0; i < natom(); i++) {
    double at_rad_i = atominfo_->atomic_radius(Z_[i]);
    SCVector3 ai(r(i));

    os << scprintf("CONECT%5d",i+1);

    for (int j=0; j < natom(); j++) {

      if (j==i) continue;

      double at_rad_j = atominfo_->atomic_radius(Z_[j]);
      SCVector3 aj(r(j));

      if (ai.dist(aj) < 1.1*(at_rad_i+at_rad_j))
          os << scprintf("%5d",j+1);
    }

    os << endl;
  }

  os << "END" << endl;
  os.flush();
}

double
Molecule::mass(int atom) const
{
  if (!mass_ || mass_[atom] == 0) {
      return atominfo_->mass(Z_[atom]);
    }
  return mass_[atom];
}

const char *
Molecule::label(int atom) const
{
  if (!labels_) return 0;
  return labels_[atom];
}

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

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