Logo Search packages:      
Sourcecode: ball version File versions

HINFile.C

// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: HINFile.C,v 1.64.28.1 2007/08/06 15:57:39 oliver Exp $
//

#include <BALL/FORMAT/HINFile.h>
#include <BALL/CONCEPT/composite.h>
#include <BALL/KERNEL/residue.h>
#include <BALL/KERNEL/system.h>
#include <BALL/KERNEL/protein.h>
#include <BALL/KERNEL/atom.h>
#include <BALL/KERNEL/PDBAtom.h>
#include <BALL/KERNEL/bond.h>
#include <BALL/KERNEL/PTE.h>

#include <stack>

namespace BALL 
{

  struct HINFileBondStruct
  {
    Size        atom1;
    Size        atom2;
    Bond::Order order;
      };
 
00029       HINFile::HINFile()
            
            :     GenericMolFile(),
                  box_(0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
                  temperature_(0.0)
      {
      }

00037       HINFile::HINFile(const String& name, File::OpenMode open_mode)
            throw(Exception::FileNotFound)
            : GenericMolFile(),
                  box_(0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
                  temperature_(0.0)
      {
            open(name, open_mode);
      }

00046       HINFile::HINFile(const HINFile& file)
            throw(Exception::FileNotFound)
            : GenericMolFile(),
                  box_(0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
                  temperature_(0.0)
      {
            if (file.getName() != "")
            {
                  try
                  {     
                        open(file.getName());
                  }
                  catch (Exception::FileNotFound&)
                  {
                  }
            }
      }

00064       HINFile::~HINFile()
            
      {
      }
      
00069       const HINFile& HINFile::operator = (const HINFile& rhs)
            throw(Exception::FileNotFound)
      {
            box_ = SimpleBox3(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
            temperature_ = 0.0;

            GenericMolFile::operator = (rhs);
            
            return *this;
      }

      

      void HINFile::writeAtom_(const Atom& atom, Size number, Size atom_offset)
      {
            getFileStream() << "atom " << number + 1 - atom_offset << " ";
            String name = atom.getName();
            if (name != "") 
            {
                  // if the atom name contains blanks or some bullshit like that, 
                  // truncate it to the first field and complain about it
                  if (name.countFields() > 1)
                  {
                        getFileStream() << name.getField(0) <<  " ";
                        Log.warn() << "HINFile::write: truncated atom name '" << name << "' to '" << name.getField(0) << "'." << std::endl;
                  }
                  else
                  {
                        getFileStream() << name.trim() << " ";
                  }
            } 
            else 
            {
                  getFileStream() << "- ";
            }

            getFileStream() << atom.getElement().getSymbol().trim() << " ";
            if ((atom.getTypeName() == "?") || (atom.getTypeName() == ""))
            {
                  getFileStream() << "**";
            } 
            else 
            {
                  getFileStream() << atom.getTypeName();    
            }
            getFileStream() << " - ";
            getFileStream() << atom.getCharge() << " ";
            getFileStream() << atom.getPosition().x << " ";
            getFileStream() << atom.getPosition().y << " ";
            getFileStream() << atom.getPosition().z << " ";

            Size number_of_bonds = 0;
            String bond_string(" ");
            const Atom* bond_partner;

            // count the valid bonds (bonds to atoms inside the system)
            for (Position i = 0; i < atom.countBonds(); ++i) 
            {
                  const Bond* bond = atom.getBond(i);
                  bond_partner = bond->getPartner(atom);

                  Size index = bond_partner->getProperty("__HINFILE_INDEX").getUnsignedInt();
                  if (index != 0)
                  {
                        number_of_bonds++;

                        bond_string += String(index);

                        switch (bond->getOrder()) 
                        {
                              case Bond::ORDER__DOUBLE:           bond_string += " d "; break;
                              case Bond::ORDER__TRIPLE:           bond_string += " t "; break;
                              case Bond::ORDER__AROMATIC:   bond_string += " a "; break;
                              // default is single bond!
                              default:                                                          bond_string += " s ";
                        }
                  }
            }
            
            // write the bonds
            getFileStream() << number_of_bonds << bond_string << std::endl;

            // HyperChem uses A/ps, as does BALL. So, no conversion is needed.
            getFileStream() << "vel " << number + 1 - atom_offset << " " 
                                                 << atom.getVelocity().x << " " 
                                                 << atom.getVelocity().y << " "
                                                 << atom.getVelocity().z << std::endl;
      }

00158       bool HINFile::write(const Molecule& molecule)
            throw(File::CannotWrite)
      {
            System S;
            S.insert(*(Molecule*)molecule.create(true));
            return write(S);
      }
      
00166       bool HINFile::write(const System& system)
            throw(File::CannotWrite)
      {
            if (!isOpen() || getOpenMode() != std::ios::out)
            {
                  throw(File::CannotWrite(__FILE__, __LINE__, name_));
            }

            // the atom_vector contains the atoms in the order of
            // the atom iterator
            vector<const Atom*> atom_vector;

            // create a vector containing pointers to the atoms
            AtomConstIterator atom_it;          
            for (atom_it = system.beginAtom(); +atom_it; ++atom_it)
            {
                  atom_vector.push_back(&(*atom_it));
                  (const_cast<Atom&>(*atom_it)).setProperty("__HINFILE_INDEX", (unsigned int)atom_vector.size());
            }

            // the index_vector contains the index of the connected component
            // (HyperChem molecule) it is in and initialize it to zero
            vector<Index>           index_vector(atom_vector.size(), -1);

            typedef list<Size> Component;
            typedef     vector<Component> ComponentVector;

            // now calculate all connected components in the graph
            // formed by atoms and bonds of the system
            // each of these connected components represents a molecule
            // for the new HyperChem file
                  
            // index of the current connected component
            Size current_index = 0;

            // the index of the atom to start a new component
            Size start_index = 0;

            // search components until all atoms have been considered
            while (start_index < atom_vector.size())
            {
                  while ((start_index < atom_vector.size()) && (index_vector[start_index] >= 0))
                                     
                  {
                        start_index++;
                  }
                              
                  if (start_index < atom_vector.size())
                  {
                        // create a stack containing all atoms to be axamined for this component
                        std::stack<Size> atom_stack;
                        
                        // our start atom is the first to be considered and is marked, too
                        atom_stack.push(start_index);
                        index_vector[start_index] = (Index)current_index;

                        // never examine this atom again as start atom
                        start_index++;

                        // calculate the connected component
                        while (atom_stack.size() > 0)
                        {
                              // check all bonds of this atom
                              // and remove it from the stack
                              const Atom& current_atom = *atom_vector[atom_stack.top()];
                              atom_stack.pop();

                              Atom::BondConstIterator bond_it = current_atom.beginBond();
                              for (; +bond_it; ++bond_it) 
                              {
                                    // add the atom if it is not marked yet
                                    // ignore all bonds to atoms outside the system 
                                    // (these atoms have not been marked, so getProperty will return 0)
                                    Size atom_index = bond_it->getPartner(current_atom)
                                                                                                ->getProperty("__HINFILE_INDEX").getUnsignedInt();
                                    if ((atom_index != 0) && (index_vector[atom_index - 1] == -1))
                                    {
                                          // remember this atom in the stack
                                          // and mark it in the index_vector
                                          atom_stack.push(atom_index - 1);
                                          index_vector[atom_index - 1] = (Index)current_index;
                                    }
                              }
                        }
                        
                        // done with this component, increase the component counter
                        current_index++;
                  }
            }

            // create an empty vector for the connected components
            // 
            ComponentVector components(current_index);

            // now extract lists of all connected component 
            // in the order of the atom_vector to keep the
            // order of residues
            for (Size i = 0; i < atom_vector.size(); i++)
            {
                  // remember the index of the atom
                  components[index_vector[i]].push_back(i);
                  
                  // and set the atom's HINFILE_INDEX properly
                  // (i.e. to the index in the right connected component
                  (const_cast<Atom*>(atom_vector[i]))->setProperty("__HINFILE_INDEX", 
                                                                                                      (unsigned int)components[index_vector[i]].size());
            }

            // write some default header
            getFileStream() << "; HyperChem file created by BALL" << std::endl;
            getFileStream() << ";" << std::endl;
            getFileStream() << "forcefield AMBER" << std::endl;

            // ?????: 
            // insert system temperature here
            getFileStream() << "sys 0" << std::endl;

            // ?????:
            // insert the periodic box size (if any)
            
            Size atom_count = 0;
            Size atom_offset = 0;

            // for each connected component: write a molecule
            for (Size j = 0; j < current_index; j++)
            {
                  // renumber all atoms relative to the first atom of
                  // each connected component to obtain a numbering starting
                  // with 1 for each molecule
                  atom_offset = atom_count;

                  // try to find a name for the molecule
                  const Molecule* mol = atom_vector[components[j].front()]->getMolecule();
                  String name = "-";
                  if (mol != 0)
                  {
                        name = mol->getName();
                        name.trim();
                  }

                  // Write the molecule identifier    
                  // For recent versions of HyperChem, the name has to be enclosed in double quotes.
                  if (name != "")
                  {
                        // Make sure the name does not contain double quotes.
                        while (name.substitute("\"", "") != String::EndPos) {};
                        getFileStream() << "mol " << j + 1 << " \"" << name << "\"" << std::endl;
                  }
                  else
                  {
                        getFileStream() << "mol " << j + 1 << std::endl;
                  }
                  
                  // now iterate over all atoms and write them
                  const Residue* current_residue = 0;

                  // the residues start at zero in each molecule as do the atoms
                  Size res_count = 0;                 
                  Component::const_iterator comp_it = components[j].begin();
                  for (; comp_it != components[j].end(); comp_it++)
                  {
                        // counter for the residues
                        const Atom* this_atom = atom_vector[*comp_it];
                              
                        const Residue* this_residue = this_atom->Composite::getAncestor(RTTI::getDefault<Residue>());
                        if (this_residue != current_residue)
                        {
                              if (current_residue != 0)
                              {
                                    // write and "endres" tag
                                    getFileStream() << "endres " << res_count << std::endl;
                              }

                              current_residue = this_residue;

                              if (this_residue != 0)
                              {
                                    res_count++;

                                    // write a new residue tag
                                    name = this_residue->getName();
                                    name.trim();
                                    if (name == "")
                                    {
                                          name = "-";
                                    }
                                    
                                    getFileStream() << "res " << res_count << " ";
                                    getFileStream() << name << " ";
                                    
                                    name = this_residue->getID();
                                    name.trim();
                                    if (name == "")
                                    {
                                          // If the ID is not set, it defaults to the
                                          // current residue number.
                                          name = String(res_count - 1);
                                    }

                                    getFileStream() << name << " - ";
                                    
                                    // write the chain name
                                    const Chain*      chain = this_residue->getChain();
                                    if ((chain != 0)) 
                                    {
                                          name = chain->getName();
                                          name.trim();
                                          
                                          if (name == "")
                                          {
                                                name = "-";
                                          }
                                    } 
                                    else 
                                    {
                                          name = "-";
                                    }

                                    getFileStream() << "-" << std::endl;
                              }
                        }

                        // write the atom
                        writeAtom_(*this_atom, atom_count++, atom_offset);
                  }

                  // if the last atom was inside a residue, write the 
                  // endres tag
                  if (current_residue != 0)
                  {
                        // write and "endres" tag
                        getFileStream() << "endres " << res_count << std::endl;
                  }

                  // write endmol keyword
                  getFileStream() << "endmol " << j + 1 << std::endl;
            }

            // clear the atom properties
            for (atom_it = system.beginAtom(); +atom_it; ++atom_it)
            {
                  (const_cast<Atom&>(*atom_it)).clearProperty("__HINFILE_INDEX");
            }

            return true;
      }

00413       bool HINFile::read(System& system)
            throw(Exception::ParseError)
      {
            return GenericMolFile::read(system);
      }
      
00419       Molecule* HINFile::read()
            throw(Exception::ParseError)
      {
            if (!isValid())
            {
                  Log.error() << "trying to read form invalid HINFile '" << getName() << "'" << std::endl;
                  return false;
            }
            
            // we define some states for our simple parser machine
            // legal transitions are:
            // START->IN_MOLECULE, IN_MOLECULE->IN_RESIDUE
            // IN_RESIDUE->IN_MOLECULE, IN_MOLECULE->START
            enum State 
            {
                  START = 0,
                  IN_MOLECULE,
                  IN_RESIDUE
            };
            
            // define a macro to print an error message for the file (only once!)
#           define ERROR(msg)\
                        throw Exception::ParseError(__FILE__, __LINE__, getLine(), String("while reading line " + String(getLineNumber()) + " of '" + getName() + "': " + (msg)));
                  
            // the current state - we don't want to insert an atom if
            // we don't have created the molecule/residue to insert it into!
            State state = START;

            // 
            // All <mol>s that contain <res>idues are inserted
            // as single chains into this protein
    // 
            Residue*    residue = 0;
            Molecule*   molecule = 0;     
            Fragment*   fragment = 0;
            Protein*  protein = 0;
            Chain*    chain = 0;

            // initial size: 100 atoms, all set to NULL pointer, 100 bonds
            std::vector<Atom*> atom_vector(100);
            Position last_atom = 0;
            static std::vector<struct HINFileBondStruct> bond_vector;
            bond_vector.clear();

            String tag;

            try
            {
                  while (readLine() && (getLine() != "") && good()) 
                  {
                        // ignore comment lines
                        if (getLine()[0] == ';' || getLine() == "") 
                        {
                              continue;
                        }

                        // determine the hyperchem tag: always the first word in a line
                        tag = getLine().getField(0);

                        if (tag == "atom")
                        {
                              if (state == IN_RESIDUE || state == IN_MOLECULE) 
                              {
                                    Size number_of_fields = getLine().countFields();
                                    if (number_of_fields < 11) 
                                    {
                                          ERROR(String("illegal <atom> line: at least 10 arguments are required for <atom>!"))
                                    }

                                    Atom* atom;
                                    if (state == IN_RESIDUE) 
                                    {
                                          PDBAtom* prot_atom = new PDBAtom;
                                          atom = RTTI::castTo<Atom>(*prot_atom);
                                          residue->insert(*prot_atom);

                                          // check the atom flags, whether this is a PDB HETATM: 
                                          // any HETATM sets its residue to NON_STANDARD 
                                          // (compare FORMAT/PDBFile:readRecordHETATM)
                                          if (getLine().getField(5).has('h')) 
                                          {
                                                residue->clearProperty(Residue::PROPERTY__AMINO_ACID);
                                                residue->setProperty(Residue::PROPERTY__NON_STANDARD);
                                          } 
                                          else 
                                          {
                                                residue->setProperty(Residue::PROPERTY__AMINO_ACID);
                                                residue->clearProperty(Residue::PROPERTY__NON_STANDARD);
                                          }
                                    } 
                                    else 
                                    {
                                          atom = new Atom;
                                          if (molecule == 0) 
                                          {
                                                fragment->insert(*atom);
                                          } 
                                          else 
                                          {
                                                molecule->insert(*atom);
                                          }
                                    }

                                    atom->setName(getLine().getField(2));

                                    if (PTE[getLine().getField(3)] == Element::UNKNOWN)
                                    {
                                          // Treat "lone pair atoms" (Lp) written by Amber -- we just keep them
                                          // for now and remove them later.
                                          if (getLine().getField(3) != "Lp")
                                          {
                                                ERROR(String("unknown element: ") + getLine().getField(3))
                                          }
                                    }
                                    // Set the element and the atom radius (vdW radius).
                                    atom->setElement(PTE[getLine().getField(3)]);
                                    atom->setRadius(PTE[getLine().getField(3)].getVanDerWaalsRadius());

                                    if (getLine().getField(4) == "**")
                                    {
                                          atom->setTypeName("?");
                                    } 
                                    else  
                                    {
                                          atom->setTypeName(getLine().getField(4));
                                    }

                                    try 
                                    {
                                          atom->setCharge(getLine().getField(6).toFloat());
                                    }
                                    catch (Exception::InvalidFormat&)
                                    {
                                          ERROR(String("illegal charge ") + getLine().getField(6))
                                    }

                                    try 
                                    {
                                          atom->setPosition(Vector3(getLine().getField(7).toFloat(), 
                                                                                                                        getLine().getField(8).toFloat(), 
                                                                                                                        getLine().getField(9).toFloat()));
                                    }
                                    catch (Exception::InvalidFormat&)
                                    {
                                          ERROR(String("illegal position (") + getLine().getField(7) + " / " + getLine().getField(8) + " / " + getLine().getField(9) + ")")
                                    }

                                    Size number_of_atom_bonds;
                                    try 
                                    {
                                          number_of_atom_bonds = ((Size)getLine().getField(10).toInt());
                                    }
                                    catch (Exception::InvalidFormat&)
                                    {
                                          ERROR(String("illegal number of bonds: " + getLine().getField(10)))
                                    }

                                    Position atom_number;
                                    try 
                                    {
                                          atom_number = (Position)getLine().getField(1).toInt() - 1;
                                    }
                                    catch (Exception::InvalidFormat&)
                                    {
                                          ERROR(String("illegal atom number: ") + getLine().getField(1))
                                    }

                                    // Store the atom pointer in the atom_vector - we need it later to create the bonds!
                                    if (atom_number >= atom_vector.size())
                                    {
                                          atom_vector.resize(atom_vector.size() * 2);
                                    }
                                    atom_vector[atom_number] = atom;
                                    if ((atom_number <= last_atom) && (last_atom > 0))
                                    {
                                          ERROR(String("unordered atom indices: ") + getLine())
                                    }
                                    last_atom = atom_number;

                                    // now iterate over all bonds and insert them into the bond_vector 
                                    // this table will be processed afterwards to create the bonds, as most of
                                    // the bound atoms are not yet created by now
                                    if (number_of_atom_bonds > 0)
                                    {
                                          // check whether the total number of fields is consistent
                                          // with the number of bonds
                                          if (number_of_fields != (11 + 2 * number_of_atom_bonds))
                                          {
                                                // write an error message!
                                                if (number_of_fields < (11 + 2 * number_of_atom_bonds))
                                                {
                                                      ERROR(String("too few fields in atom line for ") + String(number_of_atom_bonds) + " bonds")
                                                }
                                                else 
                                                {
                                                      ERROR(String("too many fields in atom line for ") + String(number_of_atom_bonds) + " bonds")
                                                }
                                          }
                                          else 
                                          {
                                                for (Position i = 0 ; i < number_of_atom_bonds; i++) 
                                                { 
                                                      struct HINFileBondStruct bond;
                                                      bond.atom1 = atom_number;
                                                      try 
                                                      {
                                                            bond.atom2 = (Index)getLine().getField(11 + 2 * (Index)i).toInt() - 1;
                                                      }
                                                      catch (Exception::InvalidFormat&)
                                                      {
                                                            ERROR(String("illegal atom number for bond ") + String(i) + ": " + getLine().getField(11 + 2 * (Index)i))
                                                      }
                                                      Bond::Order order = Bond::ORDER__UNKNOWN;
                                                      String type_field = getLine().getField(12 + 2 * (Index)i);

                                                      if (type_field.size() == 1)
                                                      {
                                                            switch (type_field[0]) 
                                                            {
                                                                  case 's':   order = Bond::ORDER__SINGLE;        break;
                                                                  case 'd':   order = Bond::ORDER__DOUBLE;        break;
                                                                  case 't': order = Bond::ORDER__TRIPLE;          break;
                                                                  case 'a': order = Bond::ORDER__AROMATIC;  break;
                                                            }
                                                      }
                                                      bond.order = order;
                                                      bond_vector.push_back(bond);
                                                }
                                          }
                                    }
                              }     
                              else 
                              {
                                    ERROR(String("<atom> tag may appear only inside a <mol> or <res>!"))
                              }
                              continue;
                        }


                        if (tag == "vel")
                        {
                              // read the velocity of an atom
                              // since HyperChem uses the same units for velocities 
                              // (resp. A/ps) we don't need a conversion
                              Vector3 velocity;

                              try
                              {
                                    velocity.x = getLine().getField(2).toFloat();
                                    velocity.y = getLine().getField(3).toFloat();
                                    velocity.z = getLine().getField(4).toFloat();
                              }
                              catch (Exception::InvalidFormat&)
                              {
                                    ERROR(String("illegal velocity (") 
                                                                         + getLine().getField(2) + " / "
                                                                         + getLine().getField(3) + " / " 
                                                                         + getLine().getField(4) + ")")
                              }

                                    
                              // check whether the atom exists
                              Position atom_number;
                              try
                              {
                                    atom_number = (Position)getLine().getField(1).toInt() - 1;
                              }
                              catch (Exception::InvalidFormat&)
                              {
                                    ERROR(String("illegal atom number ") + getLine().getField(1))
                              }

                              if (atom_number > last_atom)
                              {
                                    ERROR(String("cannot assign velocity for atom ") + String(atom_number) + ": atom not defined!")
                              } 

                              if (atom_vector[atom_number] != 0)
                              {
                                    atom_vector[atom_number]->setVelocity(velocity);
                              } 
                              else 
                              {
                                    ERROR(String("cannot assign velocity for atom ") + String(atom_number) + ": atom not defined!")
                              }
                                    
                              continue;
                        }


                        if (tag == "res")
                        {
                              // remember where we are.
                              if (state != IN_MOLECULE) 
                              {
                                    ERROR(String("<res> tag must be inside a <mol>/<endmol>"))
                              }

                              state = IN_RESIDUE;
                              
                              // create a protein if it doesn't exist already
                              if (protein == 0)
                              {
                                    protein = new Protein;
                              }

                              // check whether we already have a chain to insert into
                              if (chain == 0)
                              {
                                    protein->insert(*(chain = new Chain));
                              }

                              // create a new residue and insert it into the chain
                              residue = new Residue;
                              chain->insert(*residue);

                              // set the residue's name
                              if (getLine().getField(2) != "-")
                              {
                                    residue->setName(getLine().getField(2));
                              }

                              // set the residue's PDB ID
                              if (getLine().getField(3) != "-")
                              {
                                    residue->setID(getLine().getField(3));
                              }

                              // set the chain name 
                              if (getLine().getField(5) != "-")
                              {
                                    chain->setName(getLine().getField(5));
                              }
                              
                              // create a fragment to insert the "loose" atoms into
                              if (fragment == 0)
                              {
                                    fragment = new Fragment;
                                    chain->AtomContainer::insert(*fragment);
                              }

                              // now check for a molecule, that might already exist
                              // and move its atoms into a fragment (to be inserted into the chain)
                              if (molecule == 0)
                              {
                                    continue;
                              }

                              AtomIterator atom_it = molecule->beginAtom();
                              vector<Atom*>     tmp_atoms;
                              for (; +atom_it; ++atom_it)
                              {
                                    tmp_atoms.push_back(&*atom_it);
                              }
                              for (Size i = 0; i < tmp_atoms.size(); ++i)
                              {
                                    fragment->insert(*tmp_atoms[i]);
                              }
                              
                              // delete the (now empty!) molecule and clear the pointer
                              delete molecule;
                              molecule = 0;
                              continue;
                        }


                        if (tag == "endres")
                        {
                              if (state != IN_RESIDUE) 
                              {
                                    ERROR(String("<endres> without <res>!"))
                              }
                              
                              state = IN_MOLECULE;

                              // reset the residue pointer to zero
                              residue = 0;
                              
                              continue;
                        }


                        if (tag == "mol")
                        {
                              if (state != START) 
                              {
                                    ERROR(String("<mol> inside <mol> or <res> definition!"))
                              }

                              state = IN_MOLECULE;

                              // create a new molecule and insert it into the system.
                              // We do not yet know, whether this contains residues.
                              // If it does, we have to convert it to a protein afterwards.
                              molecule = new Molecule;

                              if (getLine().countFields() > 2)
                              {
                                    // Determine the thrid field (the molecule name).
                                    // The name may (latest versions) be enclosed by double quotes.
                                    std::vector<String> fields;
                                    getLine().splitQuoted(fields, String::CHARACTER_CLASS__WHITESPACE, "\"");
                                    String name(fields[2]);

                                    if ((name != "") && (name != "-"))
                                    {
                                          // Remove leading/trailing blanks from the name.
                                          name.trim();
                                          
                                          // For newer versions of HyperChem, the name has to be
                                          // enclosed in double quotes.
                                          if (name[0] == '"')
                                          {
                                                name.erase(0, 1);
                                          }
                                          if (name[name.size() - 1] == '"')
                                          {
                                                name.erase(name.size() - 1, 1);
                                          }
                                          molecule->setName(name);
                                    }     
                              }
                              
                              continue;
                        }

                        if (tag == "endmol")
                        {
                              if (state != IN_MOLECULE)
                              {
                                    ERROR(String("missing <mol> or <endres> tag!"))
                              }

                              state = START;
                              
                              if (fragment != 0)
                              {
                                    if (fragment->countAtoms() == 0)
                                    {
                                          delete fragment;
                                    }
                              }
                              fragment = 0;
                              if (chain != 0)
                              {
                                    if (chain->countAtoms() == 0)
                                    {
                                          delete chain;
                                    }
                              }
                              chain = 0;


                              // Now, build all bonds
                              for (Size i = 0; i < bond_vector.size(); i++)
                              {
                                    // check whether both atoms exist
                                    if (bond_vector[i].atom1 > last_atom || 
                                                bond_vector[i].atom2 > last_atom)
                                    {
                                          // complain if one of the atoms does not exist
                                          ERROR(String("HINFile: cannot create bond from atom ") + String(bond_vector[i].atom1 + 1)
                                                            + " to atom " + String(bond_vector[i].atom2 + 1) + " of molecule " 
                                                            + getLine().getField(1) + " - non-existing atom!")

                                    } 
                                    else  
                                    {
                                          // everything all right, create the bond
                                          Bond* b = atom_vector[bond_vector[i].atom1]->createBond(*atom_vector[bond_vector[i].atom2]);
                                          b->setOrder(bond_vector[i].order);
                                          
                                          // Fix up the disulphide bridges in proteins.
                                          // Make sure we have
                                          //   -- a protein
                                          //   -- two sulfurs
                                          //   -- connecting different residues
                                          //   -- these residues are not hetero residues.
                                          PDBAtom* a1 = dynamic_cast<PDBAtom*>(const_cast<Atom*>(b->getFirstAtom()));
                                          PDBAtom* a2 = dynamic_cast<PDBAtom*>(const_cast<Atom*>(b->getSecondAtom()));
                                          if ((a1 != 0) && (a2 != 0)
                                                      && (a1->getElement() == PTE[Element::S]) && (a2->getElement() == PTE[Element::S])
                                                      && (a1->getResidue() != a2->getResidue())
                                                      && (a1->getResidue() != 0) && (a2->getResidue() != 0)
                                                      && (a1->getResidue()->hasProperty(Residue::PROPERTY__AMINO_ACID)) 
                                                      && (a2->getResidue()->hasProperty(Residue::PROPERTY__AMINO_ACID)))
                                          {
                                                a1->getResidue()->setProperty(Residue::PROPERTY__HAS_SSBOND);
                                                a2->getResidue()->setProperty(Residue::PROPERTY__HAS_SSBOND);
                                          }
                                    }
                              }

                              bond_vector.clear();
                              atom_vector.clear();
                              last_atom = 0;

                              throw(Exception::IndexOverflow(__FILE__, __LINE__));
                        }


                        if (tag == "sys")
                        {
                              // set the system's temperature
                              try
                              {
                                    temperature_ = getLine().getField(1).toFloat();
                              }
                              catch (Exception::InvalidFormat&)
                              {
                                    ERROR(String("illegal temperature ") + getLine().getField(1))
                              }

                              continue;
                        }


                        if (tag == "box")
                        {
                              // retrieve the periodic boundary
                              // we assume that the box is centered about the origin
                              // of the coordinate system
                              // The manual says the parameters are the dimensions of the box,
                              // so we have to divide by two.
                              try
                              {
                                    box_.a.x = - getLine().getField(1).toFloat() / 2.0;
                                    box_.a.y = - getLine().getField(2).toFloat() / 2.0;
                                    box_.a.z = - getLine().getField(3).toFloat() / 2.0;
                              }
                              catch (Exception::InvalidFormat&)
                              {
                                    ERROR(String("illegal box position (")
                                                                         + getLine().getField(1)  + " / "
                                                                         + getLine().getField(2)  + " / "
                                                                         + getLine().getField(3)  + " )")
                              }           

                              box_.b.x = - box_.a.x;
                              box_.b.y = - box_.a.y;
                              box_.b.z = - box_.a.z;
                              continue;
                        }


                        // ignore the irrelevant fields
                        if (tag == "forcefield" ||
                                    tag == "user1color" ||
                                    tag == "user2color" ||
                                    tag == "user3color" ||
                                    tag == "user4color" ||
                                    tag == "view"                       ||
                                    tag == "seed"                       ||
                                    tag == "mass"                       ||
                                    tag == "basisset"       ||
                                    tag == "selection"      ||
                                    tag == "endselection"          ||
                                    tag == "selectrestraint" ||
                                    tag == "selectatom"                  ||
                                    tag == "formalcharge")
                        {
                              continue;
                        }     

                        // if the tag was not recognized: complain about it
                        Log.warn() << "HINFile: unknown tag " << tag << " ignored." << std::endl;
                  }
            }
            catch (Exception::ParseError& e)
            {
                  // Delete all stray atoms. The order is important:
                  // since fragment and residue could be contained in 
                  // chain, molecule, etc., they have to be deleted first!
                  delete fragment;
                  delete residue;
                  delete chain;
                  delete molecule;
                  delete protein;
                  throw e;
            }
            catch (Exception::IndexOverflow&)
            {
            }

            // if we read a protein, return it as a molecule anyway
            if (molecule == 0)
            {
                  molecule = protein;
                  protein = 0;
            }

            // Remove the lone pairs from old AMBER HC-Files 
            if (molecule != 0)
            {
                  // a list to hold the lone pairs (for deletion)
                  list<Atom*> del_list;
                  
                  // Iterate over all atoms
                  AtomIterator it = molecule->beginAtom();
                  for (; +it; ++it)
                  {
                        if (it->getElement().getSymbol() == "Lp")
                        {
                              // store lone pair in the del_list
                              del_list.push_back(&*it);
                              
                              // sum the lone pair charge into the
                              // heavy atom it is bound to
                              if (it->countBonds() > 0)
                              {
                                    float charge = it->getCharge() / (float)it->countBonds();
                                    for (Atom::BondIterator bond_it = it->beginBond(); +bond_it; ++bond_it)
                                    {
                                          Atom* partner = bond_it->getPartner(*it);
                                          if (partner != 0)
                                          {
                                                partner->setCharge(partner->getCharge() + charge);
                                          }
                                    }
                              }

                              // remove all bonds to the lone pair
                              it->destroyBonds();
                        }
                  }
                  
                  // remove the lone pairs
                  list<Atom*>::iterator list_it = del_list.begin();
                  for (; list_it != del_list.end(); ++list_it)
                  {
                        if ((*list_it)->isAutoDeletable())
                        {
                              // delete dynamically created objects
                              delete *list_it;
                        } 
                        else 
                        {
                              // destroy static atoms
                              (*list_it)->destroy();
                        }
                  }
            }

            // return the resulting molecule
            return molecule;
      }

01066       void HINFile::initRead_()
      {
            // reset some private members
            box_.a.set(0.0);
            box_.b.set(0.0);
            temperature_ = 0.0;
      }

01074       bool HINFile::hasPeriodicBoundary() const
      {
            return (box_.a != box_.b);
      }

01079       SimpleBox3 HINFile::getPeriodicBoundary() const
      {
            return box_;
      }

01084       float HINFile::getTemperature() const
      {
            return temperature_;
      }

} // namespace BALL

Generated by  Doxygen 1.6.0   Back to index