Changes in src/molecule_dynamics.cpp [112b09:389cc8]

File:

• src/molecule_dynamics.cpp (modified) (11 diffs)

src/molecule_dynamics.cpp

@@ -13 +13 @@
 #include "element.hpp"
 #include "info.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "parser.hpp"
 #include "Plane.hpp"
+#include "ThermoStatContainer.hpp"
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
 
 /************************************* Functions for class molecule *********************************/
@@ -472 +477 @@
  * \param startstep stating initial configuration in molecule::Trajectories
  * \param endstep stating final configuration in molecule::Trajectories
+ * \param &prefix path and prefix
  * \param &config configuration structure
  * \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()
  * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
  */
-bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
+bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, std::string &prefix, config &configuration, bool MapByIdentity)
 {
   molecule *mol = NULL;
@@ -524 +530 @@
   for (int i=getAtomCount(); i--; )
     SortIndex[i] = i;
-  status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
+
+  status = MoleculePerStep->OutputConfigForListOfFragments(prefix, SortIndex);
   delete[](SortIndex);
 
@@ -537 +544 @@
  * have to transform them).
  * This adds a new MD step to the config file.
- * \param *out output stream for debugging
  * \param *file filename
  * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
- * \param delta_t time step width in atomic units
- * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
- * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential()
+ * \param offset offset in matrix file to the first force component
  * \return true - file found and parsed, false - file not found or imparsable
  * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
  */
-bool molecule::VerletForceIntegration(char *file, config &configuration)
+bool molecule::VerletForceIntegration(char *file, config &configuration, const size_t offset)
 {
   Info FunctionInfo(__func__);
@@ -556 +560 @@
   ForceMatrix Force;
 
-  CountElements();  // make sure ElementsInMolecule is up to date
-
+  const int AtomCount = getAtomCount();
   // check file
   if (input == NULL) {
@@ -568 +571 @@
       return false;
     }
-    if (Force.RowCounter[0] != getAtomCount()) {
+    if (Force.RowCounter[0] != AtomCount) {
       DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
       performCriticalExit();
@@ -575 +578 @@
     // correct Forces
     Velocity.Zero();
-    for(int i=0;i<getAtomCount();i++)
+    for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Velocity[d] += Force.Matrix[0][i][d+5];
+        Velocity[d] += Force.Matrix[0][i][d+offset];
       }
-    for(int i=0;i<getAtomCount();i++)
+    for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
+        Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
       }
     // solve a constrained potential if we are meant to
@@ -596 +599 @@
     //ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
 
-    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force);
+    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force, (const size_t) 0);
   }
   // correct velocities (rather momenta) so that center of mass remains motionless
@@ -643 +646 @@
 
   // calculate scale configuration
-  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
+  ScaleTempFactor = configuration.Thermostats->TargetTemp/ActualTemp;
 
   // differentating between the various thermostats
@@ -651 +654 @@
       break;
      case Woodcock:
-      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
+      if ((configuration.Thermostats->ScaleTempStep > 0) && ((MDSteps-1) % configuration.Thermostats->ScaleTempStep == 0)) {
         DoLog(2) && (Log() << Verbose(2) <<  "Applying Woodcock thermostat..." << endl);
         ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
@@ -684 +687 @@
       delta_alpha = 0.;
       ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
-      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
-      configuration.alpha += delta_alpha*configuration.Deltat;
-      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
+      delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.Thermostats->TargetTemp)/(configuration.Thermostats->HooverMass*Units2Electronmass);
+      configuration.Thermostats->alpha += delta_alpha*configuration.Deltat;
+      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.Thermostats->alpha << "." << endl);
       // apply updated alpha as additional force
       ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );