Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_dynamics.cpp

    r112b09 r1513a74  
    66 */
    77
    8 #include "Helpers/MemDebug.hpp"
    9 
    108#include "World.hpp"
    119#include "atom.hpp"
    1210#include "config.hpp"
    1311#include "element.hpp"
    14 #include "info.hpp"
    1512#include "log.hpp"
    1613#include "memoryallocator.hpp"
    … …  
    3128  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    3229  gsl_vector *x = gsl_vector_alloc(NDIM);
     30  atom * Runner = mol->start;
    3331  atom *Sprinter = NULL;
    3432  Vector trajectory1, trajectory2, normal, TestVector;
    3533  double Norm1, Norm2, tmp, result = 0.;
    3634
    37   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    38     if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     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++)
    3938      break;
    4039    // determine normalized trajectories direction vector (n1, n2)
    … …  
    4342    trajectory1.Normalize();
    4443    Norm1 = trajectory1.Norm();
    45     Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
    46     trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
     44    Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
     45    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
    4746    trajectory2.Normalize();
    4847    Norm2 = trajectory1.Norm();
    4948    // check whether either is zero()
    5049    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    51       tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
     50      tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
    5251    } else if (Norm1 < MYEPSILON) {
    5352      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    54       trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
     53      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
    5554      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    5655      trajectory1 -= trajectory2;   // project the part in norm direction away
    5756      tmp = trajectory1.Norm();  // remaining norm is distance
    5857    } else if (Norm2 < MYEPSILON) {
    59       Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
     58      Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
    6059      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
    6160      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    … …  
    6766  //        Log() << Verbose(0) << " and ";
    6867  //        Log() << Verbose(0) << trajectory2;
    69       tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
     68      tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
    7069  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    7170    } else { // determine distance by finding minimum distance
    72   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
     71  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
    7372  //        Log() << Verbose(0) << endl;
    7473  //        Log() << Verbose(0) << "First Trajectory: ";
    … …  
    8685        gsl_matrix_set(A, 1, i, trajectory2[i]);
    8786        gsl_matrix_set(A, 2, i, normal[i]);
    88         gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
     87        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - Runner->Trajectory.R.at(Params.startstep)[i]));
    8988      }
    9089      // solve the linear system by Householder transformations
    … …  
    9796      trajectory2.Scale(gsl_vector_get(x,1));
    9897      normal.Scale(gsl_vector_get(x,2));
    99       TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
     98      TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal
    10099                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
    101100      if (TestVector.Norm() < MYEPSILON) {
    … …  
    126125{
    127126  double result = 0.;
    128   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    129     if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[(*iter)->nr]) && (Walker->nr < (*iter)->nr)) {
     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)) {
    130131  //    atom *Sprinter = PermutationMap[Walker->nr];
    131   //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
     132  //        Log() << Verbose(0) << *Walker << " and " << *Runner << " are heading to the same target at ";
    132133  //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
    133134  //        Log() << Verbose(0) << ", penalting." << endl;
    … …  
    160161  // go through every atom
    161162  atom *Runner = NULL;
    162   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     163  atom *Walker = start;
     164  while (Walker->next != end) {
     165    Walker = Walker->next;
    163166    // first term: distance to target
    164     Runner = Params.PermutationMap[(*iter)->nr];   // find target point
    165     tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
     167    Runner = Params.PermutationMap[Walker->nr];   // find target point
     168    tmp = (Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
    166169    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    167170    result += Params.PenaltyConstants[0] * tmp;
    … …  
    169172
    170173    // second term: sum of distances to other trajectories
    171     result += SumDistanceOfTrajectories((*iter), this, Params);
     174    result += SumDistanceOfTrajectories(Walker, this, Params);
    172175
    173176    // third term: penalty for equal targets
    174     result += PenalizeEqualTargets((*iter), this, Params);
     177    result += PenalizeEqualTargets(Walker, this, Params);
    175178  }
    176179
    … …  
    186189{
    187190  stringstream zeile1, zeile2;
    188   int *DoubleList = new int[AtomCount];
    189   for(int i=0;i<AtomCount;i++)
    190     DoubleList[i] = 0;
     191  int *DoubleList = Calloc<int>(AtomCount, "PrintPermutationMap: *DoubleList");
    191192  int doubles = 0;
    192193  zeile1 << "PermutationMap: ";
    … …  
    202203  if (doubles >0)
    203204    DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
    204   delete[](DoubleList);
     205  Free(&DoubleList);
    205206//  Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl;
    206207};
    … …  
    212213void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
    213214{
    214   for (int i=mol->getAtomCount(); i--;) {
     215  for (int i=mol->AtomCount; i--;) {
    215216    Params.DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    216217    Params.DistanceList[i]->clear();
    217218  }
    218219
    219   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    220     for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
    221       Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).distance((*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
     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) );
    222228    }
    223229  }
    … …  
    231237void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
    232238{
    233   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    234     Params.StepList[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // stores the step to the next iterator that could be a possible next target
    235     Params.PermutationMap[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin()->second;   // always pick target with the smallest distance
    236     Params.DoubleList[Params.DistanceList[(*iter)->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
    237     Params.DistanceIterators[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // and remember which one we picked
    238     DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << Params.DistanceList[(*iter)->nr]->begin()->first << "." << endl);
     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);
    239247  }
    240248};
    … …  
    277285void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
    278286{
    279   molecule::const_iterator iter = mol->begin();
     287  atom *Walker = mol->start;
    280288  DistanceMap::iterator NewBase;
    281289  double Potential = fabs(mol->ConstrainedPotential(Params));
    282290
    283   if (mol->empty()) {
    284     eLog() << Verbose(1) << "Molecule is empty." << endl;
    285     return;
    286   }
    287291  while ((Potential) > Params.PenaltyConstants[2]) {
    288     PrintPermutationMap(mol->getAtomCount(), Params);
    289     iter++;
    290     if (iter == mol->end()) // round-robin at the end
    291       iter = mol->begin();
    292     if (Params.DoubleList[Params.DistanceIterators[(*iter)->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     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
    293297      continue;
    294298    // now, try finding a new one
    295     Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params);
    296   }
    297   for (int i=mol->getAtomCount(); i--;) // now each single entry in the DoubleList should be <=1
     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
    298302    if (Params.DoubleList[i] > 1) {
    299303      DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
    … …  
    334338  double Potential, OldPotential, OlderPotential;
    335339  struct EvaluatePotential Params;
    336   Params.PermutationMap = new atom *[getAtomCount()];
    337   Params.DistanceList = new DistanceMap *[getAtomCount()];
    338   Params.DistanceIterators = new DistanceMap::iterator[getAtomCount()];
    339   Params.DoubleList = new int[getAtomCount()];
    340   Params.StepList = new DistanceMap::iterator[getAtomCount()];
     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");
    341345  int round;
    342   atom *Sprinter = NULL;
     346  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
    343347  DistanceMap::iterator Rider, Strider;
    344 
    345   // set to zero
    346   for (int i=0;i<getAtomCount();i++) {
    347     Params.PermutationMap[i] = NULL;
    348     Params.DoubleList[i] = 0;
    349   }
    350348
    351349  /// Minimise the potential
    … …  
    364362  DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl);
    365363  MakeInjectivePermutation(this, Params);
    366   delete[](Params.DoubleList);
     364  Free(&Params.DoubleList);
    367365
    368366  // argument minimise the constrained potential in this injective PermutationMap
    … …  
    373371    DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl);
    374372    OlderPotential = OldPotential;
    375     molecule::const_iterator iter;
    376373    do {
    377       iter = begin();
    378       for (; iter != end(); ++iter) {
    379         PrintPermutationMap(getAtomCount(), Params);
    380         Sprinter = Params.DistanceIterators[(*iter)->nr]->second;   // store initial partner
    381         Strider = Params.DistanceIterators[(*iter)->nr];  //remember old iterator
    382         Params.DistanceIterators[(*iter)->nr] = Params.StepList[(*iter)->nr];
    383         if (Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->end()) {// stop, before we run through the list and still on
    384           Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->begin();
     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();
    385383          break;
    386384        }
    387         //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->nr]->second << "." << endl;
     385        //Log() << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
    388386        // find source of the new target
    389         molecule::const_iterator runner = begin();
    390         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)
    391           if (Params.PermutationMap[(*runner)->nr] == Params.DistanceIterators[(*iter)->nr]->second) {
    392             //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->nr] << "." << endl;
     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;
    393391            break;
    394392          }
     393          Runner = Runner->next;
    395394        }
    396         if (runner != end()) { // we found the other source
     395        if (Runner != end) { // we found the other source
    397396          // then look in its distance list for Sprinter
    398           Rider = Params.DistanceList[(*runner)->nr]->begin();
    399           for (; Rider != Params.DistanceList[(*runner)->nr]->end(); Rider++)
     397          Rider = Params.DistanceList[Runner->nr]->begin();
     398          for (; Rider != Params.DistanceList[Runner->nr]->end(); Rider++)
    400399            if (Rider->second == Sprinter)
    401400              break;
    402           if (Rider != Params.DistanceList[(*runner)->nr]->end()) { // if we have found one
    403             //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->nr] << "/" << *Rider->second << "." << endl;
     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;
    404403            // exchange both
    405             Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // put next farther distance into PermutationMap
    406             Params.PermutationMap[(*runner)->nr] = Sprinter;  // and hand the old target to its respective owner
    407             PrintPermutationMap(getAtomCount(), Params);
     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);
    408407            // calculate the new potential
    409408            //Log() << Verbose(2) << "Checking new potential ..." << endl;
    … …  
    411410            if (Potential > OldPotential) { // we made everything worse! Undo ...
    412411              //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    413               //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->nr]->second << "." << endl;
     412              //Log() << Verbose(3) << "Setting " << *Runner << "'s source to " << *Params.DistanceIterators[Runner->nr]->second << "." << endl;
    414413              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    415               Params.PermutationMap[(*runner)->nr] = Params.DistanceIterators[(*runner)->nr]->second;
     414              Params.PermutationMap[Runner->nr] = Params.DistanceIterators[Runner->nr]->second;
    416415              // Undo for Walker
    417               Params.DistanceIterators[(*iter)->nr] = Strider;  // take next farther distance target
    418               //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->nr]->second << "." << endl;
    419               Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second;
     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;
    420419            } else {
    421               Params.DistanceIterators[(*runner)->nr] = Rider;  // if successful also move the pointer in the iterator list
     420              Params.DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
    422421              DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
    423422              OldPotential = Potential;
    … …  
    429428            //Log() << Verbose(0) << endl;
    430429          } else {
    431             DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl);
     430            DoeLog(1) && (eLog()<< Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl);
    432431            exit(255);
    433432          }
    434433        } else {
    435           Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // new target has no source!
     434          Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // new target has no source!
    436435        }
    437         Params.StepList[(*iter)->nr]++; // take next farther distance target
     436        Params.StepList[Walker->nr]++; // take next farther distance target
    438437      }
    439     } while (++iter != end());
     438    } while (Walker->next != end);
    440439  } while ((OlderPotential - OldPotential) > 1e-3);
    441440  DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    … …  
    443442
    444443  /// free memory and return with evaluated potential
    445   for (int i=getAtomCount(); i--;)
     444  for (int i=AtomCount; i--;)
    446445    Params.DistanceList[i]->clear();
    447   delete[](Params.DistanceList);
    448   delete[](Params.DistanceIterators);
     446  Free(&Params.DistanceList);
     447  Free(&Params.DistanceIterators);
    449448  return ConstrainedPotential(Params);
    450449};
    … …  
    484483  // Get the Permutation Map by MinimiseConstrainedPotential
    485484  atom **PermutationMap = NULL;
    486   atom *Sprinter = NULL;
     485  atom *Walker = NULL, *Sprinter = NULL;
    487486  if (!MapByIdentity)
    488487    MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
    489488  else {
    490     PermutationMap = new atom *[getAtomCount()];
     489    PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
    491490    SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
    492491  }
    … …  
    503502    mol = World::getInstance().createMolecule();
    504503    MoleculePerStep->insert(mol);
    505     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     504    Walker = start;
     505    while (Walker->next != end) {
     506      Walker = Walker->next;
    506507      // add to molecule list
    507       Sprinter = mol->AddCopyAtom((*iter));
     508      Sprinter = mol->AddCopyAtom(Walker);
    508509      for (int n=NDIM;n--;) {
    509         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);
     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);
    510511        // add to Trajectories
    511512        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    512513        if (step < MaxSteps) {
    513           (*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);
    514           (*iter)->Trajectory.U.at(step)[n] = 0.;
    515           (*iter)->Trajectory.F.at(step)[n] = 0.;
     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.;
    516517        }
    517518      }
    … …  
    521522
    522523  // store the list to single step files
    523   int *SortIndex = new int[getAtomCount()];
    524   for (int i=getAtomCount(); i--; )
     524  int *SortIndex = Malloc<int>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
     525  for (int i=AtomCount; i--; )
    525526    SortIndex[i] = i;
    526527  status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
    527   delete[](SortIndex);
    528528
    529529  // free and return
    530   delete[](PermutationMap);
     530  Free(&PermutationMap);
    531531  delete(MoleculePerStep);
    532532  return status;
    … …  
    548548bool molecule::VerletForceIntegration(char *file, config &configuration)
    549549{
    550   Info FunctionInfo(__func__);
    551550  ifstream input(file);
    552551  string token;
    … …  
    568567      return false;
    569568    }
    570     if (Force.RowCounter[0] != getAtomCount()) {
    571       DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
     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);
    572571      performCriticalExit();
    573572      return false;
    … …  
    575574    // correct Forces
    576575    Velocity.Zero();
    577     for(int i=0;i<getAtomCount();i++)
     576    for(int i=0;i<AtomCount;i++)
    578577      for(int d=0;d<NDIM;d++) {
    579578        Velocity[d] += Force.Matrix[0][i][d+5];
    580579      }
    581     for(int i=0;i<getAtomCount();i++)
     580    for(int i=0;i<AtomCount;i++)
    582581      for(int d=0;d<NDIM;d++) {
    583         Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
     582        Force.Matrix[0][i][d+5] -= Velocity[d]/(double)AtomCount;
    584583      }
    585584    // solve a constrained potential if we are meant to
    … …  
    589588      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    590589      EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force);
    591       delete[](PermutationMap);
     590      Free(&PermutationMap);
    592591    }
    593592
    594593    // and perform Verlet integration for each atom with position, velocity and force vector
    595594    // check size of vectors
    596     //ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
    597 
    598     ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force);
     595    ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
     596
     597    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps, &configuration, &Force);
    599598  }
    600599  // correct velocities (rather momenta) so that center of mass remains motionless
    601600  Velocity.Zero();
    602601  IonMass = 0.;
    603   ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps+1, &IonMass, &Velocity );
     602  ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps, &IonMass, &Velocity );
    604603
    605604  // correct velocities (rather momenta) so that center of mass remains motionless
    606605  Velocity.Scale(1./IonMass);
    607606  ActualTemp = 0.;
    608   ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
     607  ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps, &Velocity );
    609608  Thermostats(configuration, ActualTemp, Berendsen);
    610609  MDSteps++;
    … …  
    684683      delta_alpha = 0.;
    685684      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    686       delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
     685      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
    687686      configuration.alpha += delta_alpha*configuration.Deltat;
    688687      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
Note: See TracChangeset for help on using the changeset viewer.