Changeset 1907a7 for src


Ignore:
Timestamp:
Apr 2, 2009, 4:12:54 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
ca3ccc
Parents:
d8b94a
Message:

Basic implementation of Multiple molecules.

builder.cpp:

molecules.hpp:

moleculelist.cpp:

  • replaced listofmolecules array by STL list everywhere (only smaller changes necessary)
  • new merging function: SimpleMerge, SimpleAdd, SimpleMultiMerge, SimpleMultiAdd, (EmbedMerge, ScatterMerge ... both not finished). Add does not while merge does delete the src molecules.
  • new function: Enumerate(). Output of all molecules with number of atoms and chemical formula
  • new function: NumberOfActiveMolecules(). Counts the number of molecules in the list with ActiveFlag set.
  • new function: insert(). Inserts a molecule into the list with a unique index

molecules.cpp:

  • new function: SetNameFromFilename. Takes basename of a filename and sets name accordingly.
  • new function: UnlinkAtom. Only removes atom from list, does not delete it from memory.

atom.cpp:

  • Output() also accepts specific comment instead of "# molecule nr ..."
Location:
src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • src/atom.cpp

    rd8b94a r1907a7  
    11/** \file atom.cpp
    2  * 
     2 *
    33 * Function implementations for the class atom.
    4  * 
     4 *
    55 */
    66
    77#include "molecules.hpp"
    8  
     8
    99/************************************* Functions for class atom *************************************/
    1010
    1111/** Constructor of class atom.
    1212 */
    13 atom::atom() 
     13atom::atom()
    1414{
    1515        Name = NULL;
    … …  
    3333/** Destructor of class atom.
    3434 */
    35 atom::~atom() 
     35atom::~atom()
    3636{
    3737        Free((void **)&Name, "atom::~atom: *Name");
    … …  
    5858 * \param AtomNo cardinal number among these atoms of the same element
    5959 * \param *out stream to output to
     60 * \param *comment commentary after '#' sign
    6061 */
    61 bool atom::Output(int ElementNo, int AtomNo, ofstream *out) const
     62bool atom::Output(int ElementNo, int AtomNo, ofstream *out, const char *comment) const
    6263{
    6364        if (out != NULL) {
    … …  
    6768                if (v.Norm() > MYEPSILON)
    6869                        *out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
    69                 *out << " # Number in molecule " << nr << endl;
     70                if (comment != NULL)
     71                  *out << " # " << comment << endl;
     72                else
     73                  *out << " # molecule nr " << nr << endl;
    7074                return true;
    7175        } else
    … …  
    8589};
    8690
    87 ostream & operator << (ostream &ost, atom &a) 
     91ostream & operator << (ostream &ost, atom &a)
    8892{
    8993        ost << "[" << a.Name << "|" << &a << "]";
    … …  
    9498 * \param ptr atom to compare index against
    9599 * \return true - this one's is smaller, false - not
    96  */ 
     100 */
    97101bool atom::Compare(atom &ptr)
    98102{
    … …  
    103107};
    104108
    105 bool operator < (atom &a, atom &b) 
     109bool operator < (atom &a, atom &b)
    106110{
    107111        return a.Compare(b);

  • src/builder.cpp

    rd8b94a r1907a7  
    5555#include "molecules.hpp"
    5656
    57 /********************************************** Submenu routine **************************************/
     57/********************************************* Subsubmenu routine ************************************/
    5858
    5959/** Submenu for adding atoms to the molecule.
    6060 * \param *periode periodentafel
    61  * \param *mol the molecule to add to
     61 * \param *molecule molecules with atoms
    6262 */
    6363static void AddAtoms(periodentafel *periode, molecule *mol)
    … …  
    8383
    8484        switch (choice) {
     85    default:
     86      cout << Verbose(0) << "Not a valid choice." << endl;
     87      break;
    8588                        case 'a': // absolute coordinates of atom
    86                                 cout << Verbose(0) << "Enter absolute coordinates." << endl;
    87                                 first = new atom;
    88                                 first->x.AskPosition(mol->cell_size, false);
    89                                 first->type = periode->AskElement();    // give type
    90                                 mol->AddAtom(first);    // add to molecule
     89        cout << Verbose(0) << "Enter absolute coordinates." << endl;
     90        first = new atom;
     91        first->x.AskPosition(mol->cell_size, false);
     92        first->type = periode->AskElement();    // give type
     93        mol->AddAtom(first);    // add to molecule
    9194                                break;
    9295
    9396                        case 'b': // relative coordinates of atom wrt to reference point
    94                                 first = new atom;
    95                                 valid = true;
    96                                 do {
    97                                         if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
    98                                         cout << Verbose(0) << "Enter reference coordinates." << endl;
    99                                         x.AskPosition(mol->cell_size, true);
    100                                         cout << Verbose(0) << "Enter relative coordinates." << endl;
    101                                         first->x.AskPosition(mol->cell_size, false);
    102                                         first->x.AddVector((const Vector *)&x);
    103                                         cout << Verbose(0) << "\n";
    104                                 } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    105                                 first->type = periode->AskElement();    // give type
    106                                 mol->AddAtom(first);    // add to molecule
     97        first = new atom;
     98        valid = true;
     99        do {
     100          if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
     101          cout << Verbose(0) << "Enter reference coordinates." << endl;
     102          x.AskPosition(mol->cell_size, true);
     103          cout << Verbose(0) << "Enter relative coordinates." << endl;
     104          first->x.AskPosition(mol->cell_size, false);
     105          first->x.AddVector((const Vector *)&x);
     106          cout << Verbose(0) << "\n";
     107        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     108        first->type = periode->AskElement();    // give type
     109        mol->AddAtom(first);    // add to molecule
    107110                                break;
    108111
    109112                        case 'c': // relative coordinates of atom wrt to already placed atom
    110                                 first = new atom;
    111                                 valid = true;
    112                                 do {
    113                                         if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
    114                                         second = mol->AskAtom("Enter atom number: ");
    115                                         cout << Verbose(0) << "Enter relative coordinates." << endl;
    116                                         first->x.AskPosition(mol->cell_size, false);
    117                                         for (int i=NDIM;i--;) {
    118                                                 first->x.x[i] += second->x.x[i];
    119                                         }
    120                                 } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    121                                 first->type = periode->AskElement();    // give type
    122                                 mol->AddAtom(first);    // add to molecule
     113        first = new atom;
     114        valid = true;
     115        do {
     116          if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
     117          second = mol->AskAtom("Enter atom number: ");
     118          cout << Verbose(0) << "Enter relative coordinates." << endl;
     119          first->x.AskPosition(mol->cell_size, false);
     120          for (int i=NDIM;i--;) {
     121            first->x.x[i] += second->x.x[i];
     122          }
     123        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     124        first->type = periode->AskElement();    // give type
     125        mol->AddAtom(first);    // add to molecule
     126        break;
     127
     128    case 'd': // two atoms, two angles and a distance
     129        first = new atom;
     130        valid = true;
     131        do {
     132          if (!valid) {
     133            cout << Verbose(0) << "Resulting coordinates out of cell - ";
     134            first->x.Output((ofstream *)&cout);
     135            cout << Verbose(0) << endl;
     136          }
     137          cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
     138          second = mol->AskAtom("Enter central atom: ");
     139          third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
     140          fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
     141          a = ask_value("Enter distance between central (first) and new atom: ");
     142          b = ask_value("Enter angle between new, first and second atom (degrees): ");
     143          b *= M_PI/180.;
     144          bound(&b, 0., 2.*M_PI);
     145          c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
     146          c *= M_PI/180.;
     147          bound(&c, -M_PI, M_PI);
     148          cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
     149/*
     150          second->Output(1,1,(ofstream *)&cout);
     151          third->Output(1,2,(ofstream *)&cout);
     152          fourth->Output(1,3,(ofstream *)&cout);
     153          n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
     154          x.Copyvector(&second->x);
     155          x.SubtractVector(&third->x);
     156          x.Copyvector(&fourth->x);
     157          x.SubtractVector(&third->x);
     158
     159          if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
     160            cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
     161            continue;
     162          }
     163          cout << Verbose(0) << "resulting relative coordinates: ";
     164          z.Output((ofstream *)&cout);
     165          cout << Verbose(0) << endl;
     166          */
     167          // calc axis vector
     168          x.CopyVector(&second->x);
     169          x.SubtractVector(&third->x);
     170          x.Normalize();
     171          cout << "x: ",
     172          x.Output((ofstream *)&cout);
     173          cout << endl;
     174          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
     175          cout << "z: ",
     176          z.Output((ofstream *)&cout);
     177          cout << endl;
     178          y.MakeNormalVector(&x,&z);
     179          cout << "y: ",
     180          y.Output((ofstream *)&cout);
     181          cout << endl;
     182
     183          // rotate vector around first angle
     184          first->x.CopyVector(&x);
     185          first->x.RotateVector(&z,b - M_PI);
     186          cout << "Rotated vector: ",
     187          first->x.Output((ofstream *)&cout);
     188          cout << endl;
     189          // remove the projection onto the rotation plane of the second angle
     190          n.CopyVector(&y);
     191          n.Scale(first->x.Projection(&y));
     192          cout << "N1: ",
     193          n.Output((ofstream *)&cout);
     194          cout << endl;
     195          first->x.SubtractVector(&n);
     196          cout << "Subtracted vector: ",
     197          first->x.Output((ofstream *)&cout);
     198          cout << endl;
     199          n.CopyVector(&z);
     200          n.Scale(first->x.Projection(&z));
     201          cout << "N2: ",
     202          n.Output((ofstream *)&cout);
     203          cout << endl;
     204          first->x.SubtractVector(&n);
     205          cout << "2nd subtracted vector: ",
     206          first->x.Output((ofstream *)&cout);
     207          cout << endl;
     208
     209          // rotate another vector around second angle
     210          n.CopyVector(&y);
     211          n.RotateVector(&x,c - M_PI);
     212          cout << "2nd Rotated vector: ",
     213          n.Output((ofstream *)&cout);
     214          cout << endl;
     215
     216          // add the two linear independent vectors
     217          first->x.AddVector(&n);
     218          first->x.Normalize();
     219          first->x.Scale(a);
     220          first->x.AddVector(&second->x);
     221
     222          cout << Verbose(0) << "resulting coordinates: ";
     223          first->x.Output((ofstream *)&cout);
     224          cout << Verbose(0) << endl;
     225        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     226        first->type = periode->AskElement();    // give type
     227        mol->AddAtom(first);    // add to molecule
    123228                                break;
    124229
    125                         case 'd': // two atoms, two angles and a distance
    126                                 first = new atom;
    127                                 valid = true;
    128                                 do {
    129                                         if (!valid) {
    130                                                 cout << Verbose(0) << "Resulting coordinates out of cell - ";
    131                                                 first->x.Output((ofstream *)&cout);
    132                                                 cout << Verbose(0) << endl;
    133                                         }
    134                                         cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
    135                                         second = mol->AskAtom("Enter central atom: ");
    136                                         third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
    137                                         fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
    138                                         a = ask_value("Enter distance between central (first) and new atom: ");
    139                                         b = ask_value("Enter angle between new, first and second atom (degrees): ");
    140                                         b *= M_PI/180.;
    141                                         bound(&b, 0., 2.*M_PI);
    142                                         c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
    143                                         c *= M_PI/180.;
    144                                         bound(&c, -M_PI, M_PI);
    145                                         cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
    146 /*
    147                                         second->Output(1,1,(ofstream *)&cout);
    148                                         third->Output(1,2,(ofstream *)&cout);
    149                                         fourth->Output(1,3,(ofstream *)&cout);
    150                                         n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
    151                                         x.Copyvector(&second->x);
    152                                         x.SubtractVector(&third->x);
    153                                         x.Copyvector(&fourth->x);
    154                                         x.SubtractVector(&third->x);
    155 
    156                                         if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
    157                                                 cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
    158                                                 continue;
    159                                         }
    160                                         cout << Verbose(0) << "resulting relative coordinates: ";
    161                                         z.Output((ofstream *)&cout);
    162                                         cout << Verbose(0) << endl;
    163                                         */
    164                                         // calc axis vector
    165                                         x.CopyVector(&second->x);
    166                                         x.SubtractVector(&third->x);
    167                                         x.Normalize();
    168                                         cout << "x: ",
    169                                         x.Output((ofstream *)&cout);
    170                                         cout << endl;
    171                                         z.MakeNormalVector(&second->x,&third->x,&fourth->x);
    172                                         cout << "z: ",
    173                                         z.Output((ofstream *)&cout);
    174                                         cout << endl;
    175                                         y.MakeNormalVector(&x,&z);
    176                                         cout << "y: ",
    177                                         y.Output((ofstream *)&cout);
    178                                         cout << endl;
    179 
    180                                         // rotate vector around first angle
    181                                         first->x.CopyVector(&x);
    182                                         first->x.RotateVector(&z,b - M_PI);
    183                                         cout << "Rotated vector: ",
    184                                         first->x.Output((ofstream *)&cout);
    185                                         cout << endl;
    186                                         // remove the projection onto the rotation plane of the second angle
    187                                         n.CopyVector(&y);
    188                                         n.Scale(first->x.Projection(&y));
    189                                         cout << "N1: ",
    190                                         n.Output((ofstream *)&cout);
    191                                         cout << endl;
    192                                         first->x.SubtractVector(&n);
    193                                         cout << "Subtracted vector: ",
    194                                         first->x.Output((ofstream *)&cout);
    195                                         cout << endl;
    196                                         n.CopyVector(&z);
    197                                         n.Scale(first->x.Projection(&z));
    198                                         cout << "N2: ",
    199                                         n.Output((ofstream *)&cout);
    200                                         cout << endl;
    201                                         first->x.SubtractVector(&n);
    202                                         cout << "2nd subtracted vector: ",
    203                                         first->x.Output((ofstream *)&cout);
    204                                         cout << endl;
    205 
    206                                         // rotate another vector around second angle
    207                                         n.CopyVector(&y);
    208                                         n.RotateVector(&x,c - M_PI);
    209                                         cout << "2nd Rotated vector: ",
    210                                         n.Output((ofstream *)&cout);
    211                                         cout << endl;
    212 
    213                                         // add the two linear independent vectors
    214                                         first->x.AddVector(&n);
    215                                         first->x.Normalize();
    216                                         first->x.Scale(a);
    217                                         first->x.AddVector(&second->x);
    218 
    219                                         cout << Verbose(0) << "resulting coordinates: ";
    220                                         first->x.Output((ofstream *)&cout);
    221                                         cout << Verbose(0) << endl;
    222                                 } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
    223                                 first->type = periode->AskElement();    // give type
    224                                 mol->AddAtom(first);    // add to molecule
    225                                 break;
    226 
    227230                        case 'e': // least square distance position to a set of atoms
    228                                 first = new atom;
    229                                 atoms = new (Vector*[128]);
    230                                 valid = true;
    231                                 for(int i=128;i--;)
    232                                         atoms[i] = NULL;
    233                                 int i=0, j=0;
    234                                 cout << Verbose(0) << "Now we need at least three molecules.\n";
    235                                 do {
    236                                         cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
    237                                         cin >> j;
    238                                         if (j != -1) {
    239                                                 second = mol->FindAtom(j);
    240                                                 atoms[i++] = &(second->x);
    241                                         }
    242                                 } while ((j != -1) && (i<128));
    243                                 if (i >= 2) {
    244                                         first->x.LSQdistance(atoms, i);
    245 
    246                                         first->x.Output((ofstream *)&cout);
    247                                         first->type = periode->AskElement();    // give type
    248                                         mol->AddAtom(first);    // add to molecule
    249                                 } else {
    250                                         delete first;
    251                                         cout << Verbose(0) << "Please enter at least two vectors!\n";
    252                                 }
     231        first = new atom;
     232        atoms = new (Vector*[128]);
     233        valid = true;
     234        for(int i=128;i--;)
     235          atoms[i] = NULL;
     236        int i=0, j=0;
     237        cout << Verbose(0) << "Now we need at least three molecules.\n";
     238        do {
     239          cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
     240          cin >> j;
     241          if (j != -1) {
     242            second = mol->FindAtom(j);
     243            atoms[i++] = &(second->x);
     244          }
     245        } while ((j != -1) && (i<128));
     246        if (i >= 2) {
     247          first->x.LSQdistance(atoms, i);
     248
     249          first->x.Output((ofstream *)&cout);
     250          first->type = periode->AskElement();  // give type
     251          mol->AddAtom(first);  // add to molecule
     252        } else {
     253          delete first;
     254          cout << Verbose(0) << "Please enter at least two vectors!\n";
     255        }
    253256                                break;
    254257        };
    … …  
    256259
    257260/** Submenu for centering the atoms in the molecule.
    258  * \param *mol the molecule with all the atoms
     261 * \param *mol molecule with all the atoms
    259262 */
    260263static void CenterAtoms(molecule *mol)
    … …  
    314317/** Submenu for aligning the atoms in the molecule.
    315318 * \param *periode periodentafel
    316  * \param *mol the molecule with all the atoms
     319 * \param *mol molecule with all the atoms
    317320 */
    318321static void AlignAtoms(periodentafel *periode, molecule *mol)
    … …  
    382385
    383386/** Submenu for mirroring the atoms in the molecule.
    384  * \param *mol the molecule with all the atoms
     387 * \param *mol molecule with all the atoms
    385388 */
    386389static void MirrorAtoms(molecule *mol)
    … …  
    429432
    430433/** Submenu for removing the atoms from the molecule.
    431  * \param *mol the molecule with all the atoms
     434 * \param *mol molecule with all the atoms
    432435 */
    433436static void RemoveAtoms(molecule *mol)
    … …  
    487490/** Submenu for measuring out the atoms in the molecule.
    488491 * \param *periode periodentafel
    489  * \param *mol the molecule with all the atoms
     492 * \param *mol molecule with all the atoms
    490493 */
    491494static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
    … …  
    597600
    598601/** Submenu for measuring out the atoms in the molecule.
    599  * \param *mol the molecule with all the atoms
     602 * \param *mol molecule with all the atoms
    600603 * \param *configuration configuration structure for the to be written config files of all fragments
    601604 */
    … …  
    617620};
    618621
     622/********************************************** Submenu routine **************************************/
     623
     624/** Submenu for manipulating atoms.
     625 * \param *periode periodentafel
     626 * \param *molecules list of molecules whose atoms are to be manipulated
     627 */
     628static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
     629{
     630  atom *first, *second, *third, *fourth;
     631  Vector **atoms;
     632  molecule *mol = NULL;
     633  Vector x,y,z,n; // coordinates for absolute point in cell volume
     634  double *factor; // unit factor if desired
     635  double a,b,c;
     636  double bond, min_bond;
     637  char choice;  // menu choice char
     638  bool valid;
     639
     640  cout << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl;
     641  cout << Verbose(0) << "a - add an atom" << endl;
     642  cout << Verbose(0) << "r - remove an atom" << endl;
     643  cout << Verbose(0) << "b - scale a bond between atoms" << endl;
     644  cout << Verbose(0) << "u - change an atoms element" << endl;
     645  cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
     646  cout << Verbose(0) << "all else - go back" << endl;
     647  cout << Verbose(0) << "===============================================" << endl;
     648  if (molecules->NumberOfActiveMolecules() > 0)
     649    cout << Verbose(0) << "WARNING: There is more than one molecule active! Atoms will be added to each." << endl;
     650  cout << Verbose(0) << "INPUT: ";
     651  cin >> choice;
     652
     653  switch (choice) {
     654    default:
     655      cout << Verbose(0) << "Not a valid choice." << endl;
     656      break;
     657
     658    case 'a': // add atom
     659      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     660        mol = *ListRunner;
     661        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     662        AddAtoms(periode, mol);
     663      }
     664      break;
     665
     666    case 'b': // scale a bond
     667      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     668        mol = *ListRunner;
     669        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     670        cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
     671        first = mol->AskAtom("Enter first (fixed) atom: ");
     672        second = mol->AskAtom("Enter second (shifting) atom: ");
     673        min_bond = 0.;
     674        for (int i=NDIM;i--;)
     675          min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
     676        min_bond = sqrt(min_bond);
     677        cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
     678        cout << Verbose(0) << "Enter new bond length [a.u.]: ";
     679        cin >> bond;
     680        for (int i=NDIM;i--;) {
     681          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
     682        }
     683        //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
     684        //second->Output(second->type->No, 1, (ofstream *)&cout);
     685      }
     686      break;
     687
     688    case 'c': // unit scaling of the metric
     689      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     690        mol = *ListRunner;
     691        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     692       cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
     693       cout << Verbose(0) << "Enter three factors: ";
     694       factor = new double[NDIM];
     695       cin >> factor[0];
     696       cin >> factor[1];
     697       cin >> factor[2];
     698       valid = true;
     699       mol->Scale(&factor);
     700       delete[](factor);
     701      }
     702     break;
     703
     704    case 'l': // measure distances or angles
     705      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     706        mol = *ListRunner;
     707        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     708        MeasureAtoms(periode, mol, configuration);
     709      }
     710      break;
     711
     712    case 'r': // remove atom
     713      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     714        mol = *ListRunner;
     715        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     716        RemoveAtoms(mol);
     717      }
     718      break;
     719
     720    case 'u': // change an atom's element
     721      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     722        int Z;
     723        mol = *ListRunner;
     724        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     725        first = NULL;
     726        do {
     727          cout << Verbose(0) << "Change the element of which atom: ";
     728          cin >> Z;
     729        } while ((first = mol->FindAtom(Z)) == NULL);
     730        cout << Verbose(0) << "New element by atomic number Z: ";
     731        cin >> Z;
     732        first->type = periode->FindElement(Z);
     733        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
     734      }
     735      break;
     736  }
     737};
     738
     739/** Submenu for manipulating molecules.
     740 * \param *periode periodentafel
     741 * \param *molecules list of molecule to manipulate
     742 */
     743static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
     744{
     745  atom *first, *second, *third, *fourth;
     746  Vector **atoms;
     747  Vector x,y,z,n; // coordinates for absolute point in cell volume
     748  double a,b,c;
     749  int j, axis, count, faktor;
     750  char choice;  // menu choice char
     751  bool valid;
     752  molecule *mol = NULL;
     753  element **Elements;
     754  Vector **vectors;
     755  MoleculeLeafClass *Subgraphs = NULL;
     756
     757  cout << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;
     758  cout << Verbose(0) << "c - scale by unit transformation" << endl;
     759  cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
     760  cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
     761  cout << Verbose(0) << "g - center atoms in box" << endl;
     762  cout << Verbose(0) << "i - realign molecule" << endl;
     763  cout << Verbose(0) << "m - mirror all molecules" << endl;
     764  cout << Verbose(0) << "o - create connection matrix" << endl;
     765  cout << Verbose(0) << "t - translate molecule by vector" << endl;
     766  cout << Verbose(0) << "all else - go back" << endl;
     767  cout << Verbose(0) << "===============================================" << endl;
     768  if (molecules->NumberOfActiveMolecules() > 0)
     769    cout << Verbose(0) << "WARNING: There is more than one molecule active! Atoms will be added to each." << endl;
     770  cout << Verbose(0) << "INPUT: ";
     771  cin >> choice;
     772
     773  switch (choice) {
     774    default:
     775      cout << Verbose(0) << "Not a valid choice." << endl;
     776      break;
     777
     778    case 'd': // duplicate the periodic cell along a given axis, given times
     779      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     780        mol = *ListRunner;
     781        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     782        cout << Verbose(0) << "State the axis [(+-)123]: ";
     783        cin >> axis;
     784        cout << Verbose(0) << "State the factor: ";
     785        cin >> faktor;
     786
     787        mol->CountAtoms((ofstream *)&cout); // recount atoms
     788        if (mol->AtomCount != 0) {  // if there is more than none
     789          count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
     790          Elements = new element *[count];
     791          vectors = new Vector *[count];
     792          j = 0;
     793          first = mol->start;
     794          while (first->next != mol->end) { // make a list of all atoms with coordinates and element
     795            first = first->next;
     796            Elements[j] = first->type;
     797            vectors[j] = &first->x;
     798            j++;
     799          }
     800          if (count != j)
     801            cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
     802          x.Zero();
     803          y.Zero();
     804          y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     805          for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
     806            x.AddVector(&y); // per factor one cell width further
     807            for (int k=count;k--;) { // go through every atom of the original cell
     808              first = new atom(); // create a new body
     809              first->x.CopyVector(vectors[k]);  // use coordinate of original atom
     810              first->x.AddVector(&x);     // translate the coordinates
     811              first->type = Elements[k];  // insert original element
     812              mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
     813            }
     814          }
     815          if (mol->first->next != mol->last) // if connect matrix is present already, redo it
     816            mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration->GetIsAngstroem());
     817          // free memory
     818          delete[](Elements);
     819          delete[](vectors);
     820          // correct cell size
     821          if (axis < 0) { // if sign was negative, we have to translate everything
     822            x.Zero();
     823            x.AddVector(&y);
     824            x.Scale(-(faktor-1));
     825            mol->Translate(&x);
     826          }
     827          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     828        }
     829      }
     830      break;
     831
     832    case 'f':
     833      FragmentAtoms(mol, configuration);
     834      break;
     835
     836    case 'g': // center the atoms
     837      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     838        mol = *ListRunner;
     839        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     840        CenterAtoms(mol);
     841      }
     842      break;
     843
     844    case 'i': // align all atoms
     845      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     846        mol = *ListRunner;
     847        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     848        AlignAtoms(periode, mol);
     849      }
     850      break;
     851
     852    case 'm': // mirror atoms along a given axis
     853      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     854        mol = *ListRunner;
     855        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     856        MirrorAtoms(mol);
     857      }
     858      break;
     859
     860    case 'o': // create the connection matrix
     861      {
     862        double bonddistance;
     863        clock_t start,end;
     864        cout << Verbose(0) << "What's the maximum bond distance: ";
     865        cin >> bonddistance;
     866        start = clock();
     867        mol->CreateAdjacencyList((ofstream *)&cout, bonddistance, configuration->GetIsAngstroem());
     868        mol->CreateListOfBondsPerAtom((ofstream *)&cout);
     869        end = clock();
     870        cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     871      }
     872      break;
     873
     874    case 't': // translate all atoms
     875     for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
     876        mol = *ListRunner;
     877        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     878        cout << Verbose(0) << "Enter translation vector." << endl;
     879        x.AskPosition(mol->cell_size,0);
     880        mol->Translate((const Vector *)&x);
     881     }
     882     break;
     883  }
     884  // Free all
     885  if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
     886    while (Subgraphs->next != NULL) {
     887      Subgraphs = Subgraphs->next;
     888      delete(Subgraphs->previous);
     889    }
     890    delete(Subgraphs);
     891  }
     892};
     893
     894
     895/** Submenu for creating new molecules.
     896 * \param *periode periodentafel
     897 * \param *molecules list of molecules to add to
     898 */
     899static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
     900{
     901  char choice;  // menu choice char
     902  bool valid;
     903  Vector Center;
     904  int nr, count;
     905  molecule *mol = NULL;
     906  char *molname = NULL;
     907  int length;
     908  char filename[MAXSTRINGSIZE];
     909
     910  cout << Verbose(0) << "==========Edit MOLECULES=====================" << endl;
     911  cout << Verbose(0) << "c - create new molecule" << endl;
     912  cout << Verbose(0) << "l - load molecule from xyz file" << endl;
     913  cout << Verbose(0) << "n - change molecule's name" << endl;
     914  cout << Verbose(0) << "N - give molecules filename" << endl;
     915  cout << Verbose(0) << "p - parse xyz file into molecule" << endl;
     916  cout << Verbose(0) << "r - remove a molecule" << endl;
     917  cout << Verbose(0) << "all else - go back" << endl;
     918  cout << Verbose(0) << "===============================================" << endl;
     919  cout << Verbose(0) << "INPUT: ";
     920  cin >> choice;
     921
     922  switch (choice) {
     923    default:
     924      cout << Verbose(0) << "Not a valid choice." << endl;
     925      break;
     926    case 'c':
     927      mol = new molecule(periode);
     928      molecules->insert(mol);
     929      break;
     930
     931    case 'l': // laod from XYZ file
     932      cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     933      mol = new molecule(periode);
     934      do {
     935        cout << Verbose(0) << "Enter file name: ";
     936        cin >> filename;
     937      } while (!mol->AddXYZFile(filename));
     938      mol->SetNameFromFilename(filename);
     939      // center at set box dimensions
     940      mol->CenterEdge((ofstream *)&cout, &Center);
     941      mol->cell_size[0] = Center.x[0];
     942      mol->cell_size[1] = 0;
     943      mol->cell_size[2] = Center.x[1];
     944      mol->cell_size[3] = 0;
     945      mol->cell_size[4] = 0;
     946      mol->cell_size[5] = Center.x[2];
     947      molecules->insert(mol);
     948      break;
     949
     950    case 'n':
     951      do {
     952        cout << Verbose(0) << "Enter index of molecule: ";
     953        cin >> nr;
     954        mol = molecules->ReturnIndex(nr);
     955      } while (mol != NULL);
     956      cout << Verbose(0) << "Enter name: ";
     957      cin >> filename;
     958      strcpy(mol->name, filename);
     959      break;
     960
     961    case 'N':
     962      do {
     963        cout << Verbose(0) << "Enter index of molecule: ";
     964        cin >> nr;
     965        mol = molecules->ReturnIndex(nr);
     966      } while (mol != NULL);
     967      cout << Verbose(0) << "Enter name: ";
     968      cin >> filename;
     969      mol->SetNameFromFilename(filename);
     970      break;
     971
     972    case 'p': // parse XYZ file
     973      mol = NULL;
     974      do {
     975        cout << Verbose(0) << "Enter index of molecule: ";
     976        cin >> nr;
     977        mol = molecules->ReturnIndex(nr);
     978      } while (mol == NULL);
     979      cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     980      do {
     981        cout << Verbose(0) << "Enter file name: ";
     982        cin >> filename;
     983      } while (!mol->AddXYZFile(filename));
     984      mol->SetNameFromFilename(filename);
     985      break;
     986
     987    case 'r':
     988      cout << Verbose(0) << "Enter index of molecule: ";
     989      cin >> nr;
     990      count = 1;
     991      MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin();
     992      for(; ((ListRunner != molecules->ListOfMolecules.end()) && (count < nr)); ListRunner++);
     993      mol = *ListRunner;
     994      if (count == nr) {
     995        molecules->ListOfMolecules.erase(ListRunner);
     996        delete(mol);
     997      }
     998      break;
     999  }
     1000};
     1001
     1002
     1003/** Submenu for merging molecules.
     1004 * \param *periode periodentafel
     1005 * \param *molecules list of molecules to add to
     1006 */
     1007static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
     1008{
     1009  char choice;  // menu choice char
     1010  bool valid;
     1011
     1012  cout << Verbose(0) << "===========MERGE MOLECULES=====================" << endl;
     1013  cout << Verbose(0) << "e - embedding merge of two molecules" << endl;
     1014  cout << Verbose(0) << "m - multi-merge of all molecules" << endl;
     1015  cout << Verbose(0) << "s - scatter merge of two molecules" << endl;
     1016  cout << Verbose(0) << "t - simple merge of two molecules" << endl;
     1017  cout << Verbose(0) << "all else - go back" << endl;
     1018  cout << Verbose(0) << "===============================================" << endl;
     1019  cout << Verbose(0) << "INPUT: ";
     1020  cin >> choice;
     1021
     1022  switch (choice) {
     1023    default:
     1024      cout << Verbose(0) << "Not a valid choice." << endl;
     1025      break;
     1026
     1027    case 'e':
     1028      break;
     1029
     1030    case 'm':
     1031      break;
     1032
     1033    case 's':
     1034      break;
     1035
     1036    case 't':
     1037      break;
     1038  }
     1039};
     1040
     1041
    6191042/********************************************** Test routine **************************************/
    6201043
    6211044/** Is called always as option 'T' in the menu.
     1045 * \param *molecules list of molecules
    6221046 */
    623 static void testroutine(molecule *mol)
     1047static void testroutine(MoleculeListClass *molecules)
    6241048{
    6251049        // the current test routine checks the functionality of the KeySet&Graph concept:
    6261050        // We want to have a multiindex (the KeySet) describing a unique subgraph
    627         atom *Walker = mol->start;
    628         int i, comp, counter=0;
     1051  int i, comp, counter=0;
     1052
     1053  // create a clone
     1054  molecule *mol = NULL;
     1055  if (molecules->ListOfMolecules.size() != 0) // clone
     1056    mol = (molecules->ListOfMolecules.front())->CopyMolecule();
     1057  else {
     1058    cerr << "I don't have anything to test on ... ";
     1059    return;
     1060  }
     1061  atom *Walker = mol->start;
    6291062
    6301063        // generate some KeySets
    … …  
    6861119                A++;
    6871120        }
     1121        delete(mol);
    6881122};
    6891123
    … …  
    6921126 * \param *configuration pointer to configuration structure with all the values
    6931127 * \param *periode pointer to periodentafel structure with all the elements
    694  * \param *mol pointer to molecule structure with all the atoms and coordinates
     1128 * \param *molecules list of molecules structure with all the atoms and coordinates
    6951129 */
    696 static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, molecule *mol)
     1130static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
    6971131{
    6981132        char filename[MAXSTRINGSIZE];
    6991133        ofstream output;
     1134        molecule *mol = new molecule(periode);
     1135
     1136        // merge all molecules in MoleculeListClass into this molecule
     1137        int N = molecules->ListOfMolecules.size();
     1138        int *src = new int(N);
     1139        N=0;
     1140        for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1141          src[N++] = (*ListRunner)->IndexNr;
     1142        molecules->SimpleMultiAdd(mol, src, N);
    7001143
    7011144        cout << Verbose(0) << "Storing configuration ... " << endl;
    … …  
    7611204                cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
    7621205        }
     1206        delete(mol);
    7631207};
    7641208
    … …  
    7661210 * \param argc argument count
    7671211 * \param **argv arguments array
    768  * \param *mol molecule structure
     1212 * \param *molecules list of molecules structure
    7691213 * \param *periode elements structure
    7701214 * \param configuration config file structure
    … …  
    7731217 * \return exit code (0 - successful, all else - something's wrong)
    7741218 */
    775 static int ParseCommandLineOptions(int argc, char **argv, molecule *&mol, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
     1219static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
    7761220{
    7771221        Vector x,y,z,n; // coordinates for absolute point in cell volume
    … …  
    7891233        int argptr;
    7901234        PathToDatabases = LocalPath;
     1235
     1236        // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
     1237        molecule *mol = new molecule(periode);
     1238        molecules->insert(mol);
    7911239
    7921240        if (argc > 1) { // config file specified as option
    … …  
    13271775                } while (argptr < argc);
    13281776                if (SaveFlag)
    1329                         SaveConfig(ConfigFileName, &configuration, periode, mol);
     1777                        SaveConfig(ConfigFileName, &configuration, periode, molecules);
    13301778                if ((ExitFlag >= 1)) {
    13311779                        delete(mol);
    … …  
    13481796{
    13491797        periodentafel *periode = new periodentafel; // and a period table of all elements
    1350         molecule *mol = new molecule(periode);          // first we need an empty molecule
     1798        MoleculeListClass *molecules = new MoleculeListClass;  // list of all molecules
     1799        molecule *mol = NULL;
    13511800        config configuration;
    13521801        double tmp1;
    1353         double bond, min_bond;
    13541802        atom *first, *second;
    13551803        char choice;    // menu choice char
    13561804        Vector x,y,z,n; // coordinates for absolute point in cell volume
    1357         double *factor; // unit factor if desired
    13581805        bool valid; // flag if input was valid or not
    13591806        ifstream test;
    13601807        ofstream output;
    13611808        string line;
    1362         char filename[MAXSTRINGSIZE];
    13631809        char *ConfigFileName = NULL;
    13641810        char *ElementsFileName = NULL;
    13651811        int Z;
    13661812        int j, axis, count, faktor;
    1367         clock_t start,end;
    1368 //      int *MinimumRingSize = NULL;
    1369         MoleculeLeafClass *Subgraphs = NULL;
    1370 //      class StackClass<bond *> *BackEdgeStack = NULL;
    1371         element **Elements;
    1372         Vector **vectors;
    13731813
    13741814        // =========================== PARSE COMMAND LINE OPTIONS ====================================
    1375         j = ParseCommandLineOptions(argc, argv, mol, periode, configuration, ConfigFileName, ElementsFileName);
     1815        j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName, ElementsFileName);
    13761816        if (j == 1) return 0; // just for -v and -h options
    13771817        if (j) return j;        // something went wrong
    13781818
    13791819        // General stuff
    1380         if (mol->cell_size[0] == 0.) {
    1381                 cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
    1382                 for (int i=0;i<6;i++) {
    1383                         cout << Verbose(1) << "Cell size" << i << ": ";
    1384                         cin >> mol->cell_size[i];
    1385                 }
     1820        if (molecules->ListOfMolecules.size() == 0) {
     1821    mol = new molecule(periode);
     1822    if (mol->cell_size[0] == 0.) {
     1823      cout << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
     1824      for (int i=0;i<6;i++) {
     1825        cout << Verbose(1) << "Cell size" << i << ": ";
     1826        cin >> mol->cell_size[i];
     1827      }
     1828    }
     1829    molecules->insert(mol);
    13861830        }
    13871831
    … …  
    13921836        do {
    13931837                cout << Verbose(0) << endl << endl;
    1394                 cout << Verbose(0) << "============Element list=======================" << endl;
    1395                 mol->Checkout((ofstream *)&cout);
    1396                 cout << Verbose(0) << "============Atom list==========================" << endl;
    1397                 mol->Output((ofstream *)&cout);
     1838                cout << Verbose(0) << "============Molecule list=======================" << endl;
     1839                molecules->Enumerate((ofstream *)&cout);
    13981840                cout << Verbose(0) << "============Menu===============================" << endl;
    1399                 cout << Verbose(0) << "a - add an atom" << endl;
    1400                 cout << Verbose(0) << "r - remove an atom" << endl;
    1401                 cout << Verbose(0) << "b - scale a bond between atoms" << endl;
    1402                 cout << Verbose(0) << "u - change an atoms element" << endl;
    1403                 cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
    1404                 cout << Verbose(0) << "-----------------------------------------------" << endl;
    1405                 cout << Verbose(0) << "p - Parse xyz file" << endl;
    1406                 cout << Verbose(0) << "e - edit the current configuration" << endl;
    1407                 cout << Verbose(0) << "o - create connection matrix" << endl;
    1408                 cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
    1409                 cout << Verbose(0) << "-----------------------------------------------" << endl;
    1410                 cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
    1411                 cout << Verbose(0) << "i - realign molecule" << endl;
    1412                 cout << Verbose(0) << "m - mirror all molecules" << endl;
    1413                 cout << Verbose(0) << "t - translate molecule by vector" << endl;
    1414                 cout << Verbose(0) << "c - scale by unit transformation" << endl;
    1415                 cout << Verbose(0) << "g - center atoms in box" << endl;
    1416                 cout << Verbose(0) << "-----------------------------------------------" << endl;
    1417                 cout << Verbose(0) << "s - save current setup to config file" << endl;
    1418                 cout << Verbose(0) << "T - call the current test routine" << endl;
    1419                 cout << Verbose(0) << "q - quit" << endl;
    1420                 cout << Verbose(0) << "===============================================" << endl;
    1421                 cout << Verbose(0) << "Input: ";
    1422                 cin >> choice;
     1841    cout << Verbose(0) << "a - set molecule (in)active" << endl;
     1842    cout << Verbose(0) << "e - edit new molecules" << endl;
     1843    cout << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;
     1844    cout << Verbose(0) << "M - Merge molecules" << endl;
     1845    cout << Verbose(0) << "m - manipulate atoms" << endl;
     1846    cout << Verbose(0) << "-----------------------------------------------" << endl;
     1847    cout << Verbose(0) << "c - edit the current configuration" << endl;
     1848    cout << Verbose(0) << "-----------------------------------------------" << endl;
     1849    cout << Verbose(0) << "s - save current setup to config file" << endl;
     1850    cout << Verbose(0) << "T - call the current test routine" << endl;
     1851    cout << Verbose(0) << "q - quit" << endl;
     1852    cout << Verbose(0) << "===============================================" << endl;
     1853    cout << Verbose(0) << "Input: ";
     1854    cin >> choice;
    14231855
    14241856                switch (choice) {
    1425                         default:
    1426                         case 'a': // add atom
    1427                                 AddAtoms(periode, mol);
    1428                                 choice = 'a';
    1429                                 break;
    1430 
    1431                         case 'b': // scale a bond
    1432                                 cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
    1433                                 first = mol->AskAtom("Enter first (fixed) atom: ");
    1434                                 second = mol->AskAtom("Enter second (shifting) atom: ");
    1435                                 min_bond = 0.;
    1436                                 for (int i=NDIM;i--;)
    1437                                         min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
    1438                                 min_bond = sqrt(min_bond);
    1439                                 cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
    1440                                 cout << Verbose(0) << "Enter new bond length [a.u.]: ";
    1441                                 cin >> bond;
    1442                                 for (int i=NDIM;i--;) {
    1443                                         second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
    1444                                 }
    1445                                 //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
    1446                                 //second->Output(second->type->No, 1, (ofstream *)&cout);
    1447                                 break;
    1448 
    1449                         case 'c': // unit scaling of the metric
    1450                          cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
    1451                          cout << Verbose(0) << "Enter three factors: ";
    1452                          factor = new double[NDIM];
    1453                          cin >> factor[0];
    1454                          cin >> factor[1];
    1455                          cin >> factor[2];
    1456                          valid = true;
    1457                          mol->Scale(&factor);
    1458                          delete[](factor);
     1857                        case 'a':  // (in)activate molecule
     1858        {
     1859          cout << "Enter index of molecule: ";
     1860          cin >> j;
     1861          count = 1;
     1862          MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin();
     1863          for(; ((ListRunner != molecules->ListOfMolecules.end()) && (count < j)); ListRunner++);
     1864          if (count == j)
     1865            (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
     1866        }
     1867                          break;
     1868
     1869                        case 'c': // edit each field of the configuration
     1870                         configuration.Edit();
    14591871                         break;
    14601872
    1461                         case 'd': // duplicate the periodic cell along a given axis, given times
    1462                                 cout << Verbose(0) << "State the axis [(+-)123]: ";
    1463                                 cin >> axis;
    1464                                 cout << Verbose(0) << "State the factor: ";
    1465                                 cin >> faktor;
    1466 
    1467                                 mol->CountAtoms((ofstream *)&cout);     // recount atoms
    1468                                 if (mol->AtomCount != 0) {      // if there is more than none
    1469                                         count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
    1470                                         Elements = new element *[count];
    1471                                         vectors = new Vector *[count];
    1472                                         j = 0;
    1473                                         first = mol->start;
    1474                                         while (first->next != mol->end) {       // make a list of all atoms with coordinates and element
    1475                                                 first = first->next;
    1476                                                 Elements[j] = first->type;
    1477                                                 vectors[j] = &first->x;
    1478                                                 j++;
    1479                                         }
    1480                                         if (count != j)
    1481                                                 cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
    1482                                         x.Zero();
    1483                                         y.Zero();
    1484                                         y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
    1485                                         for (int i=1;i<faktor;i++) {    // then add this list with respective translation factor times
    1486                                                 x.AddVector(&y); // per factor one cell width further
    1487                                                 for (int k=count;k--;) { // go through every atom of the original cell
    1488                                                         first = new atom(); // create a new body
    1489                                                         first->x.CopyVector(vectors[k]);        // use coordinate of original atom
    1490                                                         first->x.AddVector(&x);                 // translate the coordinates
    1491                                                         first->type = Elements[k];      // insert original element
    1492                                                         mol->AddAtom(first);                            // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
    1493                                                 }
    1494                                         }
    1495                                         if (mol->first->next != mol->last) // if connect matrix is present already, redo it
    1496                                                 mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration.GetIsAngstroem());
    1497                                         // free memory
    1498                                         delete[](Elements);
    1499                                         delete[](vectors);
    1500                                         // correct cell size
    1501                                         if (axis < 0) { // if sign was negative, we have to translate everything
    1502                                                 x.Zero();
    1503                                                 x.AddVector(&y);
    1504                                                 x.Scale(-(faktor-1));
    1505                                                 mol->Translate(&x);
    1506                                         }
    1507                                         mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    1508                                 }
    1509                                 break;
    1510 
    1511                         case 'e': // edit each field of the configuration
    1512                          configuration.Edit(mol);
    1513                          break;
    1514 
    1515                         case 'f':
    1516                                 FragmentAtoms(mol, &configuration);
    1517                                 break;
    1518 
    1519                         case 'g': // center the atoms
    1520                                 CenterAtoms(mol);
    1521                                 break;
    1522 
    1523                         case 'i': // align all atoms
    1524                                 AlignAtoms(periode, mol);
    1525                                 break;
    1526 
    1527                         case 'l': // measure distances or angles
    1528                                 MeasureAtoms(periode, mol, &configuration);
    1529                                 break;
    1530 
    1531                         case 'm': // mirror atoms along a given axis
    1532                                 MirrorAtoms(mol);
    1533                                 break;
    1534 
    1535                         case 'o': // create the connection matrix
    1536                                 {
    1537                                         cout << Verbose(0) << "What's the maximum bond distance: ";
    1538                                         cin >> tmp1;
    1539                                         start = clock();
    1540                                         mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
    1541                                         mol->CreateListOfBondsPerAtom((ofstream *)&cout);
    1542 //                                      Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
    1543 //                                      while (Subgraphs->next != NULL) {
    1544 //                                              Subgraphs = Subgraphs->next;
    1545 //                                              Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
    1546 //                                              delete(Subgraphs->previous);
    1547 //                                      }
    1548 //                                      delete(Subgraphs);              // we don't need the list here, so free everything
    1549 //                                      delete[](MinimumRingSize);
    1550 //                                      Subgraphs = NULL;
    1551                                         end = clock();
    1552                                         cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    1553                                 }
    1554                                 break;
    1555 
    1556                         case 'p': // parse and XYZ file
    1557                                 cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
    1558                                 do {
    1559                                         cout << Verbose(0) << "Enter file name: ";
    1560                                         cin >> filename;
    1561                                 } while (!mol->AddXYZFile(filename));
    1562                                 break;
     1873      case 'g': // manipulate molecules
     1874        ManipulateMolecules(periode, molecules, &configuration);
     1875        break;
     1876
     1877      case 'M':  // merge molecules
     1878        MergeMolecules(periode, molecules);
     1879        break;
     1880
     1881      case 'm': // manipulate atoms
     1882        ManipulateAtoms(periode, molecules, &configuration);
     1883        break;
     1884
     1885      case 'e': // create molecule
     1886        EditMolecules(periode, molecules);
     1887        break;
    15631888
    15641889                        case 'q': // quit
    15651890                                break;
    15661891
    1567                         case 'r': // remove atom
    1568                                 RemoveAtoms(mol);
     1892                        case 's': // save to config file
     1893                                SaveConfig(ConfigFileName, &configuration, periode, molecules);
    15691894                                break;
    15701895
    1571                         case 's': // save to config file
    1572                                 SaveConfig(ConfigFileName, &configuration, periode, mol);
     1896                        case 'T':
     1897                                testroutine(molecules);
    15731898                                break;
    15741899
    1575                         case 't': // translate all atoms
    1576                          cout << Verbose(0) << "Enter translation vector." << endl;
    1577                          x.AskPosition(mol->cell_size,0);
    1578                          mol->Translate((const Vector *)&x);
    1579                          break;
    1580 
    1581                         case 'T':
    1582                                 testroutine(mol);
    1583                                 break;
    1584 
    1585                         case 'u': // change an atom's element
    1586                                 first = NULL;
    1587                                 do {
    1588                                         cout << Verbose(0) << "Change the element of which atom: ";
    1589                                         cin >> Z;
    1590                                 } while ((first = mol->FindAtom(Z)) == NULL);
    1591                                 cout << Verbose(0) << "New element by atomic number Z: ";
    1592                                 cin >> Z;
    1593                                 first->type = periode->FindElement(Z);
    1594                                 cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
    1595                                 break;
     1900                        default:
     1901                          break;
    15961902                };
    15971903        } while (choice != 'q');
    … …  
    16031909                cout << Verbose(0) << "Saving of elements.db failed." << endl;
    16041910
    1605         // Free all
    1606         if (Subgraphs != NULL) {        // free disconnected subgraph list of DFS analysis was performed
    1607                 while (Subgraphs->next != NULL) {
    1608                         Subgraphs = Subgraphs->next;
    1609                         delete(Subgraphs->previous);
    1610                 }
    1611                 delete(Subgraphs);
    1612         }
    1613         delete(mol);
     1911        delete(molecules);
    16141912        delete(periode);
    16151913        return (0);

  • src/config.cpp

    rd8b94a r1907a7  
    11/** \file config.cpp
    2  * 
     2 *
    33 * Function implementations for the class config.
    4  * 
     4 *
    55 */
    66
    … …  
    2424        configname[0]='\0';
    2525        basis = "3-21G";
    26        
     26
    2727        FastParsing = false;
    2828        ProcPEGamma=8;
    … …  
    4242        UseAddGramSch=1;
    4343        Seed=1;
    44        
     44
    4545        MaxOuterStep=0;
    4646        Deltat=1;
    … …  
    5151        MaxPsiStep=0;
    5252        EpsWannier=1e-7;
    53        
     53
    5454        MaxMinStep=100;
    5555        RelEpsTotalEnergy=1e-7;
    … …  
    6262        InitMaxMinStopStep=1;
    6363        InitMaxMinGapStopStep=0;
    64        
     64
    6565        //BoxLength[NDIM*NDIM];
    66        
     66
    6767        ECut=128.;
    6868        MaxLevel=5;
    … …  
    7777        PsiMaxNoDown=0;
    7878        AddPsis=0;
    79        
     79
    8080        RCut=20.;
    8181        StructOpt=0;
    … …  
    101101 * for each entry of the config file structure.
    102102 */
    103 void config::Edit(molecule *mol)
     103void config::Edit()
    104104{
    105105        char choice;
    106        
     106
    107107        do {
    108108                cout << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
    … …  
    137137                cout << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
    138138                cout << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
    139                 cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
     139//              cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
    140140                cout << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
    141141                cout << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
    … …  
    157157                cout << Verbose(0) << "INPUT: ";
    158158                cin >> choice;
    159                
     159
    160160                switch (choice) {
    161161                                case 'A': // mainname
    … …  
    282282                                        cin >> config::InitMaxMinStopStep;
    283283                                        break;
    284                                
    285                                 case 'j': // BoxLength
    286                                         cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
    287                                         for (int i=0;i<6;i++) {
    288                                                 cout << Verbose(0) << "Cell size" << i << ": ";
    289                                                 cin >> mol->cell_size[i];
    290                                         }
    291                                         break;
    292                                
     284
     285//                              case 'j': // BoxLength
     286//                                      cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
     287//                                      for (int i=0;i<6;i++) {
     288//                                              cout << Verbose(0) << "Cell size" << i << ": ";
     289//                                              cin >> mol->cell_size[i];
     290//                                      }
     291//                                      break;
     292
    293293                                case 'k': // ECut
    294294                                        cout << Verbose(0) << "Old: " << config::ECut << "\t new: ";
    … …  
    364364 * \param *periode pointer to a periodentafel class with all elements
    365365 * \param *mol pointer to molecule containing all atoms of the molecule
    366  * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax 
     366 * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax
    367367 */
    368368int config::TestSyntax(char *filename, periodentafel *periode, molecule *mol)
    … …  
    370370        int test;
    371371        ifstream file(filename);
    372        
     372
    373373        // search file for keyword: ProcPEGamma (new syntax)
    374374        if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
    … …  
    441441 * \param *file input file stream being the opened config file
    442442 * \param *periode pointer to a periodentafel class with all elements
    443  * \param *mol pointer to molecule containing all atoms of the molecule 
     443 * \param *mol pointer to molecule containing all atoms of the molecule
    444444 */
    445445void config::Load(char *filename, periodentafel *periode, molecule *mol)
    … …  
    452452        RetrieveConfigPathAndName(filename);
    453453        // ParseParameters
    454        
     454
    455455        /* Oeffne Hauptparameterdatei */
    456456        int di;
    … …  
    464464        int verbose = 0;
    465465        double value[3];
    466        
     466
    467467        /* Namen einlesen */
    468468
    … …  
    495495                config::DoOutCurrent = 0;
    496496        if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
    497         if (config::DoOutCurrent > 1) config::DoOutCurrent = 1; 
     497        if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
    498498        ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
    499499        if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
    … …  
    517517                if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
    518518        }
    519        
     519
    520520        ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
    521521        if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
    … …  
    527527        if (!ParseForParameter(verbose,file,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
    528528                config::EpsWannier = 1e-8;
    529        
     529
    530530        // stop conditions
    531531        //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
    532532        ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
    533533        if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
    534        
     534
    535535        ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
    536536        ParseForParameter(verbose,file,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
    … …  
    541541        if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
    542542        if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;
    543        
     543
    544544        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
    545545        ParseForParameter(verbose,file,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
    … …  
    550550        if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
    551551        if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;
    552        
     552
    553553        // Unit cell and magnetic field
    554554        ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    … …  
    566566                config::DoFullCurrent = 0;
    567567        if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
    568         if (config::DoFullCurrent > 2) config::DoFullCurrent = 2; 
     568        if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
    569569        if (config::DoOutNICS < 0) config::DoOutNICS = 0;
    570         if (config::DoOutNICS > 2) config::DoOutNICS = 2; 
     570        if (config::DoOutNICS > 2) config::DoOutNICS = 2;
    571571        if (config::DoPerturbation == 0) {
    572572                config::DoFullCurrent = 0;
    … …  
    585585        } else {
    586586                fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    587                 exit(1); 
     587                exit(1);
    588588        }
    589589        switch (config::RiemannTensor) {
    … …  
    602602                        if (config::RiemannLevel < 2) {
    603603                                config::RiemannLevel = 2;
    604                         } 
     604                        }
    605605                        if (config::RiemannLevel > config::MaxLevel-1) {
    606606                                config::RiemannLevel = config::MaxLevel-1;
    … …  
    609609                        if (config::LevRFactor < 2) {
    610610                                config::LevRFactor = 2;
    611                         } 
     611                        }
    612612                        config::Lev0Factor = 2;
    613613                        config::RTActualUse = 2;
    … …  
    619619        } else {
    620620                fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    621                 exit(1); 
     621                exit(1);
    622622        }
    623623        switch (config::PsiType) {
    … …  
    633633                break;
    634634        }
    635        
     635
    636636        // IonsInitRead
    637        
     637
    638638        ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
    639639        ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
    … …  
    653653                        ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
    654654                        elementhash[i] = periode->FindElement(Z);
    655                         cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl; 
     655                        cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
    656656                }
    657657                int repetition = 0; // which repeated keyword shall be read
    658        
     658
    659659                map<int, atom *> AtomList[config::MaxTypes];
    660660                if (!FastParsing) {
    … …  
    675675                                                } else
    676676                                                        neues = AtomList[i][j];
    677                                                 status = (status && 
     677                                                status = (status &&
    678678                                                                                ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
    679679                                                                                ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
    … …  
    681681                                                                                ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
    682682                                                if (!status) break;
    683                
     683
    684684                                                // check size of vectors
    685685                                                if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
    … …  
    689689                                                        mol->Trajectories[neues].F.resize(repetition+10);
    690690                                                }
    691                                        
     691
    692692                                                // put into trajectories list
    693693                                                for (int d=0;d<NDIM;d++)
    694694                                                        mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
    695                                                
     695
    696696                                                // parse velocities if present
    697697                                                if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
    … …  
    703703                                                for (int d=0;d<NDIM;d++)
    704704                                                        mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
    705                        
     705
    706706                                                // parse forces if present
    707707                                                if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
    … …  
    713713                                                for (int d=0;d<NDIM;d++)
    714714                                                        mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
    715        
    716         //                                              cout << "Parsed position of step " << (repetition) << ": ("; 
     715
     716        //                                              cout << "Parsed position of step " << (repetition) << ": (";
    717717        //                                              for (int d=0;d<NDIM;d++)
    718718        //                                                      cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";                                  // next step
    … …  
    770770 * \param *file input file stream being the opened config file with old pcp syntax
    771771 * \param *periode pointer to a periodentafel class with all elements
    772  * \param *mol pointer to molecule containing all atoms of the molecule 
     772 * \param *mol pointer to molecule containing all atoms of the molecule
    773773 */
    774774void config::LoadOld(char *filename, periodentafel *periode, molecule *mol)
    … …  
    781781        RetrieveConfigPathAndName(filename);
    782782        // ParseParameters
    783        
     783
    784784        /* Oeffne Hauptparameterdatei */
    785785        int l, i, di;
    … …  
    791791        int Z, No, AtomNo, found;
    792792        int verbose = 0;
    793        
     793
    794794        /* Namen einlesen */
    795795
    … …  
    823823        ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
    824824        config::EpsWannier = 1e-8;
    825        
     825
    826826        // stop conditions
    827827        //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
    828828        ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
    829829        if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
    830        
     830
    831831        ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
    832832        ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
    … …  
    836836        if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
    837837        config::MaxMinGapStopStep = 1;
    838        
     838
    839839        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
    840840        ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
    … …  
    867867        } else {
    868868                fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
    869                 exit(1); 
     869                exit(1);
    870870        }
    871871        switch (config::RiemannTensor) {
    … …  
    884884                        if (config::RiemannLevel < 2) {
    885885                                config::RiemannLevel = 2;
    886                         } 
     886                        }
    887887                        if (config::RiemannLevel > config::MaxLevel-1) {
    888888                                config::RiemannLevel = config::MaxLevel-1;
    … …  
    891891                        if (config::LevRFactor < 2) {
    892892                                config::LevRFactor = 2;
    893                         } 
     893                        }
    894894                        config::Lev0Factor = 2;
    895895                        config::RTActualUse = 2;
    … …  
    901901        } else {
    902902                fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
    903                 exit(1); 
     903                exit(1);
    904904        }
    905905        switch (config::PsiType) {
    … …  
    915915                break;
    916916        }
    917        
     917
    918918        // IonsInitRead
    919        
     919
    920920        ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
    921921        ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
    … …  
    924924
    925925        // Routine from builder.cpp
    926        
    927        
    928         for (i=MAX_ELEMENTS;i--;) 
     926
     927
     928        for (i=MAX_ELEMENTS;i--;)
    929929                elementhash[i] = NULL;
    930930        cout << Verbose(0) << "Parsing Ions ..." << endl;
    … …  
    937937                }
    938938                if (found > 0) {
    939                         if (zeile.find("Ions_Data") == 0) 
     939                        if (zeile.find("Ions_Data") == 0)
    940940                                getline(*file,zeile,'\n'); // read next line and parse this one
    941941                        istringstream input(zeile);
    … …  
    966966                        No++;
    967967                }
    968         }       
     968        }
    969969        file->close();
    970970        delete(file);
    … …  
    974974 * \param *filename name of file
    975975 * \param *periode pointer to a periodentafel class with all elements
    976  * \param *mol pointer to molecule containing all atoms of the molecule 
     976 * \param *mol pointer to molecule containing all atoms of the molecule
    977977 */
    978978bool config::Save(const char *filename, periodentafel *periode, molecule *mol) const
    … …  
    10811081 * Note that this format cannot be parsed again.
    10821082 * \param *filename name of file (without ".in" suffix!)
    1083  * \param *mol pointer to molecule containing all atoms of the molecule 
     1083 * \param *mol pointer to molecule containing all atoms of the molecule
    10841084 */
    10851085bool config::SaveMPQC(const char *filename, molecule *mol) const
    … …  
    10921092        ofstream *output = NULL;
    10931093        stringstream *fname = NULL;
    1094        
     1094
    10951095        // first without hessian
    10961096        fname = new stringstream;
    … …  
    11911191        delete(output);
    11921192        delete(fname);
    1193        
     1193
    11941194        return true;
    11951195};
    … …  
    12031203 * \param file file to be parsed
    12041204 * \param name Name of value in file (at least 3 chars!)
    1205  * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning 
     1205 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
    12061206 *                              (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
    12071207 *                               best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
    … …  
    12231223        dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char), "config::ParseForParameter: *free_dummy");
    12241224
    1225         //fprintf(stderr,"Parsing for %s\n",name);     
     1225        //fprintf(stderr,"Parsing for %s\n",name);
    12261226        if (repetition == 0)
    12271227                //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
    … …  
    12441244                                        file->clear();
    12451245                                        file->seekg(file_position, ios::beg);   // rewind to start position
    1246                                         Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");                                   
     1246                                        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    12471247                                        return 0;
    12481248                                }
    … …  
    12501250                        line++;
    12511251                } while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
    1252                
     1252
    12531253                // C++ getline removes newline at end, thus re-add
    12541254                if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
    … …  
    12761276                                //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
    12771277                                //Free((void **)&free_dummy);
    1278                                 //Error(FileOpenParams, NULL);                 
     1278                                //Error(FileOpenParams, NULL);
    12791279                        } else {
    12801280                                //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
    … …  
    12851285                if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
    12861286                        found++; // found the parameter!
    1287                         //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);         
    1288                        
     1287                        //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);
     1288
    12891289                        if (found == repetition) {
    12901290                                for (i=0;i<xth;i++) { // i = rows
    … …  
    13251325                                                dummy1[j+1] = '\0';
    13261326                                        }
    1327        
     1327
    13281328                                        int start = (type >= grid) ? 0 : yth-1 ;
    13291329                                        for (j=start;j<yth;j++) { // j = columns
    … …  
    13481348                                                                                //return 0;
    13491349                                                                                exit(255);
    1350                                                                                 //Error(FileOpenParams, NULL);                 
     1350                                                                                //Error(FileOpenParams, NULL);
    13511351                                                                        } else {
    13521352                                                                                //if (!sequential)
    … …  
    13631363                                                                // found comment, skipping rest of line
    13641364                                                                //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
    1365                                                                 if (!sequential) { // here we need it! 
     1365                                                                if (!sequential) { // here we need it!
    13661366                                                                        file->seekg(file_position, ios::beg);   // rewind to start position
    13671367                                                                }
    … …  
    13781378                                                                                //value += sizeof(int);
    13791379                                                                        break;
    1380                                                                 case(row_double): 
     1380                                                                case(row_double):
    13811381                                                                case(grid):
    13821382                                                                case(lower_trigrid):
    … …  
    14191419                        }
    14201420                }
    1421         } 
     1421        }
    14221422        if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
    14231423        Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
    … …  
    14271427        }
    14281428        //fprintf(stderr, "End of Parsing\n\n");
    1429        
     1429
    14301430        return (found); // true if found, false if not
    14311431}

  • src/helpers.cpp

    rd8b94a r1907a7  
    170170        }
    171171        // allocate string
    172         returnstring = (char *) Malloc(sizeof(char)*(order+2), "MoleculeListClass::CreateFragmentNumberForOutput: *returnstring");
     172        returnstring = (char *) Malloc(sizeof(char)*(order+2), "FixedDigitNumber: *returnstring");
    173173        // terminate    and fill string array from end backward
    174174        returnstring[order] = '\0';

  • src/moleculelist.cpp

    rd8b94a r1907a7  
    11/** \file MoleculeListClass.cpp
    2  * 
     2 *
    33 * Function implementations for the class MoleculeListClass.
    4  * 
     4 *
    55 */
    66
    … …  
    1313MoleculeListClass::MoleculeListClass()
    1414{
    15 };
    16 
    17 /** constructor for MoleculeListClass.
    18  * \param NumMolecules number of molecules to allocate for
    19  * \param NumAtoms number of atoms to allocate for
    20  */
    21 MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
    22 {
    23         ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
    24         for (int i=NumMolecules;i--;)
    25                 ListOfMolecules[i] = NULL;
    26         NumberOfMolecules = NumMolecules;
    27         NumberOfTopAtoms = NumAtoms;
    28 };
    29 
     15  // empty lists
     16  ListOfMolecules.clear();
     17  MaxIndex = 1;
     18};
    3019
    3120/** Destructor for MoleculeListClass.
    … …  
    3322MoleculeListClass::~MoleculeListClass()
    3423{
    35         cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
    36         for (int i=NumberOfMolecules;i--;) {
    37                 if (ListOfMolecules[i] != NULL) { // if NULL don't free
    38                         cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
    39                         delete(ListOfMolecules[i]);
    40                         ListOfMolecules[i] = NULL;
    41                 }
    42         }
    43         cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
    44         Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
     24  cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
     25  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     26    cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl;
     27    delete (*ListRunner);
     28  }
     29  cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
     30  ListOfMolecules.clear(); // empty list
     31};
     32
     33/** Insert a new molecule into the list and set its number.
     34 * \param *mol molecule to add to list.
     35 * \return true - add successful
     36 */
     37bool MoleculeListClass::insert(molecule *mol)
     38{
     39  mol->IndexNr = MaxIndex++;
     40  ListOfMolecules.push_back(mol);
    4541};
    4642
    … …  
    5248int MolCompare(const void *a, const void *b)
    5349{
    54         int *aList = NULL, *bList = NULL;
    55         int Count, Counter, aCounter, bCounter;
    56         int flag;
    57         atom *aWalker = NULL;
    58         atom *bWalker = NULL;
    59        
    60         // sort each atom list and put the numbers into a list, then go through
    61         //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
    62         if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
    63                 return -1;
    64         } else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
    65                 return +1;
    66                 else {
    67                         Count = (**(molecule **)a).AtomCount;
    68                         aList = new int[Count];
    69                         bList = new int[Count];
    70        
    71                         // fill the lists
    72                         aWalker = (**(molecule **)a).start;
    73                         bWalker = (**(molecule **)b).start;
    74                         Counter = 0;
    75                         aCounter = 0;
    76                         bCounter = 0;
    77                         while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
    78                                 aWalker = aWalker->next;
    79                                 bWalker = bWalker->next;
    80                                 if (aWalker->GetTrueFather() == NULL)
    81                                         aList[Counter] = Count + (aCounter++);
    82                                 else
    83                                         aList[Counter] = aWalker->GetTrueFather()->nr;
    84                                 if (bWalker->GetTrueFather() == NULL)
    85                                         bList[Counter] = Count + (bCounter++);
    86                                 else
    87                                         bList[Counter] = bWalker->GetTrueFather()->nr;
    88                                 Counter++;
    89                         }
    90                         // check if AtomCount was for real
    91                         flag = 0;
    92                         if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
    93                                 flag = -1;
    94                         } else {
    95                                 if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
    96                                         flag = 1;
    97                         }
    98                         if (flag == 0) {
    99                                 // sort the lists
    100                                 gsl_heapsort(aList,Count, sizeof(int), CompareDoubles);
    101                                 gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
    102                                 // compare the lists
    103                                
    104                                 flag = 0;
    105                                 for(int i=0;i<Count;i++) {
    106                                         if (aList[i] < bList[i]) {
    107                                                 flag = -1;
    108                                         } else {
    109                                                 if (aList[i] > bList[i])
    110                                                         flag = 1;
    111                                         }
    112                                         if (flag != 0)
    113                                                 break;
    114                                 }
    115                         }
    116                         delete[](aList);
    117                         delete[](bList);
    118                         return flag;
    119                 }
    120         }
    121         return  -1;
     50  int *aList = NULL, *bList = NULL;
     51  int Count, Counter, aCounter, bCounter;
     52  int flag;
     53  atom *aWalker = NULL;
     54  atom *bWalker = NULL;
     55
     56  // sort each atom list and put the numbers into a list, then go through
     57  //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
     58  if ((**(molecule **) a).AtomCount < (**(molecule **) b).AtomCount) {
     59    return -1;
     60  } else {
     61    if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount)
     62      return +1;
     63    else {
     64      Count = (**(molecule **) a).AtomCount;
     65      aList = new int[Count];
     66      bList = new int[Count];
     67
     68      // fill the lists
     69      aWalker = (**(molecule **) a).start;
     70      bWalker = (**(molecule **) b).start;
     71      Counter = 0;
     72      aCounter = 0;
     73      bCounter = 0;
     74      while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
     75        aWalker = aWalker->next;
     76        bWalker = bWalker->next;
     77        if (aWalker->GetTrueFather() == NULL)
     78          aList[Counter] = Count + (aCounter++);
     79        else
     80          aList[Counter] = aWalker->GetTrueFather()->nr;
     81        if (bWalker->GetTrueFather() == NULL)
     82          bList[Counter] = Count + (bCounter++);
     83        else
     84          bList[Counter] = bWalker->GetTrueFather()->nr;
     85        Counter++;
     86      }
     87      // check if AtomCount was for real
     88      flag = 0;
     89      if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
     90        flag = -1;
     91      } else {
     92        if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end))
     93          flag = 1;
     94      }
     95      if (flag == 0) {
     96        // sort the lists
     97        gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
     98        gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
     99        // compare the lists
     100
     101        flag = 0;
     102        for (int i = 0; i < Count; i++) {
     103          if (aList[i] < bList[i]) {
     104            flag = -1;
     105          } else {
     106            if (aList[i] > bList[i])
     107              flag = 1;
     108          }
     109          if (flag != 0)
     110            break;
     111        }
     112      }
     113      delete[] (aList);
     114      delete[] (bList);
     115      return flag;
     116    }
     117  }
     118  return -1;
     119};
     120
     121/** Output of a list of all molecules.
     122 * \param *out output stream
     123 */
     124void MoleculeListClass::Enumerate(ofstream *out)
     125{
     126  int i=1;
     127  element* Elemental = NULL;
     128  atom *Walker = NULL;
     129  int Counts[MAX_ELEMENTS];
     130
     131  // header
     132  *out << "Index\tName\tNo.Atoms\tformula" << endl;
     133  cout << Verbose(0) << "-----------------------------------------------" << endl;
     134  if (ListOfMolecules.size() == 0)
     135    *out << "\tNone" << endl;
     136  else {
     137    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     138      // reset element counts
     139      for (int j = 0; j<MAX_ELEMENTS;j++)
     140        Counts[j] = 0;
     141      // count atoms per element
     142      Walker = (*ListRunner)->start;
     143      while (Walker->next != (*ListRunner)->end) {
     144        Walker = Walker->next;
     145        Counts[Walker->type->Z]++;
     146      }
     147      // output Index, Name, number of atoms, chemical formula
     148      *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
     149      Elemental = (*ListRunner)->elemente->end;
     150      while(Elemental != (*ListRunner)->elemente->start) {
     151        Elemental = Elemental->previous;
     152        if (Counts[Elemental->Z] != 0)
     153          *out << Elemental->symbol << Counts[Elemental->Z];
     154      }
     155      *out << endl;
     156    }
     157  }
     158};
     159
     160/** Returns the molecule with the given index \a index.
     161 * \param index index of the desired molecule
     162 * \return pointer to molecule structure, NULL if not found
     163 */
     164molecule * MoleculeListClass::ReturnIndex(int index)
     165{
     166  int count = 1;
     167  MoleculeList::iterator ListRunner = ListOfMolecules.begin();
     168  for(; ((ListRunner != ListOfMolecules.end()) && (count < index)); ListRunner++);
     169  if (count == index)
     170    return (*ListRunner);
     171  else
     172    return NULL;
     173};
     174
     175/** Simple merge of two molecules into one.
     176 * \param *mol destination molecule
     177 * \param *srcmol source molecule
     178 * \return true - merge successful, false - merge failed (probably due to non-existant indices
     179 */
     180bool MoleculeListClass::SimpleMerge(molecule *mol, molecule *srcmol)
     181{
     182  if (srcmol == NULL)
     183    return false;
     184
     185  // put all molecules of src into mol
     186  atom *Walker = srcmol->start;
     187  atom *NextAtom = Walker->next;
     188  while (NextAtom != srcmol->end) {
     189    Walker = NextAtom;
     190    NextAtom = Walker->next;
     191    srcmol->UnlinkAtom(Walker);
     192    mol->AddAtom(Walker);
     193  }
     194
     195  // remove src
     196  ListOfMolecules.remove(srcmol);
     197  delete(srcmol);
     198  return true;
     199};
     200
     201/** Simple add of one molecules into another.
     202 * \param *mol destination molecule
     203 * \param *srcmol source molecule
     204 * \return true - merge successful, false - merge failed (probably due to non-existant indices
     205 */
     206bool MoleculeListClass::SimpleAdd(molecule *mol, molecule *srcmol)
     207{
     208  if (srcmol == NULL)
     209    return false;
     210
     211  // put all molecules of src into mol
     212  atom *Walker = srcmol->start;
     213  atom *NextAtom = Walker->next;
     214  while (NextAtom != srcmol->end) {
     215    Walker = NextAtom;
     216    NextAtom = Walker->next;
     217    Walker = mol->AddCopyAtom(Walker);
     218    Walker->father = Walker;
     219  }
     220
     221  return true;
     222};
     223
     224/** Simple merge of a given set of molecules into one.
     225 * \param *mol destination molecule
     226 * \param *src index of set of source molecule
     227 * \param N number of source molecules
     228 * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
     229 */
     230bool MoleculeListClass::SimpleMultiMerge(molecule *mol, int *src, int N)
     231{
     232  bool status = true;
     233  // check presence of all source molecules
     234  for (int i=0;i<N;i++) {
     235    molecule *srcmol = ReturnIndex(src[i]);
     236    status = status && SimpleMerge(mol, srcmol);
     237  }
     238  return status;
     239};
     240
     241/** Simple add of a given set of molecules into one.
     242 * \param *mol destination molecule
     243 * \param *src index of set of source molecule
     244 * \param N number of source molecules
     245 * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
     246 */
     247bool MoleculeListClass::SimpleMultiAdd(molecule *mol, int *src, int N)
     248{
     249  bool status = true;
     250  // check presence of all source molecules
     251  for (int i=0;i<N;i++) {
     252    molecule *srcmol = ReturnIndex(src[i]);
     253    status = status && SimpleAdd(mol, srcmol);
     254  }
     255  return status;
     256};
     257
     258/** Scatter merge of a given set of molecules into one.
     259 * Scatter merge distributes the molecules in such a manner that they don't overlap.
     260 * \param *mol destination molecule
     261 * \param *src index of set of source molecule
     262 * \param N number of source molecules
     263 * \return true - merge successful, false - merge failed (probably due to non-existant indices
     264 * \TODO find scatter center for each src molecule
     265 */
     266bool MoleculeListClass::ScatterMerge(molecule *mol, int *src, int N)
     267{
     268  // check presence of all source molecules
     269  for (int i=0;i<N;i++) {
     270    // get pointer to src molecule
     271    molecule *srcmol = ReturnIndex(src[i]);
     272    if (srcmol == NULL)
     273      return false;
     274  }
     275  // adapt each Center
     276  for (int i=0;i<N;i++) {
     277    // get pointer to src molecule
     278    molecule *srcmol = ReturnIndex(src[i]);
     279    //srcmol->Center.Zero();
     280    srcmol->Translate(&srcmol->Center);
     281  }
     282  // perform a simple multi merge
     283  SimpleMultiMerge(mol, src, N);
     284  return true;
     285};
     286
     287/** Embedding merge of a given set of molecules into one.
     288 * Embedding merge inserts one molecule into the other.
     289 * \param *mol destination molecule
     290 * \param *srcmol source molecule
     291 * \return true - merge successful, false - merge failed (probably due to non-existant indices
     292 * \TODO find embedding center
     293 */
     294bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
     295{
     296  if (srcmol == NULL)
     297    return false;
     298
     299  // calculate center for merge
     300  //srcmol->Center.Zero();
     301
     302  // perform simple merge
     303  SimpleMerge(mol, srcmol);
     304  return true;
    122305};
    123306
    … …  
    127310void MoleculeListClass::Output(ofstream *out)
    128311{
    129         *out<< Verbose(1) << "MoleculeList: ";
    130         for (int i=0;i<NumberOfMolecules;i++)
    131                 *out << ListOfMolecules[i] << "\t";
    132         *out << endl;
     312  *out << Verbose(1) << "MoleculeList: ";
     313  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
     314    *out << *ListRunner << "\t";
     315  *out << endl;
    133316};
    134317
    … …  
    142325bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
    143326{
    144         atom *Walker = NULL;
    145         atom *Runner = NULL;
    146         double ***FitConstant = NULL, **correction = NULL;
    147         int a,b;
    148         ofstream output;
    149         ifstream input;
    150         string line;
    151         stringstream zeile;
    152         double distance;
    153         char ParsedLine[1023];
    154         double tmp;
    155         char *FragmentNumber = NULL;
    156 
    157         cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
    158         // 0. parse in fit constant files that should have the same dimension as the final energy files
    159         // 0a. find dimension of matrices with constants
    160         line = path;
    161         line.append("/");
    162         line += FRAGMENTPREFIX;
    163         line += "1";
    164         line += FITCONSTANTSUFFIX;
    165         input.open(line.c_str());
    166         if (input == NULL) {
    167                 cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
    168                 return false;
    169         }
    170         a=0;
    171         b=-1; // we overcount by one
    172         while (!input.eof()) {
    173                 input.getline(ParsedLine, 1023);
    174                 zeile.str(ParsedLine);
    175                 int i=0;
    176                 while (!zeile.eof()) {
    177                         zeile >> distance;
    178                         i++;
    179                 }
    180                 if (i > a)
    181                         a = i;
    182                 b++;
    183         }
    184         cout << "I recognized " << a << " columns and " << b << " rows, ";
    185         input.close();
    186        
    187         // 0b. allocate memory for constants
    188         FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
    189         for (int k=0;k<3;k++) {
    190                 FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
    191                 for (int i=a;i--;) {
    192                         FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
    193                 }
    194         }
    195         // 0c. parse in constants
    196         for (int i=0;i<3;i++) {
    197                 line = path;
    198                 line.append("/");
    199                 line += FRAGMENTPREFIX;
    200                 sprintf(ParsedLine, "%d", i+1);
    201                 line += ParsedLine;
    202                 line += FITCONSTANTSUFFIX;
    203                 input.open(line.c_str());
    204                 if (input == NULL) {
    205                         cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
    206                         return false;
    207                 }
    208                 int k = 0,l;
    209                 while ((!input.eof()) && (k < b)) {
    210                         input.getline(ParsedLine, 1023);
    211                         //cout << "Current Line: " << ParsedLine << endl;
    212                         zeile.str(ParsedLine);
    213                         zeile.clear();
    214                         l = 0;
    215                         while ((!zeile.eof()) && (l < a)) {
    216                                 zeile >> FitConstant[i][l][k];
    217                                 //cout << FitConstant[i][l][k] << "\t";
    218                                 l++;
    219                         }
    220                         //cout << endl;
    221                         k++;
    222                 }
    223                 input.close();
    224         }
    225         for(int k=0;k<3;k++) {
    226                 cout << "Constants " << k << ":" << endl;
    227                 for (int j=0;j<b;j++) {
    228                         for (int i=0;i<a;i++) {
    229                                 cout << FitConstant[k][i][j] << "\t";
    230                         }
    231                         cout << endl;
    232                 }
    233                 cout << endl;
    234         }
    235        
    236         // 0d. allocate final correction matrix
    237         correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
    238         for (int i=a;i--;)
    239                 correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
    240                                
    241         // 1a. go through every molecule in the list
    242         for(int i=NumberOfMolecules;i--;) {
    243                 // 1b. zero final correction matrix
    244                 for (int k=a;k--;)
    245                         for (int j=b;j--;)
    246                                 correction[k][j] = 0.;
    247                 // 2. take every hydrogen that is a saturated one
    248                 Walker = ListOfMolecules[i]->start;
    249                 while (Walker->next != ListOfMolecules[i]->end) {
    250                         Walker = Walker->next;
    251                         //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
    252                         if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
    253                                 Runner = ListOfMolecules[i]->start;
    254                                 while (Runner->next != ListOfMolecules[i]->end) {
    255                                         Runner = Runner->next;
    256                                         //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
    257                                         // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
    258                                         if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {      // (hydrogens have only one bonding partner!)
    259                                                 // 4. evaluate the morse potential for each matrix component and add up
    260                                                 distance = Runner->x.Distance(&Walker->x);
    261                                                 //cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
    262                                                 for(int k=0;k<a;k++) {
    263                                                         for (int j=0;j<b;j++) {
    264                                                                 switch(k) {
    265                                                                 case 1:
    266                                                                 case 7:
    267                                                                 case 11:
    268                                                                         tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
    269                                                                         break;
    270                                                                 default:
    271                                                                         tmp = FitConstant[0][k][j] * pow( distance,     FitConstant[1][k][j]) + FitConstant[2][k][j];
    272                                                                 };
    273                                                                 correction[k][j] -= tmp;                // ground state is actually lower (disturbed by additional interaction)
    274                                                                 //cout << tmp << "\t";
    275                                                         }
    276                                                         //cout << endl;
    277                                                 }
    278                                                 //cout << endl;
    279                                         }
    280                                 }
    281                         }
    282                 }
    283                 // 5. write final matrix to file
    284                 line = path;
    285                 line.append("/");
    286                 line += FRAGMENTPREFIX;
    287                 FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
    288                 line += FragmentNumber;
    289                 delete(FragmentNumber);
    290                 line += HCORRECTIONSUFFIX;
    291                 output.open(line.c_str());
    292                 output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
    293                 for (int j=0;j<b;j++) {
    294                         for(int i=0;i<a;i++)
    295                                 output << correction[i][j] << "\t";
    296                         output << endl;
    297                 }
    298                 output.close();
    299         }
    300         line = path;
    301         line.append("/");
    302         line += HCORRECTIONSUFFIX;
    303         output.open(line.c_str());
    304         output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
    305         for (int j=0;j<b;j++) {
    306                 for(int i=0;i<a;i++)
    307                         output << 0 << "\t";
    308                 output << endl;
    309         }
    310         output.close();
    311         // 6. free memory of parsed matrices
    312         FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
    313         for (int k=0;k<3;k++) {
    314                 FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
    315                 for (int i=a;i--;) {
    316                         FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
    317                 }
    318         }
    319         cout << "done." << endl;
    320         return true;
     327  atom *Walker = NULL;
     328  atom *Runner = NULL;
     329  double ***FitConstant = NULL, **correction = NULL;
     330  int a, b;
     331  ofstream output;
     332  ifstream input;
     333  string line;
     334  stringstream zeile;
     335  double distance;
     336  char ParsedLine[1023];
     337  double tmp;
     338  char *FragmentNumber = NULL;
     339
     340  cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
     341  // 0. parse in fit constant files that should have the same dimension as the final energy files
     342  // 0a. find dimension of matrices with constants
     343  line = path;
     344  line.append("/");
     345  line += FRAGMENTPREFIX;
     346  line += "1";
     347  line += FITCONSTANTSUFFIX;
     348  input.open(line.c_str());
     349  if (input == NULL) {
     350    cerr << endl << "Unable to open " << line << ", is the directory correct?"
     351        << endl;
     352    return false;
     353  }
     354  a = 0;
     355  b = -1; // we overcount by one
     356  while (!input.eof()) {
     357    input.getline(ParsedLine, 1023);
     358    zeile.str(ParsedLine);
     359    int i = 0;
     360    while (!zeile.eof()) {
     361      zeile >> distance;
     362      i++;
     363    }
     364    if (i > a)
     365      a = i;
     366    b++;
     367  }
     368  cout << "I recognized " << a << " columns and " << b << " rows, ";
     369  input.close();
     370
     371  // 0b. allocate memory for constants
     372  FitConstant = (double ***) Malloc(sizeof(double **) * 3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
     373  for (int k = 0; k < 3; k++) {
     374    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
     375    for (int i = a; i--;) {
     376      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     377    }
     378  }
     379  // 0c. parse in constants
     380  for (int i = 0; i < 3; i++) {
     381    line = path;
     382    line.append("/");
     383    line += FRAGMENTPREFIX;
     384    sprintf(ParsedLine, "%d", i + 1);
     385    line += ParsedLine;
     386    line += FITCONSTANTSUFFIX;
     387    input.open(line.c_str());
     388    if (input == NULL) {
     389      cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
     390      return false;
     391    }
     392    int k = 0, l;
     393    while ((!input.eof()) && (k < b)) {
     394      input.getline(ParsedLine, 1023);
     395      //cout << "Current Line: " << ParsedLine << endl;
     396      zeile.str(ParsedLine);
     397      zeile.clear();
     398      l = 0;
     399      while ((!zeile.eof()) && (l < a)) {
     400        zeile >> FitConstant[i][l][k];
     401        //cout << FitConstant[i][l][k] << "\t";
     402        l++;
     403      }
     404      //cout << endl;
     405      k++;
     406    }
     407    input.close();
     408  }
     409  for (int k = 0; k < 3; k++) {
     410    cout << "Constants " << k << ":" << endl;
     411    for (int j = 0; j < b; j++) {
     412      for (int i = 0; i < a; i++) {
     413        cout << FitConstant[k][i][j] << "\t";
     414      }
     415      cout << endl;
     416    }
     417    cout << endl;
     418  }
     419
     420  // 0d. allocate final correction matrix
     421  correction = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **correction");
     422  for (int i = a; i--;)
     423    correction[i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
     424
     425  // 1a. go through every molecule in the list
     426  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     427    // 1b. zero final correction matrix
     428    for (int k = a; k--;)
     429      for (int j = b; j--;)
     430        correction[k][j] = 0.;
     431    // 2. take every hydrogen that is a saturated one
     432    Walker = (*ListRunner)->start;
     433    while (Walker->next != (*ListRunner)->end) {
     434      Walker = Walker->next;
     435      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
     436      if ((Walker->type->Z == 1) && ((Walker->father == NULL)
     437          || (Walker->father->type->Z != 1))) { // if it's a hydrogen
     438        Runner = (*ListRunner)->start;
     439        while (Runner->next != (*ListRunner)->end) {
     440          Runner = Runner->next;
     441          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
     442          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
     443          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
     444            // 4. evaluate the morse potential for each matrix component and add up
     445            distance = Runner->x.Distance(&Walker->x);
     446            //cout << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
     447            for (int k = 0; k < a; k++) {
     448              for (int j = 0; j < b; j++) {
     449                switch (k) {
     450                  case 1:
     451                  case 7:
     452                  case 11:
     453                    tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
     454                    break;
     455                  default:
     456                    tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
     457                };
     458                correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
     459                //cout << tmp << "\t";
     460              }
     461              //cout << endl;
     462            }
     463            //cout << endl;
     464          }
     465        }
     466      }
     467    }
     468    // 5. write final matrix to file
     469    line = path;
     470    line.append("/");
     471    line += FRAGMENTPREFIX;
     472    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
     473    line += FragmentNumber;
     474    delete (FragmentNumber);
     475    line += HCORRECTIONSUFFIX;
     476    output.open(line.c_str());
     477    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
     478    for (int j = 0; j < b; j++) {
     479      for (int i = 0; i < a; i++)
     480        output << correction[i][j] << "\t";
     481      output << endl;
     482    }
     483    output.close();
     484  }
     485  line = path;
     486  line.append("/");
     487  line += HCORRECTIONSUFFIX;
     488  output.open(line.c_str());
     489  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
     490  for (int j = 0; j < b; j++) {
     491    for (int i = 0; i < a; i++)
     492      output << 0 << "\t";
     493    output << endl;
     494  }
     495  output.close();
     496  // 6. free memory of parsed matrices
     497  FitConstant = (double ***) Malloc(sizeof(double **) * a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
     498  for (int k = 0; k < 3; k++) {
     499    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
     500    for (int i = a; i--;) {
     501      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
     502    }
     503  }
     504  cout << "done." << endl;
     505  return true;
    321506};
    322507
    … …  
    327512 * \return true - file written successfully, false - writing failed
    328513 */
    329 bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex)
    330 {
    331         bool status = true;
    332         ofstream ForcesFile;
    333         stringstream line;
    334         atom *Walker = NULL;
    335         element *runner = NULL;
    336 
    337         // open file for the force factors
    338         *out << Verbose(1) << "Saving   force factors ... ";
    339         line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
    340         ForcesFile.open(line.str().c_str(), ios::out);
    341         if (ForcesFile != NULL) {
    342                 //cout << Verbose(1) << "Final AtomicForcesList: ";
    343                 //output << prefix << "Forces" << endl;
    344                 for(int j=0;j<NumberOfMolecules;j++) {
    345                         //if (TEList[j] != 0) {
    346                         runner = ListOfMolecules[j]->elemente->start;
    347                         while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
    348                                 runner = runner->next;
    349                                 if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
    350                                         Walker = ListOfMolecules[j]->start;
    351                                         while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
    352                                                 Walker = Walker->next;
    353                                                 if (Walker->type->Z == runner->Z) {
    354                                                         if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
    355                                                                 //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
    356                                                                 ForcesFile <<   SortIndex[Walker->GetTrueFather()->nr] << "\t";
    357                                                                 } else  // otherwise a -1 to indicate an added saturation hydrogen
    358                                                                         ForcesFile << "-1\t";
    359                                                         }
    360                                                 }
    361                                         }
    362                         }
    363                         ForcesFile << endl;
    364                 }
    365                 ForcesFile.close();
    366                 *out << Verbose(1) << "done." << endl;
    367         } else {
    368                 status = false;
    369                 *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
    370         }
    371         ForcesFile.close();
    372        
    373         return status;
     514bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path,
     515    int *SortIndex)
     516{
     517  bool status = true;
     518  ofstream ForcesFile;
     519  stringstream line;
     520  atom *Walker = NULL;
     521  element *runner = NULL;
     522
     523  // open file for the force factors
     524  *out << Verbose(1) << "Saving force factors ... ";
     525  line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
     526  ForcesFile.open(line.str().c_str(), ios::out);
     527  if (ForcesFile != NULL) {
     528    //cout << Verbose(1) << "Final AtomicForcesList: ";
     529    //output << prefix << "Forces" << endl;
     530    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     531      runner = (*ListRunner)->elemente->start;
     532      while (runner->next != (*ListRunner)->elemente->end) { // go through every element
     533        runner = runner->next;
     534        if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms
     535          Walker = (*ListRunner)->start;
     536          while (Walker->next != (*ListRunner)->end) { // go through every atom of this element
     537            Walker = Walker->next;
     538            if (Walker->type->Z == runner->Z) {
     539              if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
     540                //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
     541                ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t";
     542              } else
     543                // otherwise a -1 to indicate an added saturation hydrogen
     544                ForcesFile << "-1\t";
     545            }
     546          }
     547        }
     548      }
     549      ForcesFile << endl;
     550    }
     551    ForcesFile.close();
     552    *out << Verbose(1) << "done." << endl;
     553  } else {
     554    status = false;
     555    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     556  }
     557  ForcesFile.close();
     558
     559  return status;
    374560};
    375561
    … …  
    380566 * \return true - success (each file was written), false - something went wrong.
    381567 */
    382 bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
    383 {
    384         ofstream outputFragment;
    385         char FragmentName[MAXSTRINGSIZE];
    386         char PathBackup[MAXSTRINGSIZE];
    387         bool result = true;
    388         bool intermediateResult = true;
    389         atom *Walker = NULL;
    390         Vector BoxDimension;
    391         char *FragmentNumber = NULL;
    392         char *path = NULL;
    393         int FragmentCounter = 0;
    394         ofstream output;
    395        
    396         // store the fragments as config and as xyz
    397         for(int i=0;i<NumberOfMolecules;i++) {
    398                 // save default path as it is changed for each fragment
    399                 path = configuration->GetDefaultPath();
    400                 if (path != NULL)
    401                         strcpy(PathBackup, path);
    402                 else
    403                         cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
    404 
    405                 // correct periodic
    406                 ListOfMolecules[i]->ScanForPeriodicCorrection(out);
    407 
    408                 // output xyz file
    409                 FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
    410                 sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    411                 outputFragment.open(FragmentName, ios::out);
    412                 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
    413                 if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
    414                         *out << " done." << endl;
    415                 else
    416                         *out << " failed." << endl;
    417                 result = result && intermediateResult;
    418                 outputFragment.close();
    419                 outputFragment.clear();
    420 
    421                 // list atoms in fragment for debugging
    422                 *out << Verbose(2) << "Contained atoms: ";
    423                 Walker = ListOfMolecules[i]->start;
    424                 while (Walker->next != ListOfMolecules[i]->end) {
    425                         Walker = Walker->next;
    426                         *out << Walker->Name << " ";
    427                 }
    428                 *out << endl;
    429                
    430                 // center on edge
    431                 ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
    432                 ListOfMolecules[i]->SetBoxDimension(&BoxDimension);     // update Box of atoms by boundary
    433                 int j = -1;
    434                 for (int k=0;k<NDIM;k++) {
    435                         j += k+1;
    436                         BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
    437                         ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
    438                 }
    439                 ListOfMolecules[i]->Translate(&BoxDimension);
    440 
    441                 // also calculate necessary orbitals
    442                 ListOfMolecules[i]->CountElements();    // this is a bugfix, atoms should should actually be added correctly to this fragment
    443                 ListOfMolecules[i]->CalculateOrbitals(*configuration);
    444                
    445                 // change path in config
    446                 //strcpy(PathBackup, configuration->configpath);
    447                 sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
    448                 configuration->SetDefaultPath(FragmentName);
    449                
    450                 // and save as config
    451                 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    452                 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
    453                 if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
    454                         *out << " done." << endl;
    455                 else
    456                         *out << " failed." << endl;
    457                 result = result && intermediateResult;
    458 
    459                 // restore old config
    460                 configuration->SetDefaultPath(PathBackup);
    461 
    462 
    463                 // and save as mpqc input file
    464                 sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    465                 *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
    466                 if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
    467                         *out << " done." << endl;
    468                 else
    469                         *out << " failed." << endl;
    470                          
    471                 result = result && intermediateResult;
    472                 //outputFragment.close();
    473                 //outputFragment.clear();
    474                 delete(FragmentNumber);
    475                 //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
    476         }
    477         cout << " done." << endl;
    478        
    479         // printing final number
    480         *out << "Final number of fragments: " << FragmentCounter << "." << endl;
    481                        
    482         return result;
    483 };
     568bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out,
     569    config *configuration, int *SortIndex)
     570{
     571  ofstream outputFragment;
     572  char FragmentName[MAXSTRINGSIZE];
     573  char PathBackup[MAXSTRINGSIZE];
     574  bool result = true;
     575  bool intermediateResult = true;
     576  atom *Walker = NULL;
     577  Vector BoxDimension;
     578  char *FragmentNumber = NULL;
     579  char *path = NULL;
     580  int FragmentCounter = 0;
     581  ofstream output;
     582
     583  // store the fragments as config and as xyz
     584  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
     585    // save default path as it is changed for each fragment
     586    path = configuration->GetDefaultPath();
     587    if (path != NULL)
     588      strcpy(PathBackup, path);
     589    else
     590      cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
     591
     592    // correct periodic
     593    (*ListRunner)->ScanForPeriodicCorrection(out);
     594
     595    // output xyz file
     596    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
     597    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     598    outputFragment.open(FragmentName, ios::out);
     599    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ...";
     600    if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
     601      *out << " done." << endl;
     602    else
     603      *out << " failed." << endl;
     604    result = result && intermediateResult;
     605    outputFragment.close();
     606    outputFragment.clear();
     607
     608    // list atoms in fragment for debugging
     609    *out << Verbose(2) << "Contained atoms: ";
     610    Walker = (*ListRunner)->start;
     611    while (Walker->next != (*ListRunner)->end) {
     612      Walker = Walker->next;
     613      *out << Walker->Name << " ";
     614    }
     615    *out << endl;
     616
     617    // center on edge
     618    (*ListRunner)->CenterEdge(out, &BoxDimension);
     619    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
     620    int j = -1;
     621    for (int k = 0; k < NDIM; k++) {
     622      j += k + 1;
     623      BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
     624      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
     625    }
     626    (*ListRunner)->Translate(&BoxDimension);
     627
     628    // also calculate necessary orbitals
     629    (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
     630    (*ListRunner)->CalculateOrbitals(*configuration);
     631
     632    // change path in config
     633    //strcpy(PathBackup, configuration->configpath);
     634    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     635    configuration->SetDefaultPath(FragmentName);
     636
     637    // and save as config
     638    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     639    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...";
     640    if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner))))
     641      *out << " done." << endl;
     642    else
     643      *out << " failed." << endl;
     644    result = result && intermediateResult;
     645
     646    // restore old config
     647    configuration->SetDefaultPath(PathBackup);
     648
     649    // and save as mpqc input file
     650    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     651    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...";
     652    if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner))))
     653      *out << " done." << endl;
     654    else
     655      *out << " failed." << endl;
     656
     657    result = result && intermediateResult;
     658    //outputFragment.close();
     659    //outputFragment.clear();
     660    delete (FragmentNumber);
     661    //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
     662  }
     663  cout << " done." << endl;
     664
     665  // printing final number
     666  *out << "Final number of fragments: " << FragmentCounter << "." << endl;
     667
     668  return result;
     669};
     670
     671/** Counts the number of molecules with the molecule::ActiveFlag set.
     672 * \return number of molecules with ActiveFlag set to true.
     673 */
     674int MoleculeListClass::NumberOfActiveMolecules()
     675{
     676  int count = 0;
     677  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
     678    count += ((*ListRunner)->ActiveFlag ? 1 : 0);
     679  return count;
     680};
     681
    484682
    485683/******************************************* Class MoleculeLeafClass ************************************************/
    … …  
    492690MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
    493691{
    494 //      if (Up != NULL)
    495 //              if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
    496 //                      Up->DownLeaf = this;
    497 //      UpLeaf = Up;
    498 //      DownLeaf = NULL;
    499         Leaf = NULL;
    500         previous = PreviousLeaf;
    501         if (previous != NULL) {
    502                 MoleculeLeafClass *Walker = previous->next;
    503                 previous->next = this;
    504                 next = Walker;
    505         } else {
    506                 next = NULL;
    507         }
     692  //    if (Up != NULL)
     693  //            if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
     694  //                    Up->DownLeaf = this;
     695  //    UpLeaf = Up;
     696  //    DownLeaf = NULL;
     697  Leaf = NULL;
     698  previous = PreviousLeaf;
     699  if (previous != NULL) {
     700    MoleculeLeafClass *Walker = previous->next;
     701    previous->next = this;
     702    next = Walker;
     703  } else {
     704    next = NULL;
     705  }
    508706};
    509707
    … …  
    512710MoleculeLeafClass::~MoleculeLeafClass()
    513711{
    514 //      if (DownLeaf != NULL) {// drop leaves further down
    515 //              MoleculeLeafClass *Walker = DownLeaf;
    516 //              MoleculeLeafClass *Next;
    517 //              do {
    518 //                      Next = Walker->NextLeaf;
    519 //                      delete(Walker);
    520 //                      Walker = Next;
    521 //              } while (Walker != NULL);
    522 //              // Last Walker sets DownLeaf automatically to NULL
    523 //      }
    524         // remove the leaf itself
    525         if (Leaf != NULL) {
    526                 delete(Leaf);
    527                 Leaf = NULL;
    528         }
    529         // remove this Leaf from level list
    530         if (previous != NULL)   
    531                 previous->next = next;
    532 //      } else { // we are first in list (connects to UpLeaf->DownLeaf)
    533 //              if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
    534 //                      NextLeaf->UpLeaf = UpLeaf;      // either null as we are top level or the upleaf of the first node
    535 //              if (UpLeaf != NULL)
    536 //                      UpLeaf->DownLeaf = NextLeaf;    // either null as we are only leaf or NextLeaf if we are just the first
    537 //      }
    538 //      UpLeaf = NULL;
    539         if (next != NULL) // are we last in list
    540                 next->previous = previous;
    541         next = NULL;
    542         previous = NULL;
     712  //    if (DownLeaf != NULL) {// drop leaves further down
     713  //            MoleculeLeafClass *Walker = DownLeaf;
     714  //            MoleculeLeafClass *Next;
     715  //            do {
     716  //                    Next = Walker->NextLeaf;
     717  //                    delete(Walker);
     718  //                    Walker = Next;
     719  //            } while (Walker != NULL);
     720  //            // Last Walker sets DownLeaf automatically to NULL
     721  //    }
     722  // remove the leaf itself
     723  if (Leaf != NULL) {
     724    delete (Leaf);
     725    Leaf = NULL;
     726  }
     727  // remove this Leaf from level list
     728  if (previous != NULL)
     729    previous->next = next;
     730  //    } else { // we are first in list (connects to UpLeaf->DownLeaf)
     731  //            if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
     732  //                    NextLeaf->UpLeaf = UpLeaf;      // either null as we are top level or the upleaf of the first node
     733  //            if (UpLeaf != NULL)
     734  //                    UpLeaf->DownLeaf = NextLeaf;    // either null as we are only leaf or NextLeaf if we are just the first
     735  //    }
     736  //    UpLeaf = NULL;
     737  if (next != NULL) // are we last in list
     738    next->previous = previous;
     739  next = NULL;
     740  previous = NULL;
    543741};
    544742
    … …  
    550748bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
    551749{
    552         return false;
     750  return false;
    553751};
    554752
    … …  
    564762bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
    565763{
    566         atom *Walker = NULL, *OtherWalker = NULL;
    567         bond *Binder = NULL;
    568         bool status = true;
    569         int AtomNo;
    570 
    571         *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
    572         // fill ListOfLocalAtoms if NULL was given
    573         if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    574                 *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    575                 return false;
    576         }
    577        
    578         if (status) {
    579                 *out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
    580                 Walker = Leaf->start;
    581                 while (Walker->next != Leaf->end) {
    582                         Walker = Walker->next;
    583                         AtomNo = Walker->GetTrueFather()->nr;   // global id of the current walker
    584                         for(int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
    585                                 Binder = reference->ListOfBondsPerAtom[AtomNo][i];
    586                                 OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr];             // local copy of current bond partner of walker
    587                                 if (OtherWalker != NULL) {
    588                                         if (OtherWalker->nr > Walker->nr)
    589                                         Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
    590                                 } else {
    591                                         *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
    592                                         status = false;
    593                                 }
    594                         }
    595                 }
    596                 Leaf->CreateListOfBondsPerAtom(out);
    597                 FragmentCounter++;
    598                 if (next != NULL)
    599                         status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
    600                 FragmentCounter--;
    601         }
    602        
    603         if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    604                 // free the index lookup list
    605                 Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
    606                 if (FragmentCounter == 0) // first fragments frees the initial pointer to list
    607                         Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
    608         }
    609         FragmentCounter--;
    610         *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
    611         return status;
     764  atom *Walker = NULL, *OtherWalker = NULL;
     765  bond *Binder = NULL;
     766  bool status = true;
     767  int AtomNo;
     768
     769  *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
     770  // fill ListOfLocalAtoms if NULL was given
     771  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
     772    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
     773    return false;
     774  }
     775
     776  if (status) {
     777    *out << Verbose(1) << "Creating adjacency list for subgraph " << this
     778        << "." << endl;
     779    Walker = Leaf->start;
     780    while (Walker->next != Leaf->end) {
     781      Walker = Walker->next;
     782      AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker
     783      for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all
     784        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
     785        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
     786        if (OtherWalker != NULL) {
     787          if (OtherWalker->nr > Walker->nr)
     788            Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
     789        } else {
     790          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
     791          status = false;
     792        }
     793      }
     794    }
     795    Leaf->CreateListOfBondsPerAtom(out);
     796    FragmentCounter++;
     797    if (next != NULL)
     798      status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
     799    FragmentCounter--;
     800  }
     801
     802  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
     803    // free the index lookup list
     804    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
     805    if (FragmentCounter == 0) // first fragments frees the initial pointer to list
     806      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
     807  }
     808  FragmentCounter--;
     809  *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
     810  return status;
    612811};
    613812
    … …  
    620819 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
    621820 */
    622 bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
    623 {
    624         atom *Walker = NULL, *Father = NULL;
    625 
    626         if (RootStack != NULL) {
    627                 // find first root candidates
    628                 if (&(RootStack[FragmentCounter]) != NULL) {
    629                         RootStack[FragmentCounter].clear();     
    630                         Walker = Leaf->start;
    631                         while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
    632                                 Walker = Walker->next;
    633                                 Father = Walker->GetTrueFather();
    634                                 if (AtomMask[Father->nr]) // apply mask
    635                 #ifdef ADDHYDROGEN
    636                                         if (Walker->type->Z != 1) // skip hydrogen
    637                 #endif
    638                                                 RootStack[FragmentCounter].push_front(Walker->nr);
    639                         }
    640                         if (next != NULL)
    641                                 next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
    642                 }       else {
    643                         *out << Verbose(1) << "Rootstack[" << FragmentCounter   << "] is NULL." << endl;
    644                         return false;
    645                 }
    646                 FragmentCounter--;
    647                 return true;
    648         } else {
    649                 *out << Verbose(1) << "Rootstack is NULL." << endl;
    650                 return false;
    651         }
     821bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out,
     822    KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
     823{
     824  atom *Walker = NULL, *Father = NULL;
     825
     826  if (RootStack != NULL) {
     827    // find first root candidates
     828    if (&(RootStack[FragmentCounter]) != NULL) {
     829      RootStack[FragmentCounter].clear();
     830      Walker = Leaf->start;
     831      while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
     832        Walker = Walker->next;
     833        Father = Walker->GetTrueFather();
     834        if (AtomMask[Father->nr]) // apply mask
     835#ifdef ADDHYDROGEN
     836          if (Walker->type->Z != 1) // skip hydrogen
     837#endif
     838          RootStack[FragmentCounter].push_front(Walker->nr);
     839      }
     840      if (next != NULL)
     841        next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
     842    } else {
     843      *out << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl;
     844      return false;
     845    }
     846    FragmentCounter--;
     847    return true;
     848  } else {
     849    *out << Verbose(1) << "Rootstack is NULL." << endl;
     850    return false;
     851  }
    652852};
    653853
    … …  
    662862bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList)
    663863{
    664         bool status = true;
    665        
    666         int Counter = Count();
    667         if (ListOfLocalAtoms == NULL) { // allocated initial pointer
    668                 // allocate and set each field to NULL
    669                 ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
    670                 if (ListOfLocalAtoms != NULL) {
    671                         for (int i=Counter;i--;)
    672                                 ListOfLocalAtoms[i] = NULL;
    673                         FreeList = FreeList && true;
    674                 } else
    675                         status = false;
    676         }
    677        
    678         if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
    679                 status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
    680                 FreeList = FreeList && true;
    681         }
    682        
    683         return status;
     864  bool status = true;
     865
     866  int Counter = Count();
     867  if (ListOfLocalAtoms == NULL) { // allocated initial pointer
     868    // allocate and set each field to NULL
     869    ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **) * Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
     870    if (ListOfLocalAtoms != NULL) {
     871      for (int i = Counter; i--;)
     872        ListOfLocalAtoms[i] = NULL;
     873      FreeList = FreeList && true;
     874    } else
     875      status = false;
     876  }
     877
     878  if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
     879    status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
     880    FreeList = FreeList && true;
     881  }
     882
     883  return status;
    684884};
    685885
    … …  
    694894 * \retuen true - success, false - failure
    695895 */
    696 bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
    697 {
    698         bool status = true;
    699         int KeySetCounter = 0;
    700        
    701         *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
    702         // fill ListOfLocalAtoms if NULL was given
    703         if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    704                 *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
    705                 return false;
    706         }
    707 
    708         // allocate fragment list
    709         if (FragmentList == NULL) {
    710                 KeySetCounter = Count();
    711                 FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
    712                 for(int i=KeySetCounter;i--;)
    713                         FragmentList[i] = NULL;
    714                 KeySetCounter = 0;
    715         }
    716        
    717         if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
    718                 // assign scanned keysets
    719                 if (FragmentList[FragmentCounter] == NULL)
    720                         FragmentList[FragmentCounter] = new Graph;
    721                 KeySet *TempSet = new KeySet;
    722                 for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
    723                         if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
    724                                 // translate keyset to local numbers
    725                                 for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
    726                                         TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
    727                                 // insert into FragmentList
    728                                 FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
    729                         }
    730                         TempSet->clear();
    731                 }
    732                 delete(TempSet);
    733                 if (KeySetCounter == 0) {// if there are no keysets, delete the list
    734                         *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
    735                         delete(FragmentList[FragmentCounter]);
    736                 } else
    737                         *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
    738                 FragmentCounter++;
    739                 if (next != NULL)
    740                         next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
    741                 FragmentCounter--;
    742         } else
    743                 *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
    744        
    745         if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    746                 // free the index lookup list
    747                 Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
    748                 if (FragmentCounter == 0) // first fragments frees the initial pointer to list
    749                         Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
    750         }
    751         *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
    752         return status;
     896bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out,
     897    molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms,
     898    Graph **&FragmentList, int &FragmentCounter, bool FreeList)
     899{
     900  bool status = true;
     901  int KeySetCounter = 0;
     902
     903  *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
     904  // fill ListOfLocalAtoms if NULL was given
     905  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
     906    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
     907    return false;
     908  }
     909
     910  // allocate fragment list
     911  if (FragmentList == NULL) {
     912    KeySetCounter = Count();
     913    FragmentList = (Graph **) Malloc(sizeof(Graph *) * KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
     914    for (int i = KeySetCounter; i--;)
     915      FragmentList[i] = NULL;
     916    KeySetCounter = 0;
     917  }
     918
     919  if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
     920    // assign scanned keysets
     921    if (FragmentList[FragmentCounter] == NULL)
     922      FragmentList[FragmentCounter] = new Graph;
     923    KeySet *TempSet = new KeySet;
     924    for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers!
     925      if (ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
     926        // translate keyset to local numbers
     927        for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
     928          TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
     929        // insert into FragmentList
     930        FragmentList[FragmentCounter]->insert(GraphPair(*TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
     931      }
     932      TempSet->clear();
     933    }
     934    delete (TempSet);
     935    if (KeySetCounter == 0) {// if there are no keysets, delete the list
     936      *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
     937      delete (FragmentList[FragmentCounter]);
     938    } else
     939      *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
     940    FragmentCounter++;
     941    if (next != NULL)
     942      next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
     943    FragmentCounter--;
     944  } else
     945    *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
     946
     947  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
     948    // free the index lookup list
     949    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
     950    if (FragmentCounter == 0) // first fragments frees the initial pointer to list
     951      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
     952  }
     953  *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
     954  return status;
    753955};
    754956
    … …  
    759961 * \param &TotalNumberOfKeySets global key set counter
    760962 * \param &TotalGraph Graph to be filled with global numbers
    761  */
    762 void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
    763 {
    764         *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
    765         KeySet *TempSet = new KeySet;
    766         if (FragmentList[FragmentCounter] != NULL) {
    767                 for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
    768                         for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
    769                                 TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
    770                         TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
    771                         TempSet->clear();
    772                 }
    773                 delete(TempSet);
    774         } else {
    775                 *out << Verbose(1) << "FragmentList is NULL." << endl;
    776         }
    777         if (next != NULL)
    778                 next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
    779         FragmentCounter--;
    780         *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
     963 */
     964void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out,
     965    Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets,
     966    Graph &TotalGraph)
     967{
     968  *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
     969  KeySet *TempSet = new KeySet;
     970  if (FragmentList[FragmentCounter] != NULL) {
     971    for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
     972      for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
     973        TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
     974      TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
     975      TempSet->clear();
     976    }
     977    delete (TempSet);
     978  } else {
     979    *out << Verbose(1) << "FragmentList is NULL." << endl;
     980  }
     981  if (next != NULL)
     982    next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
     983  FragmentCounter--;
     984  *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
    781985};
    782986
    … …  
    786990int MoleculeLeafClass::Count() const
    787991{
    788         if (next != NULL)
    789                 return next->Count()+1;
    790         else
    791                 return 1;
    792 };
     992  if (next != NULL)
     993    return next->Count() + 1;
     994  else
     995    return 1;
     996};
     997

  • src/molecules.cpp

    rd8b94a r1907a7  
    6262        cell_size[0] = cell_size[2] = cell_size[5]= 20.;
    6363        cell_size[1] = cell_size[3] = cell_size[4]= 0.;
     64        strcpy(name,"none");
    6465};
    6566
    … …  
    596597        cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
    597598        return false;
     599};
     600
     601/** Set molecule::name from the basename without suffix in the given \a *filename.
     602 * \param *filename filename
     603 */
     604void molecule::SetNameFromFilename(char *filename)
     605{
     606  int length = 0;
     607  char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
     608  char *endname = strrchr(filename, '.');
     609  if ((endname == NULL) || (endname < molname))
     610    length = strlen(molname);
     611  else
     612    length = strlen(molname) - strlen(endname);
     613  strncpy(name, molname, length);
    598614};
    599615
    … …  
    11871203};
    11881204
    1189 /** Removes atom from molecule list.
     1205/** Removes atom from molecule list and deletes it.
    11901206 * \param *pointer atom to be removed
    11911207 * \return true - succeeded, false - atom not found in list
    … …  
    12011217        Trajectories.erase(pointer);
    12021218        return remove(pointer, start, end);
     1219};
     1220
     1221/** Removes atom from molecule list, but does not delete it.
     1222 * \param *pointer atom to be removed
     1223 * \return true - succeeded, false - atom not found in list
     1224 */
     1225bool molecule::UnlinkAtom(atom *pointer)
     1226{
     1227  if (pointer == NULL)
     1228    return false;
     1229  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
     1230    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
     1231  else
     1232    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
     1233  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     1234    ElementCount--;
     1235  Trajectories.erase(pointer);
     1236  unlink(pointer);
     1237  return true;
    12031238};
    12041239
    … …  
    17271762                                                };
    17281763                                                AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
    1729                
     1764
    17301765                                        }
    17311766
    … …  
    25292564                                        cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
    25302565                                }
    2531                                 //FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
    25322566                        }
    25332567                }
    … …  
    30883122        //if (FragmentationToDo) {              // we should always store the fragments again as coordination might have changed slightly without changing bond structure
    30893123                // allocate memory for the pointer array and transmorph graphs into full molecular fragments
    3090                 BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
     3124                BondFragments = new MoleculeListClass();
    30913125                int k=0;
    30923126                for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
    30933127                        KeySet test = (*runner).first;
    30943128                        *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
    3095                         BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
     3129                        BondFragments->insert(StoreFragmentFromKeySet(out, test, configuration));
    30963130                        k++;
    30973131                }
    3098                 *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
     3132                *out << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl;
    30993133
    31003134                // ===== 9. Save fragments' configuration and keyset files et al to disk ===
    3101                 if (BondFragments->NumberOfMolecules != 0) {
     3135                if (BondFragments->ListOfMolecules.size() != 0) {
    31023136                        // create the SortIndex from BFS labels to order in the config file
    31033137                        CreateMappingLabelsToConfigSequence(out, SortIndex);
    31043138
    3105                         *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
     3139                        *out << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl;
    31063140                        if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
    31073141                                *out << Verbose(1) << "All configs written." << endl;
    … …  
    36293663};
    36303664
    3631 /** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
    3632  * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
    3633  * computer game, that winds through the connected graph representing the molecule. Color (white,
    3634  * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
    3635  * creating only unique fragments and not additional ones with vertices simply in different sequence.
    3636  * The Predecessor is always the one that came before in discovering, needed on backstepping. And
    3637  * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
    3638  * stepping.
    3639  * \param *out output stream for debugging
    3640  * \param Order number of atoms in each fragment
    3641  * \param *configuration configuration for writing config files for each fragment
    3642  * \return List of all unique fragments with \a Order atoms
    3643  */
    3644 /*
    3645 MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
    3646 {
    3647         atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
    3648         int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
    3649         int *Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
    3650         enum Shading *ColorVertexList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
    3651         enum Shading *ColorEdgeList = (enum Shading *) Malloc(sizeof(enum Shading)*BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
    3652         StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
    3653         StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
    3654         StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
    3655         MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
    3656         MoleculeListClass *FragmentList = NULL;
    3657         atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
    3658         bond *Binder = NULL;
    3659         int RunningIndex = 0, FragmentCounter = 0;
    3660 
    3661         *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
    3662 
    3663         // reset parent list
    3664         *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
    3665         for (int i=0;i<AtomCount;i++) { // reset all atom labels
    3666                 // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
    3667                 Labels[i] = -1;
    3668                 SonList[i] = NULL;
    3669                 PredecessorList[i] = NULL;
    3670                 ColorVertexList[i] = white;
    3671                 ShortestPathList[i] = -1;
    3672         }
    3673         for (int i=0;i<BondCount;i++)
    3674                 ColorEdgeList[i] = white;
    3675         RootStack->ClearStack();        // clearstack and push first atom if exists
    3676         TouchedStack->ClearStack();
    3677         Walker = start->next;
    3678         while ((Walker != end)
    3679 #ifdef ADDHYDROGEN
    3680          && (Walker->type->Z == 1)
    3681 #endif
    3682                                                                                                 ) { // search for first non-hydrogen atom
    3683                 *out << Verbose(4) << "Current Root candidate is " << Walker->Name << "." << endl;
    3684                 Walker = Walker->next;
    3685         }
    3686         if (Walker != end)
    3687                 RootStack->Push(Walker);
    3688         else
    3689                 *out << Verbose(0) << "ERROR: Could not find an appropriate Root atom!" << endl;
    3690         *out << Verbose(3) << "Root " << Walker->Name << " is on AtomStack, beginning loop through all vertices ..." << endl;
    3691 
    3692         ///// OUTER LOOP ////////////
    3693         while (!RootStack->IsEmpty()) {
    3694                 // get new root vertex from atom stack
    3695                 Root = RootStack->PopFirst();
    3696                 ShortestPathList[Root->nr] = 0;
    3697                 if (Labels[Root->nr] == -1)
    3698                         Labels[Root->nr] = RunningIndex++; // prevent it from getting again on AtomStack
    3699                 PredecessorList[Root->nr] = Root;
    3700                 TouchedStack->Push(Root);
    3701                 *out << Verbose(0) << "Root for this loop is: " << Root->Name << ".\n";
    3702 
    3703                 // clear snake stack
    3704                 SnakeStack->ClearStack();
    3705                 //SnakeStack->TestImplementation(out, start->next);
    3706 
    3707                 ///// INNER LOOP ////////////
    3708                 // Problems:
    3709                 // - what about cyclic bonds?
    3710                 Walker = Root;
    3711                 do {
    3712                         *out << Verbose(1) << "Current Walker is: " << Walker->Name;
    3713                         // initial setting of the new Walker: label, color, shortest path and put on stacks
    3714                         if (Labels[Walker->nr] == -1)   {       // give atom a unique, monotonely increasing number
    3715                                 Labels[Walker->nr] = RunningIndex++;
    3716                                 RootStack->Push(Walker);
    3717                         }
    3718                         *out << ", has label " << Labels[Walker->nr];
    3719                         if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) {     // color it if newly discovered and push on stacks (and if within reach!)
    3720                                 if ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white)) {
    3721                                         // Binder ought to be set still from last neighbour search
    3722                                         *out << ", coloring bond " << *Binder << " black";
    3723                                         ColorEdgeList[Binder->nr] = black; // mark this bond as used
    3724                                 }
    3725                                 if (ShortestPathList[Walker->nr] == -1) {
    3726                                         ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1;
    3727                                         TouchedStack->Push(Walker); // mark every atom for lists cleanup later, whose shortest path has been changed
    3728                                 }
    3729                                 if ((ShortestPathList[Walker->nr] < Order) && (ColorVertexList[Walker->nr] != darkgray)) {      // if not already on snake stack
    3730                                         SnakeStack->Push(Walker);
    3731                                         ColorVertexList[Walker->nr] = darkgray; // mark as dark gray of on snake stack
    3732                                 }
    3733                         }
    3734                         *out << ", SP of " << ShortestPathList[Walker->nr]      << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
    3735 
    3736                         // then check the stack for a newly stumbled upon fragment
    3737                         if (SnakeStack->ItemCount() == Order) { // is stack full?
    3738                                 // store the fragment if it is one and get a removal candidate
    3739                                 Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
    3740                                 // remove the candidate if one was found
    3741                                 if (Removal != NULL) {
    3742                                         *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
    3743                                         SnakeStack->RemoveItem(Removal);
    3744                                         ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
    3745                                         if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
    3746                                                 Walker = PredecessorList[Removal->nr];
    3747                                                 *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
    3748                                         }
    3749                                 }
    3750                         } else
    3751                                 Removal = NULL;
    3752 
    3753                         // finally, look for a white neighbour as the next Walker
    3754                         Binder = NULL;
    3755                         if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {    // don't look, if a new walker has been set above
    3756                                 *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
    3757                                 OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
    3758                                 if (ShortestPathList[Walker->nr] < Order) {
    3759                                         for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
    3760                                                 Binder = ListOfBondsPerAtom[Walker->nr][i];
    3761                                                 *out << Verbose(2) << "Current bond is " << *Binder << ": ";
    3762                                                 OtherAtom = Binder->GetOtherAtom(Walker);
    3763                                                 if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
    3764                                                         *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
    3765                                                         //ColorVertexList[OtherAtom->nr] = lightgray;           // mark as explored
    3766                                                 } else { // otherwise check its colour and element
    3767                                                         if (
    3768 #ifdef ADDHYDROGEN
    3769                                                         (OtherAtom->type->Z != 1) &&
    3770 #endif
    3771                                                                                 (ColorEdgeList[Binder->nr] == white)) { // skip hydrogen, look for unexplored vertices
    3772                                                                 *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
    3773                                                                 // i find it currently rather sensible to always set the predecessor in order to find one's way back
    3774                                                                 //if (PredecessorList[OtherAtom->nr] == NULL) {
    3775                                                                 PredecessorList[OtherAtom->nr] = Walker;
    3776                                                                 *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
    3777                                                                 //} else {
    3778                                                                 //      *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
    3779                                                                 //}
    3780                                                                 Walker = OtherAtom;
    3781                                                                 break;
    3782                                                         } else {
    3783                                                                 if (OtherAtom->type->Z == 1)
    3784                                                                         *out << "Links to a hydrogen atom." << endl;
    3785                                                                 else
    3786                                                                         *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
    3787                                                         }
    3788                                                 }
    3789                                         }
    3790                                 } else {        // means we have stepped beyond the horizon: Return!
    3791                                         Walker = PredecessorList[Walker->nr];
    3792                                         OtherAtom = Walker;
    3793                                         *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
    3794                                 }
    3795                                 if (Walker != OtherAtom) {      // if no white neighbours anymore, color it black
    3796                                         *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
    3797                                         ColorVertexList[Walker->nr] = black;
    3798                                         Walker = PredecessorList[Walker->nr];
    3799                                 }
    3800                         }
    3801                 } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
    3802                 *out << Verbose(2) << "Inner Looping is finished." << endl;
    3803 
    3804                 // if we reset all AtomCount atoms, we have again technically O(N^2) ...
    3805                 *out << Verbose(2) << "Resetting lists." << endl;
    3806                 Walker = NULL;
    3807                 Binder = NULL;
    3808                 while (!TouchedStack->IsEmpty()) {
    3809                         Walker = TouchedStack->PopLast();
    3810                         *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
    3811                         for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
    3812                                 ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
    3813                         PredecessorList[Walker->nr] = NULL;
    3814                         ColorVertexList[Walker->nr] = white;
    3815                         ShortestPathList[Walker->nr] = -1;
    3816                 }
    3817         }
    3818         *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
    3819 
    3820         // copy together
    3821         *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
    3822         FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
    3823         RunningIndex = 0;
    3824         while ((Leaflet != NULL) && (RunningIndex < FragmentCounter))   {
    3825                 FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
    3826                 Leaflet->Leaf = NULL; // prevent molecule from being removed
    3827                 TempLeaf = Leaflet;
    3828                 Leaflet = Leaflet->previous;
    3829                 delete(TempLeaf);
    3830         };
    3831 
    3832         // free memory and exit
    3833         Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
    3834         Free((void **)&ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
    3835         Free((void **)&Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
    3836         Free((void **)&ColorVertexList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
    3837         delete(RootStack);
    3838         delete(TouchedStack);
    3839         delete(SnakeStack);
    3840 
    3841         *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
    3842         return FragmentList;
    3843 };
    3844 */
    3845 
    38463665/** Structure containing all values in power set combination generation.
    38473666 */

  • src/molecules.hpp

    rd8b94a r1907a7  
    6868#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
    6969#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
     70
     71#define MoleculeList list <molecule *>
     72#define MoleculeListTest pair <MoleculeList::iterator, bool>
    7073
    7174/******************************** Some small functions and/or structures **********************************/
    … …  
    142145        ~atom();
    143146
    144         bool Output(int ElementNo, int AtomNo, ofstream *out) const;
     147        bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
    145148        bool OutputXYZLine(ofstream *out) const;
    146149        atom *GetTrueFather();
    … …  
    216219                int NoCyclicBonds;      //!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
    217220                double BondDistance;    //!< typical bond distance used in CreateAdjacencyList() and furtheron
     221                bool ActiveFlag;    //!< in a MoleculeListClass used to discern active from inactive molecules
     222                Vector Center;      //!< Center of molecule in a global box
     223                char name[MAXSTRINGSIZE];         //!< arbitrary name
     224                int IndexNr;        //!< index of molecule in a MoleculeListClass
    218225
    219226        molecule(periodentafel *teil);
    … …  
    223230        bool AddAtom(atom *pointer);
    224231        bool RemoveAtom(atom *pointer);
     232        bool UnlinkAtom(atom *pointer);
    225233        bool CleanupMolecule();
    226234
    … …  
    252260        Vector * DetermineCenterOfGravity(ofstream *out);
    253261        Vector * DetermineCenterOfAll(ofstream *out);
     262        void SetNameFromFilename(char *filename);
    254263        void SetBoxDimension(Vector *dim);
    255264        double * ReturnFullMatrixforSymmetric(double *cell_size);
    … …  
    326335class MoleculeListClass {
    327336        public:
    328                 molecule **ListOfMolecules;      //!< pointer list of fragment molecules to check for equality
    329                 int NumberOfMolecules;                          //!< Number of entries in \a **FragmentList and of to be returned one.
    330                 int NumberOfTopAtoms;                            //!< Number of atoms in the molecule from which all fragments originate
     337          MoleculeList ListOfMolecules; //!< List of the contained molecules
     338          int MaxIndex;
    331339
    332340        MoleculeListClass();
    333         MoleculeListClass(int Num, int NumAtoms);
    334341        ~MoleculeListClass();
    335342
    336         /// Output configs.
    337343        bool AddHydrogenCorrection(ofstream *out, char *path);
    338344        bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
     345        bool insert(molecule *mol);
     346        molecule * ReturnIndex(int index);
    339347        bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
     348        int NumberOfActiveMolecules();
     349        void Enumerate(ofstream *out);
    340350        void Output(ofstream *out);
     351
     352        // merging of molecules
     353  bool SimpleMerge(molecule *mol, molecule *srcmol);
     354  bool SimpleAdd(molecule *mol, molecule *srcmol);
     355  bool SimpleMultiMerge(molecule *mol, int *src, int N);
     356  bool SimpleMultiAdd(molecule *mol, int *src, int N);
     357  bool ScatterMerge(molecule *mol, int *src, int N);
     358  bool EmbedMerge(molecule *mol, molecule *srcmol);
    341359
    342360        private:
    … …  
    454472        bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
    455473        bool SaveMPQC(const char *filename, molecule *mol) const;
    456         void Edit(molecule *mol);
     474        void Edit();
    457475        bool GetIsAngstroem() const;
    458476        char *GetDefaultPath() const;
Note: See TracChangeset for help on using the changeset viewer.