Changeset 9dba5f


Ignore:
Timestamp:
Feb 24, 2011, 7:55:15 PM (14 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
873037
Parents:
fcac72
git-author:
Frederik Heber <heber@…> (02/10/11 22:01:26)
git-committer:
Frederik Heber <heber@…> (02/24/11 19:55:15)
Message:

PdbParser can now load and save multiple time steps.

  • WARNING: Molecules are not yet updated, it only works on simple test.
  • TEST: Parser/Pdb regression test on loading/storing multiple steps added.
  • PdbAtomInfoContainer:
    • private static knownDataKeys map from enum to strings, filled in cstor when empty by fillknownDataKeys() and cleared by befriended PdbParser's dstor.
    • operator << to easily print container.
    • getter function getDataKey.
  • PdbKey:
  • Pdbparser:
    • extended load() and added verbosity.
    • readAtomDataLine() can parse trajectories into present atoms.
    • extended save() and added verbosity.
    • new helper functions getAtomToParse() and readPdbAtomInfoContainer().
Files:
4 added
6 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified src/Parser/PdbAtomInfoContainer.cpp

    rfcac72 r9dba5f  
    2020#include "CodePatterns/MemDebug.hpp"
    2121
    22 //#include "CodePatterns/Assert.hpp"
     22#include "CodePatterns/Assert.hpp"
    2323//#include "CodePatterns/Log.hpp"
    2424#include "CodePatterns/toString.hpp"
    … …  
    2727#include "PdbAtomInfoContainer.hpp"
    2828
     29PdbAtomInfoContainer::knownDataKeysMap PdbAtomInfoContainer::knownDataKeys;
    2930
    3031PdbAtomInfoContainer::PdbAtomInfoContainer() :
    … …  
    4041  tempFactor(0.),
    4142  element(""),
    42   charge(0)
    43 {}
     43  charge(0),
     44  timestep(0)
     45{
     46  if (knownDataKeys.empty())
     47    fillknownDataKeys();
     48}
    4449
    4550PdbAtomInfoContainer::~PdbAtomInfoContainer()
    4651{}
    4752
     53void PdbAtomInfoContainer::fillknownDataKeys()
     54{
     55  knownDataKeys[PdbKey::token] = "token";
     56  knownDataKeys[PdbKey::serial] = "serial";
     57  knownDataKeys[PdbKey::name] = "name";
     58  knownDataKeys[PdbKey::altLoc] = "altLoc";
     59  knownDataKeys[PdbKey::resName] = "resName";
     60  knownDataKeys[PdbKey::chainID] = "chainID";
     61  knownDataKeys[PdbKey::resSeq] = "resSeq";
     62  knownDataKeys[PdbKey::iCode] = "iCode";
     63  knownDataKeys[PdbKey::X] = "X";
     64  knownDataKeys[PdbKey::Y] = "Y";
     65  knownDataKeys[PdbKey::Z] = "Z";
     66  knownDataKeys[PdbKey::occupancy] = "occupancy";
     67  knownDataKeys[PdbKey::tempFactor] = "tempFactor";
     68  knownDataKeys[PdbKey::element] = "element";
     69  knownDataKeys[PdbKey::charge] = "charge";
     70  knownDataKeys[PdbKey::timestep] = "timestep";
     71}
     72
    4873void PdbAtomInfoContainer::set(const PdbKey::PdbDataKey key, std::string value)
    4974{
    … …  
    93118    case PdbKey::charge :
    94119      ScanKey(charge, value);
     120      break;
     121    case PdbKey::timestep :
     122      ScanKey(timestep, value);
    95123      break;
    96124    default :
    … …  
    134162    case PdbKey::charge :
    135163      return toString(charge);
     164    case PdbKey::timestep :
     165      return toString(timestep);
    136166    default :
    137167      std::cout << "Unknown key: " << key << std::endl;
    … …  
    150180    case PdbKey::charge :
    151181      return charge;
     182    case PdbKey::timestep :
     183      return timestep;
    152184    default :
    153185      std::cout << "Unknown key or not presentable as int: " << key << std::endl;
    … …  
    175207  }
    176208}
     209
     210const std::string& PdbAtomInfoContainer::getDataKey(PdbKey::PdbDataKey key) const
     211{
     212  knownDataKeysMap::const_iterator iter;
     213  iter = knownDataKeys.find(key);
     214  ASSERT(iter != knownDataKeys.end(),
     215      "PdbAtomInfoContainer::getDataKey() - unknown key "
     216      +toString(key)+" request, did someone add a key to PdbKey?");
     217  return iter->second;
     218}
     219
     220std::ostream& operator<<(std::ostream &ost, const PdbAtomInfoContainer& m)
     221{
     222  ost << m.getDataKey(PdbKey::token) << "(" << m.get<std::string>(PdbKey::token) <<") ";
     223  ost << m.getDataKey(PdbKey::serial) << "(" << m.get<std::string>(PdbKey::serial) <<") ";
     224  ost << m.getDataKey(PdbKey::name) << "(" << m.get<std::string>(PdbKey::name) <<") ";
     225  ost << m.getDataKey(PdbKey::altLoc) << "(" << m.get<std::string>(PdbKey::altLoc) <<") ";
     226  ost << m.getDataKey(PdbKey::resName) << "(" << m.get<std::string>(PdbKey::resName) <<") ";
     227  ost << m.getDataKey(PdbKey::chainID) << "(" << m.get<std::string>(PdbKey::chainID) <<") ";
     228  ost << m.getDataKey(PdbKey::resSeq) << "(" << m.get<std::string>(PdbKey::resSeq) <<") ";
     229  ost << m.getDataKey(PdbKey::iCode) << "(" << m.get<std::string>(PdbKey::iCode) <<") ";
     230  ost << m.getDataKey(PdbKey::X) << "(" << m.get<std::string>(PdbKey::X) <<") ";
     231  ost << m.getDataKey(PdbKey::Y) << "(" << m.get<std::string>(PdbKey::Y) <<") ";
     232  ost << m.getDataKey(PdbKey::Z) << "(" << m.get<std::string>(PdbKey::Z) <<") ";
     233  ost << m.getDataKey(PdbKey::occupancy) << "(" << m.get<std::string>(PdbKey::occupancy) <<") ";
     234  ost << m.getDataKey(PdbKey::tempFactor) << "(" << m.get<std::string>(PdbKey::tempFactor) <<") ";
     235  ost << m.getDataKey(PdbKey::element) << "(" << m.get<std::string>(PdbKey::element) <<") ";
     236  ost << m.getDataKey(PdbKey::charge) << "(" << m.get<std::string>(PdbKey::charge) <<") ";
     237  ost << m.getDataKey(PdbKey::timestep) << "(" << m.get<std::string>(PdbKey::timestep) <<") ";
     238  return ost;
     239}

  • TabularUnified src/Parser/PdbAtomInfoContainer.hpp

    rfcac72 r9dba5f  
    1717#include "PdbKey.hpp"
    1818
     19class PdbParser;
    1920class Vector;
    2021
    … …  
    2627 */
    2728class PdbAtomInfoContainer {
     29  friend class PdbParser;
    2830public:
    2931  PdbAtomInfoContainer();
    … …  
    6163  }
    6264
     65  const std::string& getDataKey(PdbKey::PdbDataKey key) const;
     66
    6367private:
     68
     69  void fillknownDataKeys();
     70
     71  typedef std::map<PdbKey::PdbDataKey, std::string> knownDataKeysMap;
     72  static knownDataKeysMap knownDataKeys;
     73
    6474  std::string token;
    6575  int serial;
    … …  
    7585  std::string element;
    7686  int charge;
     87  size_t timestep;
    7788};
     89
     90std::ostream& operator<<(std::ostream &ost, const PdbAtomInfoContainer& m);
    7891
    7992template <>

  • TabularUnified src/Parser/PdbKey.hpp

    rfcac72 r9dba5f  
    2020class PdbKey {
    2121public:
    22   enum KnownTokens { Atom, Remark, Connect, Filler, EndOfFile, NoToken };
     22  enum KnownTokens { Atom, Remark, Connect, Filler, EndOfTimestep, NoToken };
    2323
    2424  enum PdbDataKey {
    … …  
    3838    tempFactor,
    3939    element,
    40     charge
     40    charge,
     41    timestep
    4142  };
    4243};

  • TabularUnified src/Parser/PdbParser.cpp

    rfcac72 r9dba5f  
    2424#include "CodePatterns/toString.hpp"
    2525#include "CodePatterns/Verbose.hpp"
     26#include "Descriptors/AtomIdDescriptor.hpp"
     27#include "World.hpp"
    2628#include "atom.hpp"
    2729#include "bond.hpp"
    … …  
    4951  knownTokens["HETATM"] = PdbKey::Atom;
    5052  knownTokens["TER"] = PdbKey::Filler;
    51   knownTokens["END"] = PdbKey::EndOfFile;
     53  knownTokens["END"] = PdbKey::EndOfTimestep;
    5254  knownTokens["CONECT"] = PdbKey::Connect;
    5355  knownTokens["REMARK"] = PdbKey::Remark;
    54   knownTokens[""] = PdbKey::EndOfFile;
     56  knownTokens[""] = PdbKey::EndOfTimestep;
    5557
    5658  // argh, why can't just PdbKey::X+(size_t)i
    … …  
    6466 */
    6567PdbParser::~PdbParser() {
     68  PdbAtomInfoContainer::knownDataKeys.clear();
    6669  additionalAtomData.clear();
    6770  atomIdMap.clear();
    … …  
    110113  atomIdMap.clear();
    111114
     115  bool NotEndOfFile = true;
    112116  molecule *newmol = World::getInstance().createMolecule();
    113117  newmol->ActiveFlag = true;
    114   bool NotEndOfFile = true;
    115118  // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
    116119  World::getInstance().getMolecules()->insert(newmol);
    117120  while (NotEndOfFile) {
    118     std::getline(*file, line, '\n');
    119     // extract first token
    120     token = getToken(line);
    121     //DoLog(1) && (Log() << Verbose(1) << " Recognized token of type : " << token << std::endl);
    122     switch (token) {
    123       case PdbKey::Atom:
    124         readAtomDataLine(line, newmol);
    125         break;
    126       case PdbKey::Remark:
    127         break;
    128       case PdbKey::Connect:
    129         readNeighbors(line);
    130         break;
    131       case PdbKey::Filler:
    132         break;
    133       case PdbKey::EndOfFile:
    134         NotEndOfFile = false;
    135         break;
    136       default:
    137         // TODO: put a throw here
    138         DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "'" << std::endl);
    139         //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+".");
    140         break;
     121    bool NotEndOfTimestep = true;
     122    while (NotEndOfTimestep) {
     123      std::getline(*file, line, '\n');
     124      // extract first token
     125      token = getToken(line);
     126      switch (token) {
     127        case PdbKey::Atom:
     128          LOG(3,"INFO: Parsing ATOM entry.");
     129          readAtomDataLine(line, newmol);
     130          break;
     131        case PdbKey::Remark:
     132          LOG(3,"INFO: Parsing REM entry.");
     133          break;
     134        case PdbKey::Connect:
     135          LOG(3,"INFO: Parsing CONECT entry.");
     136          readNeighbors(line);
     137          break;
     138        case PdbKey::Filler:
     139          LOG(3,"INFO: Stumbled upon Filler entry.");
     140          break;
     141        case PdbKey::EndOfTimestep:
     142          LOG(3,"INFO: Parsing END entry or empty line.");
     143          NotEndOfTimestep = false;
     144          break;
     145        default:
     146          // TODO: put a throw here
     147          DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "'" << std::endl);
     148          //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+".");
     149          break;
     150      }
     151      NotEndOfFile = NotEndOfFile && (file->good());
     152      linecount++;
    141153    }
    142     NotEndOfFile = NotEndOfFile && (file->good());
    143     linecount++;
    144154  }
    145155}
    … …  
    154164{
    155165  DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl);
    156   {
    157     // add initial remark
    158     *file << "REMARK created by molecuilder on ";
    159     time_t now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    160     // ctime ends in \n\0, we have to cut away the newline
    161     std::string time(ctime(&now));
    162     size_t pos = time.find('\n');
    163     if (pos != 0)
    164       *file << time.substr(0,pos);
    165     else
    166       *file << time;
    167     *file << endl;
    168   }
    169 
    170   // we distribute serials, hence clear map beforehand
     166
     167  // check for maximum number of time steps
     168  size_t max_timesteps = 0;
     169  BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms()) {
     170    if (_atom->getTrajectorySize() > max_timesteps)
     171      max_timesteps = _atom->getTrajectorySize();
     172  }
     173  LOG(2,"INFO: Found a maximum of " << max_timesteps << " time steps to store.");
     174
     175  // re-distribute serials
     176  // (new atoms might have been added)
     177  // (serials must be consistent over time steps)
    171178  atomIdMap.clear();
    172   {
    173     std::map<size_t,size_t> MolIdMap;
    174     size_t MolNo = 1;  // residue number starts at 1 in pdb
    175     for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
    176       const molecule *mol = (*atomIt)->getMolecule();
    177       if ((mol != NULL) && (MolIdMap.find(mol->getId()) == MolIdMap.end())) {
    178         MolIdMap[mol->getId()] = MolNo++;
     179  int AtomNo = 1; // serial number starts at 1 in pdb
     180  for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
     181    PdbAtomInfoContainer &atomInfo = getadditionalAtomData(*atomIt);
     182    setSerial((*atomIt)->getId(), AtomNo);
     183    atomInfo.set(PdbKey::serial, toString(AtomNo));
     184    AtomNo++;
     185  }
     186
     187  // store all time steps
     188  for (size_t step = 0; step < max_timesteps; ++step) {
     189    {
     190      // add initial remark
     191      *file << "REMARK created by molecuilder on ";
     192      time_t now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     193      // ctime ends in \n\0, we have to cut away the newline
     194      std::string time(ctime(&now));
     195      size_t pos = time.find('\n');
     196      if (pos != 0)
     197        *file << time.substr(0,pos);
     198      else
     199        *file << time;
     200      *file << ", time step " << step;
     201      *file << endl;
     202    }
     203
     204    {
     205      std::map<size_t,size_t> MolIdMap;
     206      size_t MolNo = 1;  // residue number starts at 1 in pdb
     207      for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
     208        const molecule *mol = (*atomIt)->getMolecule();
     209        if ((mol != NULL) && (MolIdMap.find(mol->getId()) == MolIdMap.end())) {
     210          MolIdMap[mol->getId()] = MolNo++;
     211        }
     212      }
     213      const size_t MaxMol = MolNo;
     214
     215      // have a count per element and per molecule (0 is for all homeless atoms)
     216      std::vector<int> **elementNo = new std::vector<int>*[MaxMol];
     217      for (size_t i = 0; i < MaxMol; ++i)
     218        elementNo[i] = new std::vector<int>(MAX_ELEMENTS,1);
     219      char name[MAXSTRINGSIZE];
     220      std::string ResidueName;
     221
     222      // write ATOMs
     223      for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
     224        PdbAtomInfoContainer &atomInfo = getadditionalAtomData(*atomIt);
     225        // gather info about residue
     226        const molecule *mol = (*atomIt)->getMolecule();
     227        if (mol == NULL) {
     228          MolNo = 0;
     229          atomInfo.set(PdbKey::resSeq, "0");
     230        } else {
     231          ASSERT(MolIdMap.find(mol->getId()) != MolIdMap.end(),
     232              "PdbParser::save() - Mol id "+toString(mol->getId())+" not present despite we set it?!");
     233          MolNo = MolIdMap[mol->getId()];
     234          atomInfo.set(PdbKey::resSeq, toString(MolIdMap[mol->getId()]));
     235          if (atomInfo.get<std::string>(PdbKey::resName) == "-")
     236            atomInfo.set(PdbKey::resName, mol->getName().substr(0,3));
     237        }
     238        // get info about atom
     239        const size_t  Z = (*atomIt)->getType()->getAtomicNumber();
     240        if (atomInfo.get<std::string>(PdbKey::name) == "-") {  // if no name set, give it a new name
     241          sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), (*elementNo[MolNo])[Z]);
     242          (*elementNo[MolNo])[Z] = ((*elementNo[MolNo])[Z]+1) % 100;   // confine to two digits
     243          atomInfo.set(PdbKey::name, name);
     244        }
     245        // set position
     246        for (size_t i=0; i<NDIM;++i) {
     247          stringstream position;
     248          position << setw(8) << fixed << setprecision(3) << (*atomIt)->getPositionAtStep(step).at(i);
     249          atomInfo.set(PositionEnumMap[i], position.str());
     250        }
     251        // change element and charge if changed
     252        if (atomInfo.get<std::string>(PdbKey::element) != (*atomIt)->getType()->getSymbol())
     253          atomInfo.set(PdbKey::element, (*atomIt)->getType()->getSymbol());
     254
     255        // finally save the line
     256        saveLine(file, atomInfo);
     257      }
     258      for (size_t i = 0; i < MaxMol; ++i)
     259        delete elementNo[i];
     260      delete elementNo;
     261
     262      // write CONECTs
     263      for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
     264        writeNeighbors(file, 4, *atomIt);
    179265      }
    180266    }
    181     const size_t MaxMol = MolNo;
    182 
    183     // have a count per element and per molecule (0 is for all homeless atoms)
    184     std::vector<int> **elementNo = new std::vector<int>*[MaxMol];
    185     for (size_t i = 0; i < MaxMol; ++i)
    186       elementNo[i] = new std::vector<int>(MAX_ELEMENTS,1);
    187     char name[MAXSTRINGSIZE];
    188     std::string ResidueName;
    189 
    190     // write ATOMs
    191     int AtomNo = 1; // serial number starts at 1 in pdb
    192     for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
    193       PdbAtomInfoContainer &atomInfo = getadditionalAtomData(*atomIt);
    194       // gather info about residue
    195       const molecule *mol = (*atomIt)->getMolecule();
    196       if (mol == NULL) {
    197         MolNo = 0;
    198         atomInfo.set(PdbKey::resSeq, "0");
    199       } else {
    200         ASSERT(MolIdMap.find(mol->getId()) != MolIdMap.end(),
    201             "PdbParser::save() - Mol id "+toString(mol->getId())+" not present despite we set it?!");
    202         MolNo = MolIdMap[mol->getId()];
    203         atomInfo.set(PdbKey::resSeq, toString(MolIdMap[mol->getId()]));
    204         if (atomInfo.get<std::string>(PdbKey::resName) == "-")
    205           atomInfo.set(PdbKey::resName, mol->getName().substr(0,3));
    206       }
    207       // get info about atom
    208       const size_t  Z = (*atomIt)->getType()->getAtomicNumber();
    209       if (atomInfo.get<std::string>(PdbKey::name) == "-") {  // if no name set, give it a new name
    210         sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), (*elementNo[MolNo])[Z]);
    211         (*elementNo[MolNo])[Z] = ((*elementNo[MolNo])[Z]+1) % 100;   // confine to two digits
    212         atomInfo.set(PdbKey::name, name);
    213       }
    214       // set position
    215       for (size_t i=0; i<NDIM;++i) {
    216         stringstream position;
    217         position << setw(8) << fixed << setprecision(3) << (*atomIt)->getPosition().at(i);
    218         atomInfo.set(PositionEnumMap[i], position.str());
    219       }
    220       // change element and charge if changed
    221       if (atomInfo.get<std::string>(PdbKey::element) != (*atomIt)->getType()->getSymbol())
    222         atomInfo.set(PdbKey::element, (*atomIt)->getType()->getSymbol());
    223       setSerial((*atomIt)->getId(), AtomNo);
    224       atomInfo.set(PdbKey::serial, toString(AtomNo));
    225 
    226       // finally save the line
    227       saveLine(file, atomInfo);
    228       AtomNo++;
    229     }
    230     for (size_t i = 0; i < MaxMol; ++i)
    231       delete elementNo[i];
    232     delete elementNo;
    233 
    234     // write CONECTs
    235     for (vector<atom *>::const_iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
    236       writeNeighbors(file, 4, *atomIt);
    237     }
    238   }
    239 
    240   // END
    241   *file << "END" << endl;
    242 }
     267    // END
     268    *file << "END" << endl;
     269  }
     270
     271}
     272
     273/** Checks whether there is an entry for the given atom's \a _id.
     274 *
     275 * @param _id atom's id we wish to check on
     276 * @return true - entry present, false - only for atom's father or no entry
     277 */
     278bool PdbParser::isPresentadditionalAtomData(unsigned int _id)
     279{
     280  return (additionalAtomData.find(_id) != additionalAtomData.end());
     281}
     282
    243283
    244284/** Either returns reference to present entry or creates new with default values.
    … …  
    376416}
    377417
     418/** Either returns present atom with given id or a newly created one.
     419 *
     420 * @param id_string
     421 * @return
     422 */
     423atom* PdbParser::getAtomToParse(std::string id_string) const
     424{
     425  // get the local ID
     426  ConvertTo<int> toInt;
     427  unsigned int AtomID = toInt(id_string);
     428  LOG(4, "INFO: Local id is "+toString(AtomID)+".");
     429  // get the atomic ID if present
     430  atom* newAtom = NULL;
     431  if (atomIdMap.count((size_t)AtomID)) {
     432    std::map<size_t, size_t>::const_iterator iter = atomIdMap.find(AtomID);
     433    AtomID = iter->second;
     434    LOG(4, "INFO: Global id present as " << AtomID << ".");
     435    // check if atom exists
     436    newAtom = World::getInstance().getAtom(AtomById(AtomID));
     437    LOG(5, "INFO: Listing all present atoms with id.");
     438    BOOST_FOREACH(atom *_atom, World::getInstance().getAllAtoms())
     439      LOG(5, "INFO: " << *_atom << " with id " << _atom->getId());
     440  }
     441  // if not exists, create
     442  if (newAtom == NULL) {
     443    newAtom = World::getInstance().createAtom();
     444    LOG(4, "INFO: No association to global id present, creating atom.");
     445  } else {
     446    LOG(4, "INFO: Existing atom found: " << *newAtom << ".");
     447  }
     448  return newAtom;
     449}
     450
     451void PdbParser::readPdbAtomInfoContainer(PdbAtomInfoContainer &atomInfo, std::string &line) const
     452{
     453  LOG(4,"INFO: Parsing token from "+line.substr(0,6)+".");
     454  atomInfo.set(PdbKey::token, line.substr(0,6));
     455  LOG(4,"INFO: Parsing serial from "+line.substr(6,5)+".");
     456  atomInfo.set(PdbKey::serial, line.substr(6,5));
     457
     458  LOG(4,"INFO: Parsing name from "+line.substr(12,4)+".");
     459  atomInfo.set(PdbKey::name, line.substr(12,4));
     460  LOG(4,"INFO: Parsing altLoc from "+line.substr(16,1)+".");
     461  atomInfo.set(PdbKey::altLoc, line.substr(16,1));
     462  LOG(4,"INFO: Parsing resName from "+line.substr(17,3)+".");
     463  atomInfo.set(PdbKey::resName, line.substr(17,3));
     464  LOG(4,"INFO: Parsing chainID from "+line.substr(21,1)+".");
     465  atomInfo.set(PdbKey::chainID, line.substr(21,1));
     466  LOG(4,"INFO: Parsing resSeq from "+line.substr(22,4)+".");
     467  atomInfo.set(PdbKey::resSeq, line.substr(22,4));
     468  LOG(4,"INFO: Parsing iCode from "+line.substr(26,1)+".");
     469  atomInfo.set(PdbKey::iCode, line.substr(26,1));
     470
     471  LOG(4,"INFO: Parsing occupancy from "+line.substr(54,6)+".");
     472  atomInfo.set(PdbKey::occupancy, line.substr(54,6));
     473  LOG(4,"INFO: Parsing tempFactor from "+line.substr(60,6)+".");
     474  atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
     475  LOG(4,"INFO: Parsing charge from "+line.substr(78,2)+".");
     476  atomInfo.set(PdbKey::charge, line.substr(78,2));
     477  LOG(4,"INFO: Parsing element from "+line.substr(76,2)+".");
     478  atomInfo.set(PdbKey::element, line.substr(76,2));
     479}
     480
    378481/** Parse an ATOM line from a PDB file.
    379482 *
    … …  
    389492void PdbParser::readAtomDataLine(std::string &line, molecule *newmol = NULL) {
    390493  vector<string>::iterator it;
    391   stringstream lineStream;
    392   atom* newAtom = World::getInstance().createAtom();
     494
     495  atom* newAtom = getAtomToParse(line.substr(6,5));
     496  LOG(3,"INFO: Parsing END entry or empty line.");
     497  bool FirstTimestep = isPresentadditionalAtomData(newAtom->getId()) ? false : true;
     498  if (FirstTimestep) {
     499    LOG(3,"INFO: Parsing new atom.");
     500  } else {
     501    LOG(3,"INFO: Parsing present atom "+toString(*newAtom)+".");
     502  }
    393503  PdbAtomInfoContainer &atomInfo = getadditionalAtomData(newAtom);
     504  LOG(4,"INFO: Information in container is "+toString(atomInfo)+".");
     505
    394506  string word;
    395507  ConvertTo<size_t> toSize_t;
    396   double tmp;
    397 
    398   lineStream << line;
    399   atomInfo.set(PdbKey::token, line.substr(0,6));
    400   atomInfo.set(PdbKey::serial, line.substr(6,5));
    401   std::pair< std::set<size_t>::const_iterator, bool> Inserter =
    402     SerialSet.insert(toSize_t(atomInfo.get<std::string>(PdbKey::serial)));
    403   ASSERT(Inserter.second,
    404       "PdbParser::readAtomDataLine() - ATOM contains entry with serial "
    405       +atomInfo.get<std::string>(PdbKey::serial)+" already present!");
    406   // assign hightest+1 instead, but then beware of CONECT entries! Another map needed!
     508
     509  // assign highest+1 instead, but then beware of CONECT entries! Another map needed!
    407510//  if (!Inserter.second) {
    408511//    const size_t id = (*SerialSet.rbegin())+1;
    … …  
    432535//      << line.substr(78,2) << std::endl);
    433536
    434   setSerial(toSize_t(atomInfo.get<std::string>(PdbKey::serial)), newAtom->getId());
    435   atomInfo.set(PdbKey::name, line.substr(12,4));
    436   atomInfo.set(PdbKey::altLoc, line.substr(16,1));
    437   atomInfo.set(PdbKey::resName, line.substr(17,3));
    438   atomInfo.set(PdbKey::chainID, line.substr(21,1));
    439   atomInfo.set(PdbKey::resSeq, line.substr(22,4));
    440   atomInfo.set(PdbKey::iCode, line.substr(26,1));
    441   PdbAtomInfoContainer::ScanKey(tmp, line.substr(30,8));
    442   newAtom->set(0, tmp);
    443   PdbAtomInfoContainer::ScanKey(tmp, line.substr(38,8));
    444   newAtom->set(1, tmp);
    445   PdbAtomInfoContainer::ScanKey(tmp, line.substr(46,8));
    446   newAtom->set(2, tmp);
    447   atomInfo.set(PdbKey::occupancy, line.substr(54,6));
    448   atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
    449   atomInfo.set(PdbKey::charge, line.substr(78,2));
    450   atomInfo.set(PdbKey::element, line.substr(76,2));
    451   const element *elem = World::getInstance().getPeriode()
    452       ->FindElement(atomInfo.get<std::string>(PdbKey::element));
    453   ASSERT(elem != NULL,
    454       "PdbParser::readAtomDataLine() - element "+atomInfo.get<std::string>(PdbKey::element)+" is unknown!");
    455   newAtom->setType(elem);
    456 
    457   if (newmol != NULL)
    458     newmol->AddAtom(newAtom);
     537  if (FirstTimestep) {
     538    // first time step
     539    // then fill info container
     540    readPdbAtomInfoContainer(atomInfo, line);
     541    // set the serial
     542    std::pair< std::set<size_t>::const_iterator, bool> Inserter =
     543      SerialSet.insert(toSize_t(atomInfo.get<std::string>(PdbKey::serial)));
     544    ASSERT(Inserter.second,
     545        "PdbParser::readAtomDataLine() - ATOM contains entry with serial "
     546        +atomInfo.get<std::string>(PdbKey::serial)+" already present!");
     547    setSerial(toSize_t(atomInfo.get<std::string>(PdbKey::serial)), newAtom->getId());
     548    // set position
     549    Vector tempVector;
     550    LOG(4,"INFO: Parsing position from ("
     551        +line.substr(30,8)+","
     552        +line.substr(38,8)+","
     553        +line.substr(46,8)+").");
     554    PdbAtomInfoContainer::ScanKey(tempVector[0], line.substr(30,8));
     555    PdbAtomInfoContainer::ScanKey(tempVector[1], line.substr(38,8));
     556    PdbAtomInfoContainer::ScanKey(tempVector[2], line.substr(46,8));
     557    newAtom->setPosition(tempVector);
     558    // set element
     559    const element *elem = World::getInstance().getPeriode()
     560        ->FindElement(atomInfo.get<std::string>(PdbKey::element));
     561    ASSERT(elem != NULL,
     562        "PdbParser::readAtomDataLine() - element "+atomInfo.get<std::string>(PdbKey::element)+" is unknown!");
     563    newAtom->setType(elem);
     564
     565    if (newmol != NULL)
     566      newmol->AddAtom(newAtom);
     567  } else {
     568    // not first time step
     569    // then parse into different container
     570    PdbAtomInfoContainer consistencyInfo;
     571    readPdbAtomInfoContainer(consistencyInfo, line);
     572    // then check additional info for consistency
     573    ASSERT(atomInfo.get<std::string>(PdbKey::token) == consistencyInfo.get<std::string>(PdbKey::token),
     574        "PdbParser::readAtomDataLine() - difference in token on multiple time step for atom with id "
     575        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     576    ASSERT(atomInfo.get<std::string>(PdbKey::name) == consistencyInfo.get<std::string>(PdbKey::name),
     577        "PdbParser::readAtomDataLine() - difference in name on multiple time step for atom with id "
     578        +atomInfo.get<std::string>(PdbKey::serial)+":"
     579        +atomInfo.get<std::string>(PdbKey::name)+"!="
     580        +consistencyInfo.get<std::string>(PdbKey::name)+".");
     581    ASSERT(atomInfo.get<std::string>(PdbKey::altLoc) == consistencyInfo.get<std::string>(PdbKey::altLoc),
     582        "PdbParser::readAtomDataLine() - difference in altLoc on multiple time step for atom with id "
     583        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     584    ASSERT(atomInfo.get<std::string>(PdbKey::resName) == consistencyInfo.get<std::string>(PdbKey::resName),
     585        "PdbParser::readAtomDataLine() - difference in resName on multiple time step for atom with id "
     586        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     587    ASSERT(atomInfo.get<std::string>(PdbKey::chainID) == consistencyInfo.get<std::string>(PdbKey::chainID),
     588        "PdbParser::readAtomDataLine() - difference in chainID on multiple time step for atom with id "
     589        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     590    ASSERT(atomInfo.get<std::string>(PdbKey::resSeq) == consistencyInfo.get<std::string>(PdbKey::resSeq),
     591        "PdbParser::readAtomDataLine() - difference in resSeq on multiple time step for atom with id "
     592        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     593    ASSERT(atomInfo.get<std::string>(PdbKey::iCode) == consistencyInfo.get<std::string>(PdbKey::iCode),
     594        "PdbParser::readAtomDataLine() - difference in iCode on multiple time step for atom with id "
     595        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     596    ASSERT(atomInfo.get<std::string>(PdbKey::occupancy) == consistencyInfo.get<std::string>(PdbKey::occupancy),
     597        "PdbParser::readAtomDataLine() - difference in occupancy on multiple time step for atom with id "
     598        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     599    ASSERT(atomInfo.get<std::string>(PdbKey::tempFactor) == consistencyInfo.get<std::string>(PdbKey::tempFactor),
     600        "PdbParser::readAtomDataLine() - difference in tempFactor on multiple time step for atom with id "
     601        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     602    ASSERT(atomInfo.get<std::string>(PdbKey::charge) == consistencyInfo.get<std::string>(PdbKey::charge),
     603        "PdbParser::readAtomDataLine() - difference in charge on multiple time step for atom with id "
     604        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     605    ASSERT(atomInfo.get<std::string>(PdbKey::element) == consistencyInfo.get<std::string>(PdbKey::element),
     606        "PdbParser::readAtomDataLine() - difference in element on multiple time step for atom with id "
     607        +atomInfo.get<std::string>(PdbKey::serial)+"!");
     608    // and parse in trajectory
     609    Vector tempVector;
     610    LOG(4,"INFO: Parsing trajectory position from ("
     611        +line.substr(30,8)+","
     612        +line.substr(38,8)+","
     613        +line.substr(46,8)+").");
     614    PdbAtomInfoContainer::ScanKey(tempVector[0], line.substr(30,8));
     615    PdbAtomInfoContainer::ScanKey(tempVector[1], line.substr(38,8));
     616    PdbAtomInfoContainer::ScanKey(tempVector[2], line.substr(46,8));
     617    // increase time step
     618    size_t timestep = atomInfo.get<size_t>(PdbKey::timestep)+1;
     619    atomInfo.set(PdbKey::timestep, toString(timestep));
     620    LOG(4,"INFO: Adding time step "+atomInfo.get<std::string>(PdbKey::timestep)+".");
     621    // and set position at new time step
     622    newAtom->setPositionAtStep(timestep, tempVector);
     623  }
     624
    459625
    460626//  printAtomInfo(newAtom);

  • TabularUnified src/Parser/PdbParser.hpp

    rfcac72 r9dba5f  
    4747
    4848  // internal getter and setter
     49  bool isPresentadditionalAtomData(unsigned int _id);
    4950  PdbAtomInfoContainer& getadditionalAtomData(atom *_atom);
    5051  size_t getSerial(const size_t atomid) const;
    5152  void setSerial(const size_t localatomid, const size_t atomid);
     53
     54  // internal helper functions
     55  atom* getAtomToParse(std::string id_string) const;
     56  void readPdbAtomInfoContainer(PdbAtomInfoContainer &atomInfo, std::string &line) const;
    5257
    5358  /**
    … …  
    7984
    8085  /**
    81    * Maps original atom IDs received from the parsed file to atom IDs in the
    82    * world.
     86   * Ascertaining uniqueness of Ids.
    8387   */
    8488  std::set<size_t> SerialSet;

  • TabularUnified tests/regression/Parser/Pdb/testsuite-parser-pdb.at

    rfcac72 r9dba5f  
    77AT_CLEANUP
    88
     9AT_SETUP([Parser - loading pdb file with multiple time steps])
     10AT_KEYWORDS([parser,pdb])
     11AT_CHECK([../../molecuilder -v 4 -i testmulti.pdb -o pdb -l ${abs_top_srcdir}/${AUTOTEST_PATH}/Parser/Pdb/pre/testmulti.pdb], 0, [ignore], [ignore])
     12AT_CHECK([file=testmulti.pdb; diff -I '.*created by molecuilder.*' $file ${abs_top_srcdir}/${AUTOTEST_PATH}/Parser/Pdb/post/$file], 0, [ignore], [ignore])
     13AT_CHECK([../../molecuilder -v 4 -i testmulti2.pdb -o pdb -l ${abs_top_srcdir}/${AUTOTEST_PATH}/Parser/Pdb/pre/testmulti2.pdb], 0, [ignore], [ignore])
     14AT_CHECK([file=testmulti2.pdb; diff -I '.*created by molecuilder.*' $file ${abs_top_srcdir}/${AUTOTEST_PATH}/Parser/Pdb/post/$file], 0, [ignore], [ignore])
     15AT_CLEANUP
     16
    917AT_SETUP([Parser - saving pdb file])
    1018AT_KEYWORDS([parser,pdb])
Note: See TracChangeset for help on using the changeset viewer.