Changeset 51cdfd for src/Dynamics/VerletForceIntegration.hpp

Timestamp:

Aug 20, 2014, 12:58:28 PM (10 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:

442495

Parents:

343c5a

git-author:

Frederik Heber <heber@…> (08/20/14 12:50:41)

git-committer:

Frederik Heber <heber@…> (08/20/14 12:58:28)

Message:

Extracted common functions from VerletForceIntegration into AtomicForceManipulator.

• AtomicForceManipulator contains e.g. parsing forces from file.
• removed performConstraintMinimization from VerletForceIntegration.
• this is preparatory for ForceAnnealing, also derived from AtomicForceManipulator.
• FIX: MoleculeVerletIntegrationAction checks whether atoms are selected and increments WorldTime only in success case.

File:

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

src/Dynamics/VerletForceIntegration.hpp

@@ -20 +20 @@
 #include "CodePatterns/Log.hpp"
 #include "CodePatterns/Verbose.hpp"
-#include "Dynamics/MinimiseConstrainedPotential.hpp"
+#include "Dynamics/AtomicForceManipulator.hpp"
 #include "Fragmentation/ForceMatrix.hpp"
 #include "Helpers/helpers.hpp"
@@ -30 +30 @@
 
 template <class T>
-class VerletForceIntegration
+class VerletForceIntegration : public AtomicForceManipulator<T>
 {
 public:
@@ -40 +40 @@
    */
   VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, bool _IsAngstroem) :
-    Deltat(_Deltat),
-    IsAngstroem(_IsAngstroem),
-    atoms(_atoms)
+    AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem)
   {}
   /** Destructor of class VerletForceIntegration.
@@ -49 +47 @@
   ~VerletForceIntegration()
   {}
-
-  /** Parses nuclear forces from file.
-   * Forces are stored in the time step \a TimeStep in the atomicForces in \a atoms.
-   * \param *file filename
-   * \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__);
-    ForceMatrix Force;
-
-    // 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 << ".");
-      return false;
-    }
-    input.close();
-    if (Force.RowCounter[0] != (int)atoms.size()) {
-      ELOG(0, "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << atoms.size() << ".");
-      return false;
-    }
-
-    addForceMatrixToAtomicForce(Force, TimeStep, 1);
-    return true;
-  }
 
   /** Performs Verlet integration.
@@ -94 +65 @@
     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-1);
+      AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(offset,NextStep-1);
 
     // solve a constrained potential if we are meant to
-    if (DoConstrainedMD)
-      performConstraintMinimization(DoConstrainedMD,NextStep-1);
+//    if (DoConstrainedMD)
+//      performConstraintMinimization(DoConstrainedMD,NextStep-1);
 
     if (NextStep > 0) {
-      for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+      for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
+          iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
         //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
-        (*iter)->VelocityVerletUpdateU((*iter)->getId(), NextStep-1, Deltat, IsAngstroem);
+        (*iter)->VelocityVerletUpdateU(
+            (*iter)->getId(),
+            NextStep-1,
+            AtomicForceManipulator<T>::Deltat,
+            AtomicForceManipulator<T>::IsAngstroem);
       }
 
       // make sum of velocities equal zero
       if (FixedCenterOfMass)
-        correctVelocitiesForFixedCenterOfMass(NextStep-1);
+        AtomicForceManipulator<T>::correctVelocitiesForFixedCenterOfMass(NextStep-1);
 
       // thermostat
-      performThermostatControl(NextStep-1);
+      AtomicForceManipulator<T>::performThermostatControl(NextStep-1);
     }
 
@@ -125 +95 @@
     // and perform Verlet integration for each atom with position, velocity and force vector
     // check size of vectors
-    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
+    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
+        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
       //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
-      (*iter)->VelocityVerletUpdateX((*iter)->getId(), NextStep, Deltat, IsAngstroem);
+      (*iter)->VelocityVerletUpdateX(
+          (*iter)->getId(),
+          NextStep,
+          AtomicForceManipulator<T>::Deltat,
+          AtomicForceManipulator<T>::IsAngstroem);
     }
 
     // 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:
-  double Deltat;
-  bool IsAngstroem;
-  AtomSetMixin<T> atoms;
 };