Changeset b8d1aeb for src/builder.cpp


Ignore:
Timestamp:
Feb 2, 2010, 11:38:06 AM (15 years ago)
Author:
Tillmann Crueger <crueger@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
7ba324
Parents:
9fe36b (diff), 2ededc2 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'MenuRefactoring' into QT4Refactoring

Conflicts:

molecuilder/src/Makefile.am
molecuilder/src/builder.cpp
molecuilder/src/unittests/Makefile.am

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    r9fe36b rb8d1aeb  
    5252using namespace std;
    5353
     54#include <cstring>
     55
    5456#include "analysis_correlation.hpp"
    5557#include "atom.hpp"
    … …  
    7476#include "Actions/MethodAction.hpp"
    7577#include "Actions/small_actions.hpp"
    76 
     78#include "version.h"
     79
     80/********************************************* Subsubmenu routine ************************************/
     81#if 0
     82/** Submenu for adding atoms to the molecule.
     83 * \param *periode periodentafel
     84 * \param *molecule molecules with atoms
     85 */
     86static void AddAtoms(periodentafel *periode, molecule *mol)
     87{
     88  atom *first, *second, *third, *fourth;
     89  Vector **atoms;
     90  Vector x,y,z,n;  // coordinates for absolute point in cell volume
     91  double a,b,c;
     92  char choice;  // menu choice char
     93  bool valid;
     94
     95  Log() << Verbose(0) << "===========ADD ATOM============================" << endl;
     96  Log() << Verbose(0) << " a - state absolute coordinates of atom" << endl;
     97  Log() << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
     98  Log() << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
     99  Log() << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
     100  Log() << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
     101  Log() << Verbose(0) << "all else - go back" << endl;
     102  Log() << Verbose(0) << "===============================================" << endl;
     103  Log() << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
     104  Log() << Verbose(0) << "INPUT: ";
     105  cin >> choice;
     106
     107  switch (choice) {
     108    default:
     109      eLog() << Verbose(2) << "Not a valid choice." << endl;
     110      break;
     111      case 'a': // absolute coordinates of atom
     112        Log() << Verbose(0) << "Enter absolute coordinates." << endl;
     113        first = new atom;
     114        first->x.AskPosition(mol->cell_size, false);
     115        first->type = periode->AskElement();  // give type
     116        mol->AddAtom(first);  // add to molecule
     117        break;
     118
     119      case 'b': // relative coordinates of atom wrt to reference point
     120        first = new atom;
     121        valid = true;
     122        do {
     123          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
     124          Log() << Verbose(0) << "Enter reference coordinates." << endl;
     125          x.AskPosition(mol->cell_size, true);
     126          Log() << Verbose(0) << "Enter relative coordinates." << endl;
     127          first->x.AskPosition(mol->cell_size, false);
     128          first->x.AddVector((const Vector *)&x);
     129          Log() << Verbose(0) << "\n";
     130        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     131        first->type = periode->AskElement();  // give type
     132        mol->AddAtom(first);  // add to molecule
     133        break;
     134
     135      case 'c': // relative coordinates of atom wrt to already placed atom
     136        first = new atom;
     137        valid = true;
     138        do {
     139          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
     140          second = mol->AskAtom("Enter atom number: ");
     141          Log() << Verbose(0) << "Enter relative coordinates." << endl;
     142          first->x.AskPosition(mol->cell_size, false);
     143          for (int i=NDIM;i--;) {
     144            first->x.x[i] += second->x.x[i];
     145          }
     146        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     147        first->type = periode->AskElement();  // give type
     148        mol->AddAtom(first);  // add to molecule
     149        break;
     150
     151    case 'd': // two atoms, two angles and a distance
     152        first = new atom;
     153        valid = true;
     154        do {
     155          if (!valid) {
     156            eLog() << Verbose(2) << "Resulting coordinates out of cell - " << first->x << endl;
     157          }
     158          Log() << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
     159          second = mol->AskAtom("Enter central atom: ");
     160          third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
     161          fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
     162          a = ask_value("Enter distance between central (first) and new atom: ");
     163          b = ask_value("Enter angle between new, first and second atom (degrees): ");
     164          b *= M_PI/180.;
     165          bound(&b, 0., 2.*M_PI);
     166          c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
     167          c *= M_PI/180.;
     168          bound(&c, -M_PI, M_PI);
     169          Log() << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
     170/*
     171          second->Output(1,1,(ofstream *)&cout);
     172          third->Output(1,2,(ofstream *)&cout);
     173          fourth->Output(1,3,(ofstream *)&cout);
     174          n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
     175          x.Copyvector(&second->x);
     176          x.SubtractVector(&third->x);
     177          x.Copyvector(&fourth->x);
     178          x.SubtractVector(&third->x);
     179
     180          if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
     181            Log() << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
     182            continue;
     183          }
     184          Log() << Verbose(0) << "resulting relative coordinates: ";
     185          z.Output();
     186          Log() << Verbose(0) << endl;
     187          */
     188          // calc axis vector
     189          x.CopyVector(&second->x);
     190          x.SubtractVector(&third->x);
     191          x.Normalize();
     192          Log() << Verbose(0) << "x: ",
     193          x.Output();
     194          Log() << Verbose(0) << endl;
     195          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
     196          Log() << Verbose(0) << "z: ",
     197          z.Output();
     198          Log() << Verbose(0) << endl;
     199          y.MakeNormalVector(&x,&z);
     200          Log() << Verbose(0) << "y: ",
     201          y.Output();
     202          Log() << Verbose(0) << endl;
     203
     204          // rotate vector around first angle
     205          first->x.CopyVector(&x);
     206          first->x.RotateVector(&z,b - M_PI);
     207          Log() << Verbose(0) << "Rotated vector: ",
     208          first->x.Output();
     209          Log() << Verbose(0) << endl;
     210          // remove the projection onto the rotation plane of the second angle
     211          n.CopyVector(&y);
     212          n.Scale(first->x.ScalarProduct(&y));
     213          Log() << Verbose(0) << "N1: ",
     214          n.Output();
     215          Log() << Verbose(0) << endl;
     216          first->x.SubtractVector(&n);
     217          Log() << Verbose(0) << "Subtracted vector: ",
     218          first->x.Output();
     219          Log() << Verbose(0) << endl;
     220          n.CopyVector(&z);
     221          n.Scale(first->x.ScalarProduct(&z));
     222          Log() << Verbose(0) << "N2: ",
     223          n.Output();
     224          Log() << Verbose(0) << endl;
     225          first->x.SubtractVector(&n);
     226          Log() << Verbose(0) << "2nd subtracted vector: ",
     227          first->x.Output();
     228          Log() << Verbose(0) << endl;
     229
     230          // rotate another vector around second angle
     231          n.CopyVector(&y);
     232          n.RotateVector(&x,c - M_PI);
     233          Log() << Verbose(0) << "2nd Rotated vector: ",
     234          n.Output();
     235          Log() << Verbose(0) << endl;
     236
     237          // add the two linear independent vectors
     238          first->x.AddVector(&n);
     239          first->x.Normalize();
     240          first->x.Scale(a);
     241          first->x.AddVector(&second->x);
     242
     243          Log() << Verbose(0) << "resulting coordinates: ";
     244          first->x.Output();
     245          Log() << Verbose(0) << endl;
     246        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
     247        first->type = periode->AskElement();  // give type
     248        mol->AddAtom(first);  // add to molecule
     249        break;
     250
     251      case 'e': // least square distance position to a set of atoms
     252        first = new atom;
     253        atoms = new (Vector*[128]);
     254        valid = true;
     255        for(int i=128;i--;)
     256          atoms[i] = NULL;
     257        int i=0, j=0;
     258        Log() << Verbose(0) << "Now we need at least three molecules.\n";
     259        do {
     260          Log() << Verbose(0) << "Enter " << i+1 << "th atom: ";
     261          cin >> j;
     262          if (j != -1) {
     263            second = mol->FindAtom(j);
     264            atoms[i++] = &(second->x);
     265          }
     266        } while ((j != -1) && (i<128));
     267        if (i >= 2) {
     268          first->x.LSQdistance((const Vector **)atoms, i);
     269
     270          first->x.Output();
     271          first->type = periode->AskElement();  // give type
     272          mol->AddAtom(first);  // add to molecule
     273        } else {
     274          delete first;
     275          Log() << Verbose(0) << "Please enter at least two vectors!\n";
     276        }
     277        break;
     278  };
     279};
     280
     281/** Submenu for centering the atoms in the molecule.
     282 * \param *mol molecule with all the atoms
     283 */
     284static void CenterAtoms(molecule *mol)
     285{
     286  Vector x, y, helper;
     287  char choice;  // menu choice char
     288
     289  Log() << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
     290  Log() << Verbose(0) << " a - on origin" << endl;
     291  Log() << Verbose(0) << " b - on center of gravity" << endl;
     292  Log() << Verbose(0) << " c - within box with additional boundary" << endl;
     293  Log() << Verbose(0) << " d - within given simulation box" << endl;
     294  Log() << Verbose(0) << "all else - go back" << endl;
     295  Log() << Verbose(0) << "===============================================" << endl;
     296  Log() << Verbose(0) << "INPUT: ";
     297  cin >> choice;
     298
     299  switch (choice) {
     300    default:
     301      Log() << Verbose(0) << "Not a valid choice." << endl;
     302      break;
     303    case 'a':
     304      Log() << Verbose(0) << "Centering atoms in config file on origin." << endl;
     305      mol->CenterOrigin();
     306      break;
     307    case 'b':
     308      Log() << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
     309      mol->CenterPeriodic();
     310      break;
     311    case 'c':
     312      Log() << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
     313      for (int i=0;i<NDIM;i++) {
     314        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
     315        cin >> y.x[i];
     316      }
     317      mol->CenterEdge(&x);  // make every coordinate positive
     318      mol->Center.AddVector(&y); // translate by boundary
     319      helper.CopyVector(&y);
     320      helper.Scale(2.);
     321      helper.AddVector(&x);
     322      mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
     323      break;
     324    case 'd':
     325      Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     326      for (int i=0;i<NDIM;i++) {
     327        Log() << Verbose(0) << "Enter axis " << i << " boundary: ";
     328        cin >> x.x[i];
     329      }
     330      // update Box of atoms by boundary
     331      mol->SetBoxDimension(&x);
     332      // center
     333      mol->CenterInBox();
     334      break;
     335  }
     336};
     337
     338/** Submenu for aligning the atoms in the molecule.
     339 * \param *periode periodentafel
     340 * \param *mol molecule with all the atoms
     341 */
     342static void AlignAtoms(periodentafel *periode, molecule *mol)
     343{
     344  atom *first, *second, *third;
     345  Vector x,n;
     346  char choice;  // menu choice char
     347
     348  Log() << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
     349  Log() << Verbose(0) << " a - state three atoms defining align plane" << endl;
     350  Log() << Verbose(0) << " b - state alignment vector" << endl;
     351  Log() << Verbose(0) << " c - state two atoms in alignment direction" << endl;
     352  Log() << Verbose(0) << " d - align automatically by least square fit" << endl;
     353  Log() << Verbose(0) << "all else - go back" << endl;
     354  Log() << Verbose(0) << "===============================================" << endl;
     355  Log() << Verbose(0) << "INPUT: ";
     356  cin >> choice;
     357
     358  switch (choice) {
     359    default:
     360    case 'a': // three atoms defining mirror plane
     361      first = mol->AskAtom("Enter first atom: ");
     362      second = mol->AskAtom("Enter second atom: ");
     363      third = mol->AskAtom("Enter third atom: ");
     364
     365      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
     366      break;
     367    case 'b': // normal vector of mirror plane
     368      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
     369      n.AskPosition(mol->cell_size,0);
     370      n.Normalize();
     371      break;
     372    case 'c': // three atoms defining mirror plane
     373      first = mol->AskAtom("Enter first atom: ");
     374      second = mol->AskAtom("Enter second atom: ");
     375
     376      n.CopyVector((const Vector *)&first->x);
     377      n.SubtractVector((const Vector *)&second->x);
     378      n.Normalize();
     379      break;
     380    case 'd':
     381      char shorthand[4];
     382      Vector a;
     383      struct lsq_params param;
     384      do {
     385        fprintf(stdout, "Enter the element of atoms to be chosen: ");
     386        fscanf(stdin, "%3s", shorthand);
     387      } while ((param.type = periode->FindElement(shorthand)) == NULL);
     388      Log() << Verbose(0) << "Element is " << param.type->name << endl;
     389      mol->GetAlignvector(&param);
     390      for (int i=NDIM;i--;) {
     391        x.x[i] = gsl_vector_get(param.x,i);
     392        n.x[i] = gsl_vector_get(param.x,i+NDIM);
     393      }
     394      gsl_vector_free(param.x);
     395      Log() << Verbose(0) << "Offset vector: ";
     396      x.Output();
     397      Log() << Verbose(0) << endl;
     398      n.Normalize();
     399      break;
     400  };
     401  Log() << Verbose(0) << "Alignment vector: ";
     402  n.Output();
     403  Log() << Verbose(0) << endl;
     404  mol->Align(&n);
     405};
     406
     407/** Submenu for mirroring the atoms in the molecule.
     408 * \param *mol molecule with all the atoms
     409 */
     410static void MirrorAtoms(molecule *mol)
     411{
     412  atom *first, *second, *third;
     413  Vector n;
     414  char choice;  // menu choice char
     415
     416  Log() << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
     417  Log() << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
     418  Log() << Verbose(0) << " b - state normal vector of mirror plane" << endl;
     419  Log() << Verbose(0) << " c - state two atoms in normal direction" << endl;
     420  Log() << Verbose(0) << "all else - go back" << endl;
     421  Log() << Verbose(0) << "===============================================" << endl;
     422  Log() << Verbose(0) << "INPUT: ";
     423  cin >> choice;
     424
     425  switch (choice) {
     426    default:
     427    case 'a': // three atoms defining mirror plane
     428      first = mol->AskAtom("Enter first atom: ");
     429      second = mol->AskAtom("Enter second atom: ");
     430      third = mol->AskAtom("Enter third atom: ");
     431
     432      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
     433      break;
     434    case 'b': // normal vector of mirror plane
     435      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
     436      n.AskPosition(mol->cell_size,0);
     437      n.Normalize();
     438      break;
     439    case 'c': // three atoms defining mirror plane
     440      first = mol->AskAtom("Enter first atom: ");
     441      second = mol->AskAtom("Enter second atom: ");
     442
     443      n.CopyVector((const Vector *)&first->x);
     444      n.SubtractVector((const Vector *)&second->x);
     445      n.Normalize();
     446      break;
     447  };
     448  Log() << Verbose(0) << "Normal vector: ";
     449  n.Output();
     450  Log() << Verbose(0) << endl;
     451  mol->Mirror((const Vector *)&n);
     452};
     453>>>>>>> MenuRefactoring:molecuilder/src/builder.cpp
     454
     455/** Submenu for removing the atoms from the molecule.
     456 * \param *mol molecule with all the atoms
     457 */
     458static void RemoveAtoms(molecule *mol)
     459{
     460  atom *first, *second;
     461  int axis;
     462  double tmp1, tmp2;
     463  char choice;  // menu choice char
     464
     465  Log() << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
     466  Log() << Verbose(0) << " a - state atom for removal by number" << endl;
     467  Log() << Verbose(0) << " b - keep only in radius around atom" << endl;
     468  Log() << Verbose(0) << " c - remove this with one axis greater value" << endl;
     469  Log() << Verbose(0) << "all else - go back" << endl;
     470  Log() << Verbose(0) << "===============================================" << endl;
     471  Log() << Verbose(0) << "INPUT: ";
     472  cin >> choice;
     473
     474  switch (choice) {
     475    default:
     476    case 'a':
     477      if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
     478        Log() << Verbose(1) << "Atom removed." << endl;
     479      else
     480        Log() << Verbose(1) << "Atom not found." << endl;
     481      break;
     482    case 'b':
     483      second = mol->AskAtom("Enter number of atom as reference point: ");
     484      Log() << Verbose(0) << "Enter radius: ";
     485      cin >> tmp1;
     486      first = mol->start;
     487      second = first->next;
     488      while(second != mol->end) {
     489        first = second;
     490        second = first->next;
     491        if (first->x.DistanceSquared((const Vector *)&second->x) > tmp1*tmp1) // distance to first above radius ...
     492          mol->RemoveAtom(first);
     493      }
     494      break;
     495    case 'c':
     496      Log() << Verbose(0) << "Which axis is it: ";
     497      cin >> axis;
     498      Log() << Verbose(0) << "Lower boundary: ";
     499      cin >> tmp1;
     500      Log() << Verbose(0) << "Upper boundary: ";
     501      cin >> tmp2;
     502      first = mol->start;
     503      second = first->next;
     504      while(second != mol->end) {
     505        first = second;
     506        second = first->next;
     507        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
     508          //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
     509          mol->RemoveAtom(first);
     510        }
     511      }
     512      break;
     513  };
     514  //mol->Output();
     515  choice = 'r';
     516};
     517
     518/** Submenu for measuring out the atoms in the molecule.
     519 * \param *periode periodentafel
     520 * \param *mol molecule with all the atoms
     521 */
     522static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
     523{
     524  atom *first, *second, *third;
     525  Vector x,y;
     526  double min[256], tmp1, tmp2, tmp3;
     527  int Z;
     528  char choice;  // menu choice char
     529
     530  Log() << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
     531  Log() << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
     532  Log() << Verbose(0) << " b - calculate bond length between two atoms" << endl;
     533  Log() << Verbose(0) << " c - calculate bond angle" << endl;
     534  Log() << Verbose(0) << " d - calculate principal axis of the system" << endl;
     535  Log() << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
     536  Log() << Verbose(0) << " f - calculate temperature from current velocity" << endl;
     537  Log() << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
     538  Log() << Verbose(0) << "all else - go back" << endl;
     539  Log() << Verbose(0) << "===============================================" << endl;
     540  Log() << Verbose(0) << "INPUT: ";
     541  cin >> choice;
     542
     543  switch(choice) {
     544    default:
     545      Log() << Verbose(1) << "Not a valid choice." << endl;
     546      break;
     547    case 'a':
     548      first = mol->AskAtom("Enter first atom: ");
     549      for (int i=MAX_ELEMENTS;i--;)
     550        min[i] = 0.;
     551
     552      second = mol->start;
     553      while ((second->next != mol->end)) {
     554        second = second->next; // advance
     555        Z = second->type->Z;
     556        tmp1 = 0.;
     557        if (first != second) {
     558          x.CopyVector((const Vector *)&first->x);
     559          x.SubtractVector((const Vector *)&second->x);
     560          tmp1 = x.Norm();
     561        }
     562        if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
     563        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
     564      }
     565      for (int i=MAX_ELEMENTS;i--;)
     566        if (min[i] != 0.) Log() << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
     567      break;
     568
     569    case 'b':
     570      first = mol->AskAtom("Enter first atom: ");
     571      second = mol->AskAtom("Enter second atom: ");
     572      for (int i=NDIM;i--;)
     573        min[i] = 0.;
     574      x.CopyVector((const Vector *)&first->x);
     575      x.SubtractVector((const Vector *)&second->x);
     576      tmp1 = x.Norm();
     577      Log() << Verbose(1) << "Distance vector is ";
     578      x.Output();
     579      Log() << Verbose(0) << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
     580      break;
     581
     582    case 'c':
     583      Log() << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
     584      first = mol->AskAtom("Enter first atom: ");
     585      second = mol->AskAtom("Enter central atom: ");
     586      third  = mol->AskAtom("Enter last atom: ");
     587      tmp1 = tmp2 = tmp3 = 0.;
     588      x.CopyVector((const Vector *)&first->x);
     589      x.SubtractVector((const Vector *)&second->x);
     590      y.CopyVector((const Vector *)&third->x);
     591      y.SubtractVector((const Vector *)&second->x);
     592      Log() << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     593      Log() << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
     594      break;
     595    case 'd':
     596      Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;
     597      Log() << Verbose(0) << "Shall we rotate? [0/1]: ";
     598      cin >> Z;
     599      if ((Z >=0) && (Z <=1))
     600        mol->PrincipalAxisSystem((bool)Z);
     601      else
     602        mol->PrincipalAxisSystem(false);
     603      break;
     604    case 'e':
     605      {
     606        Log() << Verbose(0) << "Evaluating volume of the convex envelope.";
     607        class Tesselation *TesselStruct = NULL;
     608        const LinkedCell *LCList = NULL;
     609        LCList = new LinkedCell(mol, 10.);
     610        FindConvexBorder(mol, TesselStruct, LCList, NULL);
     611        double clustervolume = VolumeOfConvexEnvelope(TesselStruct, configuration);
     612        Log() << Verbose(0) << "The tesselated surface area is " << clustervolume << "." << endl;\
     613        delete(LCList);
     614        delete(TesselStruct);
     615      }
     616      break;
     617    case 'f':
     618      mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps);
     619      break;
     620    case 'g':
     621      {
     622        char filename[255];
     623        Log() << Verbose(0) << "Please enter filename: " << endl;
     624        cin >> filename;
     625        Log() << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
     626        ofstream *output = new ofstream(filename, ios::trunc);
     627        if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
     628          Log() << Verbose(2) << "File could not be written." << endl;
     629        else
     630          Log() << Verbose(2) << "File stored." << endl;
     631        output->close();
     632        delete(output);
     633      }
     634      break;
     635  }
     636};
     637
     638/** Submenu for measuring out the atoms in the molecule.
     639 * \param *mol molecule with all the atoms
     640 * \param *configuration configuration structure for the to be written config files of all fragments
     641 */
     642static void FragmentAtoms(molecule *mol, config *configuration)
     643{
     644  int Order1;
     645  clock_t start, end;
     646
     647  Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
     648  Log() << Verbose(0) << "What's the desired bond order: ";
     649  cin >> Order1;
     650  if (mol->first->next != mol->last) {  // there are bonds
     651    start = clock();
     652    mol->FragmentMolecule(Order1, configuration);
     653    end = clock();
     654    Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     655  } else
     656    Log() << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
     657};
     658
     659/********************************************** Submenu routine **************************************/
     660
     661/** Submenu for manipulating atoms.
     662 * \param *periode periodentafel
     663 * \param *molecules list of molecules whose atoms are to be manipulated
     664 */
     665static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
     666{
     667  atom *first, *second;
     668  molecule *mol = NULL;
     669  Vector x,y,z,n; // coordinates for absolute point in cell volume
     670  double *factor; // unit factor if desired
     671  double bond, minBond;
     672  char choice;  // menu choice char
     673  bool valid;
     674
     675  Log() << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl;
     676  Log() << Verbose(0) << "a - add an atom" << endl;
     677  Log() << Verbose(0) << "r - remove an atom" << endl;
     678  Log() << Verbose(0) << "b - scale a bond between atoms" << endl;
     679  Log() << Verbose(0) << "u - change an atoms element" << endl;
     680  Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
     681  Log() << Verbose(0) << "all else - go back" << endl;
     682  Log() << Verbose(0) << "===============================================" << endl;
     683  if (molecules->NumberOfActiveMolecules() > 1)
     684    eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;
     685  Log() << Verbose(0) << "INPUT: ";
     686  cin >> choice;
     687
     688  switch (choice) {
     689    default:
     690      Log() << Verbose(0) << "Not a valid choice." << endl;
     691      break;
     692
     693    case 'a': // add atom
     694      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     695        if ((*ListRunner)->ActiveFlag) {
     696        mol = *ListRunner;
     697        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     698        AddAtoms(periode, mol);
     699      }
     700      break;
     701
     702    case 'b': // scale a bond
     703      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     704        if ((*ListRunner)->ActiveFlag) {
     705        mol = *ListRunner;
     706        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     707        Log() << Verbose(0) << "Scaling bond length between two atoms." << endl;
     708        first = mol->AskAtom("Enter first (fixed) atom: ");
     709        second = mol->AskAtom("Enter second (shifting) atom: ");
     710        minBond = 0.;
     711        for (int i=NDIM;i--;)
     712          minBond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
     713        minBond = sqrt(minBond);
     714        Log() << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << minBond << " a.u." << endl;
     715        Log() << Verbose(0) << "Enter new bond length [a.u.]: ";
     716        cin >> bond;
     717        for (int i=NDIM;i--;) {
     718          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/minBond*(minBond-bond);
     719        }
     720        //Log() << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
     721        //second->Output(second->type->No, 1);
     722      }
     723      break;
     724
     725    case 'c': // unit scaling of the metric
     726      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     727        if ((*ListRunner)->ActiveFlag) {
     728        mol = *ListRunner;
     729        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     730       Log() << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
     731       Log() << Verbose(0) << "Enter three factors: ";
     732       factor = new double[NDIM];
     733       cin >> factor[0];
     734       cin >> factor[1];
     735       cin >> factor[2];
     736       valid = true;
     737       mol->Scale((const double ** const)&factor);
     738       delete[](factor);
     739      }
     740     break;
     741
     742    case 'l': // measure distances or angles
     743      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     744        if ((*ListRunner)->ActiveFlag) {
     745        mol = *ListRunner;
     746        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     747        MeasureAtoms(periode, mol, configuration);
     748      }
     749      break;
     750
     751    case 'r': // remove atom
     752      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     753        if ((*ListRunner)->ActiveFlag) {
     754        mol = *ListRunner;
     755        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     756        RemoveAtoms(mol);
     757      }
     758      break;
     759
     760    case 'u': // change an atom's element
     761      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     762        if ((*ListRunner)->ActiveFlag) {
     763        int Z;
     764        mol = *ListRunner;
     765        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     766        first = NULL;
     767        do {
     768          Log() << Verbose(0) << "Change the element of which atom: ";
     769          cin >> Z;
     770        } while ((first = mol->FindAtom(Z)) == NULL);
     771        Log() << Verbose(0) << "New element by atomic number Z: ";
     772        cin >> Z;
     773        first->type = periode->FindElement(Z);
     774        Log() << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
     775      }
     776      break;
     777  }
     778};
     779
     780/** Submenu for manipulating molecules.
     781 * \param *periode periodentafel
     782 * \param *molecules list of molecule to manipulate
     783 */
     784static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
     785{
     786  atom *first = NULL;
     787  Vector x,y,z,n; // coordinates for absolute point in cell volume
     788  int j, axis, count, faktor;
     789  char choice;  // menu choice char
     790  molecule *mol = NULL;
     791  element **Elements;
     792  Vector **vectors;
     793  MoleculeLeafClass *Subgraphs = NULL;
     794
     795  Log() << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;
     796  Log() << Verbose(0) << "c - scale by unit transformation" << endl;
     797  Log() << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
     798  Log() << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
     799  Log() << Verbose(0) << "g - center atoms in box" << endl;
     800  Log() << Verbose(0) << "i - realign molecule" << endl;
     801  Log() << Verbose(0) << "m - mirror all molecules" << endl;
     802  Log() << Verbose(0) << "o - create connection matrix" << endl;
     803  Log() << Verbose(0) << "t - translate molecule by vector" << endl;
     804  Log() << Verbose(0) << "all else - go back" << endl;
     805  Log() << Verbose(0) << "===============================================" << endl;
     806  if (molecules->NumberOfActiveMolecules() > 1)
     807    eLog() << Verbose(2) << "There is more than one molecule active! Atoms will be added to each." << endl;
     808  Log() << Verbose(0) << "INPUT: ";
     809  cin >> choice;
     810
     811  switch (choice) {
     812    default:
     813      Log() << Verbose(0) << "Not a valid choice." << endl;
     814      break;
     815
     816    case 'd': // duplicate the periodic cell along a given axis, given times
     817      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     818        if ((*ListRunner)->ActiveFlag) {
     819        mol = *ListRunner;
     820        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     821        Log() << Verbose(0) << "State the axis [(+-)123]: ";
     822        cin >> axis;
     823        Log() << Verbose(0) << "State the factor: ";
     824        cin >> faktor;
     825
     826        mol->CountAtoms(); // recount atoms
     827        if (mol->AtomCount != 0) {  // if there is more than none
     828          count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
     829          Elements = new element *[count];
     830          vectors = new Vector *[count];
     831          j = 0;
     832          first = mol->start;
     833          while (first->next != mol->end) { // make a list of all atoms with coordinates and element
     834            first = first->next;
     835            Elements[j] = first->type;
     836            vectors[j] = &first->x;
     837            j++;
     838          }
     839          if (count != j)
     840            eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
     841          x.Zero();
     842          y.Zero();
     843          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
     844          for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
     845            x.AddVector(&y); // per factor one cell width further
     846            for (int k=count;k--;) { // go through every atom of the original cell
     847              first = new atom(); // create a new body
     848              first->x.CopyVector(vectors[k]);  // use coordinate of original atom
     849              first->x.AddVector(&x);     // translate the coordinates
     850              first->type = Elements[k];  // insert original element
     851              mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
     852            }
     853          }
     854          if (mol->first->next != mol->last) // if connect matrix is present already, redo it
     855            mol->CreateAdjacencyList(mol->BondDistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
     856          // free memory
     857          delete[](Elements);
     858          delete[](vectors);
     859          // correct cell size
     860          if (axis < 0) { // if sign was negative, we have to translate everything
     861            x.Zero();
     862            x.AddVector(&y);
     863            x.Scale(-(faktor-1));
     864            mol->Translate(&x);
     865          }
     866          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     867        }
     868      }
     869      break;
     870
     871    case 'f':
     872      FragmentAtoms(mol, configuration);
     873      break;
     874
     875    case 'g': // center the atoms
     876      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     877        if ((*ListRunner)->ActiveFlag) {
     878        mol = *ListRunner;
     879        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     880        CenterAtoms(mol);
     881      }
     882      break;
     883
     884    case 'i': // align all atoms
     885      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     886        if ((*ListRunner)->ActiveFlag) {
     887        mol = *ListRunner;
     888        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     889        AlignAtoms(periode, mol);
     890      }
     891      break;
     892
     893    case 'm': // mirror atoms along a given axis
     894      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     895        if ((*ListRunner)->ActiveFlag) {
     896        mol = *ListRunner;
     897        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     898        MirrorAtoms(mol);
     899      }
     900      break;
     901
     902    case 'o': // create the connection matrix
     903      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     904        if ((*ListRunner)->ActiveFlag) {
     905          mol = *ListRunner;
     906          double bonddistance;
     907          clock_t start,end;
     908          Log() << Verbose(0) << "What's the maximum bond distance: ";
     909          cin >> bonddistance;
     910          start = clock();
     911          mol->CreateAdjacencyList(bonddistance, configuration->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
     912          end = clock();
     913          Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
     914        }
     915      break;
     916
     917    case 't': // translate all atoms
     918      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     919        if ((*ListRunner)->ActiveFlag) {
     920        mol = *ListRunner;
     921        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
     922        Log() << Verbose(0) << "Enter translation vector." << endl;
     923        x.AskPosition(mol->cell_size,0);
     924        mol->Center.AddVector((const Vector *)&x);
     925     }
     926     break;
     927  }
     928  // Free all
     929  if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
     930    while (Subgraphs->next != NULL) {
     931      Subgraphs = Subgraphs->next;
     932      delete(Subgraphs->previous);
     933    }
     934    delete(Subgraphs);
     935  }
     936};
     937
     938
     939/** Submenu for creating new molecules.
     940 * \param *periode periodentafel
     941 * \param *molecules list of molecules to add to
     942 */
     943static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
     944{
     945  char choice;  // menu choice char
     946  Vector center;
     947  int nr, count;
     948  molecule *mol = NULL;
     949
     950  Log() << Verbose(0) << "==========EDIT MOLECULES=====================" << endl;
     951  Log() << Verbose(0) << "c - create new molecule" << endl;
     952  Log() << Verbose(0) << "l - load molecule from xyz file" << endl;
     953  Log() << Verbose(0) << "n - change molecule's name" << endl;
     954  Log() << Verbose(0) << "N - give molecules filename" << endl;
     955  Log() << Verbose(0) << "p - parse atoms in xyz file into molecule" << endl;
     956  Log() << Verbose(0) << "r - remove a molecule" << endl;
     957  Log() << Verbose(0) << "all else - go back" << endl;
     958  Log() << Verbose(0) << "===============================================" << endl;
     959  Log() << Verbose(0) << "INPUT: ";
     960  cin >> choice;
     961
     962  switch (choice) {
     963    default:
     964      Log() << Verbose(0) << "Not a valid choice." << endl;
     965      break;
     966    case 'c':
     967      mol = new molecule(periode);
     968      molecules->insert(mol);
     969      break;
     970
     971    case 'l': // load from XYZ file
     972      {
     973        char filename[MAXSTRINGSIZE];
     974        Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     975        mol = new molecule(periode);
     976        do {
     977          Log() << Verbose(0) << "Enter file name: ";
     978          cin >> filename;
     979        } while (!mol->AddXYZFile(filename));
     980        mol->SetNameFromFilename(filename);
     981        // center at set box dimensions
     982        mol->CenterEdge(&center);
     983        mol->cell_size[0] = center.x[0];
     984        mol->cell_size[1] = 0;
     985        mol->cell_size[2] = center.x[1];
     986        mol->cell_size[3] = 0;
     987        mol->cell_size[4] = 0;
     988        mol->cell_size[5] = center.x[2];
     989        molecules->insert(mol);
     990      }
     991      break;
     992
     993    case 'n':
     994      {
     995        char filename[MAXSTRINGSIZE];
     996        do {
     997          Log() << Verbose(0) << "Enter index of molecule: ";
     998          cin >> nr;
     999          mol = molecules->ReturnIndex(nr);
     1000        } while (mol == NULL);
     1001        Log() << Verbose(0) << "Enter name: ";
     1002        cin >> filename;
     1003        strcpy(mol->name, filename);
     1004      }
     1005      break;
     1006
     1007    case 'N':
     1008      {
     1009        char filename[MAXSTRINGSIZE];
     1010        do {
     1011          Log() << Verbose(0) << "Enter index of molecule: ";
     1012          cin >> nr;
     1013          mol = molecules->ReturnIndex(nr);
     1014        } while (mol == NULL);
     1015        Log() << Verbose(0) << "Enter name: ";
     1016        cin >> filename;
     1017        mol->SetNameFromFilename(filename);
     1018      }
     1019      break;
     1020
     1021    case 'p': // parse XYZ file
     1022      {
     1023        char filename[MAXSTRINGSIZE];
     1024        mol = NULL;
     1025        do {
     1026          Log() << Verbose(0) << "Enter index of molecule: ";
     1027          cin >> nr;
     1028          mol = molecules->ReturnIndex(nr);
     1029        } while (mol == NULL);
     1030        Log() << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
     1031        do {
     1032          Log() << Verbose(0) << "Enter file name: ";
     1033          cin >> filename;
     1034        } while (!mol->AddXYZFile(filename));
     1035        mol->SetNameFromFilename(filename);
     1036      }
     1037      break;
     1038
     1039    case 'r':
     1040      Log() << Verbose(0) << "Enter index of molecule: ";
     1041      cin >> nr;
     1042      count = 1;
     1043      for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1044        if (nr == (*ListRunner)->IndexNr) {
     1045          mol = *ListRunner;
     1046          molecules->ListOfMolecules.erase(ListRunner);
     1047          delete(mol);
     1048          break;
     1049        }
     1050      break;
     1051  }
     1052};
     1053
     1054
     1055/** Submenu for merging molecules.
     1056 * \param *periode periodentafel
     1057 * \param *molecules list of molecules to add to
     1058 */
     1059static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
     1060{
     1061  char choice;  // menu choice char
     1062
     1063  Log() << Verbose(0) << "===========MERGE MOLECULES=====================" << endl;
     1064  Log() << Verbose(0) << "a - simple add of one molecule to another" << endl;
     1065  Log() << Verbose(0) << "e - embedding merge of two molecules" << endl;
     1066  Log() << Verbose(0) << "m - multi-merge of all molecules" << endl;
     1067  Log() << Verbose(0) << "s - scatter merge of two molecules" << endl;
     1068  Log() << Verbose(0) << "t - simple merge of two molecules" << endl;
     1069  Log() << Verbose(0) << "all else - go back" << endl;
     1070  Log() << Verbose(0) << "===============================================" << endl;
     1071  Log() << Verbose(0) << "INPUT: ";
     1072  cin >> choice;
     1073
     1074  switch (choice) {
     1075    default:
     1076      Log() << Verbose(0) << "Not a valid choice." << endl;
     1077      break;
     1078
     1079    case 'a':
     1080      {
     1081        int src, dest;
     1082        molecule *srcmol = NULL, *destmol = NULL;
     1083        {
     1084          do {
     1085            Log() << Verbose(0) << "Enter index of destination molecule: ";
     1086            cin >> dest;
     1087            destmol = molecules->ReturnIndex(dest);
     1088          } while ((destmol == NULL) && (dest != -1));
     1089          do {
     1090            Log() << Verbose(0) << "Enter index of source molecule to add from: ";
     1091            cin >> src;
     1092            srcmol = molecules->ReturnIndex(src);
     1093          } while ((srcmol == NULL) && (src != -1));
     1094          if ((src != -1) && (dest != -1))
     1095            molecules->SimpleAdd(srcmol, destmol);
     1096        }
     1097      }
     1098      break;
     1099
     1100    case 'e':
     1101      {
     1102        int src, dest;
     1103        molecule *srcmol = NULL, *destmol = NULL;
     1104        do {
     1105          Log() << Verbose(0) << "Enter index of matrix molecule (the variable one): ";
     1106          cin >> src;
     1107          srcmol = molecules->ReturnIndex(src);
     1108        } while ((srcmol == NULL) && (src != -1));
     1109        do {
     1110          Log() << Verbose(0) << "Enter index of molecule to merge into (the fixed one): ";
     1111          cin >> dest;
     1112          destmol = molecules->ReturnIndex(dest);
     1113        } while ((destmol == NULL) && (dest != -1));
     1114        if ((src != -1) && (dest != -1))
     1115          molecules->EmbedMerge(destmol, srcmol);
     1116      }
     1117      break;
     1118
     1119    case 'm':
     1120      {
     1121        int nr;
     1122        molecule *mol = NULL;
     1123        do {
     1124          Log() << Verbose(0) << "Enter index of molecule to merge into: ";
     1125          cin >> nr;
     1126          mol = molecules->ReturnIndex(nr);
     1127        } while ((mol == NULL) && (nr != -1));
     1128        if (nr != -1) {
     1129          int N = molecules->ListOfMolecules.size()-1;
     1130          int *src = new int(N);
     1131          for(MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1132            if ((*ListRunner)->IndexNr != nr)
     1133              src[N++] = (*ListRunner)->IndexNr;       
     1134          molecules->SimpleMultiMerge(mol, src, N);
     1135          delete[](src);
     1136        }
     1137      }
     1138      break;
     1139
     1140    case 's':
     1141      Log() << Verbose(0) << "Not implemented yet." << endl;
     1142      break;
     1143
     1144    case 't':
     1145      {
     1146        int src, dest;
     1147        molecule *srcmol = NULL, *destmol = NULL;
     1148        {
     1149          do {
     1150            Log() << Verbose(0) << "Enter index of destination molecule: ";
     1151            cin >> dest;
     1152            destmol = molecules->ReturnIndex(dest);
     1153          } while ((destmol == NULL) && (dest != -1));
     1154          do {
     1155            Log() << Verbose(0) << "Enter index of source molecule to merge into: ";
     1156            cin >> src;
     1157            srcmol = molecules->ReturnIndex(src);
     1158          } while ((srcmol == NULL) && (src != -1));
     1159          if ((src != -1) && (dest != -1))
     1160            molecules->SimpleMerge(srcmol, destmol);
     1161        }
     1162      }
     1163      break;
     1164  }
     1165};
     1166
     1167/********************************************** Test routine **************************************/
     1168
     1169/** Is called always as option 'T' in the menu.
     1170 * \param *molecules list of molecules
     1171 */
     1172static void testroutine(MoleculeListClass *molecules)
     1173{
     1174  // the current test routine checks the functionality of the KeySet&Graph concept:
     1175  // We want to have a multiindex (the KeySet) describing a unique subgraph
     1176  int i, comp, counter=0;
     1177
     1178  // create a clone
     1179  molecule *mol = NULL;
     1180  if (molecules->ListOfMolecules.size() != 0) // clone
     1181    mol = (molecules->ListOfMolecules.front())->CopyMolecule();
     1182  else {
     1183    eLog() << Verbose(0) << "I don't have anything to test on ... ";
     1184    performCriticalExit();
     1185    return;
     1186  }
     1187  atom *Walker = mol->start;
     1188
     1189  // generate some KeySets
     1190  Log() << Verbose(0) << "Generating KeySets." << endl;
     1191  KeySet TestSets[mol->AtomCount+1];
     1192  i=1;
     1193  while (Walker->next != mol->end) {
     1194    Walker = Walker->next;
     1195    for (int j=0;j<i;j++) {
     1196      TestSets[j].insert(Walker->nr);
     1197    }
     1198    i++;
     1199  }
     1200  Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
     1201  KeySetTestPair test;
     1202  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
     1203  if (test.second) {
     1204    Log() << Verbose(1) << "Insertion worked?!" << endl;
     1205  } else {
     1206    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
     1207  }
     1208  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
     1209  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
     1210
     1211  // constructing Graph structure
     1212  Log() << Verbose(0) << "Generating Subgraph class." << endl;
     1213  Graph Subgraphs;
     1214
     1215  // insert KeySets into Subgraphs
     1216  Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl;
     1217  for (int j=0;j<mol->AtomCount;j++) {
     1218    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
     1219  }
     1220  Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl;
     1221  GraphTestPair test2;
     1222  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
     1223  if (test2.second) {
     1224    Log() << Verbose(1) << "Insertion worked?!" << endl;
     1225  } else {
     1226    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
     1227  }
     1228
     1229  // show graphs
     1230  Log() << Verbose(0) << "Showing Subgraph's contents, checking that it's sorted." << endl;
     1231  Graph::iterator A = Subgraphs.begin();
     1232  while (A !=  Subgraphs.end()) {
     1233    Log() << Verbose(0) << (*A).second.first << ": ";
     1234    KeySet::iterator key = (*A).first.begin();
     1235    comp = -1;
     1236    while (key != (*A).first.end()) {
     1237      if ((*key) > comp)
     1238        Log() << Verbose(0) << (*key) << " ";
     1239      else
     1240        Log() << Verbose(0) << (*key) << "! ";
     1241      comp = (*key);
     1242      key++;
     1243    }
     1244    Log() << Verbose(0) << endl;
     1245    A++;
     1246  }
     1247  delete(mol);
     1248};
     1249
     1250#endif
    771251
    781252/** Parses the command line options.
    … …  
    1031277  int argptr;
    1041278  molecule *mol = NULL;
    105   string BondGraphFileName("");
     1279  string BondGraphFileName("\n");
    1061280  int verbosity = 0;
    1071281  strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);
    … …  
    1261300            Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
    1271301            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    128             Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
     1302            Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
    1291303            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    1301304            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    1311305            Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    1321306            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
    133             Log() << Verbose(0) << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
     1307            Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
     1308            Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
    1341309            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
    1351310            Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
     1311            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
    1361312            Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
    1371313            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    2571433       mol = new molecule(periode);
    2581434       mol->ActiveFlag = true;
     1435       if (ConfigFileName != NULL)
     1436         mol->SetNameFromFilename(ConfigFileName);
    2591437       molecules->insert(mol);
     1438     }
     1439     if (configuration.BG == NULL) {
     1440       configuration.BG = new BondGraph(configuration.GetIsAngstroem());
     1441       if ((!BondGraphFileName.empty()) && (configuration.BG->LoadBondLengthTable(BondGraphFileName))) {
     1442         Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
     1443       } else {
     1444         eLog() << Verbose(1) << "Bond length table loading failed." << endl;
     1445       }
    2601446     }
    2611447
    … …  
    3591545              //argptr+=1;
    3601546              break;
     1547            case 'I':
     1548              Log() << Verbose(1) << "Dissecting molecular system into a set of disconnected subgraphs ... " << endl;
     1549              // @TODO rather do the dissection afterwards
     1550              molecules->DissectMoleculeIntoConnectedSubgraphs(periode, &configuration);
     1551              mol = NULL;
     1552              if (molecules->ListOfMolecules.size() != 0) {
     1553                for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     1554                  if ((*ListRunner)->ActiveFlag) {
     1555                    mol = *ListRunner;
     1556                    break;
     1557                  }
     1558              }
     1559              if (mol == NULL) {
     1560                mol = *(molecules->ListOfMolecules.begin());
     1561                mol->ActiveFlag = true;
     1562              }
     1563              break;
    3611564            case 'C':
    3621565              if (ExitFlag == 0) ExitFlag = 1;
    … …  
    3661569                performCriticalExit();
    3671570              } else {
    368                 SaveFlag = false;
    3691571                ofstream output(argv[argptr+1]);
    3701572                ofstream binoutput(argv[argptr+2]);
    … …  
    3861588                counter = 0;
    3871589                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    388                   Actives[counter] = (*BigFinder)->ActiveFlag;
     1590                  Actives[counter++] = (*BigFinder)->ActiveFlag;
    3891591                  (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
    3901592                }
    … …  
    3941596                int ranges[NDIM] = {1,1,1};
    3951597                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
     1598                //OutputCorrelationToSurface(&output, surfacemap);
    3961599                BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
    3971600                OutputCorrelation ( &binoutput, binmap );
    … …  
    3991602                binoutput.close();
    4001603                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
    401                   (*BigFinder)->ActiveFlag = Actives[counter];
     1604                  (*BigFinder)->ActiveFlag = Actives[counter++];
    4021605                Free(&Actives);
    4031606                delete(LCList);
    … …  
    4221625            case 'F':
    4231626              if (ExitFlag == 0) ExitFlag = 1;
    424               if (argptr+5 >=argc) {
     1627              if (argptr+6 >=argc) {
    4251628                ExitFlag = 255;
    426                 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
     1629                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
    4271630                performCriticalExit();
    4281631              } else {
    … …  
    4301633                Log() << Verbose(1) << "Filling Box with water molecules." << endl;
    4311634                // construct water molecule
    432                 molecule *filler = new molecule(periode);;
     1635                molecule *filler = new molecule(periode);
    4331636                molecule *Filling = NULL;
    4341637                atom *second = NULL, *third = NULL;
     1638//                first = new atom();
     1639//                first->type = periode->FindElement(5);
     1640//                first->x.Zero();
     1641//                filler->AddAtom(first);
    4351642                first = new atom();
    4361643                first->type = periode->FindElement(1);
    … …  
    4511658                for (int i=0;i<NDIM;i++)
    4521659                  distance[i] = atof(argv[argptr+i]);
    453                 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
     1660                Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
    4541661                if (Filling != NULL) {
     1662                  Filling->ActiveFlag = false;
    4551663                  molecules->insert(Filling);
    4561664                }
    … …  
    4991707                start = clock();
    5001708                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
    501                 FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]);
     1709                if (!FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]))
     1710                  ExitFlag = 255;
    5021711                //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
    5031712                end = clock();
    5041713                Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
    5051714                delete(LCList);
     1715                delete(T);
    5061716                argptr+=2;
    5071717              }
    … …  
    7992009                if (volume != -1)
    8002010                  ExitFlag = 255;
    801                   eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
     2011                  eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension: -u <density>" << endl;
    8022012                  performCriticalExit();
    8032013              } else {
    … …  
    9202130  Action *eraseMoleculeAction = new MethodAction("eraseMoleculeAction",boost::bind(&MoleculeListClass::eraseMolecule,molecules));
    9212131  new ActionMenuItem('r',"remove a molecule",editMoleculesMenu,eraseMoleculeAction);
     2132
    9222133}
    9232134
    … …  
    9712182
    9722183    {
     2184      cout << ESPACKVersion << endl;
     2185
     2186      setVerbosity(0);
     2187
    9732188      menuPopulaters populaters;
    9742189      populaters.MakeEditMoleculesMenu = populateEditMoleculesMenu;
    … …  
    9812196      MainWindow *mainWindow = UIFactory::get()->makeMainWindow(populaters,molecules, configuration, periode, ConfigFileName);
    9822197      mainWindow->display();
     2198
    9832199      delete mainWindow;
    9842200    }
Note: See TracChangeset for help on using the changeset viewer.