Changeset 0d1ad0 for src


Ignore:
Timestamp:
Jun 25, 2010, 9:42:28 AM (15 years ago)
Author:
Tillmann Crueger <crueger@…>
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:
04488a, 0c5eeb, 93987b
Parents:
6d574a (diff), a356f2 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'stable' into StructureRefactoring

Conflicts:

molecuilder/src/World.cpp

Location:
src
Files:
14 added
42 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/ActionRegistry.cpp

    r6d574a r0d1ad0  
    1919using namespace std;
    2020
     21/** Constructor for class ActionRegistry.
     22 */
    2123ActionRegistry::ActionRegistry()
    2224{
    2325}
    2426
     27/** Destructor for class ActionRegistry.
     28 */
    2529ActionRegistry::~ActionRegistry()
    2630{
    … …  
    3236}
    3337
     38/** Returns pointer to an action named by \a name.
     39 * \param name name of action
     40 * \return pointer to Action
     41 */
    3442Action* ActionRegistry::getActionByName(const std::string name){
    3543  map<const string,Action*>::iterator iter;
    … …  
    3947}
    4048
     49/** States whether action is present or not.
     50 * \note This iss needed as ActionRegistry::getActionByName() ASSERT()s that action is in map.
     51 * \param name name of action
     52 * \return true - Action present, false - Action absent
     53 */
    4154bool ActionRegistry::isActionByNamePresent(const std::string name){
    4255  map<const string,Action*>::iterator iter;
    … …  
    4558}
    4659
     60/** Registers an Action with the ActionRegistry.
     61 * \param *action pointer to Action.
     62 */
    4763void ActionRegistry::registerAction(Action* action){
    4864  pair<map<const string,Action*>::iterator,bool> ret;
     65  //cout << "Trying to register action with name " << action->getName() << "." << endl;
    4966  ret = actionMap.insert(pair<const string,Action*>(action->getName(),action));
    5067  ASSERT(ret.second,"Two actions with the same name added to registry");
    5168}
    5269
     70/** Unregisters an Action.
     71 * \param *action pointer to Action.
     72 */
    5373void ActionRegistry::unregisterAction(Action* action){
     74  //cout << "Unregistering action with name " << action->getName() << "." << endl;
    5475  actionMap.erase(action->getName());
    5576}
    5677
     78/** Returns an iterator pointing to the start of the map of Action's.
     79 * \return begin iterator
     80 */
    5781std::map<const std::string,Action*>::iterator ActionRegistry::getBeginIter()
    5882{
    … …  
    6084}
    6185
     86/** Returns an iterator pointing to the end of the map of Action's.
     87 * \return end iterator
     88 */
    6289std::map<const std::string,Action*>::iterator ActionRegistry::getEndIter()
    6390{
    … …  
    6592}
    6693
     94/** Returns a const iterator pointing to the start of the map of Action's.
     95 * \return constant begin iterator
     96 */
     97std::map<const std::string,Action*>::const_iterator ActionRegistry::getBeginIter() const
     98{
     99  return actionMap.begin();
     100}
     101
     102/** Returns a const iterator pointing to the end of the map of Action's.
     103 * \return constant end iterator
     104 */
     105std::map<const std::string,Action*>::const_iterator ActionRegistry::getEndIter() const
     106{
     107  return actionMap.end();
     108}
     109
     110/** Prints the contents of the ActionRegistry \a &m to \a &ost.
     111 * \param &ost output stream
     112 * \param &m reference to ActionRegistry
     113 * \return reference to the above out stream for concatenation
     114 */
     115ostream& operator<<(ostream& ost, const ActionRegistry& m)
     116{
     117  ost << "ActionRegistry contains:" << endl;
     118  for (std::map<const std::string,Action*>::const_iterator iter = m.getBeginIter(); iter != m.getEndIter(); ++iter) {
     119    ost << "\t" << iter->first << " with pointer " << iter->second << endl;
     120  }
     121  return ost;
     122};
     123
     124
     125
    67126CONSTRUCT_SINGLETON(ActionRegistry)

  • src/Actions/ActionRegistry.hpp

    r6d574a r0d1ad0  
    99#define ACTIONREGISTRY_HPP_
    1010
     11#include <iostream>
    1112#include <string>
    1213#include <map>
    … …  
    2627
    2728  std::map<const std::string,Action*>::iterator getBeginIter();
     29  std::map<const std::string,Action*>::const_iterator getBeginIter() const;
    2830  std::map<const std::string,Action*>::iterator getEndIter();
     31  std::map<const std::string,Action*>::const_iterator getEndIter() const;
    2932
    3033private:
    … …  
    3639};
    3740
     41std::ostream& operator<<(std::ostream& ost, const ActionRegistry& m);
     42
    3843#endif /* ACTIONREGISTRY_HPP_ */

  • src/Actions/CmdAction/BondLengthTableAction.cpp

    r6d574a r0d1ad0  
    99
    1010#include "Actions/CmdAction/BondLengthTableAction.hpp"
     11#include "bondgraph.hpp"
    1112#include "config.hpp"
    1213#include "log.hpp"

  • src/Actions/FragmentationAction/DepthFirstSearchAction.cpp

    r6d574a r0d1ad0  
    1010#include "Actions/FragmentationAction/DepthFirstSearchAction.hpp"
    1111#include "atom.hpp"
     12#include "bondgraph.hpp"
    1213#include "config.hpp"
    1314#include "log.hpp"

  • src/Actions/FragmentationAction/FragmentationAction.cpp

    r6d574a r0d1ad0  
    1010#include "Actions/FragmentationAction/FragmentationAction.hpp"
    1111#include "atom.hpp"
     12#include "bondgraph.hpp"
    1213#include "config.hpp"
    1314#include "log.hpp"
    … …  
    4142  double distance = -1.;
    4243  int order = 0;
     44  std::string path;
    4345  config *configuration = World::getInstance().getConfig();
    4446  int ExitFlag = 0;
    4547
    4648  cout << "pre-dialog"<< endl;
    47   dialog->queryMolecule(NAME, &mol, MapOfActions::getInstance().getDescription(NAME));
     49  dialog->queryString(NAME, &path, MapOfActions::getInstance().getDescription(NAME));
     50  dialog->queryMolecule("molecule-by-id", &mol, MapOfActions::getInstance().getDescription("molecule-by-id"));
    4851  dialog->queryDouble("distance", &distance, MapOfActions::getInstance().getDescription("distance"));
    4952  dialog->queryInt("order", &order, MapOfActions::getInstance().getDescription("order"));
    … …  
    5861    DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl);
    5962    if (mol->hasBondStructure()) {
    60       ExitFlag = mol->FragmentMolecule(order, configuration);
     63      ExitFlag = mol->FragmentMolecule(order, path);
    6164    }
    6265    World::getInstance().setExitFlag(ExitFlag);

  • src/Actions/Makefile.am

    r6d574a r0d1ad0  
    120120  WorldAction/CenterOnEdgeAction.cpp \
    121121  WorldAction/ChangeBoxAction.cpp \
     122  WorldAction/InputAction.cpp \
     123  WorldAction/OutputAction.cpp \
    122124  WorldAction/RemoveSphereOfAtomsAction.cpp \
    123125  WorldAction/RepeatBoxAction.cpp \
    … …  
    131133  WorldAction/CenterOnEdgeAction.hpp \
    132134  WorldAction/ChangeBoxAction.hpp \
     135  WorldAction/InputAction.hpp \
     136  WorldAction/OutputAction.hpp \
    133137  WorldAction/RemoveSphereOfAtomsAction.hpp \
    134138  WorldAction/RepeatBoxAction.hpp \

  • src/Actions/MapOfActions.cpp

    r6d574a r0d1ad0  
    8383  DescriptionMap["fragment-mol"] = "create for a given molecule into fragments up to given order";
    8484  DescriptionMap["help"] = "Give this help screen";
     85  DescriptionMap["input"] = "specify input files";
    8586  DescriptionMap["linear-interpolate"] = "linear interpolation in discrete steps between start and end position of a molecule";
    8687  DescriptionMap["nonconvex-envelope"] = "create the non-convex envelope for a molecule";
    8788  DescriptionMap["molecular-volume"] = "calculate the volume of a given molecule";
     89  DescriptionMap["output"] = "specify output formats";
    8890  DescriptionMap["pair-correlation"] = "pair correlation analysis between two elements, element and point or element and surface";
    8991  DescriptionMap["parse-xyz"] = "parse xyz file into World";
    … …  
    185187  TypeMap["fastparsing"] = Boolean;
    186188  TypeMap["fill-molecule"] = String;
    187   TypeMap["fragment-mol"] = Molecule;
     189  TypeMap["fragment-mol"] = String;
    188190  TypeMap["input"] = String;
    189191  TypeMap["linear-interpolate"] = String;
    190192  TypeMap["molecular-volume"] = Molecule;
    191193  TypeMap["nonconvex-envelope"] = Molecule;
     194  TypeMap["output"] = String;
    192195  TypeMap["parse-xyz"] = String;
    193196  TypeMap["pair-correlation"] = String;
    … …  
    262265  generic.insert("fragment-mol");
    263266  generic.insert("help");
    264         generic.insert("linear-interpolate");
     267  generic.insert("input");
     268  generic.insert("linear-interpolate");
    265269//  generic.insert("molecular-volume");
    266270  generic.insert("nonconvex-envelope");
     271  generic.insert("output");
    267272        generic.insert("pair-correlation");
    268 //      generic.insert("parse-xyz");
     273        generic.insert("parse-xyz");
    269274//  generic.insert("principal-axis-system");
    270275  generic.insert("remove-atom");

  • src/Actions/MoleculeAction/FillWithMoleculeAction.cpp

    r6d574a r0d1ad0  
    102102      World::getInstance().getMolecules()->insert(Filling);
    103103    }
     104    for (molecule::iterator iter = filler->begin(); !filler->empty(); iter = filler->begin()) {
     105      atom *Walker = *iter;
     106      filler->erase(iter);
     107      World::getInstance().destroyAtom(Walker);
     108    }
    104109    World::getInstance().destroyMolecule(filler);
    105110

  • src/Actions/MoleculeAction/LinearInterpolationofTrajectoriesAction.cpp

    r6d574a r0d1ad0  
    6969    if (IdMapping)
    7070      DoLog(1) && (Log() << Verbose(1) << "Using Identity for the permutation map." << endl);
    71     char outputname[MAXSTRINGSIZE];
    72     strcpy(outputname, filename.c_str());
    73     // TODO: LinearInterpolationBetweenConfiguration should use stream, not the filename directly! (better for unit test)
    74     if (!mol->LinearInterpolationBetweenConfiguration(start, end, outputname, *(World::getInstance().getConfig()), IdMapping))
    75       DoLog(2) && (Log() << Verbose(2) << "Could not store " << outputname << " files." << endl);
     71    if (!mol->LinearInterpolationBetweenConfiguration(start, end, filename, *(World::getInstance().getConfig()), IdMapping))
     72      DoLog(2) && (Log() << Verbose(2) << "Could not store " << filename << " files." << endl);
    7673    else
    7774      DoLog(2) && (Log() << Verbose(2) << "Steps created and " << filename << " files stored." << endl);

  • src/Actions/MoleculeAction/SaveAdjacencyAction.cpp

    r6d574a r0d1ad0  
    6565    World::getInstance().getConfig()->BG->ConstructBondGraph(mol);
    6666    // TODO: sollte stream nicht filename benutzen, besser fuer unit test
    67     char outputname[MAXSTRINGSIZE];
    68     strcpy(outputname, filename.c_str());
    69     mol->StoreAdjacencyToFile(NULL, outputname);
     67    mol->StoreAdjacencyToFile(filename);
    7068    delete dialog;
    7169    return Action::success;

  • src/Actions/MoleculeAction/SaveBondsAction.cpp

    r6d574a r0d1ad0  
    6464    DoLog(0) && (Log() << Verbose(0) << "Storing bonds to path " << filename << "." << endl);
    6565    World::getInstance().getConfig()->BG->ConstructBondGraph(mol);
    66     // TODO: sollte stream, nicht filenamen direkt nutzen, beser fuer unit tests
    67     char outputname[MAXSTRINGSIZE];
    68     strcpy(outputname, filename.c_str());
    69     mol->StoreBondsToFile(NULL, outputname);
     66    // TODO: sollte stream, nicht filenamen direkt nutzen, besser fuer unit tests
     67    mol->StoreBondsToFile(filename);
    7068    delete dialog;
    7169    return Action::success;

  • src/Actions/ParserAction/LoadXyzAction.cpp

    r6d574a r0d1ad0  
    88#include "Helpers/MemDebug.hpp"
    99
     10using namespace std;
     11
    1012#include "Actions/ParserAction/LoadXyzAction.hpp"
    1113#include "Parser/XyzParser.hpp"
    1214
    1315#include <iostream>
     16#include <set>
    1417#include <string>
     18#include <vector>
    1519
    16 using namespace std;
    1720
    1821#include "UIElements/UIFactory.hpp"
    … …  
    2124
    2225#include "atom.hpp"
     26#include "log.hpp"
    2327#include "molecule.hpp"
     28#include "verbose.hpp"
     29#include "World.hpp"
    2430
    2531/****** ParserLoadXyzAction *****/
    … …  
    3743//};
    3844
    39 const char ParserLoadXyzAction::NAME[] = "LoadXyz";
     45const char ParserLoadXyzAction::NAME[] = "parse-xyz";
    4046
    4147ParserLoadXyzAction::ParserLoadXyzAction() :
    … …  
    4854Action::state_ptr ParserLoadXyzAction::performCall() {
    4955  string filename;
    50   XyzParser parser;
    5156  Dialog *dialog = UIFactory::getInstance().makeDialog();
    5257
    53   dialog->queryString("filename",&filename, "Filename of the xyz file");
     58  dialog->queryString(NAME,&filename, NAME);
    5459
    5560  if(dialog->display()) {
     61    DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl);
    5662    // parse xyz file
    5763    ifstream input;
    5864    input.open(filename.c_str());
    59     if (!input.fail())
     65    if (!input.fail()) {
     66      // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include
     67      set <atom*> UniqueList;
     68      {
     69        vector<atom *> ListBefore = World::getInstance().getAllAtoms();
     70        for (vector<atom *>::iterator runner = ListBefore.begin();runner != ListBefore.end(); ++runner)
     71          UniqueList.insert(*runner);
     72      }
     73      XyzParser parser; // briefly instantiate a parser which is removed at end of focus
    6074      parser.load(&input);
     75      {
     76        vector<atom *> ListAfter = World::getInstance().getAllAtoms();
     77        pair< set<atom *>::iterator, bool > Inserter;
     78        if (UniqueList.size() != ListAfter.size()) { // only create if new atoms have been parsed
     79          MoleculeListClass *molecules = World::getInstance().getMolecules();
     80          molecule *mol= NULL;
     81          if (molecules->ListOfMolecules.empty()) {
     82            mol = World::getInstance().createMolecule();
     83            molecules->insert(mol);
     84          } else {
     85            mol = *(molecules->ListOfMolecules.begin());
     86          }
     87          for (vector<atom *>::iterator runner = ListAfter.begin(); runner != ListAfter.end(); ++runner) {
     88            Inserter = UniqueList.insert(*runner);
     89            if (Inserter.second) { // if not present, then new (just parsed) atom, add ...
     90              cout << "Adding new atom " << **runner << " to new mol." << endl;
     91              mol->AddAtom(*runner);
     92            }
     93          }
     94          mol->doCountAtoms();
     95        } else {
     96          cout << "No atoms parsed?" << endl;
     97        }
     98      }
     99    } else {
     100      DoeLog(1) && (eLog() << Verbose(1) << "Could not open file " << filename << "." << endl);
     101    }
    61102    input.close();
    62103  }

  • src/CommandLineParser.cpp

    r6d574a r0d1ad0  
    5656{
    5757  // go through all arguments
     58  cout << Verbose(1) << "By default putting input into the sequence." << endl;
     59  // TODO: This may be bad, because const string "input" is destroyed at end of function
     60  SequenceOfActions.push_back("input");
    5861  for (int i=1;i<argc;i++) {
    5962    (cout << Verbose(1) << "Checking on " << argv[i] << endl);

  • src/Legacy/oldmenu.cpp

    r6d574a r0d1ad0  
    609609{
    610610  int Order1;
     611  std::string path;
    611612  clock_t start, end;
    612613
    … …  
    614615  Log() << Verbose(0) << "What's the desired bond order: ";
    615616  cin >> Order1;
     617  DoLog(0) && (Log() << Verbose(0) << "What's the output path and prefix [e.g. /home/foo/BondFragment]: ");
     618  cin >> path;
    616619  if (mol->hasBondStructure()) {
    617620    start = clock();
    618     mol->FragmentMolecule(Order1, configuration);
     621    mol->FragmentMolecule(Order1, path);
    619622    end = clock();
    620623    Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;

  • src/Makefile.am

    r6d574a r0d1ad0  
    6363
    6464ACTIONSHEADER = \
    65   ${ANALYSISACTIONHEADER} \
    66   ${ATOMACTIONHEADER} \
    67   ${CMDACTIONHEADER} \
    68   ${FRAGMENTATIONACTIONHEADER} \
    69   ${MOLECULEACTIONHEADER} \
    70   ${PARSERACTIONHEADER} \
    71   ${TESSELATIONACTIONHEADER} \
    72   ${WORLDACTIONHEADER} \
    7365  Actions/Action.hpp \
    7466  Actions/ActionHistory.hpp \
    … …  
    8880  Parser/ChangeTracker.cpp \
    8981  Parser/FormatParser.cpp \
     82  Parser/FormatParserStorage.cpp \
     83  Parser/MpqcParser.cpp \
     84  Parser/PcpParser.cpp \
    9085  Parser/TremoloParser.cpp \
    9186  Parser/XyzParser.cpp
     87
    9288PARSERHEADER = \
    9389  Parser/ChangeTracker.hpp \
    9490  Parser/FormatParser.hpp \
     91  Parser/FormatParserStorage.hpp \
     92  Parser/MpqcParser.hpp \
     93  Parser/PcpParser.hpp \
    9594  Parser/TremoloParser.hpp \
    9695  Parser/XyzParser.hpp
    … …  
    153152  CommandLineParser.cpp \
    154153  config.cpp \
     154  ConfigFileBuffer.cpp \
    155155  element.cpp \
    156156  elements_db.cpp \
    … …  
    178178  tesselation.cpp \
    179179  tesselationhelpers.cpp \
     180  ThermoStatContainer.cpp \
    180181  triangleintersectionlist.cpp \
    181182  vector.cpp \
    … …  
    198199  CommandLineParser.hpp \
    199200  config.hpp \
     201  ConfigFileBuffer.hpp \
    200202  defs.hpp \
    201203  element.hpp \
    … …  
    220222  tesselation.hpp \
    221223  tesselationhelpers.hpp \
     224  ThermoStatContainer.hpp \
    222225  triangleintersectionlist.hpp \
    223226  verbose.hpp \
    … …  
    234237INCLUDES = -I$(top_srcdir)/src/unittests -I$(top_srcdir)/src/Actions -I$(top_srcdir)/src/UIElements
    235238
    236 noinst_LIBRARIES = libmolecuilder.a libgslwrapper.a
     239noinst_LIBRARIES = libmolecuilder.a libgslwrapper.a libparser.a
    237240bin_PROGRAMS = molecuilder joiner analyzer
    238241molecuilderdir = ${bindir}
    239242libmolecuilder_a_SOURCES = ${SOURCE} ${HEADER}
     243libparser_a_SOURCES = ${PARSERSOURCE} ${PARSERHEADER}
    240244libgslwrapper_a_SOURCES = ${LINALGSOURCE} ${LINALGHEADER}
    241245molecuilder_DATA = elements.db valence.db orbitals.db Hbonddistance.db Hbondangle.db
    242246molecuilder_LDFLAGS = $(BOOST_LDFLAGS)
    243247molecuilder_SOURCES = builder.cpp
    244 molecuilder_LDADD =  UIElements/libMolecuilderUI.a Actions/libMolecuilderActions.a libmolecuilder.a libgslwrapper.a $(BOOST_LIB) ${BOOST_THREAD_LIB} ${BOOST_PROGRAM_OPTIONS_LIB}
     248molecuilder_LDADD =  UIElements/libMolecuilderUI.a Actions/libMolecuilderActions.a libmolecuilder.a libparser.a libgslwrapper.a $(BOOST_LIB) ${BOOST_THREAD_LIB} ${BOOST_PROGRAM_OPTIONS_LIB}
    245249joiner_SOURCES = joiner.cpp datacreator.cpp parser.cpp datacreator.hpp helpers.hpp parser.hpp periodentafel.hpp
    246250joiner_LDADD = libmolecuilder.a $(BOOST_LIB) ${BOOST_THREAD_LIB}

  • src/Parser/ChangeTracker.cpp

    r6d574a r0d1ad0  
    77
    88#include "Helpers/MemDebug.hpp"
     9#include "Parser/ChangeTracker.hpp"
     10#include "Patterns/Singleton_impl.hpp"
    911
    10 #include "ChangeTracker.hpp"
    11 
    12 ChangeTracker* ChangeTracker::instance = NULL;
    1312
    1413/**
    … …  
    2726ChangeTracker::~ChangeTracker() {
    2827  World::getInstance().signOff(this);
    29 }
    30 
    31 /**
    32  * Returns the change tracker instance.
    33  *
    34  * \return this
    35  */
    36 ChangeTracker* ChangeTracker::get() {
    37   if (instance == NULL) {
    38     instance = new ChangeTracker();
    39   }
    40 
    41   return instance;
    42 }
    43 
    44 /**
    45  * Destroys the change tracker instance. Be careful, the change tracker is a
    46  * singleton and destruction might lead to a loss of consistency.
    47  */
    48 void ChangeTracker::destroy() {
    49   delete instance;
    50   instance = NULL;
    5128}
    5229
    … …  
    7653  }
    7754}
     55
     56CONSTRUCT_SINGLETON(ChangeTracker)

  • src/Parser/ChangeTracker.hpp

    r6d574a r0d1ad0  
    1818 * changes to it.
    1919 */
    20 class ChangeTracker : public Observable {
     20class ChangeTracker : public Singleton<ChangeTracker>, public Observable {
     21  friend class Singleton<ChangeTracker>;
    2122public:
    2223  void saveStatus();
    … …  
    3132  bool isConsistent;
    3233  static ChangeTracker* instance;
     34
     35  // private constructor and destructor due to singleton
    3336  ChangeTracker();
    3437  ~ChangeTracker();

  • src/Parser/FormatParser.cpp

    r6d574a r0d1ad0  
    1919  Observer("FormatParser")
    2020{
    21   ChangeTracker::get()->signOn(this);
     21  ChangeTracker::getInstance().signOn(this);
    2222  saveStream = NULL;
    2323}
    … …  
    2727 */
    2828FormatParser::~FormatParser() {
    29   ChangeTracker::get()->signOff(this);
     29  ChangeTracker::getInstance().signOff(this);
    3030}
    3131

  • src/Parser/TremoloParser.cpp

    r6d574a r0d1ad0  
    4949  knownKeys["GrpTypeNo"] = TremoloKey::GrpTypeNo;
    5050  knownKeys["torsion"] = TremoloKey::torsion;
     51
     52  // default behavior: use all possible keys on output
     53  for (std::map<std::string, TremoloKey::atomDataKey>::iterator iter = knownKeys.begin(); iter != knownKeys.end(); ++iter)
     54    usedFields.push_back(iter->first);
    5155}
    5256
    … …  
    193197
    194198  lineStream << line.substr(offset);
     199  usedFields.clear();
    195200  while (lineStream.good()) {
    196201    lineStream >> keyword;

  • src/Parser/TremoloParser.hpp

    r6d574a r0d1ad0  
    1010
    1111#include <string>
    12 #include "FormatParser.hpp"
     12#include "Parser/FormatParser.hpp"
    1313
    1414/**

  • src/Parser/XyzParser.cpp

    r6d574a r0d1ad0  
    2727 */
    2828XyzParser::~XyzParser() {
    29   delete(&comment);
    3029}
    3130
    … …  
    5857 */
    5958void XyzParser::save(ostream* file) {
     59  if (comment == "") {
     60    time_t now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     61    comment = "\tCreated by molecuilder on ";
     62    // ctime ends in \n\0, we have to cut away the newline
     63    std::string time(ctime(&now));
     64    size_t pos = time.find('\n');
     65    if (pos != 0)
     66      comment += time.substr(0,pos);
     67    else
     68      comment += time;
     69  }
    6070  *file << World::getInstance().numAtoms() << endl << comment << endl;
    6171
    6272  vector<atom*> atoms = World::getInstance().getAllAtoms();
    6373  for(vector<atom*>::iterator it = atoms.begin(); it != atoms.end(); it++) {
    64     *file << fixed << (*it)->getType()->symbol << "\t" << (*it)->x[0] << "\t" << (*it)->x[1] << "\t" << (*it)->x[2] << endl;
     74    *file << noshowpoint << (*it)->getType()->symbol << "\t" << (*it)->x[0] << "\t" << (*it)->x[1] << "\t" << (*it)->x[2] << endl;
    6575  }
    6676}

  • src/Parser/XyzParser.hpp

    r6d574a r0d1ad0  
    1010
    1111#include <string>
    12 #include "FormatParser.hpp"
     12#include "Parser/FormatParser.hpp"
    1313
    1414/**

  • src/UIElements/CommandLineUI/CommandLineWindow.cpp

    r6d574a r0d1ad0  
    4747#include "Actions/WorldAction/CenterOnEdgeAction.hpp"
    4848#include "Actions/WorldAction/ChangeBoxAction.hpp"
     49#include "Actions/WorldAction/InputAction.hpp"
     50#include "Actions/WorldAction/OutputAction.hpp"
    4951#include "Actions/WorldAction/RemoveSphereOfAtomsAction.hpp"
    5052#include "Actions/WorldAction/RepeatBoxAction.hpp"
    … …  
    8183
    8284void CommandLineWindow::display() {
     85  //cout << ActionRegistry::getInstance();
     86
    8387  // go through all possible actions
    8488  for (std::list<std::string>::iterator CommandRunner = CommandLineParser::getInstance().SequenceOfActions.begin(); CommandRunner != CommandLineParser::getInstance().SequenceOfActions.end(); ++CommandRunner) {
    85     cout << "Checking presence of " << *CommandRunner << endl;
    86     if (ActionRegistry::getInstance().isActionByNamePresent(*CommandRunner))
     89    cout << "Checking presence of " << *CommandRunner << ": ";
     90    if (ActionRegistry::getInstance().isActionByNamePresent(*CommandRunner)) {
     91      cout << "calling " << *CommandRunner << endl;
    8792      ActionRegistry::getInstance().getActionByName(*CommandRunner)->call();
     93    } else {
     94      cout << "absent." << endl;
     95    }
    8896  }
    8997}
    … …  
    152160  new WorldCenterOnEdgeAction();
    153161  new WorldChangeBoxAction();
     162  new WorldInputAction();
     163  new WorldOutputAction();
    154164  new WorldRemoveSphereOfAtomsAction();
    155165  new WorldRepeatBoxAction();

  • src/World.cpp

    r6d574a r0d1ad0  
    1414#include "molecule.hpp"
    1515#include "periodentafel.hpp"
     16#include "ThermoStatContainer.hpp"
    1617#include "Descriptors/AtomDescriptor.hpp"
    1718#include "Descriptors/AtomDescriptor_impl.hpp"
    … …  
    9091  defaultName = name;
    9192};
     93
     94class ThermoStatContainer * World::getThermostats()
     95{
     96  return Thermostats;
     97}
     98
    9299
    93100int World::getExitFlag() {
    … …  
    282289    periode(new periodentafel),
    283290    configuration(new config),
     291    Thermostats(new ThermoStatContainer),
    284292    ExitFlag(0),
    285293    atoms(),
    … …  
    307315  delete periode;
    308316  delete configuration;
     317  delete Thermostats;
    309318  MoleculeSet::iterator molIter;
    310319  for(molIter=molecules.begin();molIter!=molecules.end();++molIter){

  • src/World.hpp

    r6d574a r0d1ad0  
    3030
    3131// forward declarations
    32 class config;
    33 class periodentafel;
    34 class MoleculeListClass;
    3532class atom;
    36 class molecule;
    3733class AtomDescriptor;
    3834class AtomDescriptor_impl;
     35template<typename T> class AtomsCalculation;
     36class config;
     37class ManipulateAtomsProcess;
     38class molecule;
    3939class MoleculeDescriptor;
    4040class MoleculeDescriptor_impl;
    41 class ManipulateAtomsProcess;
    42 template<typename T>
    43 class AtomsCalculation;
     41class MoleculeListClass;
     42class periodentafel;
     43class ThermoStatContainer;
    4444
    4545/****************************************** forward declarations *****************************/
    … …  
    142142  void setDefaultName(std::string name);
    143143
     144  /**
     145   * get pointer to World's ThermoStatContainer
     146   */
     147  ThermoStatContainer * getThermostats();
     148
    144149  /*
    145150   * get the ExitFlag
    … …  
    254259  static double *cell_size;
    255260  std::string defaultName;
     261  class ThermoStatContainer *Thermostats;
    256262  int ExitFlag;
    257263public:

  • src/atom.cpp

    r6d574a r0d1ad0  
    159159  * \return true - \a *out present, false - \a *out is NULL
    160160 */
    161 bool atom::OutputArrayIndexed(ofstream * const out, const int *ElementNo, int *AtomNo, const char *comment) const
     161bool atom::OutputArrayIndexed(ostream * const out, const int *ElementNo, int *AtomNo, const char *comment) const
    162162{
    163163  AtomNo[type->Z]++;  // increment number
    … …  
    236236 * \param *AtomNo pointer to atom counter that is increased by one
    237237 */
    238 void atom::OutputMPQCLine(ofstream * const out, const Vector *center, int *AtomNo = NULL) const
     238void atom::OutputMPQCLine(ostream * const out, const Vector *center, int *AtomNo = NULL) const
    239239{
    240240  *out << "\t\t" << type->symbol << " [ " << x[0]-center->at(0) << "\t" << x[1]-center->at(1) << "\t" << x[2]-center->at(2) << " ]" << endl;

  • src/atom.hpp

    r6d574a r0d1ad0  
    5151
    5252  bool OutputIndexed(ofstream * const out, const int ElementNo, const int AtomNo, const char *comment = NULL) const;
    53   bool OutputArrayIndexed(ofstream * const out, const int *ElementNo, int *AtomNo, const char *comment = NULL) const;
     53  bool OutputArrayIndexed(ostream * const out, const int *ElementNo, int *AtomNo, const char *comment = NULL) const;
    5454  bool OutputXYZLine(ofstream *out) const;
    5555  bool OutputTrajectory(ofstream * const out, const int *ElementNo, int *AtomNo, const int step) const;
    5656  bool OutputTrajectoryXYZ(ofstream * const out, const int step) const;
    57   void OutputMPQCLine(ofstream * const out, const Vector *center, int *AtomNo) const;
     57  void OutputMPQCLine(ostream * const out, const Vector *center, int *AtomNo) const;
    5858
    5959  void InitComponentNr();

  • src/atom_trajectoryparticle.cpp

    r6d574a r0d1ad0  
    1515#include "log.hpp"
    1616#include "parser.hpp"
     17#include "ThermoStatContainer.hpp"
    1718#include "verbose.hpp"
    1819
    … …  
    197198void TrajectoryParticle::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
    198199{
    199   double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
     200  double sigma  = sqrt(configuration->Thermostats->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
    200201  Vector &U = Trajectory.U.at(Step);
    201202  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    202203    // throw a dice to determine whether it gets hit by a heat bath particle
    203     if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
     204    if (((((rand()/(double)RAND_MAX))*configuration->Thermostats->TempFrequency) < 1.)) {
    204205      DoLog(3) && (Log() << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ");
    205206      // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
    … …  
    225226  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    226227    for (int d=0; d<NDIM; d++) {
    227       U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
     228      U[d] *= sqrt(1+(configuration->Deltat/configuration->Thermostats->TempFrequency)*(ScaleTempFactor-1));
    228229      *ekin += 0.5*type->mass * U[d]*U[d];
    229230    }
    … …  
    255256  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
    256257    for (int d=0; d<NDIM; d++) {
    257         U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
     258        U[d] += configuration->Deltat/type->mass * (configuration->Thermostats->alpha * (U[d] * type->mass));
    258259        *ekin += (0.5*type->mass) * U[d]*U[d];
    259260      }

  • src/builder.cpp

    r6d574a r0d1ad0  
    11/** \file builder.cpp
    22 *
    3  * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
    4  * The output is the complete configuration file for PCP for direct use.
    5  * Features:
    6  * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
    7  * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
     3 *  date: Jan 1, 2007
     4 *  author: heber
    85 *
    96 */
    107
    11 /*! \mainpage Molecuilder - a molecular set builder
     8/*! \mainpage MoleCuilder - a molecular set builder
    129 *
    13  * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
     10 * This introductory shall briefly make acquainted with the program, helping in installing and a first run.
    1411 *
    1512 * \section about About the Program
    1613 *
    17  *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
    18  *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
    19  *  already constructed atoms.
     14 *  MoleCuilder is a program, written entirely in C++, that enables the construction of a coordinate set for the
     15 *  atoms making up an molecule. It allows for both building of simple molecules by adding atom-wise giving bond
     16 *  angles and distances or absolute coordinates, but also using them as templates. Regions can be specified and
     17 *  ordered to be filled with a molecule in a certain manner. Greater conglomerations of molecules can be tesselated
     18 *  and recognized as a region themselves to be subsequently surrounded by other (surface solvated) molecules.
     19 *  In the end, MoleCuilder allows the construction of arbitrary nano structures, whether they be crystalline or
     20 *  amorphic in nature.
    2021 *
    21  *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
    22  *  molecular dynamics implementation.
    2322 *
    2423 * \section install Installation
    … …  
    3231 *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
    3332 *  -# make doxygen-doc: Creates these html pages out of the documented source
     33 *  -# make check: Runs an extensive set of unit tests and a testsuite which also gives a good overview on the set of
     34 *                 functions.
    3435 *
    3536 * \section run Running
    … …  
    3738 *  The program can be executed by running: ./molecuilder
    3839 *
    39  *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
    40  *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
    41  *  later re-execution.
     40 *  MoleCuilder has three interfaces at your disposal:
     41 *  -# Textmenu: A simple interactive console-based menu, where awaits your choices and inputs in order to set atoms
     42 *               as you like
     43 *  -# CommandLineUI: Every command can also be chained up as a sequence of actions on the command line to be executed
     44 *               with any user interaction.
     45 *  -# GraphicalUI: A graphical user interface that also display the molecular structure being built and lots of other
     46 *               informations to ease the construction of bigger geometries.
    4247 *
    43  * \section ref References
    44  *
    45  *  For the special configuration file format, see the documentation of pcp.
     48 *  The supported output formats right now are:
     49 *  -# mpqc: Configuration files of the Massively Parallel Quantum Chemistry package (Sandia labs)
     50 *  -# pcp: Configuration files of the Parallel Car-Parrinello program (Institute for Numerical Simulation)
     51 *  -# tremolo: Configuration files of TREMOLO (Institute for Numerical Simulation)
     52 *  -# xyz: the most basic format for the 3d arrangement of atoms consisting of a list of element and 3 coordinates.
    4653 *
    4754 */
    … …  
    4956#include "Helpers/MemDebug.hpp"
    5057
    51 #include <boost/bind.hpp>
    52 
    53 using namespace std;
    54 
    55 #include <cstring>
    56 #include <cstdlib>
    57 
    58 #include "analysis_bonds.hpp"
    59 #include "analysis_correlation.hpp"
    60 #include "atom.hpp"
    61 #include "bond.hpp"
    6258#include "bondgraph.hpp"
    63 #include "boundary.hpp"
    6459#include "CommandLineParser.hpp"
    6560#include "config.hpp"
    66 #include "element.hpp"
    67 #include "ellipsoid.hpp"
    68 #include "helpers.hpp"
    69 #include "leastsquaremin.hpp"
    70 #include "linkedcell.hpp"
    7161#include "log.hpp"
    72 #include "memoryusageobserver.hpp"
    73 #include "molecule.hpp"
    74 #include "periodentafel.hpp"
     62#include "verbose.hpp"
     63#include "World.hpp"
     64
     65#include "Actions/ActionRegistry.hpp"
     66#include "Actions/ActionHistory.hpp"
     67#include "Actions/MapOfActions.hpp"
     68
     69#include "Parser/ChangeTracker.hpp"
     70#include "Parser/FormatParserStorage.hpp"
     71
    7572#include "UIElements/UIFactory.hpp"
    7673#include "UIElements/TextUI/TextUIFactory.hpp"
    … …  
    7875#include "UIElements/MainWindow.hpp"
    7976#include "UIElements/Dialog.hpp"
    80 #include "Menu/ActionMenuItem.hpp"
    81 #include "Actions/ActionRegistry.hpp"
    82 #include "Actions/ActionHistory.hpp"
    83 #include "Actions/MapOfActions.hpp"
    84 #include "Actions/MethodAction.hpp"
    85 #include "Actions/MoleculeAction/ChangeNameAction.hpp"
    86 #include "World.hpp"
     77
    8778#include "version.h"
    88 #include "World.hpp"
    8979
    90 
    91 /********************************************* Subsubmenu routine ************************************/
    92 #if 0
    93 /** Submenu for adding atoms to the molecule.
    94  * \param *periode periodentafel
    95  * \param *molecule molecules with atoms
    96  */
    97 static void AddAtoms(periodentafel *periode, molecule *mol)
    98 {
    99   atom *first, *second, *third, *fourth;
    100   Vector **atoms;
    101   Vector x,y,z,n;  // coordinates for absolute point in cell volume
    102   double a,b,c;
    103   char choice;  // menu choice char
    104   bool valid;
    105 
    106   cout << Verbose(0) << "===========ADD ATOM============================" << endl;
    107   cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
    108   cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
    109   cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
    110   cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
    111   cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
    112   cout << Verbose(0) << "all else - go back" << endl;
    113   cout << Verbose(0) << "===============================================" << endl;
    114   cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
    115   cout << Verbose(0) << "INPUT: ";
    116   cin >> choice;
    117 
    118   switch (choice) {
    119     default:
    120       DoeLog(2) && (eLog()<< Verbose(2) << "Not a valid choice." << endl);
    121       break;
    122       case 'a': // absolute coordinates of atom
    123         cout << Verbose(0) << "Enter absolute coordinates." << endl;
    124         first = new atom;
    125         first->x.AskPosition(World::getInstance().getDomain(), false);
    126         first->type = periode->AskElement();  // give type
    127         mol->AddAtom(first);  // add to molecule
    128         break;
    129 
    130       case 'b': // relative coordinates of atom wrt to reference point
    131         first = new atom;
    132         valid = true;
    133         do {
    134           if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
    135           cout << Verbose(0) << "Enter reference coordinates." << endl;
    136           x.AskPosition(World::getInstance().getDomain(), true);
    137           cout << Verbose(0) << "Enter relative coordinates." << endl;
    138           first->x.AskPosition(World::getInstance().getDomain(), false);
    139           first->x.AddVector((const Vector *)&x);
    140           cout << Verbose(0) << "\n";
    141         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    142         first->type = periode->AskElement();  // give type
    143         mol->AddAtom(first);  // add to molecule
    144         break;
    145 
    146       case 'c': // relative coordinates of atom wrt to already placed atom
    147         first = new atom;
    148         valid = true;
    149         do {
    150           if (!valid) DoeLog(2) && (eLog()<< Verbose(2) << "Resulting position out of cell." << endl);
    151           second = mol->AskAtom("Enter atom number: ");
    152           DoLog(0) && (Log() << Verbose(0) << "Enter relative coordinates." << endl);
    153           first->x.AskPosition(World::getInstance().getDomain(), false);
    154           for (int i=NDIM;i--;) {
    155             first->x.x[i] += second->x.x[i];
    156           }
    157         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    158         first->type = periode->AskElement();  // give type
    159         mol->AddAtom(first);  // add to molecule
    160         break;
    161 
    162     case 'd': // two atoms, two angles and a distance
    163         first = new atom;
    164         valid = true;
    165         do {
    166           if (!valid) {
    167             DoeLog(2) && (eLog()<< Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl);
    168           }
    169           cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
    170           second = mol->AskAtom("Enter central atom: ");
    171           third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
    172           fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
    173           a = ask_value("Enter distance between central (first) and new atom: ");
    174           b = ask_value("Enter angle between new, first and second atom (degrees): ");
    175           b *= M_PI/180.;
    176           bound(&b, 0., 2.*M_PI);
    177           c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
    178           c *= M_PI/180.;
    179           bound(&c, -M_PI, M_PI);
    180           cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
    181 /*
    182           second->Output(1,1,(ofstream *)&cout);
    183           third->Output(1,2,(ofstream *)&cout);
    184           fourth->Output(1,3,(ofstream *)&cout);
    185           n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
    186           x.Copyvector(&second->x);
    187           x.SubtractVector(&third->x);
    188           x.Copyvector(&fourth->x);
    189           x.SubtractVector(&third->x);
    190 
    191           if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
    192          coutg() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
    193             continue;
    194           }
    195           DoLog(0) && (Log() << Verbose(0) << "resulting relative coordinates: ");
    196           z.Output();
    197           DoLog(0) && (Log() << Verbose(0) << endl);
    198           */
    199           // calc axis vector
    200           x.CopyVector(&second->x);
    201           x.SubtractVector(&third->x);
    202           x.Normalize();
    203           Log() << Verbose(0) << "x: ",
    204           x.Output();
    205           DoLog(0) && (Log() << Verbose(0) << endl);
    206           z.MakeNormalVector(&second->x,&third->x,&fourth->x);
    207           Log() << Verbose(0) << "z: ",
    208           z.Output();
    209           DoLog(0) && (Log() << Verbose(0) << endl);
    210           y.MakeNormalVector(&x,&z);
    211           Log() << Verbose(0) << "y: ",
    212           y.Output();
    213           DoLog(0) && (Log() << Verbose(0) << endl);
    214 
    215           // rotate vector around first angle
    216           first->x.CopyVector(&x);
    217           first->x.RotateVector(&z,b - M_PI);
    218           Log() << Verbose(0) << "Rotated vector: ",
    219           first->x.Output();
    220           DoLog(0) && (Log() << Verbose(0) << endl);
    221           // remove the projection onto the rotation plane of the second angle
    222           n.CopyVector(&y);
    223           n.Scale(first->x.ScalarProduct(&y));
    224           Log() << Verbose(0) << "N1: ",
    225           n.Output();
    226           DoLog(0) && (Log() << Verbose(0) << endl);
    227           first->x.SubtractVector(&n);
    228           Log() << Verbose(0) << "Subtracted vector: ",
    229           first->x.Output();
    230           DoLog(0) && (Log() << Verbose(0) << endl);
    231           n.CopyVector(&z);
    232           n.Scale(first->x.ScalarProduct(&z));
    233           Log() << Verbose(0) << "N2: ",
    234           n.Output();
    235           DoLog(0) && (Log() << Verbose(0) << endl);
    236           first->x.SubtractVector(&n);
    237           Log() << Verbose(0) << "2nd subtracted vector: ",
    238           first->x.Output();
    239           DoLog(0) && (Log() << Verbose(0) << endl);
    240 
    241           // rotate another vector around second angle
    242           n.CopyVector(&y);
    243           n.RotateVector(&x,c - M_PI);
    244           Log() << Verbose(0) << "2nd Rotated vector: ",
    245           n.Output();
    246           DoLog(0) && (Log() << Verbose(0) << endl);
    247 
    248           // add the two linear independent vectors
    249           first->x.AddVector(&n);
    250           first->x.Normalize();
    251           first->x.Scale(a);
    252           first->x.AddVector(&second->x);
    253 
    254           DoLog(0) && (Log() << Verbose(0) << "resulting coordinates: ");
    255           first->x.Output();
    256           DoLog(0) && (Log() << Verbose(0) << endl);
    257         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    258         first->type = periode->AskElement();  // give type
    259         mol->AddAtom(first);  // add to molecule
    260         break;
    261 
    262       case 'e': // least square distance position to a set of atoms
    263         first = new atom;
    264         atoms = new (Vector*[128]);
    265         valid = true;
    266         for(int i=128;i--;)
    267           atoms[i] = NULL;
    268         int i=0, j=0;
    269         cout << Verbose(0) << "Now we need at least three molecules.\n";
    270         do {
    271           cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
    272           cin >> j;
    273           if (j != -1) {
    274             second = mol->FindAtom(j);
    275             atoms[i++] = &(second->x);
    276           }
    277         } while ((j != -1) && (i<128));
    278         if (i >= 2) {
    279           first->x.LSQdistance((const Vector **)atoms, i);
    280           first->x.Output();
    281           first->type = periode->AskElement();  // give type
    282           mol->AddAtom(first);  // add to molecule
    283         } else {
    284           delete first;
    285           cout << Verbose(0) << "Please enter at least two vectors!\n";
    286         }
    287         break;
    288   };
    289 };
    290 
    291 /** Submenu for centering the atoms in the molecule.
    292  * \param *mol molecule with all the atoms
    293  */
    294 static void CenterAtoms(molecule *mol)
    295 {
    296   Vector x, y, helper;
    297   char choice;  // menu choice char
    298 
    299   cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
    300   cout << Verbose(0) << " a - on origin" << endl;
    301   cout << Verbose(0) << " b - on center of gravity" << endl;
    302   cout << Verbose(0) << " c - within box with additional boundary" << endl;
    303   cout << Verbose(0) << " d - within given simulation box" << endl;
    304   cout << Verbose(0) << "all else - go back" << endl;
    305   cout << Verbose(0) << "===============================================" << endl;
    306   cout << Verbose(0) << "INPUT: ";
    307   cin >> choice;
    308 
    309   switch (choice) {
    310     default:
    311       cout << Verbose(0) << "Not a valid choice." << endl;
    312       break;
    313     case 'a':
    314       cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
    315       mol->CenterOrigin();
    316       break;
    317     case 'b':
    318       cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
    319       mol->CenterPeriodic();
    320       break;
    321     case 'c':
    322       cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
    323       for (int i=0;i<NDIM;i++) {
    324         cout << Verbose(0) << "Enter axis " << i << " boundary: ";
    325         cin >> y.x[i];
    326       }
    327       mol->CenterEdge(&x);  // make every coordinate positive
    328       mol->Center.AddVector(&y); // translate by boundary
    329       helper.CopyVector(&y);
    330       helper.Scale(2.);
    331       helper.AddVector(&x);
    332       mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
    333       break;
    334     case 'd':
    335       cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
    336       for (int i=0;i<NDIM;i++) {
    337         cout << Verbose(0) << "Enter axis " << i << " boundary: ";
    338         cin >> x.x[i];
    339       }
    340       // update Box of atoms by boundary
    341       mol->SetBoxDimension(&x);
    342       // center
    343       mol->CenterInBox();
    344       break;
    345   }
    346 };
    347 
    348 /** Submenu for aligning the atoms in the molecule.
    349  * \param *periode periodentafel
    350  * \param *mol molecule with all the atoms
    351  */
    352 static void AlignAtoms(periodentafel *periode, molecule *mol)
    353 {
    354   atom *first, *second, *third;
    355   Vector x,n;
    356   char choice;  // menu choice char
    357 
    358   cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
    359   cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
    360   cout << Verbose(0) << " b - state alignment vector" << endl;
    361   cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
    362   cout << Verbose(0) << " d - align automatically by least square fit" << endl;
    363   cout << Verbose(0) << "all else - go back" << endl;
    364   cout << Verbose(0) << "===============================================" << endl;
    365   cout << Verbose(0) << "INPUT: ";
    366   cin >> choice;
    367 
    368   switch (choice) {
    369     default:
    370     case 'a': // three atoms defining mirror plane
    371       first = mol->AskAtom("Enter first atom: ");
    372       second = mol->AskAtom("Enter second atom: ");
    373       third = mol->AskAtom("Enter third atom: ");
    374 
    375       n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
    376       break;
    377     case 'b': // normal vector of mirror plane
    378       cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    379       n.AskPosition(World::getInstance().getDomain(),0);
    380       n.Normalize();
    381       break;
    382     case 'c': // three atoms defining mirror plane
    383       first = mol->AskAtom("Enter first atom: ");
    384       second = mol->AskAtom("Enter second atom: ");
    385 
    386       n.CopyVector((const Vector *)&first->x);
    387       n.SubtractVector((const Vector *)&second->x);
    388       n.Normalize();
    389       break;
    390     case 'd':
    391       char shorthand[4];
    392       Vector a;
    393       struct lsq_params param;
    394       do {
    395         fprintf(stdout, "Enter the element of atoms to be chosen: ");
    396         fscanf(stdin, "%3s", shorthand);
    397       } while ((param.type = periode->FindElement(shorthand)) == NULL);
    398       cout << Verbose(0) << "Element is " << param.type->name << endl;
    399       mol->GetAlignvector(&param);
    400       for (int i=NDIM;i--;) {
    401         x.x[i] = gsl_vector_get(param.x,i);
    402         n.x[i] = gsl_vector_get(param.x,i+NDIM);
    403       }
    404       gsl_vector_free(param.x);
    405       cout << Verbose(0) << "Offset vector: ";
    406       x.Output();
    407       DoLog(0) && (Log() << Verbose(0) << endl);
    408       n.Normalize();
    409       break;
    410   };
    411   DoLog(0) && (Log() << Verbose(0) << "Alignment vector: ");
    412   n.Output();
    413   DoLog(0) && (Log() << Verbose(0) << endl);
    414   mol->Align(&n);
    415 };
    416 
    417 /** Submenu for mirroring the atoms in the molecule.
    418  * \param *mol molecule with all the atoms
    419  */
    420 static void MirrorAtoms(molecule *mol)
    421 {
    422   atom *first, *second, *third;
    423   Vector n;
    424   char choice;  // menu choice char
    425 
    426   DoLog(0) && (Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl);
    427   DoLog(0) && (Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl);
    428   DoLog(0) && (Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl);
    429   DoLog(0) && (Log() << Verbose(0) << " c - state two atoms in normal direction" << endl);
    430   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    431   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    432   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    433   cin >> choice;
    434 
    435   switch (choice) {
    436     default:
    437     case 'a': // three atoms defining mirror plane
    438       first = mol->AskAtom("Enter first atom: ");
    439       second = mol->AskAtom("Enter second atom: ");
    440       third = mol->AskAtom("Enter third atom: ");
    441 
    442       n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
    443       break;
    444     case 'b': // normal vector of mirror plane
    445       DoLog(0) && (Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl);
    446       n.AskPosition(World::getInstance().getDomain(),0);
    447       n.Normalize();
    448       break;
    449     case 'c': // three atoms defining mirror plane
    450       first = mol->AskAtom("Enter first atom: ");
    451       second = mol->AskAtom("Enter second atom: ");
    452 
    453       n.CopyVector((const Vector *)&first->x);
    454       n.SubtractVector((const Vector *)&second->x);
    455       n.Normalize();
    456       break;
    457   };
    458   DoLog(0) && (Log() << Verbose(0) << "Normal vector: ");
    459   n.Output();
    460   DoLog(0) && (Log() << Verbose(0) << endl);
    461   mol->Mirror((const Vector *)&n);
    462 };
    463 
    464 /** Submenu for removing the atoms from the molecule.
    465  * \param *mol molecule with all the atoms
    466  */
    467 static void RemoveAtoms(molecule *mol)
    468 {
    469   atom *first, *second;
    470   int axis;
    471   double tmp1, tmp2;
    472   char choice;  // menu choice char
    473 
    474   DoLog(0) && (Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl);
    475   DoLog(0) && (Log() << Verbose(0) << " a - state atom for removal by number" << endl);
    476   DoLog(0) && (Log() << Verbose(0) << " b - keep only in radius around atom" << endl);
    477   DoLog(0) && (Log() << Verbose(0) << " c - remove this with one axis greater value" << endl);
    478   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    479   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    480   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    481   cin >> choice;
    482 
    483   switch (choice) {
    484     default:
    485     case 'a':
    486       if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
    487         DoLog(1) && (Log() << Verbose(1) << "Atom removed." << endl);
    488       else
    489         DoLog(1) && (Log() << Verbose(1) << "Atom not found." << endl);
    490       break;
    491     case 'b':
    492       second = mol->AskAtom("Enter number of atom as reference point: ");
    493       DoLog(0) && (Log() << Verbose(0) << "Enter radius: ");
    494       cin >> tmp1;
    495       first = mol->start;
    496       second = first->next;
    497       while(second != mol->end) {
    498         first = second;
    499         second = first->next;
    500         if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
    501           mol->RemoveAtom(first);
    502       }
    503       break;
    504     case 'c':
    505       DoLog(0) && (Log() << Verbose(0) << "Which axis is it: ");
    506       cin >> axis;
    507       DoLog(0) && (Log() << Verbose(0) << "Lower boundary: ");
    508       cin >> tmp1;
    509       DoLog(0) && (Log() << Verbose(0) << "Upper boundary: ");
    510       cin >> tmp2;
    511       first = mol->start;
    512       second = first->next;
    513       while(second != mol->end) {
    514         first = second;
    515         second = first->next;
    516         if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
    517           //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
    518           mol->RemoveAtom(first);
    519         }
    520       }
    521       break;
    522   };
    523   //mol->Output();
    524   choice = 'r';
    525 };
    526 
    527 /** Submenu for measuring out the atoms in the molecule.
    528  * \param *periode periodentafel
    529  * \param *mol molecule with all the atoms
    530  */
    531 static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
    532 {
    533   atom *first, *second, *third;
    534   Vector x,y;
    535   double min[256], tmp1, tmp2, tmp3;
    536   int Z;
    537   char choice;  // menu choice char
    538 
    539   DoLog(0) && (Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl);
    540   DoLog(0) && (Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl);
    541   DoLog(0) && (Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl);
    542   DoLog(0) && (Log() << Verbose(0) << " c - calculate bond angle" << endl);
    543   DoLog(0) && (Log() << Verbose(0) << " d - calculate principal axis of the system" << endl);
    544   DoLog(0) && (Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl);
    545   DoLog(0) && (Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl);
    546   DoLog(0) && (Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl);
    547   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    548   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    549   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    550   cin >> choice;
    551 
    552   switch(choice) {
    553     default:
    554       DoLog(1) && (Log() << Verbose(1) << "Not a valid choice." << endl);
    555       break;
    556     case 'a':
    557       first = mol->AskAtom("Enter first atom: ");
    558       for (int i=MAX_ELEMENTS;i--;)
    559         min[i] = 0.;
    560 
    561       second = mol->start;
    562       while ((second->next != mol->end)) {
    563         second = second->next; // advance
    564         Z = second->type->Z;
    565         tmp1 = 0.;
    566         if (first != second) {
    567           x.CopyVector((const Vector *)&first->x);
    568           x.SubtractVector((const Vector *)&second->x);
    569           tmp1 = x.Norm();
    570         }
    571         if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
    572         //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
    573       }
    574       for (int i=MAX_ELEMENTS;i--;)
    575         if (min[i] != 0.) Log() << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
    576       break;
    577 
    578     case 'b':
    579       first = mol->AskAtom("Enter first atom: ");
    580       second = mol->AskAtom("Enter second atom: ");
    581       for (int i=NDIM;i--;)
    582         min[i] = 0.;
    583       x.CopyVector((const Vector *)&first->x);
    584       x.SubtractVector((const Vector *)&second->x);
    585       tmp1 = x.Norm();
    586       DoLog(1) && (Log() << Verbose(1) << "Distance vector is ");
    587       x.Output();
    588       DoLog(0) && (Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl);
    589       break;
    590 
    591     case 'c':
    592       DoLog(0) && (Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl);
    593       first = mol->AskAtom("Enter first atom: ");
    594       second = mol->AskAtom("Enter central atom: ");
    595       third  = mol->AskAtom("Enter last atom: ");
    596       tmp1 = tmp2 = tmp3 = 0.;
    597       x.CopyVector((const Vector *)&first->x);
    598       x.SubtractVector((const Vector *)&second->x);
    599       y.CopyVector((const Vector *)&third->x);
    600       y.SubtractVector((const Vector *)&second->x);
    601       DoLog(0) && (Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ");
    602       DoLog(0) && (Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl);
    603       break;
    604     case 'd':
    605       DoLog(0) && (Log() << Verbose(0) << "Evaluating prinicipal axis." << endl);
    606       DoLog(0) && (Log() << Verbose(0) << "Shall we rotate? [0/1]: ");
    607       cin >> Z;
    608       if ((Z >=0) && (Z <=1))
    609         mol->PrincipalAxisSystem((bool)Z);
    610       else
    611         mol->PrincipalAxisSystem(false);
    612       break;
    613     case 'e':
    614       {
    615         DoLog(0) && (Log() << Verbose(0) << "Evaluating volume of the convex envelope.");
    616         class Tesselation *TesselStruct = NULL;
    617         const LinkedCell *LCList = NULL;
    618         LCList = new LinkedCell(mol, 10.);
    619         FindConvexBorder(mol, TesselStruct, LCList, NULL);
    620         double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
    621         DoLog(0) && (Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl);\
    622         delete(LCList);
    623         delete(TesselStruct);
    624       }
    625       break;
    626     case 'f':
    627       mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps);
    628       break;
    629     case 'g':
    630       {
    631         char filename[255];
    632         DoLog(0) && (Log() << Verbose(0) << "Please enter filename: " << endl);
    633         cin >> filename;
    634         DoLog(1) && (Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl);
    635         ofstream *output = new ofstream(filename, ios::trunc);
    636         if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
    637           DoLog(2) && (Log() << Verbose(2) << "File could not be written." << endl);
    638         else
    639           DoLog(2) && (Log() << Verbose(2) << "File stored." << endl);
    640         output->close();
    641         delete(output);
    642       }
    643       break;
    644   }
    645 };
    646 
    647 /** Submenu for measuring out the atoms in the molecule.
    648  * \param *mol molecule with all the atoms
    649  * \param *configuration configuration structure for the to be written config files of all fragments
    650  */
    651 static void FragmentAtoms(molecule *mol, config *configuration)
    652 {
    653   int Order1;
    654   clock_t start, end;
    655 
    656   DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl);
    657   DoLog(0) && (Log() << Verbose(0) << "What's the desired bond order: ");
    658   cin >> Order1;
    659   if (mol->first->next != mol->last) {  // there are bonds
    660     start = clock();
    661     mol->FragmentMolecule(Order1, configuration);
    662     end = clock();
    663     DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl);
    664   } else
    665     DoLog(0) && (Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl);
    666 };
    667 
    668 /********************************************** Submenu routine **************************************/
    669 
    670 /** Submenu for manipulating atoms.
    671  * \param *periode periodentafel
    672  * \param *molecules list of molecules whose atoms are to be manipulated
    673  */
    674 static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
    675 {
    676   atom *first, *second, *third;
    677   molecule *mol = NULL;
    678   Vector x,y,z,n; // coordinates for absolute point in cell volume
    679   double *factor; // unit factor if desired
    680   double bond, minBond;
    681   char choice;  // menu choice char
    682   bool valid;
    683 
    684   DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl);
    685   DoLog(0) && (Log() << Verbose(0) << "a - add an atom" << endl);
    686   DoLog(0) && (Log() << Verbose(0) << "r - remove an atom" << endl);
    687   DoLog(0) && (Log() << Verbose(0) << "b - scale a bond between atoms" << endl);
    688   DoLog(0) && (Log() << Verbose(0) << "t - turn an atom round another bond" << endl);
    689   DoLog(0) && (Log() << Verbose(0) << "u - change an atoms element" << endl);
    690   DoLog(0) && (Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl);
    691   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    692   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    693   if (molecules->NumberOfActiveMolecules() > 1)
    694     DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
    695   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    696   cin >> choice;
    697 
    698   switch (choice) {
    699     default:
    700       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    701       break;
    702 
    703     case 'a': // add atom
    704       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    705         if ((*ListRunner)->ActiveFlag) {
    706         mol = *ListRunner;
    707         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    708         AddAtoms(periode, mol);
    709       }
    710       break;
    711 
    712     case 'b': // scale a bond
    713       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    714         if ((*ListRunner)->ActiveFlag) {
    715         mol = *ListRunner;
    716         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    717         DoLog(0) && (Log() << Verbose(0) << "Scaling bond length between two atoms." << endl);
    718         first = mol->AskAtom("Enter first (fixed) atom: ");
    719         second = mol->AskAtom("Enter second (shifting) atom: ");
    720         minBond = 0.;
    721         for (int i=NDIM;i--;)
    722           minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
    723         minBond = sqrt(minBond);
    724         DoLog(0) && (Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl);
    725         DoLog(0) && (Log() << Verbose(0) << "Enter new bond length [a.u.]: ");
    726         cin >> bond;
    727         for (int i=NDIM;i--;) {
    728           second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
    729         }
    730         //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
    731         //second->Output(second->type->No, 1);
    732       }
    733       break;
    734 
    735     case 'c': // unit scaling of the metric
    736       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    737         if ((*ListRunner)->ActiveFlag) {
    738         mol = *ListRunner;
    739         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    740        DoLog(0) && (Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl);
    741        DoLog(0) && (Log() << Verbose(0) << "Enter three factors: ");
    742        factor = new double[NDIM];
    743        cin >> factor[0];
    744        cin >> factor[1];
    745        cin >> factor[2];
    746        valid = true;
    747        mol->Scale((const double ** const)&factor);
    748        delete[](factor);
    749       }
    750      break;
    751 
    752     case 'l': // measure distances or angles
    753       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    754         if ((*ListRunner)->ActiveFlag) {
    755         mol = *ListRunner;
    756         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    757         MeasureAtoms(periode, mol, configuration);
    758       }
    759       break;
    760 
    761     case 'r': // remove atom
    762       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    763         if ((*ListRunner)->ActiveFlag) {
    764         mol = *ListRunner;
    765         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    766         RemoveAtoms(mol);
    767       }
    768       break;
    769 
    770     case 't': // turn/rotate atom
    771       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    772         if ((*ListRunner)->ActiveFlag) {
    773           mol = *ListRunner;
    774           DoLog(0) && (Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl);
    775           first = mol->AskAtom("Enter turning atom: ");
    776           second = mol->AskAtom("Enter central atom: ");
    777           third  = mol->AskAtom("Enter bond atom: ");
    778           cout << Verbose(0) << "Enter new angle in degrees: ";
    779           double tmp = 0.;
    780           cin >> tmp;
    781           // calculate old angle
    782           x.CopyVector((const Vector *)&first->x);
    783           x.SubtractVector((const Vector *)&second->x);
    784           y.CopyVector((const Vector *)&third->x);
    785           y.SubtractVector((const Vector *)&second->x);
    786           double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
    787           cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
    788           cout << Verbose(0) << alpha << " degrees" << endl;
    789           // rotate
    790           z.MakeNormalVector(&x,&y);
    791           x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
    792           x.AddVector(&second->x);
    793           first->x.CopyVector(&x);
    794           // check new angle
    795           x.CopyVector((const Vector *)&first->x);
    796           x.SubtractVector((const Vector *)&second->x);
    797           alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
    798           cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
    799           cout << Verbose(0) << alpha << " degrees" << endl;
    800         }
    801       break;
    802 
    803     case 'u': // change an atom's element
    804       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    805         if ((*ListRunner)->ActiveFlag) {
    806         int Z;
    807         mol = *ListRunner;
    808         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    809         first = NULL;
    810         do {
    811           DoLog(0) && (Log() << Verbose(0) << "Change the element of which atom: ");
    812           cin >> Z;
    813         } while ((first = mol->FindAtom(Z)) == NULL);
    814         DoLog(0) && (Log() << Verbose(0) << "New element by atomic number Z: ");
    815         cin >> Z;
    816         first->type = periode->FindElement(Z);
    817         DoLog(0) && (Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl);
    818       }
    819       break;
    820   }
    821 };
    822 
    823 /** Submenu for manipulating molecules.
    824  * \param *periode periodentafel
    825  * \param *molecules list of molecule to manipulate
    826  */
    827 static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
    828 {
    829   atom *first = NULL;
    830   Vector x,y,z,n; // coordinates for absolute point in cell volume
    831   int j, axis, count, faktor;
    832   char choice;  // menu choice char
    833   molecule *mol = NULL;
    834   element **Elements;
    835   Vector **vectors;
    836   MoleculeLeafClass *Subgraphs = NULL;
    837 
    838   DoLog(0) && (Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl);
    839   DoLog(0) && (Log() << Verbose(0) << "c - scale by unit transformation" << endl);
    840   DoLog(0) && (Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl);
    841   DoLog(0) && (Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl);
    842   DoLog(0) && (Log() << Verbose(0) << "g - center atoms in box" << endl);
    843   DoLog(0) && (Log() << Verbose(0) << "i - realign molecule" << endl);
    844   DoLog(0) && (Log() << Verbose(0) << "m - mirror all molecules" << endl);
    845   DoLog(0) && (Log() << Verbose(0) << "o - create connection matrix" << endl);
    846   DoLog(0) && (Log() << Verbose(0) << "t - translate molecule by vector" << endl);
    847   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    848   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    849   if (molecules->NumberOfActiveMolecules() > 1)
    850     DoeLog(2) && (eLog()<< Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl);
    851   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    852   cin >> choice;
    853 
    854   switch (choice) {
    855     default:
    856       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    857       break;
    858 
    859     case 'd': // duplicate the periodic cell along a given axis, given times
    860       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    861         if ((*ListRunner)->ActiveFlag) {
    862         mol = *ListRunner;
    863         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    864         DoLog(0) && (Log() << Verbose(0) << "State the axis [(+-)123]: ");
    865         cin >> axis;
    866         DoLog(0) && (Log() << Verbose(0) << "State the factor: ");
    867         cin >> faktor;
    868 
    869         mol->CountAtoms(); // recount atoms
    870         if (mol->getAtomCount() != 0) {  // if there is more than none
    871           count = mol->getAtomCount();  // is changed becausing of adding, thus has to be stored away beforehand
    872           Elements = new element *[count];
    873           vectors = new Vector *[count];
    874           j = 0;
    875           first = mol->start;
    876           while (first->next != mol->end) { // make a list of all atoms with coordinates and element
    877             first = first->next;
    878             Elements[j] = first->type;
    879             vectors[j] = &first->x;
    880             j++;
    881           }
    882           if (count != j)
    883             DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
    884           x.Zero();
    885           y.Zero();
    886           y.x[abs(axis)-1] = World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
    887           for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    888             x.AddVector(&y); // per factor one cell width further
    889             for (int k=count;k--;) { // go through every atom of the original cell
    890               first = new atom(); // create a new body
    891               first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    892               first->x.AddVector(&x);     // translate the coordinates
    893               first->type = Elements[k];  // insert original element
    894               mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
    895             }
    896           }
    897           if (mol->first->next != mol->last) // if connect matrix is present already, redo it
    898             mol->CreateAdjacencyList(mol->BondDistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
    899           // free memory
    900           delete[](Elements);
    901           delete[](vectors);
    902           // correct cell size
    903           if (axis < 0) { // if sign was negative, we have to translate everything
    904             x.Zero();
    905             x.AddVector(&y);
    906             x.Scale(-(faktor-1));
    907             mol->Translate(&x);
    908           }
    909           World::getInstance().getDomain()[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    910         }
    911       }
    912       break;
    913 
    914     case 'f':
    915       FragmentAtoms(mol, configuration);
    916       break;
    917 
    918     case 'g': // center the atoms
    919       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    920         if ((*ListRunner)->ActiveFlag) {
    921         mol = *ListRunner;
    922         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    923         CenterAtoms(mol);
    924       }
    925       break;
    926 
    927     case 'i': // align all atoms
    928       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    929         if ((*ListRunner)->ActiveFlag) {
    930         mol = *ListRunner;
    931         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    932         AlignAtoms(periode, mol);
    933       }
    934       break;
    935 
    936     case 'm': // mirror atoms along a given axis
    937       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    938         if ((*ListRunner)->ActiveFlag) {
    939         mol = *ListRunner;
    940         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    941         MirrorAtoms(mol);
    942       }
    943       break;
    944 
    945     case 'o': // create the connection matrix
    946       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    947         if ((*ListRunner)->ActiveFlag) {
    948           mol = *ListRunner;
    949           double bonddistance;
    950           clock_t start,end;
    951           DoLog(0) && (Log() << Verbose(0) << "What's the maximum bond distance: ");
    952           cin >> bonddistance;
    953           start = clock();
    954           mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
    955           end = clock();
    956           DoLog(0) && (Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl);
    957         }
    958       break;
    959 
    960     case 't': // translate all atoms
    961       for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    962         if ((*ListRunner)->ActiveFlag) {
    963         mol = *ListRunner;
    964         DoLog(0) && (Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl);
    965         DoLog(0) && (Log() << Verbose(0) << "Enter translation vector." << endl);
    966         x.AskPosition(World::getInstance().getDomain(),0);
    967         mol->Center.AddVector((const Vector *)&x);
    968      }
    969      break;
    970   }
    971   // Free all
    972   if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
    973     while (Subgraphs->next != NULL) {
    974       Subgraphs = Subgraphs->next;
    975       delete(Subgraphs->previous);
    976     }
    977     delete(Subgraphs);
    978   }
    979 };
    980 
    981 
    982 /** Submenu for creating new molecules.
    983  * \param *periode periodentafel
    984  * \param *molecules list of molecules to add to
    985  */
    986 static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
    987 {
    988   char choice;  // menu choice char
    989   Vector center;
    990   int nr, count;
    991   molecule *mol = NULL;
    992 
    993   DoLog(0) && (Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl);
    994   DoLog(0) && (Log() << Verbose(0) << "c - create new molecule" << endl);
    995   DoLog(0) && (Log() << Verbose(0) << "l - load molecule from xyz file" << endl);
    996   DoLog(0) && (Log() << Verbose(0) << "n - change molecule's name" << endl);
    997   DoLog(0) && (Log() << Verbose(0) << "N - give molecules filename" << endl);
    998   DoLog(0) && (Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl);
    999   DoLog(0) && (Log() << Verbose(0) << "r - remove a molecule" << endl);
    1000   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    1001   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    1002   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    1003   cin >> choice;
    1004 
    1005   switch (choice) {
    1006     default:
    1007       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    1008       break;
    1009     case 'c':
    1010       mol = World::getInstance().createMolecule();
    1011       molecules->insert(mol);
    1012       break;
    1013 
    1014     case 'l': // load from XYZ file
    1015       {
    1016         char filename[MAXSTRINGSIZE];
    1017         DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl);
    1018         mol = World::getInstance().createMolecule();
    1019         do {
    1020           DoLog(0) && (Log() << Verbose(0) << "Enter file name: ");
    1021           cin >> filename;
    1022         } while (!mol->AddXYZFile(filename));
    1023         mol->SetNameFromFilename(filename);
    1024         // center at set box dimensions
    1025         mol->CenterEdge(&center);
    1026         double * const cell_size = World::getInstance().getDomain();
    1027         cell_size[0] = center.x[0];
    1028         cell_size[1] = 0;
    1029         cell_size[2] = center.x[1];
    1030         cell_size[3] = 0;
    1031         cell_size[4] = 0;
    1032         cell_size[5] = center.x[2];
    1033         molecules->insert(mol);
    1034       }
    1035       break;
    1036 
    1037     case 'n':
    1038       {
    1039         char filename[MAXSTRINGSIZE];
    1040         do {
    1041           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1042           cin >> nr;
    1043           mol = molecules->ReturnIndex(nr);
    1044         } while (mol == NULL);
    1045         DoLog(0) && (Log() << Verbose(0) << "Enter name: ");
    1046         cin >> filename;
    1047         strcpy(mol->name, filename);
    1048       }
    1049       break;
    1050 
    1051     case 'N':
    1052       {
    1053         char filename[MAXSTRINGSIZE];
    1054         do {
    1055           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1056           cin >> nr;
    1057           mol = molecules->ReturnIndex(nr);
    1058         } while (mol == NULL);
    1059         DoLog(0) && (Log() << Verbose(0) << "Enter name: ");
    1060         cin >> filename;
    1061         mol->SetNameFromFilename(filename);
    1062       }
    1063       break;
    1064 
    1065     case 'p': // parse XYZ file
    1066       {
    1067         char filename[MAXSTRINGSIZE];
    1068         mol = NULL;
    1069         do {
    1070           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1071           cin >> nr;
    1072           mol = molecules->ReturnIndex(nr);
    1073         } while (mol == NULL);
    1074         DoLog(0) && (Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl);
    1075         do {
    1076           DoLog(0) && (Log() << Verbose(0) << "Enter file name: ");
    1077           cin >> filename;
    1078         } while (!mol->AddXYZFile(filename));
    1079         mol->SetNameFromFilename(filename);
    1080       }
    1081       break;
    1082 
    1083     case 'r':
    1084       DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule: ");
    1085       cin >> nr;
    1086       count = 1;
    1087       for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    1088         if (nr == (*ListRunner)->IndexNr) {
    1089           mol = *ListRunner;
    1090           molecules->ListOfMolecules.erase(ListRunner);
    1091           delete(mol);
    1092           break;
    1093         }
    1094       break;
    1095   }
    1096 };
    1097 
    1098 
    1099 /** Submenu for merging molecules.
    1100  * \param *periode periodentafel
    1101  * \param *molecules list of molecules to add to
    1102  */
    1103 static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
    1104 {
    1105   char choice;  // menu choice char
    1106 
    1107   DoLog(0) && (Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl);
    1108   DoLog(0) && (Log() << Verbose(0) << "a - simple add of one molecule to another" << endl);
    1109   DoLog(0) && (Log() << Verbose(0) << "b - count the number of bonds of two elements" << endl);
    1110   DoLog(0) && (Log() << Verbose(0) << "B - count the number of bonds of three elements " << endl);
    1111   DoLog(0) && (Log() << Verbose(0) << "e - embedding merge of two molecules" << endl);
    1112   DoLog(0) && (Log() << Verbose(0) << "h - count the number of hydrogen bonds" << endl);
    1113   DoLog(0) && (Log() << Verbose(0) << "b - count the number of hydrogen bonds" << endl);
    1114   DoLog(0) && (Log() << Verbose(0) << "m - multi-merge of all molecules" << endl);
    1115   DoLog(0) && (Log() << Verbose(0) << "s - scatter merge of two molecules" << endl);
    1116   DoLog(0) && (Log() << Verbose(0) << "t - simple merge of two molecules" << endl);
    1117   DoLog(0) && (Log() << Verbose(0) << "all else - go back" << endl);
    1118   DoLog(0) && (Log() << Verbose(0) << "===============================================" << endl);
    1119   DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    1120   cin >> choice;
    1121 
    1122   switch (choice) {
    1123     default:
    1124       DoLog(0) && (Log() << Verbose(0) << "Not a valid choice." << endl);
    1125       break;
    1126 
    1127     case 'a':
    1128       {
    1129         int src, dest;
    1130         molecule *srcmol = NULL, *destmol = NULL;
    1131         {
    1132           do {
    1133             DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: ");
    1134             cin >> dest;
    1135             destmol = molecules->ReturnIndex(dest);
    1136           } while ((destmol == NULL) && (dest != -1));
    1137           do {
    1138             DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to add from: ");
    1139             cin >> src;
    1140             srcmol = molecules->ReturnIndex(src);
    1141           } while ((srcmol == NULL) && (src != -1));
    1142           if ((src != -1) && (dest != -1))
    1143             molecules->SimpleAdd(srcmol, destmol);
    1144         }
    1145       }
    1146       break;
    1147 
    1148     case 'b':
    1149       {
    1150         const int nr = 2;
    1151         char *names[nr] = {"first", "second"};
    1152         int Z[nr];
    1153         element *elements[nr];
    1154         for (int i=0;i<nr;i++) {
    1155           Z[i] = 0;
    1156           do {
    1157             cout << "Enter " << names[i] << " element: ";
    1158             cin >> Z[i];
    1159           } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
    1160           elements[i] = periode->FindElement(Z[i]);
    1161         }
    1162         const int count = CountBondsOfTwo(molecules, elements[0], elements[1]);
    1163         cout << endl << "There are " << count << " ";
    1164         for (int i=0;i<nr;i++) {
    1165           if (i==0)
    1166             cout << elements[i]->symbol;
    1167           else
    1168             cout << "-" << elements[i]->symbol;
    1169         }
    1170         cout << " bonds." << endl;
    1171       }
    1172     break;
    1173 
    1174     case 'B':
    1175       {
    1176         const int nr = 3;
    1177         char *names[nr] = {"first", "second", "third"};
    1178         int Z[nr];
    1179         element *elements[nr];
    1180         for (int i=0;i<nr;i++) {
    1181           Z[i] = 0;
    1182           do {
    1183             cout << "Enter " << names[i] << " element: ";
    1184             cin >> Z[i];
    1185           } while ((Z[i] <= 0) && (Z[i] > MAX_ELEMENTS));
    1186           elements[i] = periode->FindElement(Z[i]);
    1187         }
    1188         const int count = CountBondsOfThree(molecules, elements[0], elements[1], elements[2]);
    1189         cout << endl << "There are " << count << " ";
    1190         for (int i=0;i<nr;i++) {
    1191           if (i==0)
    1192             cout << elements[i]->symbol;
    1193           else
    1194             cout << "-" << elements[i]->symbol;
    1195         }
    1196         cout << " bonds." << endl;
    1197       }
    1198     break;
    1199 
    1200     case 'e':
    1201       {
    1202         int src, dest;
    1203         molecule *srcmol = NULL, *destmol = NULL;
    1204         do {
    1205           DoLog(0) && (Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ");
    1206           cin >> src;
    1207           srcmol = molecules->ReturnIndex(src);
    1208         } while ((srcmol == NULL) && (src != -1));
    1209         do {
    1210           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ");
    1211           cin >> dest;
    1212           destmol = molecules->ReturnIndex(dest);
    1213         } while ((destmol == NULL) && (dest != -1));
    1214         if ((src != -1) && (dest != -1))
    1215           molecules->EmbedMerge(destmol, srcmol);
    1216       }
    1217       break;
    1218 
    1219     case 'h':
    1220       {
    1221         int Z;
    1222         cout << "Please enter interface element: ";
    1223         cin >> Z;
    1224         element * const InterfaceElement = periode->FindElement(Z);
    1225         cout << endl << "There are " << CountHydrogenBridgeBonds(molecules, InterfaceElement) << " hydrogen bridges with connections to " << (InterfaceElement != 0 ? InterfaceElement->name : "None") << "." << endl;
    1226       }
    1227       break;
    1228 
    1229     case 'm':
    1230       {
    1231         int nr;
    1232         molecule *mol = NULL;
    1233         do {
    1234           DoLog(0) && (Log() << Verbose(0) << "Enter index of molecule to merge into: ");
    1235           cin >> nr;
    1236           mol = molecules->ReturnIndex(nr);
    1237         } while ((mol == NULL) && (nr != -1));
    1238         if (nr != -1) {
    1239           int N = molecules->ListOfMolecules.size()-1;
    1240           int *src = new int(N);
    1241           for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    1242             if ((*ListRunner)->IndexNr != nr)
    1243               src[N++] = (*ListRunner)->IndexNr;       
    1244           molecules->SimpleMultiMerge(mol, src, N);
    1245           delete[](src);
    1246         }
    1247       }
    1248       break;
    1249 
    1250     case 's':
    1251       DoLog(0) && (Log() << Verbose(0) << "Not implemented yet." << endl);
    1252       break;
    1253 
    1254     case 't':
    1255       {
    1256         int src, dest;
    1257         molecule *srcmol = NULL, *destmol = NULL;
    1258         {
    1259           do {
    1260             DoLog(0) && (Log() << Verbose(0) << "Enter index of destination molecule: ");
    1261             cin >> dest;
    1262             destmol = molecules->ReturnIndex(dest);
    1263           } while ((destmol == NULL) && (dest != -1));
    1264           do {
    1265             DoLog(0) && (Log() << Verbose(0) << "Enter index of source molecule to merge into: ");
    1266             cin >> src;
    1267             srcmol = molecules->ReturnIndex(src);
    1268           } while ((srcmol == NULL) && (src != -1));
    1269           if ((src != -1) && (dest != -1))
    1270             molecules->SimpleMerge(srcmol, destmol);
    1271         }
    1272       }
    1273       break;
    1274   }
    1275 };
    1276 
    1277 /********************************************** Test routine **************************************/
    1278 
    1279 /** Is called always as option 'T' in the menu.
    1280  * \param *molecules list of molecules
    1281  */
    1282 static void testroutine(MoleculeListClass *molecules)
    1283 {
    1284   // the current test routine checks the functionality of the KeySet&Graph concept:
    1285   // We want to have a multiindex (the KeySet) describing a unique subgraph
    1286   int i, comp, counter=0;
    1287 
    1288   // create a clone
    1289   molecule *mol = NULL;
    1290   if (molecules->ListOfMolecules.size() != 0) // clone
    1291     mol = (molecules->ListOfMolecules.front())->CopyMolecule();
    1292   else {
    1293     DoeLog(0) && (eLog()<< Verbose(0) << "I don't have anything to test on ... ");
    1294     performCriticalExit();
    1295     return;
    1296   }
    1297   atom *Walker = mol->start;
    1298 
    1299   // generate some KeySets
    1300   DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl);
    1301   KeySet TestSets[mol->getAtomCount()+1];
    1302   i=1;
    1303   while (Walker->next != mol->end) {
    1304     Walker = Walker->next;
    1305     for (int j=0;j<i;j++) {
    1306       TestSets[j].insert(Walker->nr);
    1307     }
    1308     i++;
    1309   }
    1310   DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl);
    1311   KeySetTestPair test;
    1312   test = TestSets[mol->getAtomCount()-1].insert(Walker->nr);
    1313   if (test.second) {
    1314     DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
    1315   } else {
    1316     DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl);
    1317   }
    1318   TestSets[mol->getAtomCount()].insert(mol->end->previous->nr);
    1319   TestSets[mol->getAtomCount()].insert(mol->end->previous->previous->previous->nr);
    1320 
    1321   // constructing Graph structure
    1322   DoLog(0) && (Log() << Verbose(0) << "Generating Subgraph class." << endl);
    1323   Graph Subgraphs;
    1324 
    1325   // insert KeySets into Subgraphs
    1326   DoLog(0) && (Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl);
    1327   for (int j=0;j<mol->getAtomCount();j++) {
    1328     Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
    1329   }
    1330   DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl);
    1331   GraphTestPair test2;
    1332   test2 = Subgraphs.insert(GraphPair (TestSets[mol->getAtomCount()],pair<int, double>(counter++, 1.)));
    1333   if (test2.second) {
    1334     DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
    1335   } else {
    1336     DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl);
    1337   }
    1338 
    1339   // show graphs
    1340   DoLog(0) && (Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl);
    1341   Graph::iterator A = Subgraphs.begin();
    1342   while (A !=  Subgraphs.end()) {
    1343     DoLog(0) && (Log() << Verbose(0) << (*A).second.first << ": ");
    1344     KeySet::iterator key = (*A).first.begin();
    1345     comp = -1;
    1346     while (key != (*A).first.end()) {
    1347       if ((*key) > comp)
    1348         DoLog(0) && (Log() << Verbose(0) << (*key) << " ");
    1349       else
    1350         DoLog(0) && (Log() << Verbose(0) << (*key) << "! ");
    1351       comp = (*key);
    1352       key++;
    1353     }
    1354     DoLog(0) && (Log() << Verbose(0) << endl);
    1355     A++;
    1356   }
    1357   delete(mol);
    1358 };
    1359 
    1360 
    1361 /** Tries given filename or standard on saving the config file.
    1362  * \param *ConfigFileName name of file
    1363  * \param *configuration pointer to configuration structure with all the values
    1364  * \param *periode pointer to periodentafel structure with all the elements
    1365  * \param *molecules list of molecules structure with all the atoms and coordinates
    1366  */
    1367 static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
    1368 {
    1369   char filename[MAXSTRINGSIZE];
    1370   ofstream output;
    1371   molecule *mol = World::getInstance().createMolecule();
    1372   mol->SetNameFromFilename(ConfigFileName);
    1373 
    1374   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    1375     DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
    1376   }
    1377 
    1378 
    1379   // first save as PDB data
    1380   if (ConfigFileName != NULL)
    1381     strcpy(filename, ConfigFileName);
    1382   if (output == NULL)
    1383     strcpy(filename,"main_pcp_linux");
    1384   DoLog(0) && (Log() << Verbose(0) << "Saving as pdb input ");
    1385   if (configuration->SavePDB(filename, molecules))
    1386     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    1387   else
    1388     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1389 
    1390   // then save as tremolo data file
    1391   if (ConfigFileName != NULL)
    1392     strcpy(filename, ConfigFileName);
    1393   if (output == NULL)
    1394     strcpy(filename,"main_pcp_linux");
    1395   DoLog(0) && (Log() << Verbose(0) << "Saving as tremolo data input ");
    1396   if (configuration->SaveTREMOLO(filename, molecules))
    1397     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    1398   else
    1399     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1400 
    1401   // translate each to its center and merge all molecules in MoleculeListClass into this molecule
    1402   int N = molecules->ListOfMolecules.size();
    1403   int *src = new int[N];
    1404   N=0;
    1405   for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
    1406     src[N++] = (*ListRunner)->IndexNr;
    1407     (*ListRunner)->Translate(&(*ListRunner)->Center);
    1408   }
    1409   molecules->SimpleMultiAdd(mol, src, N);
    1410   delete[](src);
    1411 
    1412   // ... and translate back
    1413   for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
    1414     (*ListRunner)->Center.Scale(-1.);
    1415     (*ListRunner)->Translate(&(*ListRunner)->Center);
    1416     (*ListRunner)->Center.Scale(-1.);
    1417   }
    1418 
    1419   DoLog(0) && (Log() << Verbose(0) << "Storing configuration ... " << endl);
    1420   // get correct valence orbitals
    1421   mol->CalculateOrbitals(*configuration);
    1422   configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
    1423   if (ConfigFileName != NULL) { // test the file name
    1424     strcpy(filename, ConfigFileName);
    1425     output.open(filename, ios::trunc);
    1426   } else if (strlen(configuration->configname) != 0) {
    1427     strcpy(filename, configuration->configname);
    1428     output.open(configuration->configname, ios::trunc);
    1429     } else {
    1430       strcpy(filename, DEFAULTCONFIG);
    1431       output.open(DEFAULTCONFIG, ios::trunc);
    1432     }
    1433   output.close();
    1434   output.clear();
    1435   DoLog(0) && (Log() << Verbose(0) << "Saving of config file ");
    1436   if (configuration->Save(filename, periode, mol))
    1437     DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
    1438   else
    1439     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1440 
    1441   // and save to xyz file
    1442   if (ConfigFileName != NULL) {
    1443     strcpy(filename, ConfigFileName);
    1444     strcat(filename, ".xyz");
    1445     output.open(filename, ios::trunc);
    1446   }
    1447   if (output == NULL) {
    1448     strcpy(filename,"main_pcp_linux");
    1449     strcat(filename, ".xyz");
    1450     output.open(filename, ios::trunc);
    1451   }
    1452   DoLog(0) && (Log() << Verbose(0) << "Saving of XYZ file ");
    1453   if (mol->MDSteps <= 1) {
    1454     if (mol->OutputXYZ(&output))
    1455       DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
    1456     else
    1457       DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1458   } else {
    1459     if (mol->OutputTrajectoriesXYZ(&output))
    1460       DoLog(0) && (Log() << Verbose(0) << "successful." << endl);
    1461     else
    1462       DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1463   }
    1464   output.close();
    1465   output.clear();
    1466 
    1467   // and save as MPQC configuration
    1468   if (ConfigFileName != NULL)
    1469     strcpy(filename, ConfigFileName);
    1470   if (output == NULL)
    1471     strcpy(filename,"main_pcp_linux");
    1472   DoLog(0) && (Log() << Verbose(0) << "Saving as mpqc input ");
    1473   if (configuration->SaveMPQC(filename, mol))
    1474     DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    1475   else
    1476     DoLog(0) && (Log() << Verbose(0) << "failed." << endl);
    1477 
    1478   if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
    1479     DoeLog(2) && (eLog()<< Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl);
    1480   }
    1481 
    1482   World::getInstance().destroyMolecule(mol);
    1483 };
    1484 
    1485 #endif
    1486 
    1487 /** Parses the command line options.
    1488  * Note that this function is from now on transitional. All commands that are not passed
    1489  * here are handled by CommandLineParser and the actions of CommandLineUIFactory.
    1490  * \param argc argument count
    1491  * \param **argv arguments array
    1492  * \param *molecules list of molecules structure
    1493  * \param *periode elements structure
    1494  * \param configuration config file structure
    1495  * \param *ConfigFileName pointer to config file name in **argv
    1496  * \param *PathToDatabases pointer to db's path in **argv
    1497  * \param &ArgcList list of arguments that we do not parse here
    1498  * \return exit code (0 - successful, all else - something's wrong)
    1499  */
    1500 static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,
    1501                                    config& configuration, char **ConfigFileName, set<int> &ArgcList)
    1502 {
    1503   Vector x,y,z,n;  // coordinates for absolute point in cell volume
    1504   ifstream test;
    1505   ofstream output;
    1506   string line;
    1507   bool SaveFlag = false;
    1508   int ExitFlag = 0;
    1509   int j;
    1510   double volume = 0.;
    1511   enum ConfigStatus configPresent = absent;
    1512   int argptr;
    1513   molecule *mol = NULL;
    1514   string BondGraphFileName("\n");
    1515 
    1516   if (argc > 1) { // config file specified as option
    1517     // 1. : Parse options that just set variables or print help
    1518     argptr = 1;
    1519     do {
    1520       if (argv[argptr][0] == '-') {
    1521         DoLog(0) && (Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n");
    1522         argptr++;
    1523         switch(argv[argptr-1][1]) {
    1524           case 'h':
    1525           case 'H':
    1526           case '?':
    1527             ArgcList.insert(argptr-1);
    1528             return(1);
    1529             break;
    1530           case 'v':
    1531             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1532               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying verbosity: -v <level>" << endl);
    1533               performCriticalExit();
    1534             } else {
    1535               setVerbosity(atoi(argv[argptr]));
    1536               ArgcList.insert(argptr-1);
    1537               ArgcList.insert(argptr);
    1538               argptr++;
    1539             }
    1540             break;
    1541           case 'V':
    1542             ArgcList.insert(argptr-1);
    1543             return(1);
    1544             break;
    1545           case 'B':
    1546             if (ExitFlag == 0) ExitFlag = 1;
    1547             if ((argptr+5 >= argc)) {
    1548               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for setting Box: -B <xx> <<xy> <<xz> <yy> <yz> <zz>" << endl);
    1549               performCriticalExit();
    1550             } else {
    1551               ArgcList.insert(argptr-1);
    1552               ArgcList.insert(argptr);
    1553               ArgcList.insert(argptr+1);
    1554               ArgcList.insert(argptr+2);
    1555               ArgcList.insert(argptr+3);
    1556               ArgcList.insert(argptr+4);
    1557               ArgcList.insert(argptr+5);
    1558               argptr+=6;
    1559             }
    1560             break;
    1561           case 'e':
    1562             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1563               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl);
    1564               performCriticalExit();
    1565             } else {
    1566               ArgcList.insert(argptr-1);
    1567               ArgcList.insert(argptr);
    1568               argptr+=1;
    1569             }
    1570             break;
    1571           case 'g':
    1572             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1573               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl);
    1574               performCriticalExit();
    1575             } else {
    1576               ArgcList.insert(argptr-1);
    1577               ArgcList.insert(argptr);
    1578               argptr+=1;
    1579             }
    1580             break;
    1581           case 'M':
    1582             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1583               ExitFlag = 255;
    1584               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -M <basis name>" << endl);
    1585               performCriticalExit();
    1586             } else {
    1587               ArgcList.insert(argptr-1);
    1588               ArgcList.insert(argptr);
    1589               argptr+=1;
    1590             }
    1591             break;
    1592           case 'n':
    1593             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1594               ExitFlag = 255;
    1595               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting fast-parsing: -n <0/1>" << endl);
    1596               performCriticalExit();
    1597             } else {
    1598               ArgcList.insert(argptr-1);
    1599               ArgcList.insert(argptr);
    1600               argptr+=1;
    1601             }
    1602             break;
    1603           case 'X':
    1604             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1605               ExitFlag = 255;
    1606               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting default molecule name: -X <name>" << endl);
    1607               performCriticalExit();
    1608             } else {
    1609               ArgcList.insert(argptr-1);
    1610               ArgcList.insert(argptr);
    1611               argptr+=1;
    1612             }
    1613             break;
    1614           default:   // no match? Step on
    1615             argptr++;
    1616             break;
    1617         }
    1618       } else
    1619         argptr++;
    1620     } while (argptr < argc);
    1621 
    1622     // 3b. Find config file name and parse if possible, also BondGraphFileName
    1623     if (argv[1][0] != '-') {
    1624       // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
    1625       DoLog(0) && (Log() << Verbose(0) << "Config file given." << endl);
    1626       test.open(argv[1], ios::in);
    1627       if (test == NULL) {
    1628         //return (1);
    1629         output.open(argv[1], ios::out);
    1630         if (output == NULL) {
    1631           DoLog(1) && (Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl);
    1632           configPresent = absent;
    1633         } else {
    1634           DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl);
    1635           strcpy(*ConfigFileName, argv[1]);
    1636           configPresent = empty;
    1637           output.close();
    1638         }
    1639       } else {
    1640         test.close();
    1641         strcpy(*ConfigFileName, argv[1]);
    1642         DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... ");
    1643         switch (configuration.TestSyntax(*ConfigFileName, periode)) {
    1644           case 1:
    1645             DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl);
    1646             configuration.Load(*ConfigFileName, BondGraphFileName, periode, molecules);
    1647             configPresent = present;
    1648             break;
    1649           case 0:
    1650             DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl);
    1651             configuration.LoadOld(*ConfigFileName, BondGraphFileName, periode, molecules);
    1652             configPresent = present;
    1653             break;
    1654           default:
    1655             DoLog(0) && (Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl);
    1656             configPresent = empty;
    1657        }
    1658       }
    1659     } else
    1660       configPresent = absent;
    1661      // set mol to first active molecule
    1662      if (molecules->ListOfMolecules.size() != 0) {
    1663        for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    1664          if ((*ListRunner)->ActiveFlag) {
    1665            mol = *ListRunner;
    1666            break;
    1667          }
    1668      }
    1669      if (mol == NULL) {
    1670        mol = World::getInstance().createMolecule();
    1671        mol->ActiveFlag = true;
    1672        if (*ConfigFileName != NULL)
    1673          mol->SetNameFromFilename(*ConfigFileName);
    1674        molecules->insert(mol);
    1675      }
    1676 
    1677     // 4. parse again through options, now for those depending on elements db and config presence
    1678     argptr = 1;
    1679     do {
    1680       DoLog(0) && (Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl);
    1681       if (argv[argptr][0] == '-') {
    1682         argptr++;
    1683         if ((configPresent == present) || (configPresent == empty)) {
    1684           switch(argv[argptr-1][1]) {
    1685             case 'p':
    1686               if (ExitFlag == 0) ExitFlag = 1;
    1687               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1688                 ExitFlag = 255;
    1689                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl);
    1690                 performCriticalExit();
    1691               } else {
    1692                 SaveFlag = true;
    1693                 DoLog(1) && (Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl);
    1694                 if (!mol->AddXYZFile(argv[argptr]))
    1695                   DoLog(2) && (Log() << Verbose(2) << "File not found." << endl);
    1696                 else {
    1697                   DoLog(2) && (Log() << Verbose(2) << "File found and parsed." << endl);
    1698                   configPresent = present;
    1699                 }
    1700               }
    1701               break;
    1702             case 'a':
    1703               if (ExitFlag == 0) ExitFlag = 1;
    1704               if ((argptr+4 >= argc) || (argv[argptr][0] == '-')) {
    1705                 ExitFlag = 255;
    1706                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough arguments for adding atom: -a <Z> --position <x> <y> <z>" << endl);
    1707                 performCriticalExit();
    1708               } else {
    1709                 ArgcList.insert(argptr-1);
    1710                 ArgcList.insert(argptr);
    1711                 ArgcList.insert(argptr+1);
    1712                 ArgcList.insert(argptr+2);
    1713                 ArgcList.insert(argptr+3);
    1714                 ArgcList.insert(argptr+4);
    1715                 argptr+=5;
    1716               }
    1717               break;
    1718             default:   // no match? Don't step on (this is done in next switch's default)
    1719               break;
    1720           }
    1721         }
    1722         if (configPresent == present) {
    1723           switch(argv[argptr-1][1]) {
    1724             case 'D':
    1725               if (ExitFlag == 0) ExitFlag = 1;
    1726               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1727                 ExitFlag = 255;
    1728                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for depth-first-search analysis: -D <max. bond distance>" << endl);
    1729                 performCriticalExit();
    1730               } else {
    1731                 ArgcList.insert(argptr-1);
    1732                 ArgcList.insert(argptr);
    1733                 argptr+=1;
    1734               }
    1735               break;
    1736             case 'I':
    1737               DoLog(1) && (Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl);
    1738               ArgcList.insert(argptr-1);
    1739               argptr+=0;
    1740               break;
    1741             case 'C':
    1742               {
    1743                 if (ExitFlag == 0) ExitFlag = 1;
    1744                 if ((argptr+11 >= argc)) {
    1745                   ExitFlag = 255;
    1746                   DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C[p] <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl);
    1747                   performCriticalExit();
    1748                 } else {
    1749                   switch(argv[argptr][0]) {
    1750                     case 'E':
    1751                       ArgcList.insert(argptr-1);
    1752                       ArgcList.insert(argptr);
    1753                       ArgcList.insert(argptr+1);
    1754                       ArgcList.insert(argptr+2);
    1755                       ArgcList.insert(argptr+3);
    1756                       ArgcList.insert(argptr+4);
    1757                       ArgcList.insert(argptr+5);
    1758                       ArgcList.insert(argptr+6);
    1759                       ArgcList.insert(argptr+7);
    1760                       ArgcList.insert(argptr+8);
    1761                       ArgcList.insert(argptr+9);
    1762                       ArgcList.insert(argptr+10);
    1763                       ArgcList.insert(argptr+11);
    1764                       argptr+=12;
    1765                       break;
    1766 
    1767                     case 'P':
    1768                       ArgcList.insert(argptr-1);
    1769                       ArgcList.insert(argptr);
    1770                       ArgcList.insert(argptr+1);
    1771                       ArgcList.insert(argptr+2);
    1772                       ArgcList.insert(argptr+3);
    1773                       ArgcList.insert(argptr+4);
    1774                       ArgcList.insert(argptr+5);
    1775                       ArgcList.insert(argptr+6);
    1776                       ArgcList.insert(argptr+7);
    1777                       ArgcList.insert(argptr+8);
    1778                       ArgcList.insert(argptr+9);
    1779                       ArgcList.insert(argptr+10);
    1780                       ArgcList.insert(argptr+11);
    1781                       ArgcList.insert(argptr+12);
    1782                       ArgcList.insert(argptr+13);
    1783                       ArgcList.insert(argptr+14);
    1784                       argptr+=15;
    1785                       break;
    1786 
    1787                     case 'S':
    1788                       ArgcList.insert(argptr-1);
    1789                       ArgcList.insert(argptr);
    1790                       ArgcList.insert(argptr+1);
    1791                       ArgcList.insert(argptr+2);
    1792                       ArgcList.insert(argptr+3);
    1793                       ArgcList.insert(argptr+4);
    1794                       ArgcList.insert(argptr+5);
    1795                       ArgcList.insert(argptr+6);
    1796                       ArgcList.insert(argptr+7);
    1797                       ArgcList.insert(argptr+8);
    1798                       ArgcList.insert(argptr+9);
    1799                       ArgcList.insert(argptr+10);
    1800                       ArgcList.insert(argptr+11);
    1801                       ArgcList.insert(argptr+12);
    1802                       ArgcList.insert(argptr+13);
    1803                       ArgcList.insert(argptr+14);
    1804                       argptr+=15;
    1805                       break;
    1806 
    1807                     default:
    1808                       ExitFlag = 255;
    1809                       DoeLog(0) && (eLog()<< Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl);
    1810                       performCriticalExit();
    1811                       break;
    1812                   }
    1813                 }
    1814                 break;
    1815               }
    1816             case 'E':
    1817               if (ExitFlag == 0) ExitFlag = 1;
    1818               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr]))) {
    1819                 ExitFlag = 255;
    1820                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> --element <Z>" << endl);
    1821                 performCriticalExit();
    1822               } else {
    1823                 ArgcList.insert(argptr-1);
    1824                 ArgcList.insert(argptr);
    1825                 ArgcList.insert(argptr+1);
    1826                 ArgcList.insert(argptr+2);
    1827                 argptr+=3;
    1828               }
    1829               break;
    1830             case 'F':
    1831               if (ExitFlag == 0) ExitFlag = 1;
    1832               if ((argptr+12 >= argc) || (argv[argptr][0] == '-')) {
    1833                 ExitFlag = 255;
    1834                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling with molecule: -F <filler xyz file> --MaxDistance <distance or -1> --distances <x> <y> <z>  --lengths <surface> <randatm> <randmol> --DoRotate <0/1>" << endl);
    1835                 performCriticalExit();
    1836               } else {
    1837                 ArgcList.insert(argptr-1);
    1838                 ArgcList.insert(argptr);
    1839                 ArgcList.insert(argptr+1);
    1840                 ArgcList.insert(argptr+2);
    1841                 ArgcList.insert(argptr+3);
    1842                 ArgcList.insert(argptr+4);
    1843                 ArgcList.insert(argptr+5);
    1844                 ArgcList.insert(argptr+6);
    1845                 ArgcList.insert(argptr+7);
    1846                 ArgcList.insert(argptr+8);
    1847                 ArgcList.insert(argptr+9);
    1848                 ArgcList.insert(argptr+10);
    1849                 ArgcList.insert(argptr+11);
    1850                 ArgcList.insert(argptr+12);
    1851                 argptr+=13;
    1852               }
    1853               break;
    1854             case 'A':
    1855               if (ExitFlag == 0) ExitFlag = 1;
    1856               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1857                 ExitFlag =255;
    1858                 DoeLog(0) && (eLog()<< Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile> --molecule-by-id <molecule_id>" << endl);
    1859                 performCriticalExit();
    1860               } else {
    1861                 ArgcList.insert(argptr-1);
    1862                 ArgcList.insert(argptr);
    1863                 ArgcList.insert(argptr+1);
    1864                 ArgcList.insert(argptr+2);
    1865                 argptr+=3;
    1866               }
    1867               break;
    1868 
    1869             case 'J':
    1870               if (ExitFlag == 0) ExitFlag = 1;
    1871               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1872                 ExitFlag =255;
    1873                 DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of adjacency file: -J <path> --molecule-by-id <molecule_id>" << endl);
    1874                 performCriticalExit();
    1875               } else {
    1876                 ArgcList.insert(argptr-1);
    1877                 ArgcList.insert(argptr);
    1878                 ArgcList.insert(argptr+1);
    1879                 ArgcList.insert(argptr+2);
    1880                 argptr+=3;
    1881               }
    1882               break;
    1883 
    1884             case 'j':
    1885               if (ExitFlag == 0) ExitFlag = 1;
    1886               if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1887                 ExitFlag =255;
    1888                 DoeLog(0) && (eLog()<< Verbose(0) << "Missing path of bonds file: -j <path> --molecule-by-id <molecule_id>" << endl);
    1889                 performCriticalExit();
    1890               } else {
    1891                 ArgcList.insert(argptr-1);
    1892                 ArgcList.insert(argptr);
    1893                 ArgcList.insert(argptr+1);
    1894                 ArgcList.insert(argptr+2);
    1895                 argptr+=3;
    1896               }
    1897               break;
    1898 
    1899             case 'N':
    1900               if (ExitFlag == 0) ExitFlag = 1;
    1901               if ((argptr+4 >= argc) || (argv[argptr][0] == '-')){
    1902                 ExitFlag = 255;
    1903                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -N <molecule_id> --sphere-radius <radius> --nonconvex-file <output prefix>" << endl);
    1904                 performCriticalExit();
    1905               } else {
    1906                 ArgcList.insert(argptr-1);
    1907                 ArgcList.insert(argptr);
    1908                 ArgcList.insert(argptr+1);
    1909                 ArgcList.insert(argptr+2);
    1910                 ArgcList.insert(argptr+3);
    1911                 ArgcList.insert(argptr+4);
    1912                 argptr+=5;
    1913               }
    1914               break;
    1915             case 'S':
    1916               if (ExitFlag == 0) ExitFlag = 1;
    1917               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1918                 ExitFlag = 255;
    1919                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file> --molecule-by-id 0" << endl);
    1920                 performCriticalExit();
    1921               } else {
    1922                 ArgcList.insert(argptr-1);
    1923                 ArgcList.insert(argptr);
    1924                 ArgcList.insert(argptr+1);
    1925                 ArgcList.insert(argptr+2);
    1926                 argptr+=3;
    1927               }
    1928               break;
    1929             case 'L':
    1930               if (ExitFlag == 0) ExitFlag = 1;
    1931               if ((argptr+8 >= argc) || (argv[argptr][0] == '-')) {
    1932                 ExitFlag = 255;
    1933                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for linear interpolation: -L <prefix> --start-step <step0> --end-step <step1> --molecule-by-id 0 --id-mapping <0/1>" << endl);
    1934                 performCriticalExit();
    1935               } else {
    1936                 ArgcList.insert(argptr-1);
    1937                 ArgcList.insert(argptr);
    1938                 ArgcList.insert(argptr+1);
    1939                 ArgcList.insert(argptr+2);
    1940                 ArgcList.insert(argptr+3);
    1941                 ArgcList.insert(argptr+4);
    1942                 ArgcList.insert(argptr+5);
    1943                 ArgcList.insert(argptr+6);
    1944                 ArgcList.insert(argptr+7);
    1945                 ArgcList.insert(argptr+8);
    1946                 argptr+=9;
    1947               }
    1948               break;
    1949             case 'P':
    1950               if (ExitFlag == 0) ExitFlag = 1;
    1951               if ((argptr+2 >= argc) || (argv[argptr][0] == '-')) {
    1952                 ExitFlag = 255;
    1953                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file> --molecule-by-id <molecule_id>" << endl);
    1954                 performCriticalExit();
    1955               } else {
    1956                 ArgcList.insert(argptr-1);
    1957                 ArgcList.insert(argptr);
    1958                 ArgcList.insert(argptr+1);
    1959                 ArgcList.insert(argptr+2);
    1960                 argptr+=3;
    1961               }
    1962               break;
    1963             case 'R':
    1964               if (ExitFlag == 0) ExitFlag = 1;
    1965               if ((argptr+4 >= argc) || (argv[argptr][0] == '-'))  {
    1966                 ExitFlag = 255;
    1967                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <distance> --position <x> <y> <z>" << endl);
    1968                 performCriticalExit();
    1969               } else {
    1970                 ArgcList.insert(argptr-1);
    1971                 ArgcList.insert(argptr);
    1972                 ArgcList.insert(argptr+1);
    1973                 ArgcList.insert(argptr+2);
    1974                 ArgcList.insert(argptr+3);
    1975                 ArgcList.insert(argptr+4);
    1976                 argptr+=5;
    1977               }
    1978               break;
    1979             case 't':
    1980               if (ExitFlag == 0) ExitFlag = 1;
    1981               if ((argptr+4 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    1982                 ExitFlag = 255;
    1983                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z> --molecule-by-id <molecule_id> --periodic <0/1>" << endl);
    1984                 performCriticalExit();
    1985               } else {
    1986                 ArgcList.insert(argptr-1);
    1987                 ArgcList.insert(argptr);
    1988                 ArgcList.insert(argptr+1);
    1989                 ArgcList.insert(argptr+2);
    1990                 ArgcList.insert(argptr+3);
    1991                 ArgcList.insert(argptr+4);
    1992                 ArgcList.insert(argptr+5);
    1993                 ArgcList.insert(argptr+6);
    1994                 argptr+=7;
    1995               }
    1996               break;
    1997             case 's':
    1998               if (ExitFlag == 0) ExitFlag = 1;
    1999               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    2000                 ExitFlag = 255;
    2001                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl);
    2002                 performCriticalExit();
    2003               } else {
    2004                 ArgcList.insert(argptr-1);
    2005                 ArgcList.insert(argptr);
    2006                 ArgcList.insert(argptr+1);
    2007                 ArgcList.insert(argptr+2);
    2008                 argptr+=3;
    2009               }
    2010               break;
    2011             case 'b':
    2012               if (ExitFlag == 0) ExitFlag = 1;
    2013               if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
    2014                 ExitFlag = 255;
    2015                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
    2016                 performCriticalExit();
    2017               } else {
    2018                 ArgcList.insert(argptr-1);
    2019                 ArgcList.insert(argptr);
    2020                 ArgcList.insert(argptr+1);
    2021                 ArgcList.insert(argptr+2);
    2022                 ArgcList.insert(argptr+3);
    2023                 ArgcList.insert(argptr+4);
    2024                 ArgcList.insert(argptr+5);
    2025                 argptr+=6;
    2026               }
    2027               break;
    2028             case 'B':
    2029               if (ExitFlag == 0) ExitFlag = 1;
    2030               if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
    2031                 ExitFlag = 255;
    2032                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
    2033                 performCriticalExit();
    2034               } else {
    2035                 ArgcList.insert(argptr-1);
    2036                 ArgcList.insert(argptr);
    2037                 ArgcList.insert(argptr+1);
    2038                 ArgcList.insert(argptr+2);
    2039                 ArgcList.insert(argptr+3);
    2040                 ArgcList.insert(argptr+4);
    2041                 ArgcList.insert(argptr+5);
    2042                 argptr+=6;
    2043               }
    2044               break;
    2045             case 'c':
    2046               if (ExitFlag == 0) ExitFlag = 1;
    2047               if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    2048                 ExitFlag = 255;
    2049                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl);
    2050                 performCriticalExit();
    2051               } else {
    2052                 ArgcList.insert(argptr-1);
    2053                 ArgcList.insert(argptr);
    2054                 ArgcList.insert(argptr+1);
    2055                 ArgcList.insert(argptr+2);
    2056                 argptr+=3;
    2057               }
    2058               break;
    2059             case 'O':
    2060               if (ExitFlag == 0) ExitFlag = 1;
    2061               ArgcList.insert(argptr-1);
    2062               argptr+=0;
    2063               break;
    2064             case 'r':
    2065               if (ExitFlag == 0) ExitFlag = 1;
    2066               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])))  {
    2067                 ExitFlag = 255;
    2068                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl);
    2069                 performCriticalExit();
    2070               } else {
    2071                 ArgcList.insert(argptr-1);
    2072                 ArgcList.insert(argptr);
    2073                 argptr+=1;
    2074               }
    2075               break;
    2076             case 'f':
    2077               if (ExitFlag == 0) ExitFlag = 1;
    2078               if ((argptr+1 >= argc) || (argv[argptr][0] == '-')) {
    2079                 ExitFlag = 255;
    2080                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl);
    2081                 performCriticalExit();
    2082               } else {
    2083                 ArgcList.insert(argptr-1);
    2084                 ArgcList.insert(argptr);
    2085                 ArgcList.insert(argptr+1);
    2086                 ArgcList.insert(argptr+2);
    2087                 ArgcList.insert(argptr+3);
    2088                 ArgcList.insert(argptr+4);
    2089                 argptr+=5;
    2090               }
    2091               break;
    2092             case 'm':
    2093               if (ExitFlag == 0) ExitFlag = 1;
    2094               j = atoi(argv[argptr++]);
    2095               if ((j<0) || (j>1)) {
    2096                 DoeLog(1) && (eLog()<< Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl);
    2097                 j = 0;
    2098               }
    2099               if (j) {
    2100                 SaveFlag = true;
    2101                 DoLog(0) && (Log() << Verbose(0) << "Converting to prinicipal axis system." << endl);
    2102                 mol->PrincipalAxisSystem((bool)j);
    2103               } else
    2104                 ArgcList.insert(argptr-1);
    2105                 argptr+=0;
    2106               break;
    2107             case 'o':
    2108               if (ExitFlag == 0) ExitFlag = 1;
    2109               if ((argptr+4 >= argc) || (argv[argptr][0] == '-')){
    2110                 ExitFlag = 255;
    2111                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <molecule_id> --output-file <output file> --output-file <binned output file>" << endl);
    2112                 performCriticalExit();
    2113               } else {
    2114                 ArgcList.insert(argptr-1);
    2115                 ArgcList.insert(argptr);
    2116                 ArgcList.insert(argptr+1);
    2117                 ArgcList.insert(argptr+2);
    2118                 ArgcList.insert(argptr+3);
    2119                 ArgcList.insert(argptr+4);
    2120                 argptr+=5;
    2121               }
    2122               break;
    2123             case 'U':
    2124               if (ExitFlag == 0) ExitFlag = 1;
    2125               if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
    2126                 ExitFlag = 255;
    2127                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl);
    2128                 performCriticalExit();
    2129               } else {
    2130                 volume = atof(argv[argptr++]);
    2131                 DoLog(0) && (Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl);
    2132               }
    2133             case 'u':
    2134               if (ExitFlag == 0) ExitFlag = 1;
    2135               if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
    2136                 if (volume != -1)
    2137                   ExitFlag = 255;
    2138                   DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl);
    2139                   performCriticalExit();
    2140               } else {
    2141                 ArgcList.insert(argptr-1);
    2142                 ArgcList.insert(argptr);
    2143                 argptr+=1;
    2144               }
    2145               break;
    2146             case 'd':
    2147               if (ExitFlag == 0) ExitFlag = 1;
    2148               if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    2149                 ExitFlag = 255;
    2150                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl);
    2151                 performCriticalExit();
    2152               } else {
    2153                 ArgcList.insert(argptr-1);
    2154                 ArgcList.insert(argptr);
    2155                 ArgcList.insert(argptr+1);
    2156                 ArgcList.insert(argptr+2);
    2157                 argptr+=3;
    2158               }
    2159               break;
    2160             default:   // no match? Step on
    2161               if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
    2162                 argptr++;
    2163               break;
    2164           }
    2165         }
    2166       } else argptr++;
    2167     } while (argptr < argc);
    2168     if (SaveFlag)
    2169       configuration.SaveAll(*ConfigFileName, periode, molecules);
    2170   } else {  // no arguments, hence scan the elements db
    2171     if (periode->LoadPeriodentafel(configuration.databasepath))
    2172       DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl);
    2173     else
    2174       DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl);
    2175     configuration.RetrieveConfigPathAndName("main_pcp_linux");
    2176   }
    2177   return(ExitFlag);
    2178 };
    217980
    218081/********************************************** Main routine **************************************/
    218182
    218283void cleanUp(){
     84  FormatParserStorage::purgeInstance();
     85  ChangeTracker::purgeInstance();
    218386  World::purgeInstance();
    218487  logger::purgeInstance();
    … …  
    2197100int main(int argc, char **argv)
    2198101{
    2199     // while we are non interactive, we want to abort from asserts
    2200     //ASSERT_DO(Assert::Abort);
    2201     Vector x, y, z, n;
    2202     ifstream test;
    2203     ofstream output;
    2204     string line;
    2205     char **Arguments = NULL;
    2206     int ArgcSize = 0;
    2207     int ExitFlag = 0;
    2208     bool ArgumentsCopied = false;
    2209     char *ConfigFileName = new char[MAXSTRINGSIZE];
     102  // while we are non interactive, we want to abort from asserts
     103  //ASSERT_DO(Assert::Abort);
     104  string line;
     105  char **Arguments = NULL;
     106  int ArgcSize = 0;
     107  int ExitFlag = 0;
     108  bool ArgumentsCopied = false;
     109  std::string BondGraphFileName("\n");
     110  FormatParserStorage::getInstance().addMpqc();
     111  FormatParserStorage::getInstance().addPcp();
     112  FormatParserStorage::getInstance().addXyz();
    2210113
    2211     // print version check whether arguments are present at all
    2212     cout << ESPACKVersion << endl;
    2213     if (argc < 2) {
    2214       cout << "Obtain help with " << argv[0] << " -h." << endl;
    2215       cleanUp();
    2216       Memory::getState();
    2217       return(1);
    2218     }
     114  // print version check whether arguments are present at all
     115  cout << ESPACKVersion << endl;
     116  if (argc < 2) {
     117    cout << "Obtain help with " << argv[0] << " -h." << endl;
     118    cleanUp();
     119    Memory::getState();
     120    return(1);
     121  }
    2219122
    2220123
    2221     setVerbosity(0);
    2222     // need to init the history before any action is created
    2223     ActionHistory::init();
     124  setVerbosity(0);
     125  // need to init the history before any action is created
     126  ActionHistory::init();
    2224127
    2225     // In the interactive mode, we can leave the user the choice in case of error
    2226     ASSERT_DO(Assert::Ask);
     128  // In the interactive mode, we can leave the user the choice in case of error
     129  ASSERT_DO(Assert::Ask);
    2227130
    2228     // from this moment on, we need to be sure to deeinitialize in the correct order
    2229     // this is handled by the cleanup function
    2230     atexit(cleanUp);
     131  // from this moment on, we need to be sure to deeinitialize in the correct order
     132  // this is handled by the cleanup function
     133  atexit(cleanUp);
    2231134
    2232     // Parse command line options and if present create respective UI
    2233     {
    2234       set<int> ArgcList;
    2235       ArgcList.insert(0); // push back program!
    2236       ArgcList.insert(1); // push back config file name
    2237       // handle arguments by ParseCommandLineOptions()
    2238       ExitFlag = ParseCommandLineOptions(argc,argv,World::getInstance().getMolecules(),World::getInstance().getPeriode(),*World::getInstance().getConfig(), &ConfigFileName, ArgcList);
    2239       World::getInstance().setExitFlag(ExitFlag);
    2240       // copy all remaining arguments to a new argv
    2241       Arguments = new char *[ArgcList.size()];
    2242       cout << "The following arguments are handled by CommandLineParser: ";
    2243       for (set<int>::iterator ArgcRunner = ArgcList.begin(); ArgcRunner != ArgcList.end(); ++ArgcRunner) {
    2244         Arguments[ArgcSize] = new char[strlen(argv[*ArgcRunner])+2];
    2245         strcpy(Arguments[ArgcSize], argv[*ArgcRunner]);
    2246         cout << " " << argv[*ArgcRunner];
    2247         ArgcSize++;
    2248       }
    2249       cout << endl;
    2250       ArgumentsCopied = true;
    2251       // handle remaining arguments by CommandLineParser
    2252       MapOfActions::getInstance().AddOptionsToParser();
    2253       map <std::string, std::string> ShortFormToActionMap = MapOfActions::getInstance().getShortFormToActionMap();
    2254       CommandLineParser::getInstance().Run(ArgcSize,Arguments, ShortFormToActionMap);
    2255       if (!CommandLineParser::getInstance().isEmpty()) {
    2256         DoLog(0) && (Log() << Verbose(0) << "Setting UI to CommandLine." << endl);
    2257         UIFactory::registerFactory(new CommandLineUIFactory::description());
    2258         UIFactory::makeUserInterface("CommandLine");
     135  // Parse command line options and if present create respective UI
     136  {
     137    // construct bond graph
     138    if (World::getInstance().getConfig()->BG == NULL) {
     139      World::getInstance().getConfig()->BG = new BondGraph(World::getInstance().getConfig()->GetIsAngstroem());
     140      if (World::getInstance().getConfig()->BG->LoadBondLengthTable(BondGraphFileName)) {
     141        DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
    2259142      } else {
    2260         DoLog(0) && (Log() << Verbose(0) << "Setting UI to Text." << endl);
    2261         UIFactory::registerFactory(new TextUIFactory::description());
    2262         UIFactory::makeUserInterface("Text");
     143        DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
    2263144      }
    2264145    }
     146    // handle remaining arguments by CommandLineParser
     147    MapOfActions::getInstance().AddOptionsToParser();
     148    map <std::string, std::string> ShortFormToActionMap = MapOfActions::getInstance().getShortFormToActionMap();
     149    CommandLineParser::getInstance().Run(argc,argv, ShortFormToActionMap);
     150    if (!CommandLineParser::getInstance().isEmpty()) {
     151      DoLog(0) && (Log() << Verbose(0) << "Setting UI to CommandLine." << endl);
     152      UIFactory::registerFactory(new CommandLineUIFactory::description());
     153      UIFactory::makeUserInterface("CommandLine");
     154    } else {
     155      DoLog(0) && (Log() << Verbose(0) << "Setting UI to Text." << endl);
     156      UIFactory::registerFactory(new TextUIFactory::description());
     157      UIFactory::makeUserInterface("Text");
     158    }
     159  }
    2265160
    2266     {
    2267       MainWindow *mainWindow = UIFactory::getInstance().makeMainWindow();
    2268       mainWindow->display();
    2269       delete mainWindow;
    2270     }
     161  {
     162    MainWindow *mainWindow = UIFactory::getInstance().makeMainWindow();
     163    mainWindow->display();
     164    delete mainWindow;
     165  }
    2271166
    2272     Log() << Verbose(0) << "Saving to " << ConfigFileName << "." << endl;
    2273     World::getInstance().getConfig()->SaveAll(ConfigFileName, World::getInstance().getPeriode(), World::getInstance().getMolecules());
     167  FormatParserStorage::getInstance().SaveAll();
     168  ChangeTracker::getInstance().saveStatus();
    2274169
    2275170  // free the new argv
    … …  
    2279174    delete[](Arguments);
    2280175  }
    2281   delete[](ConfigFileName);
     176  //delete[](ConfigFileName);
    2282177
    2283178  ExitFlag = World::getInstance().getExitFlag();

  • src/config.cpp

    r6d574a r0d1ad0  
    1010#include <cstring>
    1111
    12 #include "World.hpp"
    1312#include "atom.hpp"
    1413#include "bond.hpp"
     14#include "bondgraph.hpp"
    1515#include "config.hpp"
     16#include "ConfigFileBuffer.hpp"
    1617#include "element.hpp"
    1718#include "helpers.hpp"
    … …  
    2324#include "molecule.hpp"
    2425#include "periodentafel.hpp"
     26#include "ThermoStatContainer.hpp"
    2527#include "World.hpp"
    26 
    27 /******************************** Functions for class ConfigFileBuffer **********************/
    28 
    29 /** Structure containing compare function for Ion_Type sorting.
    30  */
    31 struct IonTypeCompare {
    32   bool operator()(const char* s1, const char *s2) const {
    33     char number1[8];
    34     char number2[8];
    35     const char *dummy1, *dummy2;
    36     //Log() << Verbose(0) << s1 << "  " << s2 << endl;
    37     dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
    38     dummy2 = strchr(dummy1, '_');
    39     strncpy(number1, dummy1, dummy2-dummy1); // copy the number
    40     number1[dummy2-dummy1]='\0';
    41     dummy1 = strchr(s2, '_')+sizeof(char)*5;  // go just after "Ion_Type"
    42     dummy2 = strchr(dummy1, '_');
    43     strncpy(number2, dummy1, dummy2-dummy1); // copy the number
    44     number2[dummy2-dummy1]='\0';
    45     if (atoi(number1) != atoi(number2))
    46       return (atoi(number1) < atoi(number2));
    47     else {
    48       dummy1 = strchr(s1, '_')+sizeof(char);
    49       dummy1 = strchr(dummy1, '_')+sizeof(char);
    50       dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
    51       strncpy(number1, dummy1, dummy2-dummy1); // copy the number
    52       number1[dummy2-dummy1]='\0';
    53       dummy1 = strchr(s2, '_')+sizeof(char);
    54       dummy1 = strchr(dummy1, '_')+sizeof(char);
    55       dummy2 = strchr(dummy1, ' ') < strchr(dummy1, '\t') ? strchr(dummy1, ' ') : strchr(dummy1, '\t');
    56       strncpy(number2, dummy1, dummy2-dummy1); // copy the number
    57       number2[dummy2-dummy1]='\0';
    58       return (atoi(number1) < atoi(number2));
    59     }
    60   }
    61 };
    62 
    63 /** Constructor for ConfigFileBuffer class.
    64  */
    65 ConfigFileBuffer::ConfigFileBuffer() : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
    66 {
    67 };
    68 
    69 /** Constructor for ConfigFileBuffer class with filename to be parsed.
    70  * \param *filename file name
    71  */
    72 ConfigFileBuffer::ConfigFileBuffer(const char * const filename) : buffer(NULL), LineMapping(NULL), CurrentLine(0), NoLines(0)
    73 {
    74   ifstream *file = NULL;
    75   char line[MAXSTRINGSIZE];
    76 
    77   // prescan number of lines
    78   file= new ifstream(filename);
    79   if (file == NULL) {
    80     DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
    81     return;
    82   }
    83   NoLines = 0; // we're overcounting by one
    84   long file_position = file->tellg(); // mark current position
    85   do {
    86     file->getline(line, 256);
    87     NoLines++;
    88   } while (!file->eof());
    89   file->clear();
    90   file->seekg(file_position, ios::beg);
    91   DoLog(1) && (Log() << Verbose(1) << NoLines-1 << " lines were recognized." << endl);
    92 
    93   // allocate buffer's 1st dimension
    94   if (buffer != NULL) {
    95     DoeLog(1) && (eLog()<< Verbose(1) << "FileBuffer->buffer is not NULL!" << endl);
    96     return;
    97   } else
    98     buffer = new char *[NoLines];
    99 
    100   // scan each line and put into buffer
    101   int lines=0;
    102   int i;
    103   do {
    104     buffer[lines] = new char[MAXSTRINGSIZE];
    105     file->getline(buffer[lines], MAXSTRINGSIZE-1);
    106     i = strlen(buffer[lines]);
    107     buffer[lines][i] = '\n';
    108     buffer[lines][i+1] = '\0';
    109     lines++;
    110   } while((!file->eof()) && (lines < NoLines));
    111   DoLog(1) && (Log() << Verbose(1) << lines-1 << " lines were read into the buffer." << endl);
    112 
    113   // close and exit
    114   file->close();
    115   file->clear();
    116   delete(file);
    117 }
    118 
    119 /** Destructor for ConfigFileBuffer class.
    120  */
    121 ConfigFileBuffer::~ConfigFileBuffer()
    122 {
    123   for(int i=0;i<NoLines;++i)
    124     delete[](buffer[i]);
    125   delete[](buffer);
    126   delete[](LineMapping);
    127 }
    128 
    129 
    130 /** Create trivial mapping.
    131  */
    132 void ConfigFileBuffer::InitMapping()
    133 {
    134   LineMapping = new int[NoLines];
    135   for (int i=0;i<NoLines;i++)
    136     LineMapping[i] = i;
    137 }
    138 
    139 /** Creates a mapping for the \a *FileBuffer's lines containing the Ion_Type keyword such that they are sorted.
    140  * \a *map on return contains a list of NoAtom entries such that going through the list, yields indices to the
    141  * lines in \a *FileBuffer in a sorted manner of the Ion_Type?_? keywords. We assume that ConfigFileBuffer::CurrentLine
    142  * points to first Ion_Type entry.
    143  * \param *FileBuffer pointer to buffer structure
    144  * \param NoAtoms of subsequent lines to look at
    145  */
    146 void ConfigFileBuffer::MapIonTypesInBuffer(const int NoAtoms)
    147 {
    148   map<const char *, int, IonTypeCompare> IonTypeLineMap;
    149   if (LineMapping == NULL) {
    150     DoeLog(0) && (eLog()<< Verbose(0) << "map pointer is NULL: " << LineMapping << endl);
    151     performCriticalExit();
    152     return;
    153   }
    154 
    155   // put all into hashed map
    156   for (int i=0; i<NoAtoms; ++i) {
    157     IonTypeLineMap.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i));
    158   }
    159 
    160   // fill map
    161   int nr=0;
    162   for (map<const char *, int, IonTypeCompare>::iterator runner = IonTypeLineMap.begin(); runner != IonTypeLineMap.end(); ++runner) {
    163     if (CurrentLine+nr < NoLines)
    164       LineMapping[CurrentLine+(nr++)] = runner->second;
    165     else {
    166       DoeLog(0) && (eLog()<< Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl);
    167       performCriticalExit();
    168     }
    169   }
    170 }
    17128
    17229/************************************* Functions for class config ***************************/
    … …  
    17431/** Constructor for config file class.
    17532 */
    176 config::config() : BG(NULL), PsiType(0), MaxPsiDouble(0), PsiMaxNoUp(0), PsiMaxNoDown(0), MaxMinStopStep(1), InitMaxMinStopStep(1), ProcPEGamma(8), ProcPEPsi(1), configpath(NULL),
    177     configname(NULL), FastParsing(false), Deltat(0.01), basis(""), databasepath(NULL), DoConstrainedMD(0), MaxOuterStep(0), Thermostat(4), ThermostatImplemented(NULL),
    178     ThermostatNames(NULL), TempFrequency(2.5), alpha(0.), HooverMass(0.), TargetTemp(0.00095004455), ScaleTempStep(25),  mainname(NULL), defaultpath(NULL), pseudopotpath(NULL),
     33config::config() : BG(NULL), Thermostats(0), PsiType(0), MaxPsiDouble(0), PsiMaxNoUp(0), PsiMaxNoDown(0), MaxMinStopStep(1), InitMaxMinStopStep(1), ProcPEGamma(8), ProcPEPsi(1),
     34    configname(NULL), FastParsing(false), Deltat(0.01), basis(""), databasepath(NULL), DoConstrainedMD(0), MaxOuterStep(0), mainname(NULL), defaultpath(NULL), pseudopotpath(NULL),
    17935    DoOutVis(0), DoOutMes(1), DoOutNICS(0), DoOutOrbitals(0), DoOutCurrent(0), DoFullCurrent(0), DoPerturbation(0), DoWannier(0), CommonWannier(0), SawtoothStart(0.01),
    18036    VectorPlane(0), VectorCut(0.), UseAddGramSch(1), Seed(1), OutVisStep(10), OutSrcStep(5), MaxPsiStep(0), EpsWannier(1e-7), MaxMinStep(100), RelEpsTotalEnergy(1e-7),
    … …  
    18642  pseudopotpath = new char[MAXSTRINGSIZE];
    18743  databasepath = new char[MAXSTRINGSIZE];
    188   configpath = new char[MAXSTRINGSIZE];
    18944  configname = new char[MAXSTRINGSIZE];
     45  Thermostats = new ThermoStatContainer();
    19046  strcpy(mainname,"pcp");
    19147  strcpy(defaultpath,"not specified");
    19248  strcpy(pseudopotpath,"not specified");
    193   configpath[0]='\0';
    19449  configname[0]='\0';
    19550  basis = "3-21G";
    196 
    197   InitThermostats();
    19851};
    19952
    … …  
    20659  delete[](pseudopotpath);
    20760  delete[](databasepath);
    208   delete[](configpath);
    20961  delete[](configname);
    210   delete[](ThermostatImplemented);
    211   for (int j=0;j<MaxThermostats;j++)
    212     delete[](ThermostatNames[j]);
    213   delete[](ThermostatNames);
     62  if (Thermostats != NULL)
     63    delete(Thermostats);
    21464
    21565  if (BG != NULL)
    21666    delete(BG);
    21767};
    218 
    219 /** Initialises variables in class config for Thermostats.
    220  */
    221 void config::InitThermostats()
    222 {
    223   ThermostatImplemented = new int[MaxThermostats];
    224   ThermostatNames = new char *[MaxThermostats];
    225   for (int j=0;j<MaxThermostats;j++)
    226     ThermostatNames[j] = new char[12];
    227 
    228   strcpy(ThermostatNames[0],"None");
    229   ThermostatImplemented[0] = 1;
    230   strcpy(ThermostatNames[1],"Woodcock");
    231   ThermostatImplemented[1] = 1;
    232   strcpy(ThermostatNames[2],"Gaussian");
    233   ThermostatImplemented[2] = 1;
    234   strcpy(ThermostatNames[3],"Langevin");
    235   ThermostatImplemented[3] = 1;
    236   strcpy(ThermostatNames[4],"Berendsen");
    237   ThermostatImplemented[4] = 1;
    238   strcpy(ThermostatNames[5],"NoseHoover");
    239   ThermostatImplemented[5] = 1;
    240 };
    241 
    242 /** Readin of Thermostat related values from parameter file.
    243  * \param *fb file buffer containing the config file
    244  */
    245 void config::ParseThermostats(class ConfigFileBuffer * const fb)
    246 {
    247   char * const thermo = new char[12];
    248   const int verbose = 0;
    249 
    250   // read desired Thermostat from file along with needed additional parameters
    251   if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
    252     if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
    253       if (ThermostatImplemented[0] == 1) {
    254         Thermostat = None;
    255       } else {
    256         DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    257         Thermostat = None;
    258       }
    259     } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
    260       if (ThermostatImplemented[1] == 1) {
    261         Thermostat = Woodcock;
    262         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
    263       } else {
    264         DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    265         Thermostat = None;
    266       }
    267     } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
    268       if (ThermostatImplemented[2] == 1) {
    269         Thermostat = Gaussian;
    270         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
    271       } else {
    272         DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    273         Thermostat = None;
    274       }
    275     } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
    276       if (ThermostatImplemented[3] == 1) {
    277         Thermostat = Langevin;
    278         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
    279         if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
    280           DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl);
    281         } else {
    282           alpha = 1.;
    283         }
    284       } else {
    285         DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    286         Thermostat = None;
    287       }
    288     } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
    289       if (ThermostatImplemented[4] == 1) {
    290         Thermostat = Berendsen;
    291         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
    292       } else {
    293         DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    294         Thermostat = None;
    295       }
    296     } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
    297       if (ThermostatImplemented[5] == 1) {
    298         Thermostat = NoseHoover;
    299         ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
    300         alpha = 0.;
    301       } else {
    302         DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    303         Thermostat = None;
    304       }
    305     } else {
    306       DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
    307       Thermostat = None;
    308     }
    309   } else {
    310     if ((MaxOuterStep > 0) && (TargetTemp != 0))
    311       DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    312     Thermostat = None;
    313   }
    314   delete[](thermo);
    315 };
    316 
    31768
    31869/** Displays menu for editing each entry of the config file.
    … …  
    629380};
    630381
    631 /** Retrieves the path in the given config file name.
    632  * \param filename config file string
    633  */
    634 void config::RetrieveConfigPathAndName(const string filename)
    635 {
    636   char *ptr = NULL;
    637   char *buffer = new char[MAXSTRINGSIZE];
    638   strncpy(buffer, filename.c_str(), MAXSTRINGSIZE);
    639   int last = -1;
    640   for(last=MAXSTRINGSIZE;last--;) {
    641     if (buffer[last] == '/')
    642       break;
    643   }
    644   if (last == -1) { // no path in front, set to local directory.
    645     strcpy(configpath, "./");
    646     ptr = buffer;
    647   } else {
    648     strncpy(configpath, buffer, last+1);
    649     ptr = &buffer[last+1];
    650     if (last < 254)
    651       configpath[last+1]='\0';
    652   }
    653   strcpy(configname, ptr);
    654   DoLog(0) && (Log() << Verbose(0) << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl);
    655   delete[](buffer);
    656 };
    657 
    658 /** Initializes ConfigFileBuffer from a file.
    659  * \param *file input file stream being the opened config file
    660  * \param *FileBuffer pointer to FileBuffer on return, should point to NULL
    661  */
    662 void PrepareFileBuffer(const char * const filename, struct ConfigFileBuffer *&FileBuffer)
    663 {
    664   if (FileBuffer != NULL) {
    665     DoeLog(2) && (eLog()<< Verbose(2) << "deleting present FileBuffer in PrepareFileBuffer()." << endl);
    666     delete(FileBuffer);
    667   }
    668   FileBuffer = new ConfigFileBuffer(filename);
    669 
    670   FileBuffer->InitMapping();
    671 };
    672 
    673382/** Loads a molecule from a ConfigFileBuffer.
    674383 * \param *mol molecule to load
    … …  
    864573  file->close();
    865574  delete(file);
    866   RetrieveConfigPathAndName(filename);
    867575
    868576  // ParseParameterFile
    869   struct ConfigFileBuffer *FileBuffer = NULL;
    870   PrepareFileBuffer(filename,FileBuffer);
     577  class ConfigFileBuffer *FileBuffer = new ConfigFileBuffer(filename);
    871578
    872579  /* Oeffne Hauptparameterdatei */
    … …  
    877584  int verbose = 0;
    878585 
     586  //TODO: This is actually sensible?: if (MaxOuterStep > 0)
    879587  ParseThermostats(FileBuffer);
    880588 
    … …  
    941649  ParseForParameter(verbose,FileBuffer,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    942650  ParseForParameter(verbose,FileBuffer,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    943   ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
     651  ParseForParameter(verbose,FileBuffer,"TargetTemp", 0, 1, 1, double_type, &(Thermostats->TargetTemp), 1, optional);
    944652  //ParseForParameter(verbose,FileBuffer,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
    945653  if (!ParseForParameter(verbose,FileBuffer,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
    … …  
    1101809    return;
    1102810  }
    1103   RetrieveConfigPathAndName(filename);
    1104811  // ParseParameters
    1105812
    … …  
    1150857  ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
    1151858  ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
    1152   ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
    1153   ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
     859  ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(Thermostats->TargetTemp), 1, optional);
     860  ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(Thermostats->ScaleTempStep), 1, optional);
    1154861  config::EpsWannier = 1e-8;
    1155862
    … …  
    13391046    *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
    13401047    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
    1341     *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
    1342     switch(Thermostat) {
     1048    *output << "Thermostat\t" << Thermostats->ThermostatNames[Thermostats->Thermostat] << "\t";
     1049    switch(Thermostats->Thermostat) {
    13431050      default:
    13441051      case None:
    13451052        break;
    13461053      case Woodcock:
    1347         *output << ScaleTempStep;
     1054        *output << Thermostats->ScaleTempStep;
    13481055        break;
    13491056      case Gaussian:
    1350         *output << ScaleTempStep;
     1057        *output << Thermostats->ScaleTempStep;
    13511058        break;
    13521059      case Langevin:
    1353         *output << TempFrequency << "\t" << alpha;
     1060        *output << Thermostats->TempFrequency << "\t" << Thermostats->alpha;
    13541061        break;
    13551062      case Berendsen:
    1356         *output << TempFrequency;
     1063        *output << Thermostats->TempFrequency;
    13571064        break;
    13581065      case NoseHoover:
    1359         *output << HooverMass;
     1066        *output << Thermostats->HooverMass;
    13601067        break;
    13611068    };
    … …  
    13721079    *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
    13731080    *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
    1374     *output << "TargetTemp\t" << config::TargetTemp << "\t# Target temperature" << endl;
     1081    *output << "TargetTemp\t" << Thermostats->TargetTemp << "\t# Target temperature" << endl;
    13751082    *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
    13761083    *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
    … …  
    14831190    // output of atoms
    14841191    AtomNo = 0;
    1485     mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
     1192    mol->ActOnAllAtoms( &atom::OutputMPQCLine, (ostream * const) output, (const Vector *)center, &AtomNo );
    14861193    delete(center);
    14871194    *output << "\t}" << endl;
    … …  
    15251232    // output of atoms
    15261233    AtomNo = 0;
    1527     mol->ActOnAllAtoms( &atom::OutputMPQCLine, output, (const Vector *)center, &AtomNo );
     1234    mol->ActOnAllAtoms( &atom::OutputMPQCLine, (ostream * const) output, (const Vector *)center, &AtomNo );
    15281235    delete(center);
    15291236    *output << "\t}" << endl;
    … …  
    17861493  molecule *mol = NULL;
    17871494
    1788   if (!strcmp(configpath, GetDefaultPath())) {
    1789     eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
    1790   }
    1791 
    1792 
    17931495  // first save as PDB data
    17941496  if (ConfigFileName != NULL)
    … …  
    18271529    mol->doCountAtoms();
    18281530    mol->CountElements();
    1829     mol->CalculateOrbitals(*this);
     1531    //mol->CalculateOrbitals(*this);
    18301532    delete[](src);
    18311533  } else {
    … …  
    18331535      mol = *(molecules->ListOfMolecules.begin());
    18341536      mol->doCountAtoms();
    1835       mol->CalculateOrbitals(*this);
     1537      //mol->CalculateOrbitals(*this);
    18361538    } else {
    18371539      DoeLog(1) && (eLog() << Verbose(1) << "There are no molecules to save!" << endl);
    … …  
    18951597  else
    18961598    Log() << Verbose(0) << "\t... failed." << endl;
    1897 
    1898   if (!strcmp(configpath, GetDefaultPath())) {
    1899     eLog() << Verbose(2) << "config is found under different path then stated in config file::defaultpath!" << endl;
    1900   }
    19011599
    19021600  // don't destroy molecule as it contains all our atoms
    … …  
    23462044  return (found); // true if found, false if not
    23472045}
     2046
     2047/** Reading of Thermostat related values from parameter file.
     2048 * \param *fb file buffer containing the config file
     2049 */
     2050void config::ParseThermostats(class ConfigFileBuffer * const fb)
     2051{
     2052  char * const thermo = new char[12];
     2053  const int verbose = 0;
     2054
     2055  // read desired Thermostat from file along with needed additional parameters
     2056  if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
     2057    if (strcmp(thermo, Thermostats->ThermostatNames[0]) == 0) { // None
     2058      if (Thermostats->ThermostatImplemented[0] == 1) {
     2059        Thermostats->Thermostat = None;
     2060      } else {
     2061        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
     2062        Thermostats->Thermostat = None;
     2063      }
     2064    } else if (strcmp(thermo, Thermostats->ThermostatNames[1]) == 0) { // Woodcock
     2065      if (Thermostats->ThermostatImplemented[1] == 1) {
     2066        Thermostats->Thermostat = Woodcock;
     2067        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read scaling frequency
     2068      } else {
     2069        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
     2070        Thermostats->Thermostat = None;
     2071      }
     2072    } else if (strcmp(thermo, Thermostats->ThermostatNames[2]) == 0) { // Gaussian
     2073      if (Thermostats->ThermostatImplemented[2] == 1) {
     2074        Thermostats->Thermostat = Gaussian;
     2075        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &Thermostats->ScaleTempStep, 1, critical); // read collision rate
     2076      } else {
     2077        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
     2078        Thermostats->Thermostat = None;
     2079      }
     2080    } else if (strcmp(thermo, Thermostats->ThermostatNames[3]) == 0) { // Langevin
     2081      if (Thermostats->ThermostatImplemented[3] == 1) {
     2082        Thermostats->Thermostat = Langevin;
     2083        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read gamma
     2084        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &Thermostats->alpha, 1, optional)) {
     2085          DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << Thermostats->alpha << "." << endl);
     2086        } else {
     2087          Thermostats->alpha = 1.;
     2088        }
     2089      } else {
     2090        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
     2091        Thermostats->Thermostat = None;
     2092      }
     2093    } else if (strcmp(thermo, Thermostats->ThermostatNames[4]) == 0) { // Berendsen
     2094      if (Thermostats->ThermostatImplemented[4] == 1) {
     2095        Thermostats->Thermostat = Berendsen;
     2096        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->TempFrequency, 1, critical); // read \tau_T
     2097      } else {
     2098        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
     2099        Thermostats->Thermostat = None;
     2100      }
     2101    } else if (strcmp(thermo, Thermostats->ThermostatNames[5]) == 0) { // Nose-Hoover
     2102      if (Thermostats->ThermostatImplemented[5] == 1) {
     2103        Thermostats->Thermostat = NoseHoover;
     2104        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &Thermostats->HooverMass, 1, critical); // read Hoovermass
     2105        Thermostats->alpha = 0.;
     2106      } else {
     2107        DoLog(1) && (Log() << Verbose(1) << "Warning: " << Thermostats->ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
     2108        Thermostats->Thermostat = None;
     2109      }
     2110    } else {
     2111      DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
     2112      Thermostats->Thermostat = None;
     2113    }
     2114  } else {
     2115    if ((Thermostats->TargetTemp != 0))
     2116      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
     2117    Thermostats->Thermostat = None;
     2118  }
     2119  delete[](thermo);
     2120};
     2121

  • src/config.hpp

    r6d574a r0d1ad0  
    2020#include <string>
    2121
    22 #include "bondgraph.hpp"
    23 
    2422/****************************************** forward declarations *****************************/
    2523
     24class BondGraph;
     25class ConfigFileBuffer;
    2626class molecule;
    2727class MoleculeListClass;
    2828class periodentafel;
     29class ThermoStatContainer;
    2930
    3031/********************************************** declarations *******************************/
    31 
    32 class ConfigFileBuffer {
    33   public:
    34     char **buffer;
    35     int *LineMapping;
    36     int CurrentLine;
    37     int NoLines;
    38 
    39     ConfigFileBuffer();
    40     ConfigFileBuffer(const char * const filename);
    41     ~ConfigFileBuffer();
    42 
    43     void InitMapping();
    44     void MapIonTypesInBuffer(const int NoAtoms);
    45 };
    4632
    4733/** The config file.
    … …  
    5137  public:
    5238    class BondGraph *BG;
     39    class ThermoStatContainer *Thermostats;
    5340
    5441    int PsiType;
    … …  
    6047    int ProcPEGamma;
    6148    int ProcPEPsi;
    62     char *configpath;
    6349    char *configname;
    6450    bool FastParsing;
    … …  
    7056    int DoConstrainedMD;
    7157    int MaxOuterStep;
    72     int Thermostat;
    73     int *ThermostatImplemented;
    74     char **ThermostatNames;
    75     double TempFrequency;
    76     double alpha;
    77     double HooverMass;
    78     double TargetTemp;
    79     int ScaleTempStep;
    8058
    8159  private:
    … …  
    138116  void Load(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList);
    139117  void LoadOld(const char * const filename, const string &BondGraphFileName, const periodentafel * const periode, MoleculeListClass * const &MolList);
    140   void RetrieveConfigPathAndName(const string filename);
    141118  bool Save(const char * const filename, const periodentafel * const periode, molecule * const mol) const;
    142119  bool SaveMPQC(const char * const filename, const molecule * const mol) const;
    … …  
    152129  char *GetDefaultPath() const;
    153130  void SetDefaultPath(const char * const path);
    154   void InitThermostats();
    155131  void ParseThermostats(class ConfigFileBuffer * const fb);
    156132};

  • src/molecule.cpp

    r6d574a r0d1ad0  
    7979void molecule::setName(const std::string _name){
    8080  OBSERVE;
     81  cout << "Set name of molecule " << getId() << " to " << _name << endl;
    8182  strncpy(name,_name.c_str(),MAXSTRINGSIZE);
    8283}
    … …  
    739740  else
    740741    length = strlen(molname) - strlen(endname);
     742  cout << "Set name of molecule " << getId() << " to " << molname << endl;
    741743  strncpy(name, molname, length);
    742744  name[length]='\0';
    … …  
    880882        ElementNo[i] = current++;
    881883    }
    882     ActOnAllAtoms( &atom::OutputArrayIndexed, output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
     884    ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
    883885    return true;
    884886  }
    … …  
    10031005  for(int i=MAX_ELEMENTS;i--;)
    10041006    ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
    1005 };
    1006 
    1007 
    1008 /** Counts necessary number of valence electrons and returns number and SpinType.
    1009  * \param configuration containing everything
    1010  */
    1011 void molecule::CalculateOrbitals(class config &configuration)
    1012 {
    1013   configuration.MaxPsiDouble = configuration.PsiMaxNoDown = configuration.PsiMaxNoUp = configuration.PsiType = 0;
    1014   for(int i=MAX_ELEMENTS;i--;) {
    1015     if (ElementsInMolecule[i] != 0) {
    1016       //Log() << Verbose(0) << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
    1017       configuration.MaxPsiDouble += ElementsInMolecule[i]*((int)elemente->FindElement(i)->Valence);
    1018     }
    1019   }
    1020   configuration.PsiMaxNoDown = configuration.MaxPsiDouble/2 + (configuration.MaxPsiDouble % 2);
    1021   configuration.PsiMaxNoUp = configuration.MaxPsiDouble/2;
    1022   configuration.MaxPsiDouble /= 2;
    1023   configuration.PsiType = (configuration.PsiMaxNoDown == configuration.PsiMaxNoUp) ? 0 : 1;
    1024   if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2) && ((configuration.PsiMaxNoDown != 1) || (configuration.PsiMaxNoUp != 0))) {
    1025     configuration.ProcPEGamma /= 2;
    1026     configuration.ProcPEPsi *= 2;
    1027   } else {
    1028     configuration.ProcPEGamma *= configuration.ProcPEPsi;
    1029     configuration.ProcPEPsi = 1;
    1030   }
    1031   cout << configuration.PsiMaxNoDown << ">" << configuration.PsiMaxNoUp << endl;
    1032   if (configuration.PsiMaxNoDown > configuration.PsiMaxNoUp) {
    1033     configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.PsiMaxNoDown;
    1034     cout << configuration.PsiMaxNoDown << " " << configuration.InitMaxMinStopStep << endl;
    1035   } else {
    1036     configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.PsiMaxNoUp;
    1037     cout << configuration.PsiMaxNoUp << " " << configuration.InitMaxMinStopStep << endl;
    1038   }
    10391007};
    10401008

  • src/molecule.hpp

    r6d574a r0d1ad0  
    8080  double *PenaltyConstants;   //!<  penalty constant in front of each term
    8181};
    82 
    83 #define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
    84 enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
    85 
    8682
    8783/** The complete molecule.
    … …  
    265261  /// Count and change present atoms' coordination.
    266262  void CountElements();
    267   void CalculateOrbitals(class config &configuration);
    268263  bool CenterInBox();
    269264  bool BoundInBox();
    … …  
    292287  double MinimiseConstrainedPotential(atom **&permutation, int startstep, int endstep, bool IsAngstroem);
    293288  void EvaluateConstrainedForces(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
    294   bool LinearInterpolationBetweenConfiguration(int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity);
     289  bool LinearInterpolationBetweenConfiguration(int startstep, int endstep, std::string &prefix, config &configuration, bool MapByIdentity);
    295290       
    296291  bool CheckBounds(const Vector *x) const;
    … …  
    324319
    325320  /// Fragment molecule by two different approaches:
    326   int FragmentMolecule(int Order, config *configuration);
    327   bool CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
    328   bool StoreBondsToFile(char *path, char *filename);
    329   bool StoreAdjacencyToFile(char *path, char *filename);
    330   bool CheckAdjacencyFileAgainstMolecule(char *path, atom **ListOfAtoms);
    331   bool ParseOrderAtSiteFromFile(char *path);
    332   bool StoreOrderAtSiteFile(char *path);
    333   bool StoreForcesFile(MoleculeListClass *BondFragments, char *path, int *SortIndex);
     321  int FragmentMolecule(int Order, std::string &prefix);
     322  bool CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, std::string path = "");
     323  bool StoreBondsToFile(std::string &filename, std::string path = "");
     324  bool StoreAdjacencyToFile(std::string &filename, std::string path = "");
     325  bool CheckAdjacencyFileAgainstMolecule(std::string &path, atom **ListOfAtoms);
     326  bool ParseOrderAtSiteFromFile(std::string &path);
     327  bool StoreOrderAtSiteFile(std::string &path);
     328  bool StoreForcesFile(MoleculeListClass *BondFragments, std::string &path, int *SortIndex);
    334329  bool CreateMappingLabelsToConfigSequence(int *&SortIndex);
    335330  bool CreateFatherLookupTable(atom **&LookupTable, int count = 0);
    … …  
    380375  ~MoleculeListClass();
    381376
    382   bool AddHydrogenCorrection(char *path);
    383   bool StoreForcesFile(char *path, int *SortIndex);
     377  bool AddHydrogenCorrection(std::string &path);
     378  bool StoreForcesFile(std::string &path, int *SortIndex);
    384379  void insert(molecule *mol);
    385380  void erase(molecule *mol);
    386381  molecule * ReturnIndex(int index);
    387   bool OutputConfigForListOfFragments(config *configuration, int *SortIndex);
     382  bool OutputConfigForListOfFragments(std::string &prefix, int *SortIndex);
    388383  int NumberOfActiveMolecules();
    389384  void Enumerate(ostream *out);

  • src/molecule_dynamics.cpp

    r6d574a r0d1ad0  
    1818#include "parser.hpp"
    1919#include "Plane.hpp"
     20#include "ThermoStatContainer.hpp"
    2021
    2122/************************************* Functions for class molecule *********************************/
    … …  
    472473 * \param startstep stating initial configuration in molecule::Trajectories
    473474 * \param endstep stating final configuration in molecule::Trajectories
     475 * \param &prefix path and prefix
    474476 * \param &config configuration structure
    475477 * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
    476478 * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
    477479 */
    478 bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
     480bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, std::string &prefix, config &configuration, bool MapByIdentity)
    479481{
    480482  molecule *mol = NULL;
    … …  
    524526  for (int i=getAtomCount(); i--; )
    525527    SortIndex[i] = i;
    526   status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
     528
     529  status = MoleculePerStep->OutputConfigForListOfFragments(prefix, SortIndex);
    527530  delete[](SortIndex);
    528531
    … …  
    643646
    644647  // calculate scale configuration
    645   ScaleTempFactor = configuration.TargetTemp/ActualTemp;
     648  ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
    646649
    647650  // differentating between the various thermostats
    … …  
    651654      break;
    652655     case Woodcock:
    653       if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
     656      if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
    654657        DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
    655658        ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
    … …  
    684687      delta_alpha = 0.;
    685688      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    686       delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
    687       configuration.alpha += delta_alpha*configuration.Deltat;
    688       DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
     689      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
     690      configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
     691      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
    689692      // apply updated alpha as additional force
    690693      ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );

  • src/molecule_fragmentation.cpp

    r6d574a r0d1ad0  
    8282 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
    8383 * Finally, the temporary graph is inserted into the given \a FragmentList for return.
    84  * \param *out output stream for debugging
    85  * \param *path path to file
     84 * \param &path path to file
    8685 * \param *FragmentList empty, filled on return
    8786 * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
    8887 */
    89 bool ParseKeySetFile(char *path, Graph *&FragmentList)
     88bool ParseKeySetFile(std::string &path, Graph *&FragmentList)
    9089{
    9190  bool status = true;
    … …  
    9493  GraphTestPair testGraphInsert;
    9594  int NumberOfFragments = 0;
    96   char filename[MAXSTRINGSIZE];
     95  string filename;
    9796
    9897  if (FragmentList == NULL) { // check list pointer
    … …  
    102101  // 1st pass: open file and read
    103102  DoLog(1) && (Log() << Verbose(1) << "Parsing the KeySet file ... " << endl);
    104   sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, KEYSETFILE);
    105   InputFile.open(filename);
    106   if (InputFile != NULL) {
     103  filename = path + KEYSETFILE;
     104  InputFile.open(filename.c_str());
     105  if (InputFile.good()) {
    107106    // each line represents a new fragment
    108107    char buffer[MAXSTRINGSIZE];
    … …  
    181180
    182181/** Stores key sets to file.
    183  * \param *out output stream for debugging
    184182 * \param KeySetList Graph with Keysets
    185  * \param *path path to file
     183 * \param &path path to file
    186184 * \return true - file written successfully, false - writing failed
    187185 */
    188 bool StoreKeySetFile(Graph &KeySetList, char *path)
    189 {
    190   ofstream output;
     186bool StoreKeySetFile(Graph &KeySetList, std::string &path)
     187{
    191188  bool status =  true;
    192   string line;
     189  string line = path + KEYSETFILE;
     190  ofstream output(line.c_str());
    193191
    194192  // open KeySet file
    195   line = path;
    196   line.append("/");
    197   line += FRAGMENTPREFIX;
    198   line += KEYSETFILE;
    199   output.open(line.c_str(), ios::out);
    200193  DoLog(1) && (Log() << Verbose(1) << "Saving key sets of the total graph ... ");
    201   if(output != NULL) {
     194  if(output.good()) {
    202195    for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
    203196      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
    … …  
    302295
    303296/** Scans the adaptive order file and insert (index, value) into map.
    304  * \param *out output stream for debugging
    305  * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
     297 * \param &path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
    306298 * \param &IndexedKeySetList list to find key set for a given index \a No
    307299 * \return adaptive criteria list from file
    308300 */
    309 map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(char *path, map<int,KeySet> &IndexKeySetList)
     301map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(std::string &path, map<int,KeySet> &IndexKeySetList)
    310302{
    311303  map<int, pair<double,int> > *AdaptiveCriteriaList = new map<int, pair<double,int> >;
    … …  
    313305  double Value = 0.;
    314306  char buffer[MAXSTRINGSIZE];
    315   sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
    316   ifstream InputFile(buffer, ios::in);
     307  string filename = path + ENERGYPERFRAGMENT;
     308  ifstream InputFile(filename.c_str());
     309
     310  if (InputFile.fail()) {
     311    DoeLog(1) && (eLog() << Verbose(1) << "Cannot find file " << filename << "." << endl);
     312    return AdaptiveCriteriaList;
     313  }
    317314
    318315  if (CountLinesinFile(InputFile) > 0) {
    … …  
    419416
    420417/** Checks whether the OrderAtSite is still below \a Order at some site.
    421  * \param *out output stream for debugging
    422418 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
    423419 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
    424420 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
    425421 * \param *MinimumRingSize array of max. possible order to avoid loops
    426  * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
     422 * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
    427423 * \return true - needs further fragmentation, false - does not need fragmentation
    428424 */
    429 bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
     425bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, std::string path)
    430426{
    431427  bool status = false;
    … …  
    585581 * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
    586582 * subgraph in the MoleculeListClass.
    587  * \param *out output stream for debugging
    588583 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
    589  * \param *configuration configuration for writing config files for each fragment
     584 * \param &prefix path and prefix of the bond order configs to be written
    590585 * \return 1 - continue, 2 - stop (no fragmentation occured)
    591586 */
    592 int molecule::FragmentMolecule(int Order, config *configuration)
     587int molecule::FragmentMolecule(int Order, std::string &prefix)
    593588{
    594589  MoleculeListClass *BondFragments = NULL;
    … …  
    624619
    625620  // === compare it with adjacency file ===
    626   FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(configuration->configpath, ListOfAtoms);
     621  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(prefix, ListOfAtoms);
    627622  delete[](ListOfAtoms);
    628623
    … …  
    658653
    659654  // ===== 3. if structure still valid, parse key set file and others =====
    660   FragmentationToDo = FragmentationToDo && ParseKeySetFile(configuration->configpath, ParsedFragmentList);
     655  FragmentationToDo = FragmentationToDo && ParseKeySetFile(prefix, ParsedFragmentList);
    661656
    662657  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
    663   FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(configuration->configpath);
     658  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(prefix);
    664659
    665660  // =================================== Begin of FRAGMENTATION ===============================
    … …  
    672667  AtomMask[getAtomCount()] = false;
    673668  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
    674   while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
     669  while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, prefix))) {
    675670    FragmentationToDo = FragmentationToDo || CheckOrder;
    676671    AtomMask[getAtomCount()] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    … …  
    727722    KeySet test = (*runner).first;
    728723    DoLog(0) && (Log() << Verbose(0) << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl);
    729     BondFragments->insert(StoreFragmentFromKeySet(test, configuration));
     724    BondFragments->insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
    730725    k++;
    731726  }
    … …  
    739734
    740735    DoLog(1) && (Log() << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl);
    741     if (BondFragments->OutputConfigForListOfFragments(configuration, SortIndex))
     736    if (BondFragments->OutputConfigForListOfFragments(prefix, SortIndex))
    742737      DoLog(1) && (Log() << Verbose(1) << "All configs written." << endl);
    743738    else
    … …  
    745740
    746741    // store force index reference file
    747     BondFragments->StoreForcesFile(configuration->configpath, SortIndex);
     742    BondFragments->StoreForcesFile(prefix, SortIndex);
    748743
    749744    // store keysets file
    750     StoreKeySetFile(TotalGraph, configuration->configpath);
     745    StoreKeySetFile(TotalGraph, prefix);
    751746
    752747    {
    753748      // store Adjacency file
    754       char filename[MAXSTRINGSIZE];
    755       strcpy(filename, FRAGMENTPREFIX);
    756       strcat(filename, ADJACENCYFILE);
    757       StoreAdjacencyToFile(configuration->configpath, filename);
     749      std::string filename = prefix + ADJACENCYFILE;
     750      StoreAdjacencyToFile(filename);
    758751    }
    759752
    760753    // store Hydrogen saturation correction file
    761     BondFragments->AddHydrogenCorrection(configuration->configpath);
     754    BondFragments->AddHydrogenCorrection(prefix);
    762755
    763756    // store adaptive orders into file
    764     StoreOrderAtSiteFile(configuration->configpath);
     757    StoreOrderAtSiteFile(prefix);
    765758
    766759    // restore orbital and Stop values
    767     CalculateOrbitals(*configuration);
     760    //CalculateOrbitals(*configuration);
    768761
    769762    // free memory for bond part
    … …  
    782775/** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
    783776 * Atoms not present in the file get "-1".
    784  * \param *out output stream for debugging
    785  * \param *path path to file ORDERATSITEFILE
     777 * \param &path path to file ORDERATSITEFILE
    786778 * \return true - file writable, false - not writable
    787779 */
    788 bool molecule::StoreOrderAtSiteFile(char *path)
    789 {
    790   stringstream line;
     780bool molecule::StoreOrderAtSiteFile(std::string &path)
     781{
     782  string line;
    791783  ofstream file;
    792784
    793   line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
    794   file.open(line.str().c_str());
     785  line = path + ORDERATSITEFILE;
     786  file.open(line.c_str());
    795787  DoLog(1) && (Log() << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl);
    796   if (file != NULL) {
     788  if (file.good()) {
    797789    ActOnAllAtoms( &atom::OutputOrder, &file );
    798790    file.close();
    … …  
    800792    return true;
    801793  } else {
    802     DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line.str() << "." << endl);
     794    DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line << "." << endl);
    803795    return false;
    804796  }
    … …  
    807799/** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
    808800 * Atoms not present in the file get "0".
    809  * \param *out output stream for debugging
    810  * \param *path path to file ORDERATSITEFILEe
     801 * \param &path path to file ORDERATSITEFILEe
    811802 * \return true - file found and scanned, false - file not found
    812803 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
    813804 */
    814 bool molecule::ParseOrderAtSiteFromFile(char *path)
     805bool molecule::ParseOrderAtSiteFromFile(std::string &path)
    815806{
    816807  unsigned char *OrderArray = new unsigned char[getAtomCount()];
    … …  
    818809  bool status;
    819810  int AtomNr, value;
    820   stringstream line;
     811  string line;
    821812  ifstream file;
    822813
    … …  
    827818
    828819  DoLog(1) && (Log() << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl);
    829   line << path << "/" << FRAGMENTPREFIX << ORDERATSITEFILE;
    830   file.open(line.str().c_str());
    831   if (file != NULL) {
     820  line = path + ORDERATSITEFILE;
     821  file.open(line.c_str());
     822  if (file.good()) {
    832823    while (!file.eof()) { // parse from file
    833824      AtomNr = -1;
    … …  
    850841    status = true;
    851842  } else {
    852     DoLog(1) && (Log() << Verbose(1) << "\t ... failed to open file " << line.str() << "." << endl);
     843    DoLog(1) && (Log() << Verbose(1) << "\t ... failed to open file " << line << "." << endl);
    853844    status = false;
    854845  }

  • src/molecule_graph.cpp

    r6d574a r0d1ad0  
    1212#include "bondgraph.hpp"
    1313#include "config.hpp"
     14#include "defs.hpp"
    1415#include "element.hpp"
    1516#include "helpers.hpp"
    … …  
    10241025/** Storing the bond structure of a molecule to file.
    10251026 * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
    1026  * \param *path path to file
    1027  * \param *filename name of file
     1027 * \param &filename name of file
     1028 * \param path path to file, defaults to empty
    10281029 * \return true - file written successfully, false - writing failed
    10291030 */
    1030 bool molecule::StoreAdjacencyToFile(char *path, char *filename)
     1031bool molecule::StoreAdjacencyToFile(std::string &filename, std::string path)
    10311032{
    10321033  ofstream AdjacencyFile;
    1033   stringstream line;
     1034  string line;
    10341035  bool status = true;
    10351036
    1036   if (path != NULL)
    1037     line << path << "/" << filename;
     1037  if (path != "")
     1038    line = path + "/" + filename;
    10381039  else
    1039     line << filename;
    1040   AdjacencyFile.open(line.str().c_str(), ios::out);
     1040    line = filename;
     1041  AdjacencyFile.open(line.c_str(), ios::out);
    10411042  DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
    1042   if (AdjacencyFile != NULL) {
     1043  if (AdjacencyFile.good()) {
    10431044    AdjacencyFile << "m\tn" << endl;
    10441045    ActOnAllAtoms(&atom::OutputAdjacency, &AdjacencyFile);
    … …  
    10461047    DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
    10471048  } else {
    1048     DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line.str() << "." << endl);
     1049    DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
    10491050    status = false;
    10501051  }
    … …  
    10561057/** Storing the bond structure of a molecule to file.
    10571058 * Simply stores Atom::nr and then the Atom::nr of all bond partners, one per line.
    1058  * \param *path path to file
    1059  * \param *filename name of file
     1059 * \param &filename name of file
     1060 * \param path path to file, defaults to empty
    10601061 * \return true - file written successfully, false - writing failed
    10611062 */
    1062 bool molecule::StoreBondsToFile(char *path, char *filename)
     1063bool molecule::StoreBondsToFile(std::string &filename, std::string path)
    10631064{
    10641065  ofstream BondFile;
    1065   stringstream line;
     1066  string line;
    10661067  bool status = true;
    10671068
    1068   if (path != NULL)
    1069     line << path << "/" << filename;
     1069  if (path != "")
     1070    line = path + "/" + filename;
    10701071  else
    1071     line << filename;
    1072   BondFile.open(line.str().c_str(), ios::out);
     1072    line = filename;
     1073  BondFile.open(line.c_str(), ios::out);
    10731074  DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
    1074   if (BondFile != NULL) {
     1075  if (BondFile.good()) {
    10751076    BondFile << "m\tn" << endl;
    10761077    ActOnAllAtoms(&atom::OutputBonds, &BondFile);
    … …  
    10781079    DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
    10791080  } else {
    1080     DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line.str() << "." << endl);
     1081    DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
    10811082    status = false;
    10821083  }
    … …  
    10861087;
    10871088
    1088 bool CheckAdjacencyFileAgainstMolecule_Init(char *path, ifstream &File, int *&CurrentBonds)
    1089 {
    1090   stringstream filename;
    1091   filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
    1092   File.open(filename.str().c_str(), ios::out);
     1089bool CheckAdjacencyFileAgainstMolecule_Init(std::string &path, ifstream &File, int *&CurrentBonds)
     1090{
     1091  string filename;
     1092  filename = path + ADJACENCYFILE;
     1093  File.open(filename.c_str(), ios::out);
    10931094  DoLog(1) && (Log() << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... " << endl);
    1094   if (File == NULL)
     1095  if (File.fail())
    10951096    return false;
    10961097
    … …  
    11461147 * \return true - structure is equal, false - not equivalence
    11471148 */
    1148 bool molecule::CheckAdjacencyFileAgainstMolecule(char *path, atom **ListOfAtoms)
     1149bool molecule::CheckAdjacencyFileAgainstMolecule(std::string &path, atom **ListOfAtoms)
    11491150{
    11501151  ifstream File;

  • src/moleculelist.cpp

    r6d574a r0d1ad0  
    1212#include "atom.hpp"
    1313#include "bond.hpp"
     14#include "bondgraph.hpp"
    1415#include "boundary.hpp"
    1516#include "config.hpp"
    … …  
    379380 * bonded to the same atom, then we add for this pair a correction term constructed from a Morse
    380381 * potential function fit to QM calculations with respecting to the interatomic hydrogen distance.
    381  * \param *out output stream for debugging
    382  * \param *path path to file
    383  */
    384 bool MoleculeListClass::AddHydrogenCorrection(char *path)
     382 * \param &path path to file
     383 */
     384bool MoleculeListClass::AddHydrogenCorrection(std::string &path)
    385385{
    386386  bond *Binder = NULL;
    … …  
    400400  // 0a. find dimension of matrices with constants
    401401  line = path;
    402   line.append("/");
    403   line += FRAGMENTPREFIX;
    404402  line += "1";
    405403  line += FITCONSTANTSUFFIX;
    406404  input.open(line.c_str());
    407   if (input == NULL) {
     405  if (input.fail()) {
    408406    DoLog(1) && (Log() << Verbose(1) << endl << "Unable to open " << line << ", is the directory correct?" << endl);
    409407    return false;
    … …  
    569567
    570568/** Store force indices, i.e. the connection between the nuclear index in the total molecule config and the respective atom in fragment config.
    571  * \param *out output stream for debugging
    572  * \param *path path to file
     569 * \param &path path to file
    573570 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
    574571 * \return true - file written successfully, false - writing failed
    575572 */
    576 bool MoleculeListClass::StoreForcesFile(char *path,
    577     int *SortIndex)
     573bool MoleculeListClass::StoreForcesFile(std::string &path, int *SortIndex)
    578574{
    579575  bool status = true;
    580   ofstream ForcesFile;
    581   stringstream line;
     576  string filename(path);
     577  filename += FORCESFILE;
     578  ofstream ForcesFile(filename.c_str());
    582579  periodentafel *periode=World::getInstance().getPeriode();
    583580
    584581  // open file for the force factors
    585582  DoLog(1) && (Log() << Verbose(1) << "Saving  force factors ... ");
    586   line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
    587   ForcesFile.open(line.str().c_str(), ios::out);
    588   if (ForcesFile != NULL) {
     583  if (!ForcesFile.fail()) {
    589584    //Log() << Verbose(1) << "Final AtomicForcesList: ";
    590585    //output << prefix << "Forces" << endl;
    … …  
    611606  } else {
    612607    status = false;
    613     DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line.str() << "." << endl);
     608    DoLog(1) && (Log() << Verbose(1) << "failed to open file " << filename << "." << endl);
    614609  }
    615610  ForcesFile.close();
    … …  
    620615/** Writes a config file for each molecule in the given \a **FragmentList.
    621616 * \param *out output stream for debugging
    622  * \param *configuration standard configuration to attach atoms in fragment molecule to.
     617 * \param &prefix path and prefix to the fragment config files
    623618 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
    624619 * \return true - success (each file was written), false - something went wrong.
    625620 */
    626 bool MoleculeListClass::OutputConfigForListOfFragments(config *configuration, int *SortIndex)
     621bool MoleculeListClass::OutputConfigForListOfFragments(std::string &prefix, int *SortIndex)
    627622{
    628623  ofstream outputFragment;
    629   char FragmentName[MAXSTRINGSIZE];
     624  std::string FragmentName;
    630625  char PathBackup[MAXSTRINGSIZE];
    631626  bool result = true;
    … …  
    645640  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    646641    // save default path as it is changed for each fragment
    647     path = configuration->GetDefaultPath();
     642    path = World::getInstance().getConfig()->GetDefaultPath();
    648643    if (path != NULL)
    649644      strcpy(PathBackup, path);
    … …  
    658653    // output xyz file
    659654    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
    660     sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    661     outputFragment.open(FragmentName, ios::out);
     655    FragmentName = prefix + FragmentNumber + ".conf.xyz";
     656    outputFragment.open(FragmentName.c_str(), ios::out);
    662657    DoLog(2) && (Log() << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ...");
    663658    if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
    … …  
    682677    for (int k = 0; k < NDIM; k++) {
    683678      j += k + 1;
    684       BoxDimension[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
     679      BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
    685680      cell_size[j] = BoxDimension[k] * 2.;
    686681    }
    … …  
    689684    // also calculate necessary orbitals
    690685    (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
    691     (*ListRunner)->CalculateOrbitals(*configuration);
     686    //(*ListRunner)->CalculateOrbitals(*World::getInstance().getConfig);
    692687
    693688    // change path in config
    694     //strcpy(PathBackup, configuration->configpath);
    695     sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
    696     configuration->SetDefaultPath(FragmentName);
     689    FragmentName = PathBackup;
     690    FragmentName += "/";
     691    FragmentName += FRAGMENTPREFIX;
     692    FragmentName += FragmentNumber;
     693    FragmentName += "/";
     694    World::getInstance().getConfig()->SetDefaultPath(FragmentName.c_str());
    697695
    698696    // and save as config
    699     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     697    FragmentName = prefix + FragmentNumber + ".conf";
    700698    DoLog(2) && (Log() << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...");
    701     if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner))))
     699    if ((intermediateResult = World::getInstance().getConfig()->Save(FragmentName.c_str(), (*ListRunner)->elemente, (*ListRunner))))
    702700      DoLog(0) && (Log() << Verbose(0) << " done." << endl);
    703701    else
    … …  
    706704
    707705    // restore old config
    708     configuration->SetDefaultPath(PathBackup);
     706    World::getInstance().getConfig()->SetDefaultPath(PathBackup);
    709707
    710708    // and save as mpqc input file
    711     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     709    FragmentName = prefix + FragmentNumber + ".conf";
    712710    DoLog(2) && (Log() << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...");
    713     if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner))))
     711    if ((intermediateResult = World::getInstance().getConfig()->SaveMPQC(FragmentName.c_str(), (*ListRunner))))
    714712      DoLog(2) && (Log() << Verbose(2) << " done." << endl);
    715713    else
    … …  
    767765
    768766  // 1. dissect the molecule into connected subgraphs
    769   if (!configuration->BG->ConstructBondGraph(mol)) {
    770     World::getInstance().destroyMolecule(mol);
    771     DoeLog(1) && (eLog()<< Verbose(1) << "There are no bonds." << endl);
     767  if (configuration->BG != NULL) {
     768    if (!configuration->BG->ConstructBondGraph(mol)) {
     769      World::getInstance().destroyMolecule(mol);
     770      DoeLog(1) && (eLog()<< Verbose(1) << "There are no bonds." << endl);
     771      return;
     772    }
     773  } else {
     774    DoeLog(1) && (eLog()<< Verbose(1) << "There is no BondGraph class present to create bonds." << endl);
    772775    return;
    773776  }

  • src/periodentafel.cpp

    r6d574a r0d1ad0  
    3232periodentafel::periodentafel()
    3333{
    34   bool status = true;
    35   status = LoadElementsDatabase(new stringstream(elementsDB,ios_base::in));
    36   ASSERT(status,  "General element initialization failed");
    37   status = LoadValenceDatabase(new stringstream(valenceDB,ios_base::in));
    38   ASSERT(status, "Valence entry of element initialization failed");
    39   status = LoadOrbitalsDatabase(new stringstream(orbitalsDB,ios_base::in));
    40   ASSERT(status, "Orbitals entry of element initialization failed");
    41   status = LoadHBondAngleDatabase(new stringstream(HbondangleDB,ios_base::in));
    42   ASSERT(status, "HBond angle entry of element initialization failed");
    43   status = LoadHBondLengthsDatabase(new stringstream(HbonddistanceDB,ios_base::in));
    44   ASSERT(status, "HBond distance entry of element initialization failed");
     34  {
     35    stringstream input(elementsDB,ios_base::in);
     36    bool status = LoadElementsDatabase(&input);
     37    ASSERT(status,  "General element initialization failed");
     38  }
     39  {
     40    stringstream input(valenceDB,ios_base::in);
     41    bool status = LoadValenceDatabase(&input);
     42    ASSERT(status, "Valence entry of element initialization failed");
     43  }
     44  {
     45    stringstream input(orbitalsDB,ios_base::in);
     46    bool status = LoadOrbitalsDatabase(&input);
     47    ASSERT(status, "Orbitals entry of element initialization failed");
     48  }
     49  {
     50    stringstream input(HbondangleDB,ios_base::in);
     51    bool status = LoadHBondAngleDatabase(&input);
     52    ASSERT(status, "HBond angle entry of element initialization failed");
     53  }
     54  {
     55    stringstream input(HbonddistanceDB,ios_base::in);
     56    bool status = LoadHBondLengthsDatabase(&input);
     57    ASSERT(status, "HBond distance entry of element initialization failed");
     58  }
    4559};
    4660
    … …  
    333347          ASSERT(Elemental != NULL, "element should be present but is not??");
    334348          *Elemental = *neues;
     349          delete(neues);
     350          neues = Elemental;
    335351        } else {
    336352          InserterTest = elements.insert(pair <atomicNumber_t,element*> (neues->getNumber(), neues));

  • src/unittests/LinkedCellUnitTest.cpp

    r6d574a r0d1ad0  
    264264  Vector tester;
    265265  LinkedCell::LinkedNodes *ListOfPoints = NULL;
    266   atom *Walker = NULL;
    267266  size_t size = 0;
    268267
    … …  
    326325  Vector tester;
    327326  LinkedCell::LinkedNodes *ListOfPoints = NULL;
    328   atom *Walker = NULL;
    329327  size_t size = 0;
    330328

  • src/unittests/Makefile.am

    r6d574a r0d1ad0  
    5252GSLLIBS = ../libgslwrapper.a $(BOOST_LIB) ${BOOST_THREAD_LIB}
    5353ALLLIBS = ../libmolecuilder.a ${GSLLIBS}
     54PARSERLIBS = ../libparser.a ${ALLLIBS}
    5455
    5556TESTSOURCES = \
    … …  
    201202
    202203ParserUnitTest_SOURCES = UnitTestMain.cpp ParserUnitTest.cpp ParserUnitTest.hpp
    203 ParserUnitTest_LDADD = ${ALLLIBS}
     204ParserUnitTest_LDADD = ${PARSERLIBS}
    204205
    205206periodentafelTest_SOURCES = UnitTestMain.cpp periodentafelTest.cpp periodentafelTest.hpp

  • src/unittests/ParserUnitTest.cpp

    r6d574a r0d1ad0  
    1212#include <cppunit/ui/text/TestRunner.h>
    1313
     14#include "Parser/MpqcParser.hpp"
     15#include "Parser/PcpParser.hpp"
     16#include "Parser/TremoloParser.hpp"
    1417#include "Parser/XyzParser.hpp"
    15 #include "Parser/TremoloParser.hpp"
    1618#include "World.hpp"
    1719#include "atom.hpp"
    … …  
    2931CPPUNIT_TEST_SUITE_REGISTRATION( ParserUnitTest );
    3032
     33static string waterPcp = "# ParallelCarParinello - main configuration file - created with molecuilder\n\
     34\n\
     35mainname\tpcp\t# programm name (for runtime files)\n\
     36defaultpath\not specified\t# where to put files during runtime\n\
     37pseudopotpath\not specified\t# where to find pseudopotentials\n\
     38\n\
     39ProcPEGamma\t8\t# for parallel computing: share constants\n\
     40ProcPEPsi\t1\t# for parallel computing: share wave functions\n\
     41DoOutVis\t0\t# Output data for OpenDX\n\
     42DoOutMes\t1\t# Output data for measurements\n\
     43DoOutOrbitals\t0\t# Output all Orbitals\n\
     44DoOutCurr\t0\t# Ouput current density for OpenDx\n\
     45DoOutNICS\t0\t# Output Nucleus independent current shieldings\n\
     46DoPerturbation\t0\t# Do perturbation calculate and determine susceptibility and shielding\n\
     47DoFullCurrent\t0\t# Do full perturbation\n\
     48DoConstrainedMD\t0\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD\n\
     49Thermostat\tBerendsen\t2.5\t# Which Thermostat and its parameters to use in MD case.\n\
     50CommonWannier\t0\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center\n\
     51SawtoothStart\t0.01\t# Absolute value for smooth transition at cell border \n\
     52VectorPlane\t0\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot\n\
     53VectorCut\t0\t# Cut plane axis value\n\
     54AddGramSch\t1\t# Additional GramSchmidtOrtogonalization to be safe\n\
     55Seed\t1\t# initial value for random seed for Psi coefficients\n\
     56\n\
     57MaxOuterStep\t0\t# number of MolecularDynamics/Structure optimization steps\n\
     58Deltat\t0.01\t# time per MD step\n\
     59OutVisStep\t10\t# Output visual data every ...th step\n\
     60OutSrcStep\t5\t# Output \"restart\" data every ..th step\n\
     61TargetTemp\t0.000950045\t# Target temperature\n\
     62MaxPsiStep\t3\t# number of Minimisation steps per state (0 - default)\n\
     63EpsWannier\t1e-07\t# tolerance value for spread minimisation of orbitals\n\
     64# Values specifying when to stop\n\
     65MaxMinStep\t100\t# Maximum number of steps\n\
     66RelEpsTotalE\t1e-07\t# relative change in total energy\n\
     67RelEpsKineticE\t1e-05\t# relative change in kinetic energy\n\
     68MaxMinStopStep\t2\t# check every ..th steps\n\
     69MaxMinGapStopStep\t1\t# check every ..th steps\n\
     70\n\
     71# Values specifying when to stop for INIT, otherwise same as above\n\
     72MaxInitMinStep\t100\t# Maximum number of steps\n\
     73InitRelEpsTotalE\t1e-05\t# relative change in total energy\n\
     74InitRelEpsKineticE\t0.0001\t# relative change in kinetic energy\n\
     75InitMaxMinStopStep\t2\t# check every ..th steps\n\
     76InitMaxMinGapStopStep\t1\t# check every ..th steps\n\
     77\n\
     78BoxLength\t# (Length of a unit cell)\n\
     7920\n\
     800\t20\n\
     810\t0\t20\n\
     82\n\
     83ECut\t128\t# energy cutoff for discretization in Hartrees\n\
     84MaxLevel\t5\t# number of different levels in the code, >=2\n\
     85Level0Factor\t2\t# factor by which node number increases from S to 0 level\n\
     86RiemannTensor\t0\t# (Use metric)\n\
     87PsiType\t0\t# 0 - doubly occupied, 1 - SpinUp,SpinDown\n\
     88MaxPsiDouble\t2\t# here: specifying both maximum number of SpinUp- and -Down-states\n\
     89PsiMaxNoUp\t2\t# here: specifying maximum number of SpinUp-states\n\
     90PsiMaxNoDown\t2\t# here: specifying maximum number of SpinDown-states\n\
     91AddPsis\t0\t# Additional unoccupied Psis for bandgap determination\n\
     92\n\
     93RCut\t20\t# R-cut for the ewald summation\n\
     94StructOpt\t0\t# Do structure optimization beforehand\n\
     95IsAngstroem\t1\t# 0 - Bohr, 1 - Angstroem\n\
     96RelativeCoord\t0\t# whether ion coordinates are relative (1) or absolute (0)\n\
     97MaxTypes\t2\t# maximum number of different ion types\n\
     98\n\
     99# Ion type data (PP = PseudoPotential, Z = atomic number)\n\
     100#Ion_TypeNr.\tAmount\tZ\tRGauss\tL_Max(PP)L_Loc(PP)IonMass\t# chemical name, symbol\n\
     101Ion_Type1\t2\t1\t1.0\t3\t3\t1.008\tHydrogen\tH\n\
     102Ion_Type2\t1\t8\t1.0\t3\t3\t15.999\tOxygen\tO\n\
     103#Ion_TypeNr._Nr.R[0]\tR[1]\tR[2]\tMoveType (0 MoveIon, 1 FixedIon)\n\
     104Ion_Type2_1\t0.000000000\t0.000000000\t0.000000000\t0 # molecule nr 0\n\
     105Ion_Type1_1\t0.758602\t0.000000000\t0.504284\t0 # molecule nr 1\n\
     106Ion_Type1_2\t0.758602\t0.000000000\t-0.504284\t0 # molecule nr 2\n";
     107static string waterMpqc ="% Created by MoleCuilder\n\
     108mpqc: (\n\
     109\tsavestate = no\n\
     110\tdo_gradient = yes\n\
     111\tmole<MBPT2>: (\n\
     112\t\tmaxiter = 200\n\
     113\t\tbasis = $:basis\n\
     114\t\tmolecule = $:molecule\n\
     115\t\treference<CLHF>: (\n\
     116\t\t\tbasis = $:basis\n\
     117\t\t\tmolecule = $:molecule\n\
     118\t\t)\n\
     119\t)\n\
     120)\n\
     121molecule<Molecule>: (\n\
     122\tunit = angstrom\n\
     123\t{ atoms geometry } = {\n\
     124\t\tO [ -0.505735\t0\t0 ]\n\
     125\t\tH [ 0.252867\t0\t0.504284 ]\n\
     126\t\tH [ 0.252867\t0\t-0.504284 ]\n\
     127\t}\n\
     128)\n\
     129basis<GaussianBasisSet>: (\n\
     130\tname = \"3-21G\"\n\
     131\tmolecule = $:molecule\n\
     132)\n";
     133static string waterXyz = "3\nH2O: water molecule\nO\t0\t0\t0\nH\t0.758602\t0\t0.504284\nH\t0.758602\t0\t-0.504284\n";
     134static string Tremolo_Atomdata1 = "# ATOMDATA\tId\tname\tType\tx=3\n";
     135static string Tremolo_Atomdata2 = "#\n#ATOMDATA Id name Type x=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
     136static string Tremolo_invalidkey = "#\n#ATOMDATA Id name foo Type x=3\n\n\n";
     137static string Tremolo_velocity = "#\n#ATOMDATA Id name Type u=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
     138static string Tremolo_neighbours = "#\n#ATOMDATA Id Type neighbors=2\n1 H 3 0\n2 H 3 0\n3 O 1 2\n";
     139static string Tremolo_improper = "#\n#ATOMDATA Id Type imprData\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
     140static string Tremolo_torsion = "#\n#ATOMDATA Id Type torsion\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
     141static string Tremolo_full = "# ATOMDATA\tx=3\tu=3\tF\tstress\tId\tneighbors=5\timprData\tGroupMeasureTypeNo\tType\textType\tname\tresName\tchainID\tresSeq\toccupancy\ttempFactor\tsegID\tCharge\tcharge\tGrpTypeNo\ttorsion\n0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t-\t0\tH\t-\t-\t-\t0\t0\t0\t0\t0\t0\t0\t0\t-\t\n";
    31142
    32143void ParserUnitTest::setUp() {
    33144  World::getInstance();
     145
     146  // we need hydrogens and oxygens in the following tests
     147  CPPUNIT_ASSERT(World::getInstance().getPeriode()->FindElement(1) != NULL);
     148  CPPUNIT_ASSERT(World::getInstance().getPeriode()->FindElement(8) != NULL);
    34149}
    35150
    36151void ParserUnitTest::tearDown() {
     152  ChangeTracker::purgeInstance();
    37153  World::purgeInstance();
    38154}
    … …  
    43159  cout << "Testing the XYZ parser." << endl;
    44160  XyzParser* testParser = new XyzParser();
    45   string waterXyz = "3\nH2O: water molecule\nO\t0.000000\t0.000000\t0.000000\nH\t0.758602\t0.000000\t0.504284\nH\t0.758602\t0.000000\t-0.504284\n";
    46161  stringstream input;
    47162  input << waterXyz;
    … …  
    62177  TremoloParser* testParser = new TremoloParser();
    63178  stringstream input, output;
    64   string waterTremolo;
    65179
    66180  // Atomdata beginning with "# ATOMDATA"
    67   waterTremolo = "# ATOMDATA\tId\tname\tType\tx=3\n";
    68   input << waterTremolo;
    69   testParser->load(&input);
    70   testParser->save(&output);
    71   CPPUNIT_ASSERT(waterTremolo == output.str());
     181  input << Tremolo_Atomdata1;
     182  testParser->load(&input);
     183  testParser->save(&output);
     184  CPPUNIT_ASSERT(Tremolo_Atomdata1 == output.str());
    72185  input.clear();
    73186  output.clear();
    74187
    75188  // Atomdata beginning with "#ATOMDATA"
    76   waterTremolo = "#\n#ATOMDATA Id name Type x=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
    77   input << waterTremolo;
     189  input << Tremolo_Atomdata2;
    78190  testParser->load(&input);
    79191  testParser->save(&output);
    … …  
    83195
    84196  // Invalid key in Atomdata line
    85   waterTremolo = "#\n#ATOMDATA Id name foo Type x=3\n\n\n";
    86   input << waterTremolo;
     197  input << Tremolo_invalidkey;
    87198  testParser->load(&input);
    88199  //TODO: proove invalidity
    … …  
    93204  TremoloParser* testParser = new TremoloParser();
    94205  stringstream input;
    95   string waterTremolo;
    96206
    97207  // One simple data line
    98   waterTremolo = "#\n#ATOMDATA Id name Type x=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
    99   input << waterTremolo;
     208  input << Tremolo_Atomdata2;
    100209  testParser->load(&input);
    101210  CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->x[0] == 3.0);
    … …  
    106215  TremoloParser* testParser = new TremoloParser();
    107216  stringstream input;
    108   string waterTremolo;
    109217
    110218  // One simple data line
    111   waterTremolo = "#\n#ATOMDATA Id name Type u=3\n1 hydrogen H 3.0 4.5 0.1\n\n";
    112   input << waterTremolo;
     219  input << Tremolo_velocity;
    113220  testParser->load(&input);
    114221  CPPUNIT_ASSERT(World::getInstance().getAtom(AtomByType(1))->v[0] == 3.0);
    … …  
    119226  TremoloParser* testParser = new TremoloParser();
    120227  stringstream input;
    121   string waterTremolo;
    122228
    123229  // Neighbor data
    124   waterTremolo = "#\n#ATOMDATA Id Type neighbors=2\n1 H 3 0\n2 H 3 0\n3 O 1 2\n";
    125   input << waterTremolo;
     230  input << Tremolo_neighbours;
    126231  testParser->load(&input);
    127232
    … …  
    135240  TremoloParser* testParser = new TremoloParser();
    136241  stringstream input, output;
    137   string waterTremolo;
    138242
    139243  // Neighbor data
    140   waterTremolo = "#\n#ATOMDATA Id Type imprData\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
    141   input << waterTremolo;
     244  input << Tremolo_improper;
    142245  testParser->load(&input);
    143246  testParser->save(&output);
    … …  
    151254  TremoloParser* testParser = new TremoloParser();
    152255  stringstream input, output;
    153   string waterTremolo;
    154256
    155257  // Neighbor data
    156   waterTremolo = "#\n#ATOMDATA Id Type torsion\n8 H 9-10\n9 H 10-8,8-10\n10 O -\n";
    157   input << waterTremolo;
     258  input << Tremolo_torsion;
    158259  testParser->load(&input);
    159260  testParser->save(&output);
    … …  
    173274  testParser->setFieldsForSave("x=3 u=3 F stress Id neighbors=5 imprData GroupMeasureTypeNo Type extType name resName chainID resSeq occupancy tempFactor segID Charge charge GrpTypeNo torsion");
    174275  testParser->save(&output);
    175   CPPUNIT_ASSERT(output.str() == "# ATOMDATA\tx=3\tu=3\tF\tstress\tId\tneighbors=5\timprData\tGroupMeasureTypeNo\tType\textType\tname\tresName\tchainID\tresSeq\toccupancy\ttempFactor\tsegID\tCharge\tcharge\tGrpTypeNo\ttorsion\n0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t-\t0\tH\t-\t-\t-\t0\t0\t0\t0\t0\t0\t0\t0\t-\t\n");
     276  CPPUNIT_ASSERT(output.str() == Tremolo_full);
    176277
    177278  cout << "testing the tremolo parser is done" << endl;
    178279}
     280
     281void ParserUnitTest::readwritePcpTest() {
     282  stringstream input(waterPcp);
     283  PcpParser* testParser = new PcpParser();
     284  testParser->load(&input);
     285  input.clear();
     286
     287  CPPUNIT_ASSERT_EQUAL(3, World::getInstance().numAtoms());
     288
     289  string newWaterPcp = "";
     290  stringstream output;
     291  testParser->save(&output);
     292
     293  input << output;
     294  PcpParser* testParser2 = new PcpParser();
     295  testParser2->load(&input);
     296
     297  CPPUNIT_ASSERT_EQUAL(6, World::getInstance().numAtoms());
     298
     299  CPPUNIT_ASSERT(*testParser == *testParser2);
     300}
     301
     302void ParserUnitTest::writeMpqcTest() {
     303  // build up water molecule
     304  atom *Walker = NULL;
     305  Walker = World::getInstance().createAtom();
     306  Walker->type = World::getInstance().getPeriode()->FindElement(8);
     307  Walker->x = Vector(0,0,0);
     308  Walker = World::getInstance().createAtom();
     309  Walker->type = World::getInstance().getPeriode()->FindElement(1);
     310  Walker->x = Vector(0.758602,0,0.504284);
     311  Walker = World::getInstance().createAtom();
     312  Walker->type = World::getInstance().getPeriode()->FindElement(1);
     313  Walker->x = Vector(0.758602,0,-0.504284);
     314  CPPUNIT_ASSERT_EQUAL(3, World::getInstance().numAtoms());
     315
     316  // create two stringstreams, one stored, one created
     317  stringstream input(waterMpqc);
     318  MpqcParser* testParser = new MpqcParser();
     319  stringstream output;
     320  testParser->save(&output);
     321
     322  // compare both configs
     323  string first = input.str();
     324  string second = output.str();
     325  CPPUNIT_ASSERT(first == second);
     326}

  • src/unittests/ParserUnitTest.hpp

    r6d574a r0d1ad0  
    2222  CPPUNIT_TEST ( readAndWriteTremoloTorsionInformationTest );
    2323  CPPUNIT_TEST ( writeTremoloTest );
     24  CPPUNIT_TEST ( readwritePcpTest );
     25  CPPUNIT_TEST ( writeMpqcTest );
    2426  CPPUNIT_TEST_SUITE_END();
    2527
    … …  
    3638  void readAndWriteTremoloTorsionInformationTest();
    3739  void writeTremoloTest();
     40  void readwritePcpTest();
     41  void writeMpqcTest();
    3842};
    3943
Note: See TracChangeset for help on using the changeset viewer.