Ignore:
Timestamp:
Nov 7, 2009, 12:16:57 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Children:
1fcedd
Parents:
341850
git-author:
Frederik Heber <heber@…> (11/07/09 09:48:07)
git-committer:
Frederik Heber <heber@…> (11/07/09 12:16:57)
Message:

Added config::SavePDB() and config::SaveMPQC().

  • note: for CODICE we need to know about the different connected subgraphs created by the DFSAnalysis(). Hence, we write a pdb file which contains a resid number to discern in VMD between the molecules. Also, we need neighbour construction from a simple xyz file for TREMOLO in order to measure MSDs per water molecule.
  • new function in config.cpp: config::SavePDB() gets molecule or MoleculeListClass and writes PDB file
  • new function in config.cpp: config::SaveTREMOLO() gets molecule or MoleculeListClass and writes TREMOLO data file (with neighbours, mapped to global ids, and resid and resname)
  • new function in moleculelist.cpp: MoleculeListClass::CountAllAtoms() - counts all atoms.
  • BUGFIX: In MoleculeListClass::DissectMoleculeIntoConnectedSubgraphs() we did not shift the chained bond list from mol into the connected subgraphs. Thus, they were free'd on delete(mol) and no bonds were present afterwards. This is fixed, now we unlink() in mol and re-link() in the respective subgraph
File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/config.cpp

    r341850 r486aa5  
    55 */
    66
     7#include <stdio.h>
     8
    79#include "atom.hpp"
     10#include "bond.hpp"
    811#include "config.hpp"
    912#include "element.hpp"
    … …  
    200203    Free(&ThermostatNames[j]);
    201204  Free(&ThermostatNames);
     205
     206  if (BG != NULL)
     207    delete(BG);
    202208};
    203209
    … …  
    10551061  LoadMolecule(mol, FileBuffer, periode, FastParsing);
    10561062  mol->ActiveFlag = true;
    1057   MolList->insert(mol);
    10581063
    10591064  // 4. dissect the molecule into connected subgraphs
    1060   BG->ConstructBondGraph(mol);
    1061 
     1065  MolList->DissectMoleculeIntoConnectedSubgraphs(mol,this);
     1066
     1067  delete(mol);
    10621068  delete(FileBuffer);
    10631069};
    … …  
    14101416    delete(output);
    14111417    return result;
    1412   } else
     1418  } else {
     1419    eLog() << Verbose(1) << "Cannot open output file:" << filename << endl;
    14131420    return false;
     1421  }
    14141422};
    14151423
    … …  
    14301438    *fname << filename << ".in";
    14311439    output = new ofstream(fname->str().c_str(), ios::out);
     1440    if (output == NULL) {
     1441      eLog() << Verbose(1) << "Cannot open mpqc output file:" << fname << endl;
     1442      delete(fname);
     1443      return false;
     1444    }
    14321445    *output << "% Created by MoleCuilder" << endl;
    14331446    *output << "mpqc: (" << endl;
    … …  
    14681481    *fname << filename << ".hess.in";
    14691482    output = new ofstream(fname->str().c_str(), ios::out);
     1483    if (output == NULL) {
     1484      eLog() << Verbose(1) << "Cannot open mpqc hessian output file:" << fname << endl;
     1485      delete(fname);
     1486      return false;
     1487    }
    14701488    *output << "% Created by MoleCuilder" << endl;
    14711489    *output << "mpqc: (" << endl;
    … …  
    14991517    delete(fname);
    15001518  }
     1519
     1520  return true;
     1521};
     1522
     1523/** Stores all atoms from all molecules in a PDB input file.
     1524 * Note that this format cannot be parsed again.
     1525 * \param *filename name of file (without ".in" suffix!)
     1526 * \param *MolList pointer to MoleculeListClass containing all atoms
     1527 */
     1528bool config::SavePDB(const char * const filename, const MoleculeListClass * const MolList) const
     1529{
     1530  int AtomNo = -1;
     1531  int MolNo = 0;
     1532  atom *Walker = NULL;
     1533  FILE *f = NULL;
     1534
     1535  char name[MAXSTRINGSIZE];
     1536  strncpy(name, filename, MAXSTRINGSIZE-1);
     1537  strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
     1538  f = fopen(name, "w" );
     1539  if (f == NULL) {
     1540    eLog() << Verbose(1) << "Cannot open pdb output file:" << name << endl;
     1541    return false;
     1542  }
     1543  fprintf(f, "# Created by MoleCuilder\n");
     1544
     1545  for (MoleculeList::const_iterator Runner = MolList->ListOfMolecules.begin(); Runner != MolList->ListOfMolecules.end(); Runner++) {
     1546    Walker = (*Runner)->start;
     1547    int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
     1548    AtomNo = 0;
     1549    while (Walker->next != (*Runner)->end) {
     1550      Walker = Walker->next;
     1551      sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
     1552      elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
     1553      fprintf(f,
     1554             "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
     1555             Walker->nr,                /* atom serial number */
     1556             name,         /* atom name */
     1557             (*Runner)->name,      /* residue name */
     1558             'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
     1559             MolNo,         /* residue sequence number */
     1560             Walker->node->x[0],                 /* position X in Angstroem */
     1561             Walker->node->x[1],                 /* position Y in Angstroem */
     1562             Walker->node->x[2],                 /* position Z in Angstroem */
     1563             (double)Walker->type->Valence,         /* occupancy */
     1564             (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
     1565             "0",            /* segment identifier */
     1566             Walker->type->symbol,    /* element symbol */
     1567             "0");           /* charge */
     1568      AtomNo++;
     1569    }
     1570    Free(&elementNo);
     1571    MolNo++;
     1572  }
     1573  fclose(f);
     1574
     1575  return true;
     1576};
     1577
     1578/** Stores all atoms in a PDB input file.
     1579 * Note that this format cannot be parsed again.
     1580 * \param *filename name of file (without ".in" suffix!)
     1581 * \param *mol pointer to molecule
     1582 */
     1583bool config::SavePDB(const char * const filename, const molecule * const mol) const
     1584{
     1585  int AtomNo = -1;
     1586  atom *Walker = NULL;
     1587  FILE *f = NULL;
     1588
     1589  int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
     1590  char name[MAXSTRINGSIZE];
     1591  strncpy(name, filename, MAXSTRINGSIZE-1);
     1592  strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
     1593  f = fopen(name, "w" );
     1594  if (f == NULL) {
     1595    eLog() << Verbose(1) << "Cannot open pdb output file:" << name << endl;
     1596    Free(&elementNo);
     1597    return false;
     1598  }
     1599  fprintf(f, "# Created by MoleCuilder\n");
     1600
     1601  Walker = mol->start;
     1602  AtomNo = 0;
     1603  while (Walker->next != mol->end) {
     1604    Walker = Walker->next;
     1605    sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
     1606    elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
     1607    fprintf(f,
     1608           "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
     1609           Walker->nr,                /* atom serial number */
     1610           name,         /* atom name */
     1611           mol->name,      /* residue name */
     1612           'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
     1613           0,         /* residue sequence number */
     1614           Walker->node->x[0],                 /* position X in Angstroem */
     1615           Walker->node->x[1],                 /* position Y in Angstroem */
     1616           Walker->node->x[2],                 /* position Z in Angstroem */
     1617           (double)Walker->type->Valence,         /* occupancy */
     1618           (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
     1619           "0",            /* segment identifier */
     1620           Walker->type->symbol,    /* element symbol */
     1621           "0");           /* charge */
     1622    AtomNo++;
     1623  }
     1624  fclose(f);
     1625  Free(&elementNo);
     1626
     1627  return true;
     1628};
     1629
     1630/** Stores all atoms in a TREMOLO data input file.
     1631 * Note that this format cannot be parsed again.
     1632 * \param *filename name of file (without ".in" suffix!)
     1633 * \param *mol pointer to molecule
     1634 */
     1635bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
     1636{
     1637  atom *Walker = NULL;
     1638  ofstream *output = NULL;
     1639  stringstream * const fname = new stringstream;
     1640
     1641  *fname << filename << ".data";
     1642  output = new ofstream(fname->str().c_str(), ios::out);
     1643  if (output == NULL) {
     1644    eLog() << Verbose(1) << "Cannot open tremolo output file:" << fname << endl;
     1645    delete(fname);
     1646    return false;
     1647  }
     1648
     1649  // scan maximum number of neighbours
     1650  Walker = mol->start;
     1651  int MaxNeighbours = 0;
     1652  while (Walker->next != mol->end) {
     1653    Walker = Walker->next;
     1654    const int count = Walker->ListOfBonds.size();
     1655    if (MaxNeighbours < count)
     1656      MaxNeighbours = count;
     1657  }
     1658  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
     1659
     1660  Walker = mol->start;
     1661  while (Walker->next != mol->end) {
     1662    Walker = Walker->next;
     1663    *output << Walker->nr << "\t";
     1664    *output << Walker->Name << "\t";
     1665    *output << mol->name << "\t";
     1666    *output << 0 << "\t";
     1667    *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
     1668    *output << (double)Walker->type->Valence << "\t";
     1669    *output << Walker->type->symbol << "\t";
     1670    for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
     1671      *output << (*runner)->GetOtherAtom(Walker)->nr << "\t";
     1672    for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
     1673      *output << "-\t";
     1674    *output << endl;
     1675  }
     1676  output->flush();
     1677  output->close();
     1678  delete(output);
     1679  delete(fname);
     1680
     1681  return true;
     1682};
     1683
     1684/** Stores all atoms from all molecules in a TREMOLO data input file.
     1685 * Note that this format cannot be parsed again.
     1686 * \param *filename name of file (without ".in" suffix!)
     1687 * \param *MolList pointer to MoleculeListClass containing all atoms
     1688 */
     1689bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
     1690{
     1691  atom *Walker = NULL;
     1692  ofstream *output = NULL;
     1693  stringstream * const fname = new stringstream;
     1694
     1695  *fname << filename << ".data";
     1696  output = new ofstream(fname->str().c_str(), ios::out);
     1697  if (output == NULL) {
     1698    eLog() << Verbose(1) << "Cannot open tremolo output file:" << fname << endl;
     1699    delete(fname);
     1700    return false;
     1701  }
     1702
     1703  // scan maximum number of neighbours
     1704  int MaxNeighbours = 0;
     1705  for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     1706    Walker = (*MolWalker)->start;
     1707    while (Walker->next != (*MolWalker)->end) {
     1708      Walker = Walker->next;
     1709      const int count = Walker->ListOfBonds.size();
     1710      if (MaxNeighbours < count)
     1711        MaxNeighbours = count;
     1712    }
     1713  }
     1714  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
     1715
     1716  // create global to local id map
     1717  int **LocalNotoGlobalNoMap = Calloc<int *>(MolList->ListOfMolecules.size(), "config::SaveTREMOLO - **LocalNotoGlobalNoMap");
     1718  {
     1719    int MolCounter = 0;
     1720    int AtomNo = 0;
     1721    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     1722      LocalNotoGlobalNoMap[MolCounter] = Calloc<int>(MolList->CountAllAtoms(), "config::SaveTREMOLO - *LocalNotoGlobalNoMap[]");
     1723
     1724      (*MolWalker)->SetIndexedArrayForEachAtomTo( LocalNotoGlobalNoMap[MolCounter], &atom::nr, IncrementalAbsoluteValue, &AtomNo);
     1725
     1726      MolCounter++;
     1727    }
     1728  }
     1729
     1730  // write the file
     1731  {
     1732    int MolCounter = 0;
     1733    int AtomNo = 0;
     1734    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     1735      Walker = (*MolWalker)->start;
     1736      while (Walker->next != (*MolWalker)->end) {
     1737        Walker = Walker->next;
     1738        *output << AtomNo << "\t";
     1739        *output << Walker->Name << "\t";
     1740        *output << (*MolWalker)->name << "\t";
     1741        *output << MolCounter << "\t";
     1742        *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
     1743        *output << (double)Walker->type->Valence << "\t";
     1744        *output << Walker->type->symbol << "\t";
     1745        for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
     1746          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ] << "\t";
     1747        for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
     1748          *output << "-\t";
     1749        *output << endl;
     1750        AtomNo++;
     1751      }
     1752      MolCounter++;
     1753    }
     1754  }
     1755
     1756  // store & free
     1757  output->flush();
     1758  output->close();
     1759  delete(output);
     1760  delete(fname);
     1761  for(size_t i=0;i<MolList->ListOfMolecules.size(); i++)
     1762    Free(&LocalNotoGlobalNoMap[i]);
     1763  Free(&LocalNotoGlobalNoMap);
    15011764
    15021765  return true;
Note: See TracChangeset for help on using the changeset viewer.