Ignore:
Timestamp:
Apr 2, 2009, 4:12:54 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
3021d93
Parents:
451d7a
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 ..."
File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/builder.cpp

    r451d7a rc830e8e  
    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);
Note: See TracChangeset for help on using the changeset viewer.