source: src/Parser/PdbParser.cpp@ 16462f

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * PdbParser.cpp
 *
 *  Created on: Aug 17, 2010
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "Helpers/MemDebug.hpp"

#include "Helpers/Assert.hpp"
#include "Helpers/Log.hpp"
#include "Helpers/toString.hpp"
#include "Helpers/Verbose.hpp"
#include "World.hpp"
#include "atom.hpp"
#include "bond.hpp"
#include "element.hpp"
#include "molecule.hpp"
#include "periodentafel.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Parser/PdbParser.hpp"

#include <map>
#include <vector>

#include <iostream>
#include <iomanip>

using namespace std;

/**
 * Constructor.
 */
PdbParser::PdbParser() {
  knownTokens["ATOM"] = PdbKey::Atom;
  knownTokens["HETATM"] = PdbKey::Atom;
  knownTokens["TER"] = PdbKey::Filler;
  knownTokens["END"] = PdbKey::EndOfFile;
  knownTokens["CONECT"] = PdbKey::Connect;
  knownTokens["REMARK"] = PdbKey::Remark;
  knownTokens[""] = PdbKey::EndOfFile;

  // argh, why can't just PdbKey::X+(size_t)i
  PositionEnumMap[0] = PdbKey::X;
  PositionEnumMap[1] = PdbKey::Y;
  PositionEnumMap[2] = PdbKey::Z;
}

/**
 * Destructor.
 */
PdbParser::~PdbParser() {
  additionalAtomData.clear();
  atomIdMap.clear();
}


/** Parses the initial word of the given \a line and returns the token type.
 *
 * @param line line to scan
 * @return token type
 */
enum PdbKey::KnownTokens PdbParser::getToken(string &line)
{
  // look for first space
  const size_t space_location = line.find(' ');
  const size_t tab_location = line.find('\t');
  size_t location = space_location < tab_location ? space_location : tab_location;
  string token;
  if (location != string::npos) {
    //DoLog(1) && (Log() << Verbose(1) << "Found space at position " << space_location << std::endl);
    token = line.substr(0,space_location);
  } else {
    token = line;
  }

  //DoLog(1) && (Log() << Verbose(1) << "Token is " << token << std::endl);
  if (knownTokens.count(token) == 0)
    return PdbKey::NoToken;
  else
    return knownTokens[token];

  return PdbKey::NoToken;
}

/**
 * Loads atoms from a PDB-formatted file.
 *
 * \param PDB file
 */
void PdbParser::load(istream* file) {
  string line;
  size_t linecount  = 0;
  enum PdbKey::KnownTokens token;

  // reset atomIdMap for this file (to correctly parse CONECT entries)
  atomIdMap.clear();

  molecule *newmol = World::getInstance().createMolecule();
  newmol->ActiveFlag = true;
  bool NotEndOfFile = true;
  // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
  World::getInstance().getMolecules()->insert(newmol);
  while (NotEndOfFile) {
    std::getline(*file, line, '\n');
    // extract first token
    token = getToken(line);
    //DoLog(1) && (Log() << Verbose(1) << " Recognized token of type : " << token << std::endl);
    switch (token) {
      case PdbKey::Atom:
        readAtomDataLine(line, newmol);
        break;
      case PdbKey::Remark:
        break;
      case PdbKey::Connect:
        readNeighbors(line);
        break;
      case PdbKey::Filler:
        break;
      case PdbKey::EndOfFile:
        NotEndOfFile = false;
        break;
      default:
        // TODO: put a throw here
        DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "'" << std::endl);
        //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+".");
        break;
    }
    NotEndOfFile = NotEndOfFile && (file->good());
    linecount++;
  }
}

/**
 * Saves the World's current state into as a PDB file.
 *
 * \param file where to save the state
 */
void PdbParser::save(ostream* file) {
  DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl);

  {
    // add initial remark
    *file << "REMARK created by molecuilder on ";
    time_t now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    // ctime ends in \n\0, we have to cut away the newline
    std::string time(ctime(&now));
    size_t pos = time.find('\n');
    if (pos != 0)
      *file << time.substr(0,pos);
    else
      *file << time;
    *file << endl;
  }

  // we distribute serials, hence clear map beforehand
  atomIdMap.clear();
  {
    std::vector<atom *> AtomList = World::getInstance().getAllAtoms();
    std::vector<molecule *> MolList = World::getInstance().getAllMolecules();
    std::map<size_t,size_t> MolIdMap;
    size_t MolNo = 1;  // residue number starts at 1 in pdb
    for (std::vector<molecule *>::const_iterator iter = MolList.begin();
        iter != MolList.end();
        ++iter) {
      MolIdMap[(*iter)->getId()] = MolNo++;
    }
    const size_t MaxMol = MolNo;

    // have a count per element and per molecule (0 is for all homeless atoms)
    std::vector<int> **elementNo = new std::vector<int>*[MaxMol];
    for (size_t i = 0; i < MaxMol; ++i)
      elementNo[i] = new std::vector<int>(MAX_ELEMENTS,1);
    char name[MAXSTRINGSIZE];
    std::string ResidueName;

    // write ATOMs
    int AtomNo = 1; // serial number starts at 1 in pdb
    for (vector<atom *>::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
      PdbAtomInfoContainer &atomInfo = additionalAtomData[(*atomIt)->getId()];
      // gather info about residue
      const molecule *mol = (*atomIt)->getMolecule();
      if (mol == NULL) {
        MolNo = 0;
        atomInfo.set(PdbKey::resSeq, "0");
      } else {
        ASSERT(MolIdMap.find(mol->getId()) != MolIdMap.end(),
            "PdbParser::save() - Mol id "+toString(mol->getId())+" not present despite we set it?!");
        MolNo = MolIdMap[mol->getId()];
        atomInfo.set(PdbKey::resSeq, toString(MolIdMap[mol->getId()]));
        if (atomInfo.get<std::string>(PdbKey::resName) == "-")
          atomInfo.set(PdbKey::resName, mol->getName().substr(0,3));
      }
      // get info about atom
      const size_t  Z = (*atomIt)->getType()->getAtomicNumber();
      if (atomInfo.get<std::string>(PdbKey::name) == "-") {  // if no name set, give it a new name
        sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), (*elementNo[MolNo])[Z]);
        (*elementNo[MolNo])[Z] = ((*elementNo[MolNo])[Z]+1) % 100;   // confine to two digits
        atomInfo.set(PdbKey::name, name);
      }
      // set position
      for (size_t i=0; i<NDIM;++i) {
        stringstream position;
        position << setw(8) << fixed << setprecision(3) << (*atomIt)->getPosition().at(i);
        atomInfo.set(PositionEnumMap[i], position.str());
      }
      // change element and charge if changed
      if (atomInfo.get<std::string>(PdbKey::element) != (*atomIt)->getType()->getSymbol())
        atomInfo.set(PdbKey::element, (*atomIt)->getType()->getSymbol());
      setSerial((*atomIt)->getId(), AtomNo);
      atomInfo.set(PdbKey::serial, toString(AtomNo));

      // finally save the line
      saveLine(file, atomInfo);
      AtomNo++;
    }
    for (size_t i = 0; i < MaxMol; ++i)
      delete elementNo[i];
    delete elementNo;

    // write CONECTs
    for (vector<atom *>::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) {
      writeNeighbors(file, 4, *atomIt);
    }
  }

  // END
  *file << "END" << endl;
}

/**
 * Writes one line of PDB-formatted data to the provided stream.
 *
 * \param stream where to write the line to
 * \param *currentAtom the atom of which information should be written
 * \param AtomNo serial number of atom
 * \param *name name of atom, i.e. H01
 * \param ResidueName Name of molecule
 * \param ResidueNo number of residue
 */
void PdbParser::saveLine(
    ostream* file,
    const PdbAtomInfoContainer &atomInfo)
{
  *file << setfill(' ') << left << setw(6)
      << atomInfo.get<std::string>(PdbKey::token);
  *file << setfill(' ') << right << setw(5)
      << atomInfo.get<int>(PdbKey::serial); /* atom serial number */
  *file << " "; /* char 12 is empty */
  *file << setfill(' ') << left << setw(4)
      << atomInfo.get<std::string>(PdbKey::name);  /* atom name */
  *file << setfill(' ') << left << setw(1)
      << atomInfo.get<std::string>(PdbKey::altLoc); /* alternate location/conformation */
  *file << setfill(' ') << left << setw(3)
      << atomInfo.get<std::string>(PdbKey::resName);  /* residue name */
  *file << " "; /* char 21 is empty */
  *file << setfill(' ') << left << setw(1)
      << atomInfo.get<std::string>(PdbKey::chainID); /* chain identifier */
  *file << setfill(' ') << left << setw(4)
      << atomInfo.get<int>(PdbKey::resSeq); /* residue sequence number */
  *file << setfill(' ') << left << setw(1)
      << atomInfo.get<std::string>(PdbKey::iCode); /* iCode */
  *file << "   "; /* char 28-30 are empty */
  // have the following operate on stringstreams such that format specifiers
  // only act on these
  for (size_t i=0;i<NDIM;++i) {
    stringstream position;
    position << fixed << setprecision(3) << showpoint
        << atomInfo.get<double>(PositionEnumMap[i]);
    *file << setfill(' ') << right << setw(8) << position.str();
  }
  {
    stringstream occupancy;
    occupancy << fixed << setprecision(2) << showpoint
        << atomInfo.get<double>(PdbKey::occupancy); /* occupancy */
    *file << setfill(' ') << right << setw(6) << occupancy.str();
  }
  {
    stringstream tempFactor;
    tempFactor << fixed << setprecision(2) << showpoint
        << atomInfo.get<double>(PdbKey::tempFactor); /* temperature factor */
    *file << setfill(' ') << right << setw(6) << tempFactor.str();
  }
  *file << "          "; /* char 68-76 are empty */
  *file << setfill(' ') << right << setw(2) << atomInfo.get<std::string>(PdbKey::element); /* element */
  *file << setfill(' ') << right << setw(2) << atomInfo.get<int>(PdbKey::charge); /* charge */

  *file << endl;
}

/**
 * Writes the neighbor information of one atom to the provided stream.
 *
 * \param *file  where to write neighbor information to
 * \param MaxnumberOfNeighbors of neighbors
 * \param *currentAtom to the atom of which to take the neighbor information
 */
void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) {
  int MaxNo = MaxnumberOfNeighbors;
  if (!currentAtom->ListOfBonds.empty()) {
    for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) {
      if (MaxNo >= MaxnumberOfNeighbors) {
        *file << "CONECT";
        *file << setw(5) << getSerial(currentAtom->getId());
        MaxNo = 0;
      }
      *file << setw(5) << getSerial((*currentBond)->GetOtherAtom(currentAtom)->getId());
      MaxNo++;
      if (MaxNo == MaxnumberOfNeighbors)
        *file << "\n";
    }
    if (MaxNo != MaxnumberOfNeighbors)
      *file << "\n";
  }
}

/** Retrieves a value from PdbParser::atomIdMap.
 * \param atomid key
 * \return value
 */
size_t PdbParser::getSerial(const size_t atomid) const
{
  ASSERT(atomIdMap.find(atomid) != atomIdMap.end(), "PdbParser::getAtomId: atomid not present in Map.");
  return (atomIdMap.find(atomid)->second);
}

/** Sets an entry in PdbParser::atomIdMap.
 * \param localatomid key
 * \param atomid value
 * \return true - key not present, false - value present
 */
void PdbParser::setSerial(const size_t localatomid, const size_t atomid)
{
  pair<std::map<size_t,size_t>::iterator, bool > inserter;
//  DoLog(1) && (Log() << Verbose(1) << "PdbParser::setAtomId() - Inserting ("
//      << localatomid << " -> " << atomid << ")." << std::endl);
  inserter = atomIdMap.insert( make_pair(localatomid, atomid) );
  ASSERT(inserter.second, "PdbParser::setAtomId: atomId already present in Map.");
}

/** Parse an ATOM line from a PDB file.
 *
 * Reads one data line of a pdstatus file and interprets it according to the
 * specifications of the PDB 3.2 format: http://www.wwpdb.org/docs.html
 *
 *  A new atom is created and filled with available information, non-
 *  standard information is placed in additionalAtomData at the atom's id.
 *
 * \param line to parse as an atom
 * \param newmol molecule to add parsed atoms to
 */
void PdbParser::readAtomDataLine(std::string &line, molecule *newmol = NULL) {
  vector<string>::iterator it;
  stringstream lineStream;
  atom* newAtom = World::getInstance().createAtom();
  additionalAtomData[newAtom->getId()] = *(new PdbAtomInfoContainer);
  PdbAtomInfoContainer &atomInfo = additionalAtomData[newAtom->getId()];
  string word;
  ConvertTo<size_t> toSize_t;
  double tmp;

  lineStream << line;
  atomInfo.set(PdbKey::token, line.substr(0,6));
  atomInfo.set(PdbKey::serial, line.substr(6,5));
  std::pair< std::set<size_t>::const_iterator, bool> Inserter =
    SerialSet.insert(toSize_t(atomInfo.get<std::string>(PdbKey::serial)));
  ASSERT(Inserter.second,
      "PdbParser::readAtomDataLine() - ATOM contains entry with serial "
      +atomInfo.get<std::string>(PdbKey::serial)+" already present!");
  // assign hightest+1 instead, but then beware of CONECT entries! Another map needed!
//  if (!Inserter.second) {
//    const size_t id = (*SerialSet.rbegin())+1;
//    SerialSet.insert(id);
//    atomInfo.set(PdbKey::serial, toString(id));
//    DoeLog(2) && (eLog() << Verbose(2)
//        << "Serial " << atomInfo.get<std::string>(PdbKey::serial) << " already present, "
//        << "assigning " << toString(id) << " instead." << std::endl);
//  }

  // check whether serial exists, if so, assign next available

//  DoLog(2) && (Log() << Verbose(2) << "Split line:"
//      << line.substr(6,5) << "|"
//      << line.substr(12,4) << "|"
//      << line.substr(16,1) << "|"
//      << line.substr(17,3) << "|"
//      << line.substr(21,1) << "|"
//      << line.substr(22,4) << "|"
//      << line.substr(26,1) << "|"
//      << line.substr(30,8) << "|"
//      << line.substr(38,8) << "|"
//      << line.substr(46,8) << "|"
//      << line.substr(54,6) << "|"
//      << line.substr(60,6) << "|"
//      << line.substr(76,2) << "|"
//      << line.substr(78,2) << std::endl);

  setSerial(toSize_t(atomInfo.get<std::string>(PdbKey::serial)), newAtom->getId());
  atomInfo.set(PdbKey::name, line.substr(12,4));
  atomInfo.set(PdbKey::altLoc, line.substr(16,1));
  atomInfo.set(PdbKey::resName, line.substr(17,3));
  atomInfo.set(PdbKey::chainID, line.substr(21,1));
  atomInfo.set(PdbKey::resSeq, line.substr(22,4));
  atomInfo.set(PdbKey::iCode, line.substr(26,1));
  PdbAtomInfoContainer::ScanKey(tmp, line.substr(30,8));
  newAtom->set(0, tmp);
  PdbAtomInfoContainer::ScanKey(tmp, line.substr(38,8));
  newAtom->set(1, tmp);
  PdbAtomInfoContainer::ScanKey(tmp, line.substr(46,8));
  newAtom->set(2, tmp);
  atomInfo.set(PdbKey::occupancy, line.substr(54,6));
  atomInfo.set(PdbKey::tempFactor, line.substr(60,6));
  atomInfo.set(PdbKey::charge, line.substr(78,2));
  atomInfo.set(PdbKey::element, line.substr(76,2));
  const element *elem = World::getInstance().getPeriode()
      ->FindElement(atomInfo.get<std::string>(PdbKey::element));
  ASSERT(elem != NULL,
      "PdbParser::readAtomDataLine() - element "+atomInfo.get<std::string>(PdbKey::element)+" is unknown!");
  newAtom->setType(elem);

  if (newmol != NULL)
    newmol->AddAtom(newAtom);

//  printAtomInfo(newAtom);
}

/** Prints all PDB-specific information known about an atom.
 *
 */
void PdbParser::printAtomInfo(const atom * const newAtom) const
{
  const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(newAtom->getId()); // operator[] const does not exist

  DoLog(1) && (Log() << Verbose(1) << "We know about atom " << newAtom->getId() << ":" << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\ttoken is " << atomInfo.get<std::string>(PdbKey::token) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tserial is " << atomInfo.get<int>(PdbKey::serial) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tname is " << atomInfo.get<std::string>(PdbKey::name) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\taltLoc is " << atomInfo.get<std::string>(PdbKey::altLoc) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tresName is " << atomInfo.get<std::string>(PdbKey::resName) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tchainID is " << atomInfo.get<std::string>(PdbKey::chainID) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tresSeq is " << atomInfo.get<int>(PdbKey::resSeq) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tiCode is " << atomInfo.get<std::string>(PdbKey::iCode) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tX is " << atomInfo.get<double>(PdbKey::X) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tY is " << atomInfo.get<double>(PdbKey::Y) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tZ is " << atomInfo.get<double>(PdbKey::Z) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\toccupancy is " << atomInfo.get<double>(PdbKey::occupancy) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\ttempFactor is " << atomInfo.get<double>(PdbKey::tempFactor) << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\telement is '" << *(newAtom->getType()) << "'" << std::endl);
  DoLog(1) && (Log() << Verbose(1) << "\tcharge is " << atomInfo.get<int>(PdbKey::charge) << std::endl);
}

/**
 * Reads neighbor information for one atom from the input.
 *
 * \param line to parse as an atom
 */
void PdbParser::readNeighbors(std::string &line)
{
  const size_t length = line.length();
  std::list<size_t> ListOfNeighbors;
  ConvertTo<size_t> toSize_t;

  // obtain neighbours
  // show split line for debugging
  string output;
  ASSERT(length >=16,
      "PdbParser::readNeighbors() - CONECT entry has not enough entries: "+line+"!");
//  output = "Split line:|";
//  output += line.substr(6,5) + "|";
  const size_t id = toSize_t(line.substr(6,5));
  for (size_t index = 11; index <= 26; index+=5) {
    if (index+5 <= length) {
//      output += line.substr(index,5) + "|";
      const size_t otherid = toSize_t(line.substr(index,5));
      ListOfNeighbors.push_back(otherid);
    } else  {
      break;
    }
  }
//  DoLog(2) && (Log() << Verbose(2) << output << std::endl);

  // add neighbours
  atom *_atom = World::getInstance().getAtom(AtomById(getSerial(id)));
  for (std::list<size_t>::const_iterator iter = ListOfNeighbors.begin();
      iter != ListOfNeighbors.end();
      ++iter) {
//    DoLog(1) && (Log() << Verbose(1) << "Adding Bond (" << getAtomId(id) << "," << getAtomId(*iter) << ")" << std::endl);
    atom * const _Otheratom = World::getInstance().getAtom(AtomById(getSerial(*iter)));
    _atom->addBond(_Otheratom);
  }
}

/**
 * Replaces atom IDs read from the file by the corresponding world IDs. All IDs
 * IDs of the input string will be replaced; expected separating characters are
 * "-" and ",".
 *
 * \param string in which atom IDs should be adapted
 *
 * \return input string with modified atom IDs
 */
//string PdbParser::adaptIdDependentDataString(string data) {
//  // there might be no IDs
//  if (data == "-") {
//    return "-";
//  }
//
//  char separator;
//  int id;
//  stringstream line, result;
//  line << data;
//
//  line >> id;
//  result << atomIdMap[id];
//  while (line.good()) {
//    line >> separator >> id;
//    result << separator << atomIdMap[id];
//  }
//
//  return result.str();
//  return "";
//}


bool PdbParser::operator==(const PdbParser& b) const
{
  bool status = true;
  World::AtomComposite atoms = World::getInstance().getAllAtoms();
  for (World::AtomComposite::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    if ((additionalAtomData.find((*iter)->getId()) != additionalAtomData.end())
        && (b.additionalAtomData.find((*iter)->getId()) != b.additionalAtomData.end())) {
      const PdbAtomInfoContainer &atomInfo = additionalAtomData.at((*iter)->getId());
      const PdbAtomInfoContainer &OtheratomInfo = b.additionalAtomData.at((*iter)->getId());

      status = status && (atomInfo.get<std::string>(PdbKey::serial) == OtheratomInfo.get<std::string>(PdbKey::serial));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in serials!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::name) == OtheratomInfo.get<std::string>(PdbKey::name));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in names!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::altLoc) == OtheratomInfo.get<std::string>(PdbKey::altLoc));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in altLocs!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::resName) == OtheratomInfo.get<std::string>(PdbKey::resName));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resNames!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::chainID) == OtheratomInfo.get<std::string>(PdbKey::chainID));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in chainIDs!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::resSeq) == OtheratomInfo.get<std::string>(PdbKey::resSeq));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resSeqs!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::iCode) == OtheratomInfo.get<std::string>(PdbKey::iCode));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in iCodes!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::occupancy) == OtheratomInfo.get<std::string>(PdbKey::occupancy));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in occupancies!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::tempFactor) == OtheratomInfo.get<std::string>(PdbKey::tempFactor));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in tempFactors!" << std::endl);
      status = status && (atomInfo.get<std::string>(PdbKey::charge) == OtheratomInfo.get<std::string>(PdbKey::charge));
      if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in charges!" << std::endl);
    }
  }

  return status;
}