Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_dynamics.cpp

    r1513a74 r35b698  
    66 */
    77
     8#include "Helpers/MemDebug.hpp"
     9
    810#include "World.hpp"
    911#include "atom.hpp"
    1012#include "config.hpp"
    1113#include "element.hpp"
     14#include "info.hpp"
    1215#include "log.hpp"
    1316#include "memoryallocator.hpp"
    … …  
    1518#include "parser.hpp"
    1619#include "Plane.hpp"
     20#include "ThermoStatContainer.hpp"
    1721
    1822/************************************* Functions for class molecule *********************************/
    … …  
    2832  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    2933  gsl_vector *x = gsl_vector_alloc(NDIM);
    30   atom * Runner = mol->start;
    3134  atom *Sprinter = NULL;
    3235  Vector trajectory1, trajectory2, normal, TestVector;
    3336  double Norm1, Norm2, tmp, result = 0.;
    3437
    35   while (Runner->next != mol->end) {
    36     Runner = Runner->next;
    37     if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     38  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     39    if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
    3840      break;
    3941    // determine normalized trajectories direction vector (n1, n2)
    … …  
    4244    trajectory1.Normalize();
    4345    Norm1 = trajectory1.Norm();
    44     Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
    45     trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
     46    Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
     47    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
    4648    trajectory2.Normalize();
    4749    Norm2 = trajectory1.Norm();
    4850    // check whether either is zero()
    4951    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    50       tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
     52      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
    5153    } else if (Norm1 < MYEPSILON) {
    5254      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    53       trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
     55      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
    5456      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    5557      trajectory1 -= trajectory2;   // project the part in norm direction away
    5658      tmp = trajectory1.Norm();  // remaining norm is distance
    5759    } else if (Norm2 < MYEPSILON) {
    58       Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
     60      Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
    5961      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
    6062      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    … …  
    6668  //        Log() << Verbose(0) << " and ";
    6769  //        Log() << Verbose(0) << trajectory2;
    68       tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
     70      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
    6971  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    7072    } else { // determine distance by finding minimum distance
    71   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
     73  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
    7274  //        Log() << Verbose(0) << endl;
    7375  //        Log() << Verbose(0) << "First Trajectory: ";
    … …  
    8587        gsl_matrix_set(A, 1, i, trajectory2[i]);
    8688        gsl_matrix_set(A, 2, i, normal[i]);
    87         gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - Runner->Trajectory.R.at(Params.startstep)[i]));
     89        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
    8890      }
    8991      // solve the linear system by Householder transformations
    … …  
    9698      trajectory2.Scale(gsl_vector_get(x,1));
    9799      normal.Scale(gsl_vector_get(x,2));
    98       TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal
     100      TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
    99101                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
    100102      if (TestVector.Norm() < MYEPSILON) {
    … …  
    125127{
    126128  double result = 0.;
    127   atom * Runner = mol->start;
    128   while (Runner->next != mol->end) {
    129     Runner = Runner->next;
    130     if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
     129  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     130    if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[(*iter)->nr]) && (Walker->nr < (*iter)->nr)) {
    131131  //    atom *Sprinter = PermutationMap[Walker->nr];
    132   //        Log() << Verbose(0) << *Walker << " and " << *Runner << " are heading to the same target at ";
     132  //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
    133133  //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
    134134  //        Log() << Verbose(0) << ", penalting." << endl;
    … …  
    161161  // go through every atom
    162162  atom *Runner = NULL;
    163   atom *Walker = start;
    164   while (Walker->next != end) {
    165     Walker = Walker->next;
     163  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    166164    // first term: distance to target
    167     Runner = Params.PermutationMap[Walker->nr];   // find target point
    168     tmp = (Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
     165    Runner = Params.PermutationMap[(*iter)->nr];   // find target point
     166    tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
    169167    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    170168    result += Params.PenaltyConstants[0] * tmp;
    … …  
    172170
    173171    // second term: sum of distances to other trajectories
    174     result += SumDistanceOfTrajectories(Walker, this, Params);
     172    result += SumDistanceOfTrajectories((*iter), this, Params);
    175173
    176174    // third term: penalty for equal targets
    177     result += PenalizeEqualTargets(Walker, this, Params);
     175    result += PenalizeEqualTargets((*iter), this, Params);
    178176  }
    179177
    … …  
    189187{
    190188  stringstream zeile1, zeile2;
    191   int *DoubleList = Calloc<int>(AtomCount, "PrintPermutationMap: *DoubleList");
     189  int *DoubleList = new int[AtomCount];
     190  for(int i=0;i<AtomCount;i++)
     191    DoubleList[i] = 0;
    192192  int doubles = 0;
    193193  zeile1 << "PermutationMap: ";
    … …  
    203203  if (doubles >0)
    204204    DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
    205   Free(&DoubleList);
     205  delete[](DoubleList);
    206206//  Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl;
    207207};
    … …  
    213213void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
    214214{
    215   for (int i=mol->AtomCount; i--;) {
     215  for (int i=mol->getAtomCount(); i--;) {
    216216    Params.DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    217217    Params.DistanceList[i]->clear();
    218218  }
    219219
    220   atom *Runner = NULL;
    221   atom *Walker = mol->start;
    222   while (Walker->next != mol->end) {
    223     Walker = Walker->next;
    224     Runner = mol->start;
    225     while(Runner->next != mol->end) {
    226       Runner = Runner->next;
    227       Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)), Runner) );
     220  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     221    for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
     222      Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).distance((*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
    228223    }
    229224  }
    … …  
    237232void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
    238233{
    239   atom *Walker = mol->start;
    240   while (Walker->next != mol->end) {
    241     Walker = Walker->next;
    242     Params.StepList[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
    243     Params.PermutationMap[Walker->nr] = Params.DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
    244     Params.DoubleList[Params.DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
    245     Params.DistanceIterators[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // and remember which one we picked
    246     DoLog(2) && (Log() << Verbose(2) << *Walker << " starts with distance " << Params.DistanceList[Walker->nr]->begin()->first << "." << endl);
     234  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     235    Params.StepList[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // stores the step to the next iterator that could be a possible next target
     236    Params.PermutationMap[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin()->second;   // always pick target with the smallest distance
     237    Params.DoubleList[Params.DistanceList[(*iter)->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
     238    Params.DistanceIterators[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // and remember which one we picked
     239    DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << Params.DistanceList[(*iter)->nr]->begin()->first << "." << endl);
    247240  }
    248241};
    … …  
    285278void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
    286279{
    287   atom *Walker = mol->start;
     280  molecule::const_iterator iter = mol->begin();
    288281  DistanceMap::iterator NewBase;
    289282  double Potential = fabs(mol->ConstrainedPotential(Params));
    290283
     284  if (mol->empty()) {
     285    eLog() << Verbose(1) << "Molecule is empty." << endl;
     286    return;
     287  }
    291288  while ((Potential) > Params.PenaltyConstants[2]) {
    292     PrintPermutationMap(mol->AtomCount, Params);
    293     Walker = Walker->next;
    294     if (Walker == mol->end) // round-robin at the end
    295       Walker = mol->start->next;
    296     if (Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     289    PrintPermutationMap(mol->getAtomCount(), Params);
     290    iter++;
     291    if (iter == mol->end()) // round-robin at the end
     292      iter = mol->begin();
     293    if (Params.DoubleList[Params.DistanceIterators[(*iter)->nr]->second->nr] <= 1)  // no need to make those injective that aren't
    297294      continue;
    298295    // now, try finding a new one
    299     Potential = TryNextNearestNeighbourForInjectivePermutation(mol, Walker, Potential, Params);
    300   }
    301   for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1
     296    Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params);
     297  }
     298  for (int i=mol->getAtomCount(); i--;) // now each single entry in the DoubleList should be <=1
    302299    if (Params.DoubleList[i] > 1) {
    303300      DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
    … …  
    338335  double Potential, OldPotential, OlderPotential;
    339336  struct EvaluatePotential Params;
    340   Params.PermutationMap = Calloc<atom*>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.**PermutationMap");
    341   Params.DistanceList = Malloc<DistanceMap*>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.**DistanceList");
    342   Params.DistanceIterators = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*DistanceIterators");
    343   Params.DoubleList = Calloc<int>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*DoubleList");
    344   Params.StepList = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*StepList");
     337  Params.PermutationMap = new atom *[getAtomCount()];
     338  Params.DistanceList = new DistanceMap *[getAtomCount()];
     339  Params.DistanceIterators = new DistanceMap::iterator[getAtomCount()];
     340  Params.DoubleList = new int[getAtomCount()];
     341  Params.StepList = new DistanceMap::iterator[getAtomCount()];
    345342  int round;
    346   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
     343  atom *Sprinter = NULL;
    347344  DistanceMap::iterator Rider, Strider;
     345
     346  // set to zero
     347  for (int i=0;i<getAtomCount();i++) {
     348    Params.PermutationMap[i] = NULL;
     349    Params.DoubleList[i] = 0;
     350  }
    348351
    349352  /// Minimise the potential
    … …  
    362365  DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl);
    363366  MakeInjectivePermutation(this, Params);
    364   Free(&Params.DoubleList);
     367  delete[](Params.DoubleList);
    365368
    366369  // argument minimise the constrained potential in this injective PermutationMap
    … …  
    371374    DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl);
    372375    OlderPotential = OldPotential;
     376    molecule::const_iterator iter;
    373377    do {
    374       Walker = start;
    375       while (Walker->next != end) { // pick one
    376         Walker = Walker->next;
    377         PrintPermutationMap(AtomCount, Params);
    378         Sprinter = Params.DistanceIterators[Walker->nr]->second;   // store initial partner
    379         Strider = Params.DistanceIterators[Walker->nr];  //remember old iterator
    380         Params.DistanceIterators[Walker->nr] = Params.StepList[Walker->nr];
    381         if (Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
    382           Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->begin();
     378      iter = begin();
     379      for (; iter != end(); ++iter) {
     380        PrintPermutationMap(getAtomCount(), Params);
     381        Sprinter = Params.DistanceIterators[(*iter)->nr]->second;   // store initial partner
     382        Strider = Params.DistanceIterators[(*iter)->nr];  //remember old iterator
     383        Params.DistanceIterators[(*iter)->nr] = Params.StepList[(*iter)->nr];
     384        if (Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->end()) {// stop, before we run through the list and still on
     385          Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->begin();
    383386          break;
    384387        }
    385         //Log() << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
     388        //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->nr]->second << "." << endl;
    386389        // find source of the new target
    387         Runner = start->next;
    388         while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    389           if (Params.PermutationMap[Runner->nr] == Params.DistanceIterators[Walker->nr]->second) {
    390             //Log() << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
     390        molecule::const_iterator runner = begin();
     391        for (; runner != end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
     392          if (Params.PermutationMap[(*runner)->nr] == Params.DistanceIterators[(*iter)->nr]->second) {
     393            //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->nr] << "." << endl;
    391394            break;
    392395          }
    393           Runner = Runner->next;
    394396        }
    395         if (Runner != end) { // we found the other source
     397        if (runner != end()) { // we found the other source
    396398          // then look in its distance list for Sprinter
    397           Rider = Params.DistanceList[Runner->nr]->begin();
    398           for (; Rider != Params.DistanceList[Runner->nr]->end(); Rider++)
     399          Rider = Params.DistanceList[(*runner)->nr]->begin();
     400          for (; Rider != Params.DistanceList[(*runner)->nr]->end(); Rider++)
    399401            if (Rider->second == Sprinter)
    400402              break;
    401           if (Rider != Params.DistanceList[Runner->nr]->end()) { // if we have found one
    402             //Log() << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
     403          if (Rider != Params.DistanceList[(*runner)->nr]->end()) { // if we have found one
     404            //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->nr] << "/" << *Rider->second << "." << endl;
    403405            // exchange both
    404             Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
    405             Params.PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
    406             PrintPermutationMap(AtomCount, Params);
     406            Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // put next farther distance into PermutationMap
     407            Params.PermutationMap[(*runner)->nr] = Sprinter;  // and hand the old target to its respective owner
     408            PrintPermutationMap(getAtomCount(), Params);
    407409            // calculate the new potential
    408410            //Log() << Verbose(2) << "Checking new potential ..." << endl;
    … …  
    410412            if (Potential > OldPotential) { // we made everything worse! Undo ...
    411413              //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    412               //Log() << Verbose(3) << "Setting " << *Runner << "'s source to " << *Params.DistanceIterators[Runner->nr]->second << "." << endl;
     414              //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->nr]->second << "." << endl;
    413415              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    414               Params.PermutationMap[Runner->nr] = Params.DistanceIterators[Runner->nr]->second;
     416              Params.PermutationMap[(*runner)->nr] = Params.DistanceIterators[(*runner)->nr]->second;
    415417              // Undo for Walker
    416               Params.DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
    417               //Log() << Verbose(3) << "Setting " << *Walker << "'s source to " << *Params.DistanceIterators[Walker->nr]->second << "." << endl;
    418               Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
     418              Params.DistanceIterators[(*iter)->nr] = Strider;  // take next farther distance target
     419              //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->nr]->second << "." << endl;
     420              Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second;
    419421            } else {
    420               Params.DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
     422              Params.DistanceIterators[(*runner)->nr] = Rider;  // if successful also move the pointer in the iterator list
    421423              DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
    422424              OldPotential = Potential;
    … …  
    428430            //Log() << Verbose(0) << endl;
    429431          } else {
    430             DoeLog(1) && (eLog()<< Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl);
     432            DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl);
    431433            exit(255);
    432434          }
    433435        } else {
    434           Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // new target has no source!
     436          Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // new target has no source!
    435437        }
    436         Params.StepList[Walker->nr]++; // take next farther distance target
     438        Params.StepList[(*iter)->nr]++; // take next farther distance target
    437439      }
    438     } while (Walker->next != end);
     440    } while (++iter != end());
    439441  } while ((OlderPotential - OldPotential) > 1e-3);
    440442  DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    … …  
    442444
    443445  /// free memory and return with evaluated potential
    444   for (int i=AtomCount; i--;)
     446  for (int i=getAtomCount(); i--;)
    445447    Params.DistanceList[i]->clear();
    446   Free(&Params.DistanceList);
    447   Free(&Params.DistanceIterators);
     448  delete[](Params.DistanceList);
     449  delete[](Params.DistanceIterators);
    448450  return ConstrainedPotential(Params);
    449451};
    … …  
    471473 * \param startstep stating initial configuration in molecule::Trajectories
    472474 * \param endstep stating final configuration in molecule::Trajectories
     475 * \param &prefix path and prefix
    473476 * \param &config configuration structure
    474477 * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
    475478 * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
    476479 */
    477 bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
     480bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, std::string &prefix, config &configuration, bool MapByIdentity)
    478481{
    479482  molecule *mol = NULL;
    … …  
    483486  // Get the Permutation Map by MinimiseConstrainedPotential
    484487  atom **PermutationMap = NULL;
    485   atom *Walker = NULL, *Sprinter = NULL;
     488  atom *Sprinter = NULL;
    486489  if (!MapByIdentity)
    487490    MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
    488491  else {
    489     PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
     492    PermutationMap = new atom *[getAtomCount()];
    490493    SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
    491494  }
    … …  
    502505    mol = World::getInstance().createMolecule();
    503506    MoleculePerStep->insert(mol);
    504     Walker = start;
    505     while (Walker->next != end) {
    506       Walker = Walker->next;
     507    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    507508      // add to molecule list
    508       Sprinter = mol->AddCopyAtom(Walker);
     509      Sprinter = mol->AddCopyAtom((*iter));
    509510      for (int n=NDIM;n--;) {
    510         Sprinter->x[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
     511        Sprinter->x[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
    511512        // add to Trajectories
    512513        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    513514        if (step < MaxSteps) {
    514           Walker->Trajectory.R.at(step)[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
    515           Walker->Trajectory.U.at(step)[n] = 0.;
    516           Walker->Trajectory.F.at(step)[n] = 0.;
     515          (*iter)->Trajectory.R.at(step)[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
     516          (*iter)->Trajectory.U.at(step)[n] = 0.;
     517          (*iter)->Trajectory.F.at(step)[n] = 0.;
    517518        }
    518519      }
    … …  
    522523
    523524  // store the list to single step files
    524   int *SortIndex = Malloc<int>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
    525   for (int i=AtomCount; i--; )
     525  int *SortIndex = new int[getAtomCount()];
     526  for (int i=getAtomCount(); i--; )
    526527    SortIndex[i] = i;
    527   status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
     528
     529  status = MoleculePerStep->OutputConfigForListOfFragments(prefix, SortIndex);
     530  delete[](SortIndex);
    528531
    529532  // free and return
    530   Free(&PermutationMap);
     533  delete[](PermutationMap);
    531534  delete(MoleculePerStep);
    532535  return status;
    … …  
    548551bool molecule::VerletForceIntegration(char *file, config &configuration)
    549552{
     553  Info FunctionInfo(__func__);
    550554  ifstream input(file);
    551555  string token;
    … …  
    567571      return false;
    568572    }
    569     if (Force.RowCounter[0] != AtomCount) {
    570       DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl);
     573    if (Force.RowCounter[0] != getAtomCount()) {
     574      DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
    571575      performCriticalExit();
    572576      return false;
    … …  
    574578    // correct Forces
    575579    Velocity.Zero();
    576     for(int i=0;i<AtomCount;i++)
     580    for(int i=0;i<getAtomCount();i++)
    577581      for(int d=0;d<NDIM;d++) {
    578582        Velocity[d] += Force.Matrix[0][i][d+5];
    579583      }
    580     for(int i=0;i<AtomCount;i++)
     584    for(int i=0;i<getAtomCount();i++)
    581585      for(int d=0;d<NDIM;d++) {
    582         Force.Matrix[0][i][d+5] -= Velocity[d]/(double)AtomCount;
     586        Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
    583587      }
    584588    // solve a constrained potential if we are meant to
    … …  
    588592      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    589593      EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force);
    590       Free(&PermutationMap);
     594      delete[](PermutationMap);
    591595    }
    592596
    593597    // and perform Verlet integration for each atom with position, velocity and force vector
    594598    // check size of vectors
    595     ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
    596 
    597     ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps, &configuration, &Force);
     599    //ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
     600
     601    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force);
    598602  }
    599603  // correct velocities (rather momenta) so that center of mass remains motionless
    600604  Velocity.Zero();
    601605  IonMass = 0.;
    602   ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps, &IonMass, &Velocity );
     606  ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps+1, &IonMass, &Velocity );
    603607
    604608  // correct velocities (rather momenta) so that center of mass remains motionless
    605609  Velocity.Scale(1./IonMass);
    606610  ActualTemp = 0.;
    607   ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps, &Velocity );
     611  ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
    608612  Thermostats(configuration, ActualTemp, Berendsen);
    609613  MDSteps++;
    … …  
    642646
    643647  // calculate scale configuration
    644   ScaleTempFactor = configuration.TargetTemp/ActualTemp;
     648  ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
    645649
    646650  // differentating between the various thermostats
    … …  
    650654      break;
    651655     case Woodcock:
    652       if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
     656      if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
    653657        DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
    654658        ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
    … …  
    683687      delta_alpha = 0.;
    684688      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    685       delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
    686       configuration.alpha += delta_alpha*configuration.Deltat;
    687       DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
     689      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
     690      configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
     691      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
    688692      // apply updated alpha as additional force
    689693      ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
Note: See TracChangeset for help on using the changeset viewer.