Changeset 568be7


Ignore:
Timestamp:
Nov 7, 2009, 12:16:57 PM (15 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:
b9947d
Parents:
96c961
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
Location:
src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    r96c961 r568be7  
    12491249  molecule *mol = new molecule(periode);
    12501250
     1251  if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
     1252    eLog() << Verbose(0) << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
     1253  }
     1254
     1255
     1256  // first save as PDB data
     1257  if (ConfigFileName != NULL)
     1258    strcpy(filename, ConfigFileName);
     1259  if (output == NULL)
     1260    strcpy(filename,"main_pcp_linux");
     1261  Log() << Verbose(0) << "Saving as pdb input ";
     1262  if (configuration->SavePDB(filename, molecules))
     1263    Log() << Verbose(0) << "done." << endl;
     1264  else
     1265    Log() << Verbose(0) << "failed." << endl;
     1266
     1267  // then save as tremolo data file
     1268  if (ConfigFileName != NULL)
     1269    strcpy(filename, ConfigFileName);
     1270  if (output == NULL)
     1271    strcpy(filename,"main_pcp_linux");
     1272  Log() << Verbose(0) << "Saving as tremolo data input ";
     1273  if (configuration->SaveTREMOLO(filename, molecules))
     1274    Log() << Verbose(0) << "done." << endl;
     1275  else
     1276    Log() << Verbose(0) << "failed." << endl;
     1277
    12511278  // translate each to its center and merge all molecules in MoleculeListClass into this molecule
    12521279  int N = molecules->ListOfMolecules.size();
    … …  
    13291356    eLog() << Verbose(0) << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
    13301357  }
     1358
    13311359  delete(mol);
    13321360};

  • src/config.cpp

    r96c961 r568be7  
    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;

  • src/config.hpp

    r96c961 r568be7  
    1010
    1111using namespace std;
    12 
    13 class MoleculeListClass;
    1412
    1513/*********************************************** includes ***********************************/
    … …  
    2725
    2826class molecule;
     27class MoleculeListClass;
    2928class periodentafel;
    3029
    … …  
    142141  bool Save(const char * const filename, const periodentafel * const periode, molecule * const mol) const;
    143142  bool SaveMPQC(const char * const filename, const molecule * const mol) const;
     143  bool SavePDB(const char * const filename, const MoleculeListClass * const MolList) const;
     144  bool SavePDB(const char * const filename, const molecule * const mol) const;
     145  bool SaveTREMOLO(const char * const filename, const molecule * const mol) const;
     146  bool SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const;
     147
    144148  void Edit();
    145149  bool GetIsAngstroem() const;

  • src/molecule.hpp

    r96c961 r568be7  
    322322  void Output(ofstream *out);
    323323  void DissectMoleculeIntoConnectedSubgraphs(molecule * const mol, config * const configuration);
     324  int CountAllAtoms() const;
    324325
    325326  // merging of molecules

  • src/moleculelist.cpp

    r96c961 r568be7  
    743743{
    744744  // 1. dissect the molecule into connected subgraphs
    745   configuration ->BG->ConstructBondGraph(mol);
     745  configuration->BG->ConstructBondGraph(mol);
    746746
    747747  // 2. scan for connected subgraphs
    … …  
    780780
    781781  // 4c. relocate atoms to new molecules and remove from Leafs
    782   Walker = mol->start;
     782  Walker = NULL;
    783783  while (mol->start->next != mol->end) {
    784784    Walker = mol->start->next;
    … …  
    797797    }
    798798  }
    799   // 4d. we don't need to redo bonds, as they are connected subgraphs and still maintained their ListOfBonds
     799  // 4d. we don't need to redo bonds, as they are connected subgraphs and still maintained their ListOfBonds, but we have to remove them from first..last list
     800  bond *Binder = mol->first;
     801  while (mol->first->next != mol->last) {
     802    Binder = mol->first->next;
     803    Walker = Binder->leftatom;
     804    unlink(Binder);
     805    link(Binder,molecules[MolMap[Walker->nr]-1]->last);   // counting starts at 1
     806  }
    800807  // 4e. free Leafs
    801808  MolecularWalker = Subgraphs;
    … …  
    809816  Log() << Verbose(1) << "I scanned " << FragmentCounter << " molecules." << endl;
    810817};
     818
     819/** Count all atoms in each molecule.
     820 * \return number of atoms in the MoleculeListClass.
     821 * TODO: the inner loop should be done by some (double molecule::CountAtom()) function
     822 */
     823int MoleculeListClass::CountAllAtoms() const
     824{
     825  atom *Walker = NULL;
     826  int AtomNo = 0;
     827  for (MoleculeList::const_iterator MolWalker = ListOfMolecules.begin(); MolWalker != ListOfMolecules.end(); MolWalker++) {
     828    Walker = (*MolWalker)->start;
     829    while (Walker->next != (*MolWalker)->end) {
     830      Walker = Walker->next;
     831      AtomNo++;
     832    }
     833  }
     834  return AtomNo;
     835}
     836
    811837
    812838/******************************************* Class MoleculeLeafClass ************************************************/
Note: See TracChangeset for help on using the changeset viewer.