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

/** Write file containing only cut atoms and energies as charges. */
int Action_Pairwise::WriteCutFrame(int frameNum, Topology const& Parm, AtomMask const& CutMask, 
                                   Darray const& CutCharges,
                                   Frame const& frame, std::string const& outfilename) 
{
  if (CutMask.Nselected() != (int)CutCharges.size()) {
    mprinterr("Error: WriteCutFrame: # of charges (%u) != # mask atoms (%i)\n",
              CutCharges.size(), CutMask.Nselected());
    return 1;
  }
  Frame CutFrame(frame, CutMask);
  Topology* CutParm = Parm.modifyStateByMask( CutMask );
  if (CutParm == 0) return 1;
  // Set new charges
  for (int i = 0; i != CutParm->Natom(); i++)
    CutParm->SetAtom(i).SetCharge( CutCharges[i] );
  int err = 0;
  Trajout_Single tout;
  if (tout.PrepareTrajWrite(outfilename, "multi", CutParm, CoordinateInfo(), 1,
                            TrajectoryFile::MOL2FILE))
  {
    mprinterr("Error: Could not set up cut mol2 file %s\n", outfilename.c_str());
    err = 1;
  } else {
    tout.WriteSingle(frameNum, CutFrame);
    tout.EndTraj();
  }
  delete CutParm;
  return err;
}

/** Set up the exclusion list based on the given mask and parm.
  * \return the total number of interactions, -1 on error.
  */
int Action_Pairwise::SetupNonbondParm(AtomMask const& maskIn, Topology const& ParmIn) {
  // Check if LJ parameters present - need at least 2 atoms for it to matter.
  if (ParmIn.Natom() > 1 && !ParmIn.Nonbond().HasNonbond()) {
    mprinterr("Error: Topology does not have LJ information.\n");
    return -1;
  }

  // Determine the actual number of pairwise interactions that will be calcd.
  // This is ((N^2 - N) / 2) - SUM[ #excluded atoms]
  int N_interactions = ((maskIn.Nselected() * maskIn.Nselected()) - maskIn.Nselected()) / 2;
  for (AtomMask::const_iterator at = maskIn.begin(); at != maskIn.end(); ++at)
    N_interactions -= ParmIn[ *at ].Nexcluded();

  // DEBUG - Print total number of interactions for this parm
  mprintf("\t%i interactions for this parm.\n",N_interactions);

  // DEBUG - Print exclusion list for each atom
  /*for (unsigned int atom = 0; atom < exclusionList.size(); atom++) {
    mprintf("\t%8u:",atom + 1);
    for (std::vector<int>::iterator eat = exclusionList[atom].begin();
                                    eat != exclusionList[atom].end();
                                    eat++)
    {
      mprintf(" %i",*eat + 1);
    }
    mprintf("\n");
  }*/
  return N_interactions;
}

int DataSet_Coords_TRJ::CoordsSetup(Topology const& topIn, CoordinateInfo const& cInfoIn ) {
  if (trajinList_.empty()) {
    top_ = topIn;
    cInfo_ = cInfoIn;
  } else {
    if ( topIn.Natom() != top_.Natom() ) {
      mprinterr("Error: For TRAJ data set currently all trajectories must have same number\n"
                "Error:  of atoms: %i != %i. Recommended course of action is to\n"
                "Error:  create a trajectory where all frames have been stripped to the same\n"
                "Error:  number of atoms first.\n", topIn.Natom(), top_.Natom());
      return 1;
    }
    // Since velocity/force info is not always allocated in Frame, if one traj
    // has velocity/force info ensure that all do.
    if (cInfoIn.HasVel()) cInfo_.SetVelocity( true );
    if (cInfoIn.HasForce()) cInfo_.SetForce( true );
  }
  return 0;
}

/** Sets the temporary charge array and makes sure that we have the necessary
  * parameters in our topology to calculate nonbonded energy terms
  */
int Action_LIE::SetupParms(Topology const& ParmIn) {
  if (!ParmIn.Nonbond().HasNonbond()) {
    mprinterr("Error: Topology does not have LJ information.\n");
    return 1;
  }
  // Store the charges
  atom_charge_.clear();
  atom_charge_.reserve( ParmIn.Natom() );
  for (Topology::atom_iterator atom = ParmIn.begin();
                               atom != ParmIn.end(); ++atom)
    atom_charge_.push_back( atom->Charge() * Constants::ELECTOAMBER / sqrt(dielc_) );
  return 0;
}

Exec::RetType Exec_ParmStrip::Execute(CpptrajState& State, ArgList& argIn) {
    Topology* parm = State.DSL().GetTopByIndex( argIn );
  if (parm == 0) return CpptrajState::ERR;
  // Check if this topology has already been used to set up an input
  // trajectory, as this will break the traj read.
  bool topology_in_use = false;
  const char* fname = 0;
  for (TrajinList::trajin_it tIn = State.InputTrajList().trajin_begin();
                             tIn != State.InputTrajList().trajin_end(); ++tIn)
    if ( (*tIn)->Traj().Parm() == parm ) {
      topology_in_use = true;
      fname = (*tIn)->Traj().Filename().full();
      break;
    }
  if (!topology_in_use) {
    for (TrajinList::ensemble_it eIn = State.InputTrajList().ensemble_begin();
                                 eIn != State.InputTrajList().ensemble_end(); ++eIn)
      if ( (*eIn)->Traj().Parm() == parm ) {
        topology_in_use = true;
        fname = (*eIn)->Traj().Filename().full();
        break;
      }
  }
  if (topology_in_use) {
    mprinterr("Error: Topology '%s' has already been used to set up trajectory '%s'.\n"
              "Error:   To strip this topology use the 'strip' action.\n",
              parm->c_str(), fname);
    return CpptrajState::ERR;
  }
  AtomMask tempMask( argIn.GetMaskNext() );
  // Since want to keep atoms outside mask, invert selection
  tempMask.InvertMaskExpression();
  if (parm->SetupIntegerMask( tempMask )) return CpptrajState::ERR;
  mprintf("\tStripping atoms in mask [%s] (%i) from %s\n",tempMask.MaskString(),
           parm->Natom() - tempMask.Nselected(), parm->c_str());
  Topology* tempParm = parm->modifyStateByMask(tempMask);
  if (tempParm==0) {
    mprinterr("Error: %s: Could not strip parm.\n", argIn.Command());
    return CpptrajState::ERR;
  } else {
    // Replace parm with stripped version
    *parm = *tempParm;
    parm->Brief("Stripped parm:");
    delete tempParm;
  }
  return CpptrajState::OK;
}

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

/** Set up masks. */
int Cpptraj::MaskArray::SetupMasks(AtomMask const& maskIn, Topology const& topIn)
{
  if (type_ == BY_MOLECULE && topIn.Nmol() < 1) {
    mprintf("Warning: '%s' has no molecule information, cannot setup by molecule.\n",
             topIn.c_str());
    return 1;
  }
  masks_.clear();
  if ( maskIn.None() ) {
    mprintf("Warning: Nothing selected by mask '%s'\n", maskIn.MaskString());
    return 0;
  }
  int last = -1;
  int current = 0;
  maxAtomsPerMask_ = 0;
  sameNumAtomsPerMask_ = true;
  for (AtomMask::const_iterator atm = maskIn.begin(); atm != maskIn.end(); ++atm)
  {
    switch (type_) {
      case BY_ATOM     : current = *atm; break;
      case BY_RESIDUE  : current = topIn[*atm].ResNum(); break;
      case BY_MOLECULE : current = topIn[*atm].MolNum(); break;
    }
    if (current != last) {
      if (!masks_.empty())
        checkAtomsPerMask( masks_.back().Nselected() );
      masks_.push_back( AtomMask() );
      masks_.back().SetNatoms( topIn.Natom() );
    }
    masks_.back().AddSelectedAtom( *atm );
    last = current;
  }
  if (!masks_.empty())
    checkAtomsPerMask( masks_.back().Nselected() );

  return 0;
}

/** An atom pair list consists of 2 values for each entry, a beginning
  * index and ending index. For molecules and residues this is the first
  * and just beyond the last atom; for atoms it is just the atom itself
  * twice.
  */
Image::PairType Image::CreatePairList(Topology const& Parm, Mode modeIn,
                                       std::string const& maskExpression)
{
  PairType atomPairs;
  // Set up mask based on desired imaging mode.
  if ( modeIn == BYMOL || modeIn == BYRES ) {
    CharMask cmask( maskExpression );
    if ( Parm.SetupCharMask( cmask ) ) return atomPairs;
    cmask.MaskInfo();
    if (cmask.None()) return atomPairs;
    // Set up atom range for each entity to be imaged.
    if (modeIn == BYMOL) {
      atomPairs.reserve( Parm.Nmol()*2 );
      for (Topology::mol_iterator mol = Parm.MolStart();
                                  mol != Parm.MolEnd(); ++mol)
        CheckRange( atomPairs, cmask, mol->BeginAtom(), mol->EndAtom());
    } else { // BYRES
      atomPairs.reserve( Parm.Nres()*2 );
      for (Topology::res_iterator residue = Parm.ResStart();
                                  residue != Parm.ResEnd(); ++residue)
        CheckRange( atomPairs, cmask, residue->FirstAtom(), residue->LastAtom() );
    }
  } else { // BYATOM
    AtomMask imask( maskExpression );
    if ( Parm.SetupIntegerMask( imask ) ) return atomPairs;
    imask.MaskInfo();
    if (imask.None()) return atomPairs;
    atomPairs.reserve( Parm.Natom()*2 );
    for (AtomMask::const_iterator atom = imask.begin(); atom != imask.end(); ++atom) {
      atomPairs.push_back(  *atom    );
      atomPairs.push_back( (*atom)+1 );
    }
  }
//  mprintf("\tNumber of %ss to be imaged is %zu based on mask '%s'\n",
//           ModeString[modeIn], atomPairs.size()/2, maskIn.MaskString());
  return atomPairs;
}

//.........这里部分代码省略.........
      if (Query[qres] != qref.Parm().Res(qres).SingleCharName()) {
        mprintf("Warning: Potential residue mismatch: Query %s reference %s\n",
                Residue::ConvertResName(Query[qres]), qref.Parm().Res(qres).c_str());
      }
    }
  }
  // Build subject using coordinate from reference.
  //AtomMask sMask; // Contain atoms that should be in sTop
  Topology sTop;
  Frame sFrame;
  Iarray placeHolder; // Atom indices of placeholder residues.
  for (unsigned int sres = 0; sres != Sbjct.size(); sres++) {
    int qres = Smap[sres];
    NameType SresName( Residue::ConvertResName(Sbjct[sres]) );
    if (qres != -1) {
      Residue const& QR = qref.Parm().Res(qres);
      Residue SR(SresName, sres+1, ' ', QR.ChainID());
      if (Query[qres] == Sbjct[sres]) { // Exact match. All non-H atoms.
        for (int qat = QR.FirstAtom(); qat != QR.LastAtom(); qat++)
        {
          if (qref.Parm()[qat].Element() != Atom::HYDROGEN)
            sTop.AddTopAtom( qref.Parm()[qat], SR );
            sFrame.AddXYZ( qref.Coord().XYZ(qat) );
            //sMask.AddAtom(qat);
        }
      } else { // Partial match. Copy only backbone and CB.
        for (int qat = QR.FirstAtom(); qat != QR.LastAtom(); qat++)
        {
          if ( qref.Parm()[qat].Name().Match("N" ) ||
               qref.Parm()[qat].Name().Match("CA") ||
               qref.Parm()[qat].Name().Match("CB") ||
               qref.Parm()[qat].Name().Match("C" ) ||
               qref.Parm()[qat].Name().Match("O" ) )
          {
            sTop.AddTopAtom( qref.Parm()[qat], SR );
            sFrame.AddXYZ( qref.Coord().XYZ(qat) );
          }
        }
      }
    } else {
      // Residue in query does not exist for subject. Just put placeholder CA for now.
      Vec3 Zero(0.0);
      placeHolder.push_back( sTop.Natom() );
      sTop.AddTopAtom( Atom("CA", "C "), Residue(SresName, sres+1, ' ', ' ') );
      sFrame.AddXYZ( Zero.Dptr() );
    }
  }
  //sTop.PrintAtomInfo("*");
  mprintf("\tPlaceholder residue indices:");
  for (Iarray::const_iterator p = placeHolder.begin(); p != placeHolder.end(); ++p)
    mprintf(" %i", *p + 1);
  mprintf("\n");
  // Try to give placeholders more reasonable coordinates.
  if (!placeHolder.empty()) {
    Iarray current_indices;
    unsigned int pidx = 0;
    while (pidx < placeHolder.size()) {
      if (current_indices.empty()) {
        current_indices.push_back( placeHolder[pidx++] );
        // Search for the end of this segment
        for (; pidx != placeHolder.size(); pidx++) {
          if (placeHolder[pidx] - current_indices.back() > 1) break;
          current_indices.push_back( placeHolder[pidx] );
        }
        // DEBUG
        mprintf("\tSegment:");
        for (Iarray::const_iterator it = current_indices.begin();
                                    it != current_indices.end(); ++it)
          mprintf(" %i", *it + 1);
        // Get coordinates of residues bordering segment.
        int prev_res = sTop[current_indices.front()].ResNum() - 1;
        int next_res = sTop[current_indices.back() ].ResNum() + 1;
        mprintf(" (prev_res=%i, next_res=%i)\n", prev_res+1, next_res+1);
        Vec3 prev_crd(sFrame.XYZ(current_indices.front() - 1));
        Vec3 next_crd(sFrame.XYZ(current_indices.back()  + 1));
        prev_crd.Print("prev_crd");
        next_crd.Print("next_crd");
        Vec3 crd_step = (next_crd - prev_crd) / (double)(current_indices.size()+1);
        crd_step.Print("crd_step");
        double* xyz = sFrame.xAddress() + (current_indices.front() * 3);
        for (unsigned int i = 0; i != current_indices.size(); i++, xyz += 3) {
          prev_crd += crd_step;
          xyz[0] = prev_crd[0];
          xyz[1] = prev_crd[1];
          xyz[2] = prev_crd[2];
        }
        current_indices.clear();
      }
    }
  }
  //Topology* sTop = qref.Parm().partialModifyStateByMask( sMask );
  //if (sTop == 0) return 1;
  //Frame sFrame(qref.Coord(), sMask);
  // Write output traj
  Trajout_Single trajout;
  if (trajout.PrepareTrajWrite(outfilename, argIn, &sTop, CoordinateInfo(), 1, fmt)) return 1;
  if (trajout.WriteSingle(0, sFrame)) return 1;
  trajout.EndTraj();
  return 0;
}

/** Find potential symmetric atoms. All residues up to the last selected
  * residue are considered.
  */
int SymmetricRmsdCalc::SetupSymmRMSD(Topology const& topIn, AtomMask const& tgtMask, bool remapIn)
{
  // Allocate space for remapping selected atoms in target frame. This will
  // also put the correct masses in based on the mask.
  tgtRemap_.SetupFrameFromMask(tgtMask, topIn.Atoms());
  // Create map of original atom numbers to selected indices
  Iarray SelectedIdx( topIn.Natom(), -1 );
  int tgtIdx = 0;
  for (int originalAtom = 0; originalAtom != topIn.Natom(); ++originalAtom)
    if ( originalAtom == tgtMask[tgtIdx] )
      SelectedIdx[originalAtom] = tgtIdx++;
  if (debug_ > 0) {
    mprintf("DEBUG: Original atom -> Selected Index mapping:\n");
    for (int originalAtom = 0; originalAtom != topIn.Natom(); ++originalAtom)
      mprintf("\t%8i -> %8i\n", originalAtom + 1, SelectedIdx[originalAtom] + 1);
  }
  // Create initial 1 to 1 atom map for all selected atoms; indices in 
  // SymmetricAtomIndices will correspond to positions in AMap.
  AMap_.resize( tgtRemap_.Natom() );
  // Determine last selected residue.
  int last_res = topIn[tgtMask.back()].ResNum() + 1;
  mprintf("\tResidues up to %s will be considered for symmetry correction.\n",
          topIn.TruncResNameNum(last_res-1).c_str());
  // In each residue, determine which selected atoms are symmetric.
  SymmetricAtomIndices_.clear();
  AtomMap resmap;
  if (debug_ > 1) resmap.SetDebug(1);
  for (int res = 0; res < last_res; ++res) {
    AtomMap::AtomIndexArray residue_SymmetricGroups;
    if (resmap.SymmetricAtoms(topIn, residue_SymmetricGroups, res)) {
      mprinterr("Error: Finding symmetric atoms in residue '%s'\n",
                topIn.TruncResNameNum(res).c_str());
      return 1;
    }
    if (!residue_SymmetricGroups.empty()) {
      // Which atoms in symmetric groups are selected?
      bool resHasSelectedSymmAtoms = false;
      for (AtomMap::AtomIndexArray::const_iterator symmGroup = residue_SymmetricGroups.begin();
                                                   symmGroup != residue_SymmetricGroups.end();
                                                 ++symmGroup)
      {
        Iarray selectedAtomIndices;
        for (Iarray::const_iterator atnum = symmGroup->begin();
                                    atnum != symmGroup->end(); ++atnum)
        {
          if ( SelectedIdx[*atnum] != -1 )
            selectedAtomIndices.push_back( SelectedIdx[*atnum] ); // Store tgtMask indices
        }
        if (!selectedAtomIndices.empty()) {
          SymmetricAtomIndices_.push_back( selectedAtomIndices );
          resHasSelectedSymmAtoms = true;
        }
      }
      // If remapping and not all atoms in a residue are selected, warn user.
      // TODO: Should they just be considered even if not selected?
      if (remapIn && resHasSelectedSymmAtoms) {
        for (int atom = topIn.Res(res).FirstAtom(); atom != topIn.Res(res).LastAtom(); ++atom)
          if (SelectedIdx[atom] == -1) {
            mprintf("Warning: Not all atoms selected in residue '%s'. Re-mapped\n"
                    "Warning:   structures may appear distorted.\n", 
                    topIn.TruncResNameNum(res).c_str());
            break;
          }
      }
    }
  }
  if (debug_ > 0) {
    mprintf("DEBUG: Potential Symmetric Atom Groups:\n");
    for (AtomIndexArray::const_iterator symmatoms = SymmetricAtomIndices_.begin();
                                        symmatoms != SymmetricAtomIndices_.end();
                                        ++symmatoms)
    {
      mprintf("\t%8u) ", symmatoms - SymmetricAtomIndices_.begin());
      for (Iarray::const_iterator atom = symmatoms->begin();
                                  atom != symmatoms->end(); ++atom)
        mprintf(" %s(%i)", topIn.AtomMaskName(tgtMask[*atom]).c_str(), tgtMask[*atom] + 1);
      mprintf("\n");
    } 
  }
  return 0;
}

//.........这里部分代码省略.........
  // 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);
    int nlines = ndihedral / 2;
    if ( (ndihedral % 2) != 0) nlines++;
    for (int dihline = 0; dihline < nlines; dihline++) {
      if ( (buffer=infile.NextLine()) == 0) {
        mprinterr("Error: Reading dihedral line %i\n", dihline+1);
        return 1;
      }
      // Each line has 2 groups of 4 atom numbers
      int ndihread = 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);
      for (int dihidx=0; dihidx < ndihread; dihidx += 4)
        parmOut.AddDihedral( bondatoms[dihidx  ]-1,
                             bondatoms[dihidx+1]-1,
                             bondatoms[dihidx+2]-1,
                             bondatoms[dihidx+3]-1 );
    }
  } else
    mprintf("Warning: PSF has no dihedrals.\n");
  mprintf("\tPSF contains %i atoms, %i residues.\n", parmOut.Natom(), parmOut.Nres());

  infile.CloseFile();

  return 0;
}

int Parm_CharmmPsf::WriteParm(FileName const& fname, Topology const& parm) {
  // TODO: CMAP etc info
  CpptrajFile outfile;
  if (outfile.OpenWrite(fname)) return 1;
  // Write PSF
  outfile.Printf("PSF\n\n");
  // Write title
  std::string titleOut = parm.ParmName();
  titleOut.resize(78);
  outfile.Printf("%8i !NTITLE\n* %-78s\n\n", 1, titleOut.c_str());
  // Write NATOM section
  outfile.Printf("%8i !NATOM\n", parm.Natom());
  unsigned int idx = 1;
  // Make fake segment ids for now.
  char segid[2];
  segid[0] = 'A';
  segid[1] = '\0';
  mprintf("Warning: Assigning single letter segment IDs.\n");
  int currentMol = 0;
  bool inSolvent = false;
  for (Topology::atom_iterator atom = parm.begin(); atom != parm.end(); ++atom, ++idx) {
    int resnum = atom->ResNum();
    if (atom->MolNum() != currentMol) {
      if (!inSolvent) {
        inSolvent = parm.Mol(atom->MolNum()).IsSolvent();
        currentMol = atom->MolNum();
        segid[0]++;
      } else
        inSolvent = parm.Mol(atom->MolNum()).IsSolvent();
    }
    // TODO: Print type name for xplor-like PSF
    int typeindex = atom->TypeIndex() + 1;
    // If type begins with digit, assume charmm numbers were read as
    // type. Currently Amber types all begin with letters.
    if (isdigit(atom->Type()[0]))
      typeindex = convertToInteger( *(atom->Type()) );
    // ATOM# SEGID RES# RES ATNAME ATTYPE CHRG MASS (REST OF COLUMNS ARE LIKELY FOR CMAP AND CHEQ)
    outfile.Printf("%8i %-4s %-4i %-4s %-4s %4i %14.6G %9g  %10i\n", idx, segid,
                   parm.Res(resnum).OriginalResNum(), parm.Res(resnum).c_str(),
                   atom->c_str(), typeindex, atom->Charge(),
                   atom->Mass(), 0);
  }
  outfile.Printf("\n");
  // Write NBOND section
  outfile.Printf("%8u !NBOND: bonds\n", parm.Bonds().size() + parm.BondsH().size());
  idx = 1;
  for (BondArray::const_iterator bond = parm.BondsH().begin();
                                 bond != parm.BondsH().end(); ++bond, ++idx)
  {
    outfile.Printf("%8i%8i", bond->A1()+1, bond->A2()+1);
    if ((idx % 4)==0) outfile.Printf("\n"); 
  }
  for (BondArray::const_iterator bond = parm.Bonds().begin();
                                 bond != parm.Bonds().end(); ++bond, ++idx)
  {
    outfile.Printf("%8i%8i", bond->A1()+1, bond->A2()+1);
    if ((idx % 4)==0) outfile.Printf("\n"); 
  }
  if ((idx % 4)!=0) outfile.Printf("\n");
  outfile.Printf("\n");
  // Write NTHETA section
  outfile.Printf("%8u !NTHETA: angles\n", parm.Angles().size() + parm.AnglesH().size());
  idx = 1;
  for (AngleArray::const_iterator ang = parm.AnglesH().begin();
                                  ang != parm.AnglesH().end(); ++ang, ++idx)
  {
    outfile.Printf("%8i%8i%8i", ang->A1()+1, ang->A2()+1, ang->A3()+1);
    if ((idx % 3)==0) outfile.Printf("\n");
  }
  for (AngleArray::const_iterator ang = parm.Angles().begin();
                                  ang != parm.Angles().end(); ++ang, ++idx)
  {
    outfile.Printf("%8i%8i%8i", ang->A1()+1, ang->A2()+1, ang->A3()+1);
    if ((idx % 3)==0) outfile.Printf("\n");
  }
  if ((idx % 3)==0) outfile.Printf("\n");
  outfile.Printf("\n");
  // Write out NPHI section
  outfile.Printf("%8u !NPHI: dihedrals\n", parm.Dihedrals().size() + parm.DihedralsH().size());
  idx = 1;
  for (DihedralArray::const_iterator dih = parm.DihedralsH().begin();
                                     dih != parm.DihedralsH().end(); ++dih, ++idx)
  {
    outfile.Printf("%8i%8i%8i%8i", dih->A1()+1, dih->A2()+1, dih->A3()+1, dih->A4()+1);
    if ((idx % 2)==0) outfile.Printf("\n");
  }
  for (DihedralArray::const_iterator dih = parm.Dihedrals().begin();
                                     dih != parm.Dihedrals().end(); ++dih, ++idx)
  {
    outfile.Printf("%8i%8i%8i%8i", dih->A1()+1, dih->A2()+1, dih->A3()+1, dih->A4()+1);
    if ((idx % 2)==0) outfile.Printf("\n");
  }
  if ((idx % 2)==0) outfile.Printf("\n");
  outfile.Printf("\n");

  outfile.CloseFile();
  return 0;
}