C++ Topology::AddBond方法代码示例

/** Read file as a Tripos Mol2 file. */
int Parm_Mol2::ReadParm(FileName const& fname, Topology &parmOut) {
  Mol2File infile;
  if (infile.OpenRead(fname)) return 1;
  mprintf("    Reading Mol2 file %s as topology file.\n",infile.Filename().base());
  // Get @<TRIPOS>MOLECULE information
  if (infile.ReadMolecule()) return 1;
  parmOut.SetParmName( infile.Mol2Title(), infile.Filename() );

  // Get @<TRIPOS>ATOM information
  if (infile.ScanTo( Mol2File::ATOM)) return 1;
  double XYZ[3];
  for (int atom=0; atom < infile.Mol2Natoms(); atom++) {
    if ( infile.Mol2XYZ(XYZ) ) return 1;
    parmOut.AddTopAtom( infile.Mol2Atom(), infile.Mol2Residue(), XYZ );
  }

  // Get @<TRIPOS>BOND information [optional]
  int at1 = 0;
  int at2 = 0;
  if (infile.ScanTo(Mol2File::BOND)==0) {
    for (int bond=0; bond < infile.Mol2Nbonds(); bond++) {
      if (infile.Mol2Bond(at1, at2)) return 1;
      // mol2 atom #s start from 1
      parmOut.AddBond(at1-1, at2-1);
    }
    needsBondSearch_ = false;
  } else {
    mprintf("      Mol2 file does not contain bond information.\n");
    needsBondSearch_ = true;
  }

  // No box
  parmOut.SetParmBox( Box() );

  mprintf("    Mol2 contains %i atoms, %i residues,\n", parmOut.Natom(),parmOut.Nres());
  //mprintf("    %i bonds to H, %i other bonds.\n", parmOut.NbondsWithH,parmOut.NbondsWithoutH);

  infile.CloseFile();

  return 0;
}

/** Read file as a Tinker file. */
int Parm_Tinker::ReadParm(FileName const& fname, Topology &parmOut) {
  TinkerFile infile;
  infile.SetTinkerName( fname );
  if (infile.OpenTinker()) return 1;
  mprintf("\tReading Tinker file %s as topology file.\n",infile.Filename().base());
  // Allocate memory for coordinates.
  double* Coords = new double[ infile.TinkerNatom() * 3 ];
  std::vector<int> Bonds;
  std::vector<Atom> Atoms = infile.ReadTinkerAtoms(Coords, Bonds);
  if (Atoms.empty()) return 1;
  // Use up to first 3 chars of title as residue name.
  std::string resname;
  for (std::string::const_iterator c = infile.TinkerTitle().begin();
                                   c != infile.TinkerTitle().end(); ++c)
    resname += *c;
  if (resname.size() > 3) resname.resize(3);
  Residue tinker_res( resname, 0, ' ', ' ' );
  // Put atoms into topology
  const double* XYZ = Coords;
  for (std::vector<Atom>::const_iterator atom = Atoms.begin();
                                         atom != Atoms.end();
                                       ++atom, XYZ += 3)
    parmOut.AddTopAtom( *atom, tinker_res, XYZ );
  delete[] Coords;
  // Add bond information
  for (std::vector<int>::const_iterator bond = Bonds.begin();
                                        bond != Bonds.end(); bond += 2)
    parmOut.AddBond( *bond, *(bond+1) );
  // Try to set up residue info based on bonds.
  if (parmOut.Setup_NoResInfo()) return 1;
  // Set topology box info.
  parmOut.SetParmBox( infile.TinkerBox() );
  parmOut.SetParmName( infile.TinkerTitle(), infile.Filename() );
  mprintf("\tTinker file contains %i atoms, %i residues,\n", parmOut.Natom(),parmOut.Nres());
  //mprintf("    %i bonds to H, %i other bonds.\n", parmOut.NbondsWithH,parmOut.NbondsWithoutH);

  infile.CloseFile();

  return 0;
}

/** Open the Charmm PSF file specified by filename and set up topology data.
  * Mask selection requires natom, nres, names, resnames, resnums.
  */
int Parm_CharmmPsf::ReadParm(FileName const& fname, Topology &parmOut) {
  const size_t TAGSIZE = 10; 
  char tag[TAGSIZE];
  tag[0]='\0';

  CpptrajFile infile;
  if (infile.OpenRead(fname)) return 1;
  mprintf("    Reading Charmm PSF file %s as topology file.\n",infile.Filename().base());
  // Read the first line, should contain PSF...
  const char* buffer = 0;
  if ( (buffer=infile.NextLine()) == 0 ) return 1;
  // Advance to <ntitle> !NTITLE
  int ntitle = FindTag(tag, "!NTITLE", 7, infile); 
  // Only read in 1st title. Skip any asterisks.
  std::string psftitle;
  if (ntitle > 0) {
    buffer = infile.NextLine();
    const char* ptr = buffer;
    while (*ptr != '\0' && (*ptr == ' ' || *ptr == '*')) ++ptr;
    psftitle.assign( ptr );
  }
  parmOut.SetParmName( NoTrailingWhitespace(psftitle), infile.Filename() );
  // Advance to <natom> !NATOM
  int natom = FindTag(tag, "!NATOM", 6, infile);
  if (debug_>0) mprintf("\tPSF: !NATOM tag found, natom=%i\n", natom);
  // If no atoms, probably issue with PSF file
  if (natom < 1) {
    mprinterr("Error: No atoms in PSF file.\n");
    return 1;
  }
  // Read the next natom lines
  int psfresnum = 0;
  char psfresname[6];
  char psfname[6];
  char psftype[6];
  double psfcharge;
  double psfmass;
  for (int atom=0; atom < natom; atom++) {
    if ( (buffer=infile.NextLine()) == 0 ) {
      mprinterr("Error: ReadParmPSF(): Reading atom %i\n",atom+1);
      return 1;
    }
    // Read line
    // ATOM# SEGID RES# RES ATNAME ATTYPE CHRG MASS (REST OF COLUMNS ARE LIKELY FOR CMAP AND CHEQ)
    sscanf(buffer,"%*i %*s %i %s %s %s %lf %lf",&psfresnum, psfresname, 
           psfname, psftype, &psfcharge, &psfmass);
    parmOut.AddTopAtom( Atom( psfname, psfcharge, psfmass, psftype), 
                        Residue( psfresname, psfresnum, ' ', ' '), 0 );
  } // END loop over atoms 
  // Advance to <nbond> !NBOND
  int bondatoms[9];
  int nbond = FindTag(tag, "!NBOND", 6, infile);
  if (nbond > 0) {
    if (debug_>0) mprintf("\tPSF: !NBOND tag found, nbond=%i\n", nbond);
    int nlines = nbond / 4;
    if ( (nbond % 4) != 0) nlines++;
    for (int bondline=0; bondline < nlines; bondline++) {
      if ( (buffer=infile.NextLine()) == 0 ) {
        mprinterr("Error: ReadParmPSF(): Reading bond line %i\n",bondline+1);
        return 1;
      }
      // Each line has 4 pairs of atom numbers
      int nbondsread = sscanf(buffer,"%i %i %i %i %i %i %i %i",bondatoms,bondatoms+1,
                              bondatoms+2,bondatoms+3, bondatoms+4,bondatoms+5,
                              bondatoms+6,bondatoms+7);
      // NOTE: Charmm atom nums start from 1
      for (int bondidx=0; bondidx < nbondsread; bondidx+=2)
        parmOut.AddBond(bondatoms[bondidx]-1, bondatoms[bondidx+1]-1);
    }
  } else
    mprintf("Warning: PSF has no bonds.\n");
  // Advance to <nangles> !NTHETA
  int nangle = FindTag(tag, "!NTHETA", 7, infile);
  if (nangle > 0) {
    if (debug_>0) mprintf("\tPSF: !NTHETA tag found, nangle=%i\n", nangle);
    int nlines = nangle / 3;
    if ( (nangle % 3) != 0) nlines++;
    for (int angleline=0; angleline < nlines; angleline++) {
      if ( (buffer=infile.NextLine()) == 0) {
        mprinterr("Error: Reading angle line %i\n", angleline+1);
        return 1;
      }
      // Each line has 3 groups of 3 atom numbers
      int nanglesread = sscanf(buffer,"%i %i %i %i %i %i %i %i %i",bondatoms,bondatoms+1,
                              bondatoms+2,bondatoms+3, bondatoms+4,bondatoms+5,
                              bondatoms+6,bondatoms+7, bondatoms+8);
      for (int angleidx=0; angleidx < nanglesread; angleidx += 3)
        parmOut.AddAngle( bondatoms[angleidx  ]-1,
                          bondatoms[angleidx+1]-1,
                          bondatoms[angleidx+2]-1 );
    }
  } else
    mprintf("Warning: PSF has no angles.\n");
  // Advance to <ndihedrals> !NPHI
  int ndihedral = FindTag(tag, "!NPHI", 5, infile);
  if (ndihedral > 0) {
    if (debug_>0) mprintf("\tPSF: !NPHI tag found, ndihedral=%i\n", ndihedral);
//.........这里部分代码省略.........

/** Search for bonds between atoms in residues and atoms in adjacent residues
  * using distance-based criterion that depends on atomic elements.
  * \param top Topology to add bonds to.
  * \param frameIn Frame containing atomic coordinates.
  * \param offset Offset to add when determining if a bond is present.
  * \param debug If > 0 print extra info.
  */
int BondSearch( Topology& top, Frame const& frameIn, double offset, int debug) {
  mprintf("\tDetermining bond info from distances.\n");
  if (frameIn.empty()) {
    mprinterr("Internal Error: No coordinates set; cannot search for bonds.\n");
    return 1;
  }
# ifdef TIMER
  Timer time_total, time_within, time_between;
  time_total.Start();
  time_within.Start();
# endif
  // ----- STEP 1: Determine bonds within residues
  for (Topology::res_iterator res = top.ResStart(); res != top.ResEnd(); ++res)
  {
    int stopatom = res->LastAtom();
    // Check for bonds between each atom in the residue.
    for (int atom1 = res->FirstAtom(); atom1 != stopatom; ++atom1) {
      Atom::AtomicElementType a1Elt = top[atom1].Element();
      // If this is a hydrogen and it already has a bond, move on.
      if (a1Elt==Atom::HYDROGEN && top[atom1].Nbonds() > 0 )
        continue;
      for (int atom2 = atom1 + 1; atom2 != stopatom; ++atom2) {
        Atom::AtomicElementType a2Elt = top[atom2].Element();
        double D2 = DIST2_NoImage(frameIn.XYZ(atom1), frameIn.XYZ(atom2) );
        double cutoff2 = Atom::GetBondLength(a1Elt, a2Elt) + offset;
        cutoff2 *= cutoff2;
        if (D2 < cutoff2) {
          top.AddBond(atom1, atom2);
          // Once a bond has been made to hydrogen move on.
          if (a1Elt==Atom::HYDROGEN) break;
        }
      }
    }
  }
# ifdef TIMER
  time_within.Stop();
  time_between.Start();
# endif
  // ----- STEP 2: Determine bonds between adjacent residues
  Topology::mol_iterator nextmol = top.MolStart();
  if (top.Nmol() > 0)
    ++nextmol;
  for (Topology::res_iterator res = top.ResStart() + 1; res != top.ResEnd(); ++res)
  {
    // If molecule information is already present, check if first atom of 
    // this residue >= first atom of next molecule, which indicates this
    // residue and the previous residue are in different molecules.
    if ( (nextmol != top.MolEnd()) &&
         (res->FirstAtom() >= nextmol->BeginAtom()) )
    {
      ++nextmol;
      continue;
    }
    // If this residue is recognized as solvent, no need to check previous or
    // next residue
    if ( res->NameIsSolvent() ) {
      ++res;
      if (res == top.ResEnd()) break;
      continue;
    }
    // Get previous residue
    Topology::res_iterator previous_res = res - 1;
    // If previous residue is recognized as solvent, no need to check previous.
    if ( previous_res->NameIsSolvent() ) continue;
    // Get previous residue start atom
    int startatom = previous_res->FirstAtom();
    // Previous residue stop atom, this residue start atom
    int midatom = res->FirstAtom();
    // This residue stop atom
    int stopatom = res->LastAtom();
    // Check for bonds between adjacent residues
    for (int atom1 = startatom; atom1 != midatom; atom1++) {
      Atom::AtomicElementType a1Elt = top[atom1].Element();
      if (a1Elt==Atom::HYDROGEN) continue;
      for (int atom2 = midatom; atom2 != stopatom; atom2++) {
        Atom::AtomicElementType a2Elt = top[atom2].Element();
        if (a2Elt==Atom::HYDROGEN) continue;
        double D2 = DIST2_NoImage(frameIn.XYZ(atom1), frameIn.XYZ(atom2) );
        double cutoff2 = Atom::GetBondLength(a1Elt, a2Elt) + offset;
        cutoff2 *= cutoff2;
        if (D2 < cutoff2)
          top.AddBond(atom1, atom2);
      }
    }
  }
# ifdef TIMER
  time_between.Stop();
  time_total.Stop();
  time_within.WriteTiming(2, "Distances within residues", time_total.Total());
  time_between.WriteTiming(2, "Distances between residues", time_total.Total());
  time_total.WriteTiming(1, "Total for determining bonds via distances");
# endif
  if (debug > 0)
//.........这里部分代码省略.........