Changeset 4882d5 for src/Dynamics

Timestamp:

Dec 14, 2012, 5:39:30 PM (12 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:

5251af

Parents:

8b886f

git-author:

Frederik Heber <heber@…> (10/05/12 07:59:22)

git-committer:

Frederik Heber <heber@…> (12/14/12 17:39:30)

Message:

Rewrote VerletIntegrationAction to diminish dependence on ForceMatrix.

• atom::VelocityVerletUpdate() now works on the forces stored in AtomicForces.
• ThermoStatContainer::getActive() getter for the currently active thermostat.
• VerletForceIntegration::operator() split up into several smaller functions.

File:

• src/Dynamics/VerletForceIntegration.hpp (modified) (4 diffs)

src/Dynamics/VerletForceIntegration.hpp

@@ -23 +23 @@
 #include "Fragmentation/ForceMatrix.hpp"
 #include "Helpers/helpers.hpp"
+#include "Helpers/defs.hpp"
 #include "LinearAlgebra/Vector.hpp"
 #include "Thermostats/ThermoStatContainer.hpp"
-#include "Thermostats/Berendsen.hpp"
+#include "Thermostats/Thermostat.hpp"
 #include "World.hpp"
 
@@ -32 +33 @@
 {
 public:
-  VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, int _startstep, bool _IsAngstroem) :
+  /** Constructor of class VerletForceIntegration.
+   *
+   * \param _atoms set of atoms to integrate
+   * \param _Deltat time step width in atomic units
+   * \param _IsAngstroem whether length units are in angstroem or bohr radii
+   */
+  VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, bool _IsAngstroem) :
     Deltat(_Deltat),
     IsAngstroem(_IsAngstroem),
-    atoms(_atoms),
-    MDSteps(_startstep)
+    atoms(_atoms)
   {}
+  /** Destructor of class VerletForceIntegration.
+   *
+   */
   ~VerletForceIntegration()
   {}
 
-  /** Parses nuclear forces from file and performs Verlet integration.
-   * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
-   * have to transform them).
-   * This adds a new MD step to the config file.
+  /** Parses nuclear forces from file.
+   * Forces are stored in the time step \a TimeStep in the atomicForces in \a atoms.
    * \param *file filename
-   * \param offset offset in matrix file to the first force component
-   * \param DoConstrainedMD whether a constrained MD shall be done
-   * \param FixedCenterOfMass whether forces and velocities are correct to have fixed center of mass
-   * \return true - file found and parsed, false - file not found or imparsable
-   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
-   */
-  bool operator()(char *file, const size_t offset, int DoConstrainedMD, bool FixedCenterOfMass)
+   * \param TimeStep time step to parse forces file into
+   * \return true - file parsed, false - file not found or imparsable
+   */
+  bool parseForcesFile(const char *file, const int TimeStep)
   {
     Info FunctionInfo(__func__);
-    string token;
-    stringstream item;
-    double IonMass, ActualTemp;
     ForceMatrix Force;
 
-    const int AtomCount = atoms.size();
-    ASSERT(AtomCount != 0, "VerletForceIntegration::operator() - no atoms to integrate.");
     // parse file into ForceMatrix
     std::ifstream input(file);
     if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
       ELOG(0, "Could not parse Force Matrix file " << file << ".");
-      performCriticalExit();
       return false;
     }
     input.close();
-    if (Force.RowCounter[0] != AtomCount) {
+    if (Force.RowCounter[0] != (int)atoms.size()) {
       ELOG(0, "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << atoms.size() << ".");
-      performCriticalExit();
       return false;
     }
 
-    if (FixedCenterOfMass) {
-      Vector ForceVector;
-      // correct Forces
-      //std::cout << "Force before correction, " << Force << std::endl;
-      ForceVector.Zero();
-      for(int i=0;i<AtomCount;i++)
-        for(int d=0;d<NDIM;d++) {
-          ForceVector[d] += Force.Matrix[0][i][d+offset];
-        }
-      ForceVector.Scale(1./static_cast<double>(AtomCount));
-      //std::cout << "Force before second correction, " << Force << std::endl;
-      for(int i=0;i<AtomCount;i++)
-        for(int d=0;d<NDIM;d++) {
-          Force.Matrix[0][i][d+offset] -= ForceVector[d];
-        }
-      LOG(3, "INFO: forces correct by " << ForceVector << "each.");
-    }
+    addForceMatrixToAtomicForce(Force, TimeStep, 1);
+    return true;
+  }
+
+  /** Performs Verlet integration.
+   * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
+   * have to transform them).
+   * This adds a new MD step \f$ t + \Delta t \f$ to the config file.
+   * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param offset offset in matrix file to the first force component
+   * \param DoConstrainedMD whether a constrained MD shall be done
+   * \param FixedCenterOfMass whether forces and velocities are correct to have fixed center of mass
+   * \return true - file found and parsed, false - no atoms, file not found or imparsable
+   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
+   */
+  bool operator()(const int NextStep, const size_t offset, const int DoConstrainedMD, const bool FixedCenterOfMass)
+  {
+    Info FunctionInfo(__func__);
+
+    // check that atoms are present at all
+    if (atoms.size() == 0) {
+      ELOG(2, "VerletForceIntegration::operator() - no atoms to integrate.");
+      return false;
+    }
+
+    // make sum of forces equal zero
+    if (FixedCenterOfMass)
+      correctForceMatrixForFixedCenterOfMass(offset,NextStep);
 
     // solve a constrained potential if we are meant to
     if (DoConstrainedMD) {
-      // calculate forces and potential
-      std::map<atom *, atom*> PermutationMap;
-      MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
-      //double ConstrainedPotentialEnergy =
-      Minimiser(DoConstrainedMD, 0, IsAngstroem);
-      Minimiser.EvaluateConstrainedForces(&Force);
+      performConstraintMinimization(DoConstrainedMD,NextStep);
     }
 
@@ -109 +112 @@
     for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
       //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
-      (*iter)->VelocityVerletUpdate((*iter)->getId(), MDSteps+1, Deltat, IsAngstroem, &Force, (const size_t) 1);
-    }
-
-    if (FixedCenterOfMass) {
-      Vector Velocity;
-      // correct velocities (rather momenta) so that center of mass remains motionless
-      Velocity = atoms.totalMomentumAtStep(MDSteps+1);
-      IonMass = atoms.totalMass();
-
-      // correct velocities (rather momenta) so that center of mass remains motionless
-      Velocity *= 1./IonMass;
-      atoms.addVelocityAtStep(-1.*Velocity,MDSteps+1);
-
-      LOG(3, "INFO: Velocities corrected by " << Velocity << " each.");
-    }
+      (*iter)->VelocityVerletUpdate((*iter)->getId(), NextStep, Deltat, IsAngstroem);
+    }
+
+    // make sum of velocities equal zero
+    if (FixedCenterOfMass)
+      correctVelocitiesForFixedCenterOfMass(NextStep);
 
     // thermostat
-    ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
-    LOG(3, "INFO: Current temperature is " << ActualTemp);
-    Berendsen berendsen = Berendsen();
-    berendsen.addToContainer(World::getInstance().getThermostats());
-    double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
-    ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
-    LOG(3, "INFO: New temperature after thermostat is " << ActualTemp);
-    LOG(1, "Kinetic energy is " << ekin << ".");
-
-    // next step
-    MDSteps++;
+    performThermostatControl(NextStep);
 
     // exit
     return true;
   };
+
+private:
+  void addForceMatrixToAtomicForce(const ForceMatrix &Force, const int &TimeStep, const int offset)
+  {
+    // place forces from matrix into atoms
+    Vector tempVector;
+    size_t i=0;
+    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter,++i) {
+      for(size_t d=0;d<NDIM;d++) {
+        tempVector[d] = Force.Matrix[0][i][d+offset]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
+      }
+      tempVector += (*iter)->getAtomicForceAtStep(TimeStep);
+      (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
+    }
+  }
+
+  void correctForceMatrixForFixedCenterOfMass(const size_t offset, const int &TimeStep) {
+    Vector ForceVector;
+    // correct Forces
+    //std::cout << "Force before correction, " << Force << std::endl;
+    ForceVector.Zero();
+    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+      ForceVector += (*iter)->getAtomicForceAtStep(TimeStep);
+    }
+    ForceVector.Scale(1./(double)atoms.size());
+    //std::cout << "Force before second correction, " << Force << std::endl;
+    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+      const Vector tempVector = (*iter)->getAtomicForceAtStep(TimeStep) - ForceVector;
+      (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
+    }
+    LOG(3, "INFO: forces correct by " << ForceVector << "each.");
+  }
+
+  void correctVelocitiesForFixedCenterOfMass(const int &TimeStep) {
+    Vector Velocity;
+    double IonMass;
+    // correct velocities (rather momenta) so that center of mass remains motionless
+    Velocity = atoms.totalMomentumAtStep(TimeStep);
+    IonMass = atoms.totalMass();
+
+    // correct velocities (rather momenta) so that center of mass remains motionless
+    Velocity *= 1./IonMass;
+    atoms.addVelocityAtStep(-1.*Velocity,TimeStep);
+
+    LOG(3, "INFO: Velocities corrected by " << Velocity << " each.");
+  }
+
+  void performConstraintMinimization(const int &DoConstrainedMD, const int &TimeStep) {
+    // calculate forces and potential
+    ForceMatrix Force;
+    std::map<atom *, atom*> PermutationMap;
+    MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
+    //double ConstrainedPotentialEnergy =
+    Minimiser(DoConstrainedMD, 0, IsAngstroem);
+    Minimiser.EvaluateConstrainedForces(&Force);
+    addForceMatrixToAtomicForce(Force, TimeStep, 1);
+  }
+
+  void performThermostatControl(const int &TimeStep) {
+    double ActualTemp;
+
+    // calculate current temperature
+    ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
+    LOG(3, "INFO: Current temperature is " << ActualTemp);
+
+    // rescale to desired value
+    double ekin = World::getInstance().getThermostats()->getActive()->scaleAtoms(TimeStep,ActualTemp,atoms);
+    ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
+    LOG(3, "INFO: New temperature after thermostat is " << ActualTemp);
+    LOG(1, "Kinetic energy is " << ekin << ".");
+  }
 
 private:
@@ -146 +201 @@
   bool IsAngstroem;
   AtomSetMixin<T> atoms;
-  int MDSteps;
 };