source: src/Dynamics/VerletForceIntegration.hpp@ 321470

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
Last change on this file since 321470 was bcb593, checked in by Frederik Heber <heber@…>, 11 years ago

FIX: VerletIntegrationAction now assumes forces have just been calculated.

  • according to Wikipedia's Velocity_Verlet info, we first integrate position, then calculate forces, then integrate velocity. This assumes that forces based on next time step's position are already known. This is possible when parsed from file but not if they are calculated dynamically via fragmentation. Hence, we now integrate velocity from last time step to current, then integrate position from current time step to next. Then we are in the position to calculate forces and do this cycle again.
  • for this, VelocityVerletUpdate was split into ..X() and ..U() for independent integration of position and velocity.
  • VelocityVerletIntegration::operator() now first corrects forces, then integrates velocites, corrects them too and finally integrates positions according to the above new scheme.
  • removed option MDSteps from VerletIntegrationAction.
  • TESTFIX: Had to change regression rest on VerletIntegration since the cycle sequence has changed. It was nonsense before to have the forces already in some file respective to future position that actually first need to come out of the time integration.
  • Property mode set to 100644
File size: 7.5 KB
RevLine 
[435065]1/*
2 * VerletForceIntegration.hpp
3 *
4 * Created on: Feb 23, 2011
5 * Author: heber
6 */
7
8#ifndef VERLETFORCEINTEGRATION_HPP_
9#define VERLETFORCEINTEGRATION_HPP_
10
11// include config.h
12#ifdef HAVE_CONFIG_H
13#include <config.h>
14#endif
15
[6f0841]16#include "Atom/atom.hpp"
17#include "Atom/AtomSet.hpp"
[435065]18#include "CodePatterns/Assert.hpp"
19#include "CodePatterns/Info.hpp"
20#include "CodePatterns/Log.hpp"
21#include "CodePatterns/Verbose.hpp"
22#include "Dynamics/MinimiseConstrainedPotential.hpp"
[a9b86d]23#include "Fragmentation/ForceMatrix.hpp"
[255829]24#include "Helpers/helpers.hpp"
[4882d5]25#include "Helpers/defs.hpp"
[435065]26#include "LinearAlgebra/Vector.hpp"
[ab26c3]27#include "Thermostats/ThermoStatContainer.hpp"
[4882d5]28#include "Thermostats/Thermostat.hpp"
[435065]29#include "World.hpp"
30
31template <class T>
32class VerletForceIntegration
33{
34public:
[4882d5]35 /** Constructor of class VerletForceIntegration.
36 *
37 * \param _atoms set of atoms to integrate
38 * \param _Deltat time step width in atomic units
39 * \param _IsAngstroem whether length units are in angstroem or bohr radii
40 */
41 VerletForceIntegration(AtomSetMixin<T> &_atoms, double _Deltat, bool _IsAngstroem) :
[435065]42 Deltat(_Deltat),
43 IsAngstroem(_IsAngstroem),
[4882d5]44 atoms(_atoms)
[435065]45 {}
[4882d5]46 /** Destructor of class VerletForceIntegration.
47 *
48 */
[435065]49 ~VerletForceIntegration()
50 {}
51
[4882d5]52 /** Parses nuclear forces from file.
53 * Forces are stored in the time step \a TimeStep in the atomicForces in \a atoms.
[435065]54 * \param *file filename
[4882d5]55 * \param TimeStep time step to parse forces file into
56 * \return true - file parsed, false - file not found or imparsable
[435065]57 */
[4882d5]58 bool parseForcesFile(const char *file, const int TimeStep)
[435065]59 {
60 Info FunctionInfo(__func__);
61 ForceMatrix Force;
62
63 // parse file into ForceMatrix
64 std::ifstream input(file);
65 if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
[47d041]66 ELOG(0, "Could not parse Force Matrix file " << file << ".");
[435065]67 return false;
68 }
69 input.close();
[4882d5]70 if (Force.RowCounter[0] != (int)atoms.size()) {
[47d041]71 ELOG(0, "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << atoms.size() << ".");
[435065]72 return false;
73 }
74
[4882d5]75 addForceMatrixToAtomicForce(Force, TimeStep, 1);
76 return true;
77 }
78
79 /** Performs Verlet integration.
[bcb593]80 *
81 * We assume that forces have just been calculated. Then, we perform the velocity
82 * and the position calculation for \f$ t + \Delta t \f$, such that forces may be
83 * again calculated with respect to the new position.
84 *
[4882d5]85 * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
86 * \param offset offset in matrix file to the first force component
87 * \param DoConstrainedMD whether a constrained MD shall be done
88 * \param FixedCenterOfMass whether forces and velocities are correct to have fixed center of mass
89 * \return true - file found and parsed, false - no atoms, file not found or imparsable
90 * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
91 */
92 bool operator()(const int NextStep, const size_t offset, const int DoConstrainedMD, const bool FixedCenterOfMass)
93 {
94 Info FunctionInfo(__func__);
95
96 // check that atoms are present at all
97 if (atoms.size() == 0) {
98 ELOG(2, "VerletForceIntegration::operator() - no atoms to integrate.");
99 return false;
[435065]100 }
101
[4882d5]102 // make sum of forces equal zero
103 if (FixedCenterOfMass)
[bcb593]104 correctForceMatrixForFixedCenterOfMass(offset,NextStep-1);
[4882d5]105
[435065]106 // solve a constrained potential if we are meant to
[bcb593]107 if (DoConstrainedMD)
108 performConstraintMinimization(DoConstrainedMD,NextStep-1);
109
110 if (NextStep > 0) {
111 for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
112 //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
113 (*iter)->VelocityVerletUpdateU((*iter)->getId(), NextStep-1, Deltat, IsAngstroem);
114 }
115
116 // make sum of velocities equal zero
117 if (FixedCenterOfMass)
118 correctVelocitiesForFixedCenterOfMass(NextStep-1);
119
120 // thermostat
121 performThermostatControl(NextStep-1);
[435065]122 }
123
124 //std::cout << "Force before velocity verlet, " << Force << std::endl;
125 // and perform Verlet integration for each atom with position, velocity and force vector
126 // check size of vectors
127 for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
128 //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
[bcb593]129 (*iter)->VelocityVerletUpdateX((*iter)->getId(), NextStep, Deltat, IsAngstroem);
[435065]130 }
131
[4882d5]132 // exit
133 return true;
134 };
135
136private:
137 void addForceMatrixToAtomicForce(const ForceMatrix &Force, const int &TimeStep, const int offset)
138 {
139 // place forces from matrix into atoms
140 Vector tempVector;
141 size_t i=0;
142 for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter,++i) {
143 for(size_t d=0;d<NDIM;d++) {
144 tempVector[d] = Force.Matrix[0][i][d+offset]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
145 }
146 tempVector += (*iter)->getAtomicForceAtStep(TimeStep);
147 (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
[435065]148 }
[4882d5]149 }
[435065]150
[4882d5]151 void correctForceMatrixForFixedCenterOfMass(const size_t offset, const int &TimeStep) {
152 Vector ForceVector;
153 // correct Forces
154 //std::cout << "Force before correction, " << Force << std::endl;
155 ForceVector.Zero();
156 for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
157 ForceVector += (*iter)->getAtomicForceAtStep(TimeStep);
158 }
159 ForceVector.Scale(1./(double)atoms.size());
160 //std::cout << "Force before second correction, " << Force << std::endl;
161 for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
162 const Vector tempVector = (*iter)->getAtomicForceAtStep(TimeStep) - ForceVector;
163 (*iter)->setAtomicForceAtStep(TimeStep, tempVector);
164 }
165 LOG(3, "INFO: forces correct by " << ForceVector << "each.");
166 }
167
168 void correctVelocitiesForFixedCenterOfMass(const int &TimeStep) {
169 Vector Velocity;
170 double IonMass;
171 // correct velocities (rather momenta) so that center of mass remains motionless
172 Velocity = atoms.totalMomentumAtStep(TimeStep);
173 IonMass = atoms.totalMass();
174
175 // correct velocities (rather momenta) so that center of mass remains motionless
176 Velocity *= 1./IonMass;
177 atoms.addVelocityAtStep(-1.*Velocity,TimeStep);
178
179 LOG(3, "INFO: Velocities corrected by " << Velocity << " each.");
180 }
181
182 void performConstraintMinimization(const int &DoConstrainedMD, const int &TimeStep) {
183 // calculate forces and potential
184 ForceMatrix Force;
185 std::map<atom *, atom*> PermutationMap;
186 MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
187 //double ConstrainedPotentialEnergy =
188 Minimiser(DoConstrainedMD, 0, IsAngstroem);
189 Minimiser.EvaluateConstrainedForces(&Force);
190 addForceMatrixToAtomicForce(Force, TimeStep, 1);
191 }
192
193 void performThermostatControl(const int &TimeStep) {
194 double ActualTemp;
195
196 // calculate current temperature
197 ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
[435065]198 LOG(3, "INFO: Current temperature is " << ActualTemp);
[4882d5]199
200 // rescale to desired value
201 double ekin = World::getInstance().getThermostats()->getActive()->scaleAtoms(TimeStep,ActualTemp,atoms);
202 ActualTemp = atoms.totalTemperatureAtStep(TimeStep);
[435065]203 LOG(3, "INFO: New temperature after thermostat is " << ActualTemp);
[47d041]204 LOG(1, "Kinetic energy is " << ekin << ".");
[4882d5]205 }
[435065]206
207private:
208 double Deltat;
209 bool IsAngstroem;
210 AtomSetMixin<T> atoms;
211};
212
213#endif /* VERLETFORCEINTEGRATION_HPP_ */
Note: See TracBrowser for help on using the repository browser.