Changeset 1ca488 for src/builder.cpp


Ignore:
Timestamp:
Jan 26, 2010, 12:45:41 PM (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:
04b6f9
Parents:
8927ae (diff), 1cf5df (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 commit 'heber/Analysis_PairCorrelation' into MenuRefactoring

Conflicts:

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

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