Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_dynamics.cpp

    r920c70 r112b09  
    55 *      Author: heber
    66 */
     7
     8#include "Helpers/MemDebug.hpp"
    79
    810#include "World.hpp"
    … …  
    2931  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    3032  gsl_vector *x = gsl_vector_alloc(NDIM);
    31   atom * Runner = mol->start;
    3233  atom *Sprinter = NULL;
    3334  Vector trajectory1, trajectory2, normal, TestVector;
    3435  double Norm1, Norm2, tmp, result = 0.;
    3536
    36   while (Runner->next != mol->end) {
    37     Runner = Runner->next;
    38     if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     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++)
    3939      break;
    4040    // determine normalized trajectories direction vector (n1, n2)
    … …  
    4343    trajectory1.Normalize();
    4444    Norm1 = trajectory1.Norm();
    45     Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
    46     trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
     45    Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
     46    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
    4747    trajectory2.Normalize();
    4848    Norm2 = trajectory1.Norm();
    4949    // check whether either is zero()
    5050    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    51       tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
     51      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
    5252    } else if (Norm1 < MYEPSILON) {
    5353      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    54       trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
     54      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
    5555      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    5656      trajectory1 -= trajectory2;   // project the part in norm direction away
    5757      tmp = trajectory1.Norm();  // remaining norm is distance
    5858    } else if (Norm2 < MYEPSILON) {
    59       Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
     59      Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
    6060      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
    6161      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    … …  
    6767  //        Log() << Verbose(0) << " and ";
    6868  //        Log() << Verbose(0) << trajectory2;
    69       tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
     69      tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
    7070  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    7171    } else { // determine distance by finding minimum distance
    72   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
     72  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
    7373  //        Log() << Verbose(0) << endl;
    7474  //        Log() << Verbose(0) << "First Trajectory: ";
    … …  
    8686        gsl_matrix_set(A, 1, i, trajectory2[i]);
    8787        gsl_matrix_set(A, 2, i, normal[i]);
    88         gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - Runner->Trajectory.R.at(Params.startstep)[i]));
     88        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
    8989      }
    9090      // solve the linear system by Householder transformations
    … …  
    9797      trajectory2.Scale(gsl_vector_get(x,1));
    9898      normal.Scale(gsl_vector_get(x,2));
    99       TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal
     99      TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
    100100                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
    101101      if (TestVector.Norm() < MYEPSILON) {
    … …  
    126126{
    127127  double result = 0.;
    128   atom * Runner = mol->start;
    129   while (Runner->next != mol->end) {
    130     Runner = Runner->next;
    131     if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
     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)) {
    132130  //    atom *Sprinter = PermutationMap[Walker->nr];
    133   //        Log() << Verbose(0) << *Walker << " and " << *Runner << " are heading to the same target at ";
     131  //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
    134132  //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
    135133  //        Log() << Verbose(0) << ", penalting." << endl;
    … …  
    162160  // go through every atom
    163161  atom *Runner = NULL;
    164   atom *Walker = start;
    165   while (Walker->next != end) {
    166     Walker = Walker->next;
     162  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    167163    // first term: distance to target
    168     Runner = Params.PermutationMap[Walker->nr];   // find target point
    169     tmp = (Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
     164    Runner = Params.PermutationMap[(*iter)->nr];   // find target point
     165    tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
    170166    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    171167    result += Params.PenaltyConstants[0] * tmp;
    … …  
    173169
    174170    // second term: sum of distances to other trajectories
    175     result += SumDistanceOfTrajectories(Walker, this, Params);
     171    result += SumDistanceOfTrajectories((*iter), this, Params);
    176172
    177173    // third term: penalty for equal targets
    178     result += PenalizeEqualTargets(Walker, this, Params);
     174    result += PenalizeEqualTargets((*iter), this, Params);
    179175  }
    180176
    … …  
    216212void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
    217213{
    218   for (int i=mol->AtomCount; i--;) {
     214  for (int i=mol->getAtomCount(); i--;) {
    219215    Params.DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    220216    Params.DistanceList[i]->clear();
    221217  }
    222218
    223   atom *Runner = NULL;
    224   atom *Walker = mol->start;
    225   while (Walker->next != mol->end) {
    226     Walker = Walker->next;
    227     Runner = mol->start;
    228     while(Runner->next != mol->end) {
    229       Runner = Runner->next;
    230       Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)), Runner) );
     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)) );
    231222    }
    232223  }
    … …  
    240231void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
    241232{
    242   atom *Walker = mol->start;
    243   while (Walker->next != mol->end) {
    244     Walker = Walker->next;
    245     Params.StepList[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
    246     Params.PermutationMap[Walker->nr] = Params.DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
    247     Params.DoubleList[Params.DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
    248     Params.DistanceIterators[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // and remember which one we picked
    249     DoLog(2) && (Log() << Verbose(2) << *Walker << " starts with distance " << Params.DistanceList[Walker->nr]->begin()->first << "." << endl);
     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);
    250239  }
    251240};
    … …  
    288277void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
    289278{
    290   atom *Walker = mol->start;
     279  molecule::const_iterator iter = mol->begin();
    291280  DistanceMap::iterator NewBase;
    292281  double Potential = fabs(mol->ConstrainedPotential(Params));
    293282
     283  if (mol->empty()) {
     284    eLog() << Verbose(1) << "Molecule is empty." << endl;
     285    return;
     286  }
    294287  while ((Potential) > Params.PenaltyConstants[2]) {
    295     PrintPermutationMap(mol->AtomCount, Params);
    296     Walker = Walker->next;
    297     if (Walker == mol->end) // round-robin at the end
    298       Walker = mol->start->next;
    299     if (Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     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
    300293      continue;
    301294    // now, try finding a new one
    302     Potential = TryNextNearestNeighbourForInjectivePermutation(mol, Walker, Potential, Params);
    303   }
    304   for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1
     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
    305298    if (Params.DoubleList[i] > 1) {
    306299      DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
    … …  
    341334  double Potential, OldPotential, OlderPotential;
    342335  struct EvaluatePotential Params;
    343   Params.PermutationMap = new atom *[AtomCount];
    344   Params.DistanceList = new DistanceMap *[AtomCount];
    345   Params.DistanceIterators = new DistanceMap::iterator[AtomCount];
    346   Params.DoubleList = new int[AtomCount];
    347   Params.StepList = new DistanceMap::iterator[AtomCount];
     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()];
    348341  int round;
    349   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
     342  atom *Sprinter = NULL;
    350343  DistanceMap::iterator Rider, Strider;
    351344
    352345  // set to zero
    353   for (int i=0;i<AtomCount;i++) {
     346  for (int i=0;i<getAtomCount();i++) {
    354347    Params.PermutationMap[i] = NULL;
    355348    Params.DoubleList[i] = 0;
    … …  
    380373    DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl);
    381374    OlderPotential = OldPotential;
     375    molecule::const_iterator iter;
    382376    do {
    383       Walker = start;
    384       while (Walker->next != end) { // pick one
    385         Walker = Walker->next;
    386         PrintPermutationMap(AtomCount, Params);
    387         Sprinter = Params.DistanceIterators[Walker->nr]->second;   // store initial partner
    388         Strider = Params.DistanceIterators[Walker->nr];  //remember old iterator
    389         Params.DistanceIterators[Walker->nr] = Params.StepList[Walker->nr];
    390         if (Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
    391           Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->begin();
     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();
    392385          break;
    393386        }
    394         //Log() << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
     387        //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->nr]->second << "." << endl;
    395388        // find source of the new target
    396         Runner = start->next;
    397         while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    398           if (Params.PermutationMap[Runner->nr] == Params.DistanceIterators[Walker->nr]->second) {
    399             //Log() << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
     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;
    400393            break;
    401394          }
    402           Runner = Runner->next;
    403395        }
    404         if (Runner != end) { // we found the other source
     396        if (runner != end()) { // we found the other source
    405397          // then look in its distance list for Sprinter
    406           Rider = Params.DistanceList[Runner->nr]->begin();
    407           for (; Rider != Params.DistanceList[Runner->nr]->end(); Rider++)
     398          Rider = Params.DistanceList[(*runner)->nr]->begin();
     399          for (; Rider != Params.DistanceList[(*runner)->nr]->end(); Rider++)
    408400            if (Rider->second == Sprinter)
    409401              break;
    410           if (Rider != Params.DistanceList[Runner->nr]->end()) { // if we have found one
    411             //Log() << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
     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;
    412404            // exchange both
    413             Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
    414             Params.PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
    415             PrintPermutationMap(AtomCount, Params);
     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);
    416408            // calculate the new potential
    417409            //Log() << Verbose(2) << "Checking new potential ..." << endl;
    … …  
    419411            if (Potential > OldPotential) { // we made everything worse! Undo ...
    420412              //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    421               //Log() << Verbose(3) << "Setting " << *Runner << "'s source to " << *Params.DistanceIterators[Runner->nr]->second << "." << endl;
     413              //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->nr]->second << "." << endl;
    422414              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    423               Params.PermutationMap[Runner->nr] = Params.DistanceIterators[Runner->nr]->second;
     415              Params.PermutationMap[(*runner)->nr] = Params.DistanceIterators[(*runner)->nr]->second;
    424416              // Undo for Walker
    425               Params.DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
    426               //Log() << Verbose(3) << "Setting " << *Walker << "'s source to " << *Params.DistanceIterators[Walker->nr]->second << "." << endl;
    427               Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
     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;
    428420            } else {
    429               Params.DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
     421              Params.DistanceIterators[(*runner)->nr] = Rider;  // if successful also move the pointer in the iterator list
    430422              DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
    431423              OldPotential = Potential;
    … …  
    437429            //Log() << Verbose(0) << endl;
    438430          } else {
    439             DoeLog(1) && (eLog()<< Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl);
     431            DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl);
    440432            exit(255);
    441433          }
    442434        } else {
    443           Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // new target has no source!
     435          Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // new target has no source!
    444436        }
    445         Params.StepList[Walker->nr]++; // take next farther distance target
     437        Params.StepList[(*iter)->nr]++; // take next farther distance target
    446438      }
    447     } while (Walker->next != end);
     439    } while (++iter != end());
    448440  } while ((OlderPotential - OldPotential) > 1e-3);
    449441  DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    … …  
    451443
    452444  /// free memory and return with evaluated potential
    453   for (int i=AtomCount; i--;)
     445  for (int i=getAtomCount(); i--;)
    454446    Params.DistanceList[i]->clear();
    455447  delete[](Params.DistanceList);
    … …  
    492484  // Get the Permutation Map by MinimiseConstrainedPotential
    493485  atom **PermutationMap = NULL;
    494   atom *Walker = NULL, *Sprinter = NULL;
     486  atom *Sprinter = NULL;
    495487  if (!MapByIdentity)
    496488    MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
    497489  else {
    498     PermutationMap = new atom *[AtomCount];
     490    PermutationMap = new atom *[getAtomCount()];
    499491    SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
    500492  }
    … …  
    511503    mol = World::getInstance().createMolecule();
    512504    MoleculePerStep->insert(mol);
    513     Walker = start;
    514     while (Walker->next != end) {
    515       Walker = Walker->next;
     505    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    516506      // add to molecule list
    517       Sprinter = mol->AddCopyAtom(Walker);
     507      Sprinter = mol->AddCopyAtom((*iter));
    518508      for (int n=NDIM;n--;) {
    519         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);
     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);
    520510        // add to Trajectories
    521511        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    522512        if (step < MaxSteps) {
    523           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);
    524           Walker->Trajectory.U.at(step)[n] = 0.;
    525           Walker->Trajectory.F.at(step)[n] = 0.;
     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.;
    526516        }
    527517      }
    … …  
    531521
    532522  // store the list to single step files
    533   int *SortIndex = new int[AtomCount];
    534   for (int i=AtomCount; i--; )
     523  int *SortIndex = new int[getAtomCount()];
     524  for (int i=getAtomCount(); i--; )
    535525    SortIndex[i] = i;
    536526  status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
    … …  
    578568      return false;
    579569    }
    580     if (Force.RowCounter[0] != AtomCount) {
    581       DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl);
     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);
    582572      performCriticalExit();
    583573      return false;
    … …  
    585575    // correct Forces
    586576    Velocity.Zero();
    587     for(int i=0;i<AtomCount;i++)
     577    for(int i=0;i<getAtomCount();i++)
    588578      for(int d=0;d<NDIM;d++) {
    589579        Velocity[d] += Force.Matrix[0][i][d+5];
    590580      }
    591     for(int i=0;i<AtomCount;i++)
     581    for(int i=0;i<getAtomCount();i++)
    592582      for(int d=0;d<NDIM;d++) {
    593         Force.Matrix[0][i][d+5] -= Velocity[d]/(double)AtomCount;
     583        Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
    594584      }
    595585    // solve a constrained potential if we are meant to
    … …  
    694684      delta_alpha = 0.;
    695685      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    696       delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
     686      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
    697687      configuration.alpha += delta_alpha*configuration.Deltat;
    698688      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
Note: See TracChangeset for help on using the changeset viewer.