Changeset 4a7776a

Timestamp:

Oct 12, 2009, 10:30:02 AM (15 years ago)

Author:

Frederik Heber <heber@…>

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:

5034e1

Parents:

ccd9f5

Message:

Complete refactoring of molecule_dynamics.cpp

• new small functions: SumDistanceOfTrajectories(), PenalizeEqualTargets(), (), CreateInitialLists(), TryNextNearestNeighbourForInjectivePermutation(), MakeInjectivePermutation()
• Thermostat routines (case blocks) have been shifted to small functions in atom class (and soon to be extracted from there into its own class: ThermostatedParticle ?)
• old functions are a lot smaller now: MinimiseConstrainedPotential(), molecule::EvaluateConstrainedForces(), VerletForceIntegration(), also using ActOnAllAtoms
• atom class has a lot of new small functions:
  • trajectory: ResizeTrajectory(), CopyStepOnStep(), VelocityVerletUpdate(), VelocityVerletUpdate()
  • constraint potential: CorrectVelocity(), Thermostat_...()

Location:

src

Files:

• atom.cpp (modified) (3 diffs)
• atom.hpp (modified) (3 diffs)
• molecule.cpp (modified) (1 diff)
• molecule.hpp (modified) (1 diff)
• molecule_dynamics.cpp (modified) (7 diffs)
• molecule_template.hpp (modified) (1 diff)

src/atom.cpp

@@ -7 +7 @@
 #include "atom.hpp"
 #include "bond.hpp"
+#include "config.hpp"
 #include "element.hpp"
 #include "memoryallocator.hpp"
@@ -265 +266 @@
 };
 
+/** Extends the trajectory STL vector to the new size.
+ * Does nothing if \a MaxSteps is smaller than current size.
+ * \param MaxSteps
+ */
+void atom::ResizeTrajectory(int MaxSteps)
+{
+  if (Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
+    //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
+    Trajectory.R.resize(MaxSteps+1);
+    Trajectory.U.resize(MaxSteps+1);
+    Trajectory.F.resize(MaxSteps+1);
+  }
+};
+
+/** Copies a given trajectory step \a src onto another \a dest
+ * \param dest index of destination step
+ * \param src index of source step
+ */
+void atom::CopyStepOnStep(int dest, int src)
+{
+  if (dest == src)  // self assignment check
+    return;
+
+  for (int n=NDIM;n--;) {
+    Trajectory.R.at(dest).x[n] = Trajectory.R.at(src).x[n];
+    Trajectory.U.at(dest).x[n] = Trajectory.U.at(src).x[n];
+    Trajectory.F.at(dest).x[n] = Trajectory.F.at(src).x[n];
+  }
+};
+
+/** Performs a velocity verlet update of the trajectory.
+ * Parameters are according to those in configuration class.
+ * \param NextStep index of sequential step to set
+ * \param *configuration pointer to configuration with parameters
+ * \param *Force matrix with forces
+ */
+void atom::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
+{
+  //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
+  for (int d=0; d<NDIM; d++) {
+    Trajectory.F.at(NextStep).x[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    Trajectory.R.at(NextStep).x[d] = Trajectory.R.at(NextStep-1).x[d];
+    Trajectory.R.at(NextStep).x[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+    Trajectory.R.at(NextStep).x[d] += 0.5*configuration->Deltat*configuration->Deltat*(Trajectory.F.at(NextStep).x[d]/type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+  }
+  // Update U
+  for (int d=0; d<NDIM; d++) {
+    Trajectory.U.at(NextStep).x[d] = Trajectory.U.at(NextStep-1).x[d];
+    Trajectory.U.at(NextStep).x[d] += configuration->Deltat * (Trajectory.F.at(NextStep).x[d]+Trajectory.F.at(NextStep-1).x[d]/type->mass); // v = F/m * t
+  }
+  // Update R (and F)
+//      out << "Integrated position&velocity of step " << (NextStep) << ": (";
+//      for (int d=0;d<NDIM;d++)
+//        out << Trajectory.R.at(NextStep).x[d] << " ";          // next step
+//      out << ")\t(";
+//      for (int d=0;d<NDIM;d++)
+//        cout << Trajectory.U.at(NextStep).x[d] << " ";          // next step
+//      out << ")" << endl;
+};
+
+/** Sums up mass and kinetics.
+ * \param Step step to sum for
+ * \param *TotalMass pointer to total mass sum
+ * \param *TotalVelocity pointer to tota velocity sum
+ */
+void atom::SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity )
+{
+  *TotalMass += type->mass;  // sum up total mass
+  for(int d=0;d<NDIM;d++) {
+    TotalVelocity->x[d] += Trajectory.U.at(Step).x[d]*type->mass;
+  }
+};
+
 /** Returns squared distance to a given vector.
  * \param origin vector to calculate distance to
@@ -312 +386 @@
     Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
 };
+
+/** Correct velocity against the summed \a CoGVelocity for \a step.
+ * \param *ActualTemp sum up actual temperature meanwhile
+ * \param Step MD step in atom::Tracjetory
+ * \param *CoGVelocity remnant velocity (i.e. vector sum of all atom velocities)
+ */
+void atom::CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity)
+{
+  for(int d=0;d<NDIM;d++) {
+    Trajectory.U.at(Step).x[d] -= CoGVelocity->x[d];
+    *ActualTemp += 0.5 * type->mass * Trajectory.U.at(Step).x[d] * Trajectory.U.at(Step).x[d];
+  }
+};
+
+/** Scales velocity of atom according to Woodcock thermostat.
+ * \param ScaleTempFactor factor to scale the velocities with (i.e. sqrt of energy scale factor)
+ * \param Step MD step to scale
+ * \param *ekin sum of kinetic energy
+ */
+void atom::Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] *= ScaleTempFactor;
+      *ekin += 0.5*type->mass * U[d]*U[d];
+    }
+};
+
+/** Scales velocity of atom according to Gaussian thermostat.
+ * \param Step MD step to scale
+ * \param *G
+ * \param *E
+ */
+void atom::Thermostat_Gaussian_init(int Step, double *G, double *E)
+{
+  double *U = Trajectory.U.at(Step).x;
+  double *F = Trajectory.F.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      *G += U[d] * F[d];
+      *E += U[d]*U[d]*type->mass;
+    }
+};
+
+/** Determines scale factors according to Gaussian thermostat.
+ * \param Step MD step to scale
+ * \param GE G over E ratio
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and TargetTemp
+ */
+void atom::Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] += configuration->Deltat/type->mass * ( (G_over_E) * (U[d]*type->mass) );
+      *ekin += type->mass * U[d]*U[d];
+    }
+};
+
+/** Scales velocity of atom according to Langevin thermostat.
+ * \param Step MD step to scale
+ * \param *r random number generator
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and TargetTemp
+ */
+void atom::Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration)
+{
+  double sigma  = sqrt(configuration->TargetTemp/type->mass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    // throw a dice to determine whether it gets hit by a heat bath particle
+    if (((((rand()/(double)RAND_MAX))*configuration->TempFrequency) < 1.)) {
+      cout << Verbose(3) << "Particle " << *this << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
+      // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
+      for (int d=0; d<NDIM; d++) {
+        U[d] = gsl_ran_gaussian (r, sigma);
+      }
+      cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
+    }
+    for (int d=0; d<NDIM; d++)
+      *ekin += 0.5*type->mass * U[d]*U[d];
+  }
+};
+
+/** Scales velocity of atom according to Berendsen thermostat.
+ * \param Step MD step to scale
+ * \param ScaleTempFactor factor to scale energy (not velocity!) with
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and Deltat
+ */
+void atom::Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      U[d] *= sqrt(1+(configuration->Deltat/configuration->TempFrequency)*(ScaleTempFactor-1));
+      *ekin += 0.5*type->mass * U[d]*U[d];
+    }
+  }
+};
+
+/** Initializes current run of NoseHoover thermostat.
+ * \param Step MD step to scale
+ * \param *delta_alpha additional sum of kinetic energy on return
+ */
+void atom::Thermostat_NoseHoover_init(int Step, double *delta_alpha)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+      *delta_alpha += U[d]*U[d]*type->mass;
+    }
+  }
+};
+
+/** Initializes current run of NoseHoover thermostat.
+ * \param Step MD step to scale
+ * \param *ekin sum of kinetic energy
+ * \param *configuration configuration class with TempFrequency and Deltat
+ */
+void atom::Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration)
+{
+  double *U = Trajectory.U.at(Step).x;
+  if (FixedIon == 0) { // even FixedIon moves, only not by other's forces
+    for (int d=0; d<NDIM; d++) {
+        U[d] += configuration->Deltat/type->mass * (configuration->alpha * (U[d] * type->mass));
+        *ekin += (0.5*type->mass) * U[d]*U[d];
+      }
+  }
+};

src/atom.hpp

@@ -21 +21 @@
 #include <vector>
 
+#include <gsl/gsl_randist.h>
+
 #include "tesselation.hpp"
 
@@ -26 +28 @@
 
 class bond;
+class config;
 class element;
 class ForceMatrix;
@@ -80 +83 @@
   bool Compare(const atom &ptr);
 
+  // trajectory stuff
+  void ResizeTrajectory(int MaxSteps);
+  void CopyStepOnStep(int dest, int src);
+  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
+  void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity );
+
   double DistanceToVector(Vector &origin);
   double DistanceSquaredToVector(Vector &origin);
 
+  bool IsInParallelepiped(Vector offset, double *parallelepiped);
+
+  // constraint potential and dynamics stuff
   void AddKineticToTemperature(double *temperature, int step) const;
   void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
+  void CorrectVelocity(double *ActualTemp, int Step, Vector *CoGVelocity);
 
-  bool IsInParallelepiped(Vector offset, double *parallelepiped);
+  // thermostats
+  void Thermostat_Woodcock(double ScaleTempFactor, int Step, double *ekin);
+  void Thermostat_Gaussian_init(int Step, double *G, double *E);
+  void Thermostat_Gaussian_least_constraint(int Step, double G_over_E, double *ekin, config *configuration);
+  void Thermostat_Langevin(int Step, gsl_rng * r, double *ekin, config *configuration);
+  void Thermostat_Berendsen(int Step, double ScaleTempFactor, double *ekin, config *configuration);
+  void Thermostat_NoseHoover_init(int Step, double *delta_alpha);
+  void Thermostat_NoseHoover_scale(int Step, double *ekin, config *configuration);
+
 
   ostream & operator << (ostream &ost);

src/molecule.cpp

@@ -1140 +1140 @@
   return true;
 };
+
+void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int TesselPoint::*index)
+{
+  atom *Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    array[(Walker->*index)] = Walker;
+  }
+};

src/molecule.hpp

@@ -176 +176 @@
 
   // templates for allowing global manipulation of an array with one entry per atom
+  void SetIndexedArrayForEachAtomTo ( atom **array, int TesselPoint::* index);
   template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T));
   template <typename T> void SetIndexedArrayForEachAtomTo ( T *array, int TesselPoint::* index, void (*Setor)(T *, T), T t);

src/molecule_dynamics.cpp

@@ -497 +497 @@
   else {
     PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      PermutationMap[Walker->nr] = Walker;   // always pick target with the smallest distance
-    }
+    SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
   }
 
   // check whether we have sufficient space in Trajectories for each atom
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if (Walker->Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
-      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
-      Walker->Trajectory.R.resize(MaxSteps+1);
-      Walker->Trajectory.U.resize(MaxSteps+1);
-      Walker->Trajectory.F.resize(MaxSteps+1);
-    }
-  }
+  ActOnAllAtoms( &atom::ResizeTrajectory, MaxSteps );
   // push endstep to last one
-  Walker = start;
-  while (Walker->next != end) { // remove the endstep (is now the very last one)
-    Walker = Walker->next;
-    for (int n=NDIM;n--;) {
-      Walker->Trajectory.R.at(MaxSteps).x[n] = Walker->Trajectory.R.at(endstep).x[n];
-      Walker->Trajectory.U.at(MaxSteps).x[n] = Walker->Trajectory.U.at(endstep).x[n];
-      Walker->Trajectory.F.at(MaxSteps).x[n] = Walker->Trajectory.F.at(endstep).x[n];
-    }
-  }
+  ActOnAllAtoms( &atom::CopyStepOnStep, MaxSteps, endstep );
   endstep = MaxSteps;
 
@@ -578 +557 @@
 bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
 {
-  atom *walker = NULL;
   ifstream input(file);
   string token;
   stringstream item;
-  double IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
+  double IonMass, ConstrainedPotentialEnergy, ActualTemp;
+  Vector Velocity;
   ForceMatrix Force;
 
@@ -601 +580 @@
     }
     // correct Forces
-    for(int d=0;d<NDIM;d++)
-      Vector[d] = 0.;
+    Velocity.Zero();
     for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Vector[d] += Force.Matrix[0][i][d+5];
+        Velocity.x[d] += Force.Matrix[0][i][d+5];
       }
     for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
+        Force.Matrix[0][i][d+5] -= Velocity.x[d]/(double)AtomCount;
       }
     // solve a constrained potential if we are meant to
@@ -621 +599 @@
 
     // and perform Verlet integration for each atom with position, velocity and force vector
-    walker = start;
-    while (walker->next != end) { // go through every atom of this element
-      walker = walker->next;
-      //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-      // check size of vectors
-      if (walker->Trajectory.R.size() <= (unsigned int)(MDSteps)) {
-        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-        walker->Trajectory.R.resize(MDSteps+10);
-        walker->Trajectory.U.resize(MDSteps+10);
-        walker->Trajectory.F.resize(MDSteps+10);
-      }
-
-      // Update R (and F)
-      for (int d=0; d<NDIM; d++) {
-        walker->Trajectory.F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-        walker->Trajectory.R.at(MDSteps).x[d] = walker->Trajectory.R.at(MDSteps-1).x[d];
-        walker->Trajectory.R.at(MDSteps).x[d] += configuration.Deltat*(walker->Trajectory.U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-        walker->Trajectory.R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(walker->Trajectory.F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
-      }
-      // Update U
-      for (int d=0; d<NDIM; d++) {
-        walker->Trajectory.U.at(MDSteps).x[d] = walker->Trajectory.U.at(MDSteps-1).x[d];
-        walker->Trajectory.U.at(MDSteps).x[d] += configuration.Deltat * (walker->Trajectory.F.at(MDSteps).x[d]+walker->Trajectory.F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
-      }
-//      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
-//      for (int d=0;d<NDIM;d++)
-//        out << walker->Trajectory.R.at(MDSteps).x[d] << " ";          // next step
-//      out << ")\t(";
-//      for (int d=0;d<NDIM;d++)
-//        cout << walker->Trajectory.U.at(MDSteps).x[d] << " ";          // next step
-//      out << ")" << endl;
-            // next atom
-    }
+    // check size of vectors
+    ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
+
+    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps, &configuration, &Force);
   }
   // correct velocities (rather momenta) so that center of mass remains motionless
-  for(int d=0;d<NDIM;d++)
-    Vector[d] = 0.;
+  Velocity.Zero();
   IonMass = 0.;
-  walker = start;
-  while (walker->next != end) { // go through every atom
-    walker = walker->next;
-    IonMass += walker->type->mass;  // sum up total mass
-    for(int d=0;d<NDIM;d++) {
-      Vector[d] += walker->Trajectory.U.at(MDSteps).x[d]*walker->type->mass;
-    }
-  }
+  ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps, &IonMass, &Velocity );
+
   // correct velocities (rather momenta) so that center of mass remains motionless
-  for(int d=0;d<NDIM;d++)
-    Vector[d] /= IonMass;
+  Velocity.Scale(1./IonMass);
   ActualTemp = 0.;
-  walker = start;
-  while (walker->next != end) { // go through every atom of this element
-    walker = walker->next;
-    for(int d=0;d<NDIM;d++) {
-      walker->Trajectory.U.at(MDSteps).x[d] -= Vector[d];
-      ActualTemp += 0.5 * walker->type->mass * walker->Trajectory.U.at(MDSteps).x[d] * walker->Trajectory.U.at(MDSteps).x[d];
-    }
-  }
+  ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps, &Velocity );
   Thermostats(configuration, ActualTemp, Berendsen);
   MDSteps++;
-
 
   // exit
@@ -712 +645 @@
   double delta_alpha = 0.;
   double ScaleTempFactor;
-  double sigma;
-  double IonMass;
-  int d;
   gsl_rng * r;
   const gsl_rng_type * T;
-  double *U = NULL, *F = NULL, FConstraint[NDIM];
-  atom *walker = NULL;
 
   // calculate scale configuration
@@ -731 +659 @@
       if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
         cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
-        walker = start;
-        while (walker->next != end) { // go through every atom of this element
-          walker = walker->next;
-          IonMass = walker->type->mass;
-          U = walker->Trajectory.U.at(MDSteps).x;
-          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-            for (d=0; d<NDIM; d++) {
-              U[d] *= sqrt(ScaleTempFactor);
-              ekin += 0.5*IonMass * U[d]*U[d];
-            }
-        }
+        ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
       }
       break;
      case Gaussian:
       cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            G += U[d] * F[d];
-            E += U[d]*U[d]*IonMass;
-          }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
+
       cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            FConstraint[d] = (G/E) * (U[d]*IonMass);
-            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
-            ekin += IonMass * U[d]*U[d];
-          }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
+
       break;
      case Langevin:
@@ -780 +677 @@
       r = gsl_rng_alloc (T);
       // Go through each ion
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          // throw a dice to determine whether it gets hit by a heat bath particle
-          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) {
-            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
-            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
-            for (d=0; d<NDIM; d++) {
-              U[d] = gsl_ran_gaussian (r, sigma);
-            }
-            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
-          }
-          for (d=0; d<NDIM; d++)
-            ekin += 0.5*IonMass * U[d]*U[d];
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
       break;
+
      case Berendsen:
       cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
-            ekin += 0.5*IonMass * U[d]*U[d];
-          }
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
       break;
+
      case NoseHoover:
       cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
       // dynamically evolve alpha (the additional degree of freedom)
       delta_alpha = 0.;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            delta_alpha += U[d]*U[d]*IonMass;
-          }
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
       delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
       configuration.alpha += delta_alpha*configuration.Deltat;
       cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
       // apply updated alpha as additional force
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
-              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
-              ekin += (0.5*IonMass) * U[d]*U[d];
-            }
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
       break;
   }

src/molecule_template.hpp

@@ -436 +436 @@
 };
 
-// ===================== Accesssing arrays indexed by some integer for each atom ======================
+// ===================== Accessing arrays indexed by some integer for each atom ======================
 
 // for atom ints