Changeset 435065 for src

Timestamp:

Mar 2, 2011, 9:48:20 PM (14 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:

341f22

Parents:

eb33c4

Message:

Rewrote VerletForceIntegration into a functor in Dynamics/.

• we now have the regression test working with checking against an integrated file which has not been checked though (just by eye and logged output to make sense, not against other code).

Location:

src

Files:

• Actions/MoleculeAction/VerletIntegrationAction.cpp (modified) (2 diffs)
• Actions/MoleculeAction/VerletIntegrationAction.def (modified) (1 diff)
• Dynamics/VerletForceIntegration.hpp (added)
• Makefile.am (modified) (1 diff)
• atom_atominfo.cpp (modified) (2 diffs)
• atom_atominfo.hpp (modified) (2 diffs)
• molecule.hpp (modified) (1 diff)
• molecule_dynamics.cpp (modified) (1 diff)

src/Actions/MoleculeAction/VerletIntegrationAction.cpp

@@ -20 +20 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include "atom.hpp"
+#include "AtomSet.hpp"
 #include "CodePatterns/Log.hpp"
+#include "CodePatterns/Verbose.hpp"
+#include "Dynamics/VerletForceIntegration.hpp"
 #include "molecule.hpp"
-#include "CodePatterns/Verbose.hpp"
 #include "World.hpp"
 
+#include <vector>
 #include <iostream>
 #include <fstream>
@@ -38 +42 @@
 /** =========== define the function ====================== */
 Action::state_ptr MoleculeVerletIntegrationAction::performCall() {
-  molecule *mol = NULL;
-
   // obtain information
   getParametersfromValueStorage();
 
-  for (World::MoleculeSelectionIterator iter = World::getInstance().beginMoleculeSelection(); iter != World::getInstance().endMoleculeSelection(); ++iter) {
-    mol = iter->second;
-    DoLog(1) && (Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl);
-    // TODO: sollte besser stream nutzen, nicht filename direkt (es sei denn, ist prefix), besser fuer unit test
-    char outputname[MAXSTRINGSIZE];
-    strcpy(outputname, params.forcesfile.string().c_str());
-    if (!mol->VerletForceIntegration(outputname, *(World::getInstance().getConfig()), 0))
+
+  DoLog(1) && (Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl);
+  // TODO: sollte besser stream nutzen, nicht filename direkt (es sei denn, ist prefix), besser fuer unit test
+  char outputname[MAXSTRINGSIZE];
+  strcpy(outputname, params.forcesfile.string().c_str());
+  AtomSetMixin<std::vector<atom *> > set(World::getInstance().getSelectedAtoms());
+  for (int step = 0; step < params.MDSteps; ++step) {
+    VerletForceIntegration<std::vector<atom *> > Verlet(set, params.Deltat, step, false);
+    if (!Verlet(outputname, 1, 0, params.FixedCenterOfMass))
       DoLog(2) && (Log() << Verbose(2) << "File " << params.forcesfile << " not found." << endl);
     else
       DoLog(2) && (Log() << Verbose(2) << "File " << params.forcesfile << " found and parsed." << endl);
   }
+
   return Action::success;
 }

src/Actions/MoleculeAction/VerletIntegrationAction.def

@@ -12 +12 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NODEFAULT) for each (undefined) default value
-#define paramtypes (boost::filesystem::path)
-#define paramtokens ("verlet-integration")
-#define paramdescriptions ("perform verlet integration of a given force file")
+#define paramtypes (boost::filesystem::path)(double)(int)(bool)
+#define paramtokens ("verlet-integration")("deltat")("MDSteps")("keep-fixed-CenterOfMass")
+#define paramdescriptions ("perform verlet integration of a given force file")("time step width")("number of MDSteps to integrate")("whether forces and velocities shall be corrected such that center of mass remains at rest")
 #undef paramdefaults
-#define paramreferences (forcesfile)
+#define paramreferences (forcesfile)(Deltat)(MDSteps)(FixedCenterOfMass)
 
 #undef statetypes

src/Makefile.am

@@ -125 +125 @@
         Dynamics/LinearInterpolationBetweenSteps.hpp \
         Dynamics/MinimiseConstrainedPotential.hpp \
-        Dynamics/OutputTemperature.hpp
+        Dynamics/OutputTemperature.hpp \
+        Dynamics/VerletForceIntegration.hpp
 
 GRAPHSOURCE = \

src/atom_atominfo.cpp

@@ -30 +30 @@
 #include "atom_atominfo.hpp"
 
+#include <iomanip>
+
 /** Constructor of class AtomInfo.
  */
@@ -492 +494 @@
  * Parameters are according to those in configuration class.
  * \param NextStep index of sequential step to set
- * \param *configuration pointer to configuration with parameters
+ * \param Deltat time step width
+ * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
  * \param *Force matrix with forces
  */
-void AtomInfo::VelocityVerletUpdate(int nr, const unsigned int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
-{
+void AtomInfo::VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem, ForceMatrix *Force, const size_t offset)
+{
+  LOG(2, "INFO: Particle that currently " << *this);
+  LOG(2, "INFO: Integrating with mass=" << getMass() << " and Deltat="
+      << Deltat << " at NextStep=" << NextStep);
   // update force
   // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
   Vector tempVector;
-  for (int d=0; d<NDIM; d++)
-    tempVector[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+  for (int d=0; d<NDIM; d++) {
+    ASSERT(Force->RowCounter[0] > nr,
+        "AtomInfo::VelocityVerletUpdate() - "+toString(nr)
+        +" out of bounds [0,"+toString(Force->RowCounter[0])
+        +"] access on force matrix.");
+    ASSERT(Force->ColumnCounter[0] > d+offset,
+        "AtomInfo::VelocityVerletUpdate() - "+toString(d+offset)
+        +" out of bounds [0,"+toString(Force->ColumnCounter[0])
+        +"] access on force matrix.");
+    tempVector[d] = -Force->Matrix[0][nr][d+offset]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
+  }
+  LOG(3, "INFO: Force at step " << NextStep << " set to " << tempVector);
   setAtomicForceAtStep(NextStep, tempVector);
 
   // update position
   tempVector = getPositionAtStep(NextStep-1);
-  tempVector += configuration->Deltat*(getAtomicVelocityAtStep(NextStep-1));     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-  tempVector += 0.5*configuration->Deltat*configuration->Deltat*(getAtomicForceAtStep(NextStep))*(1./getMass());     // F = m * a and s =
+  LOG(4, "INFO: initial position from last step " << setprecision(4) << tempVector);
+  tempVector += Deltat*(getAtomicVelocityAtStep(NextStep-1));     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
+  LOG(4, "INFO: position with velocity " << getAtomicVelocityAtStep(NextStep-1) << " from last step " << tempVector);
+  tempVector += 0.5*Deltat*Deltat*(getAtomicForceAtStep(NextStep))*(1./getMass());     // F = m * a and s =
+  LOG(4, "INFO: position with force " << getAtomicForceAtStep(NextStep) << " from last step " << tempVector);
   setPositionAtStep(NextStep, tempVector);
+  LOG(3, "INFO: Position at step " << NextStep << " set to " << tempVector);
 
   // Update U
   tempVector = getAtomicVelocityAtStep(NextStep-1);
-  tempVector += configuration->Deltat * (getAtomicForceAtStep(NextStep)+getAtomicForceAtStep(NextStep-1))*(1./getMass()); // v = F/m * t
+  LOG(4, "INFO: initial velocity from last step " << tempVector);
+  tempVector += Deltat * (getAtomicForceAtStep(NextStep)+getAtomicForceAtStep(NextStep-1))*(1./getMass()); // v = F/m * t
+  LOG(4, "INFO: Velocity with force from last " << getAtomicForceAtStep(NextStep-1)
+      << " and present " << getAtomicForceAtStep(NextStep) << " step " << tempVector);
   setAtomicVelocityAtStep(NextStep, tempVector);
-
-  // some info for debugging
-  DoLog(2) && (Log() << Verbose(2)
-      << "Integrated position&velocity of step " << (NextStep) << ": ("
-      << getPositionAtStep(NextStep) << ")\t("
-      << getAtomicVelocityAtStep(NextStep) << ")" << std::endl);
+  LOG(3, "INFO: Velocity at step " << NextStep << " set to " << tempVector);
 };
 

src/atom_atominfo.hpp

@@ -27 +27 @@
 
 class AtomInfo;
-class config;
 class element;
 class ForceMatrix;
@@ -242 +241 @@
   size_t getTrajectorySize() const;
   void CopyStepOnStep(const unsigned int dest, const unsigned int src);
-  void VelocityVerletUpdate(int nr, const unsigned int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
+  void VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem, ForceMatrix *Force, const size_t offset);
   double getKineticEnergy(const unsigned int step) const;
   Vector getMomentum(const unsigned int step) const;

src/molecule.hpp

@@ -176 +176 @@
   void SetBoxDimension(Vector *dim);
   bool ScanForPeriodicCorrection();
-  bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
   double VolumeOfConvexEnvelope(bool IsAngstroem);
   RealSpaceMatrix getInertiaTensor() const;

src/molecule_dynamics.cpp

@@ -39 +39 @@
 
 
-/** 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.
- * \param *file filename
- * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
- * \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, const size_t offset)
-{
-  Info FunctionInfo(__func__);
-  string token;
-  stringstream item;
-  double IonMass, ConstrainedPotentialEnergy, ActualTemp;
-  Vector Velocity;
-  ForceMatrix Force;
-
-  const int AtomCount = getAtomCount();
-  // parse file into ForceMatrix
-  std::ifstream input(file);
-  if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
-    performCriticalExit();
-    return false;
-  }
-  input.close();
-  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();
-    return false;
-  }
-  // correct Forces
-  Velocity.Zero();
-  for(int i=0;i<AtomCount;i++)
-    for(int d=0;d<NDIM;d++) {
-      Velocity[d] += Force.Matrix[0][i][d+offset];
-    }
-  for(int i=0;i<AtomCount;i++)
-    for(int d=0;d<NDIM;d++) {
-      Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
-    }
-  // solve a constrained potential if we are meant to
-  if (configuration.DoConstrainedMD) {
-    // calculate forces and potential
-    std::map<atom*, atom *> PermutationMap;
-    MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
-    //double ConstrainedPotentialEnergy =
-    Minimiser(configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
-    Minimiser.EvaluateConstrainedForces(&Force);
-  }
-
-  // and perform Verlet integration for each atom with position, velocity and force vector
-  // check size of vectors
-  BOOST_FOREACH(atom *_atom, atoms) {
-    _atom->VelocityVerletUpdate(_atom->getNr(), MDSteps+1, &configuration, &Force, (const size_t) 0);
-  }
-
-  // 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);
-  ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1);
-  Berendsen berendsen = Berendsen();
-  berendsen.addToContainer(World::getInstance().getThermostats());
-  double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms);
-  DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
-  MDSteps++;
-
-  // exit
-  return true;
-};