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 |
|
---|
16 | #include "Atom/atom.hpp"
|
---|
17 | #include "Atom/AtomSet.hpp"
|
---|
18 | #include "CodePatterns/Assert.hpp"
|
---|
19 | #include "CodePatterns/Info.hpp"
|
---|
20 | #include "CodePatterns/Log.hpp"
|
---|
21 | #include "CodePatterns/Verbose.hpp"
|
---|
22 | #include "Dynamics/AtomicForceManipulator.hpp"
|
---|
23 | #include "Fragmentation/ForceMatrix.hpp"
|
---|
24 | #include "Helpers/helpers.hpp"
|
---|
25 | #include "Helpers/defs.hpp"
|
---|
26 | #include "LinearAlgebra/Vector.hpp"
|
---|
27 | #include "Thermostats/ThermoStatContainer.hpp"
|
---|
28 | #include "Thermostats/Thermostat.hpp"
|
---|
29 | #include "World.hpp"
|
---|
30 |
|
---|
31 | template <class T>
|
---|
32 | class VerletForceIntegration : public AtomicForceManipulator<T>
|
---|
33 | {
|
---|
34 | public:
|
---|
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) :
|
---|
42 | AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem)
|
---|
43 | {}
|
---|
44 | /** Destructor of class VerletForceIntegration.
|
---|
45 | *
|
---|
46 | */
|
---|
47 | ~VerletForceIntegration()
|
---|
48 | {}
|
---|
49 |
|
---|
50 | /** Performs Verlet integration.
|
---|
51 | *
|
---|
52 | * We assume that forces have just been calculated. Then, we perform the velocity
|
---|
53 | * and the position calculation for \f$ t + \Delta t \f$, such that forces may be
|
---|
54 | * again calculated with respect to the new position.
|
---|
55 | *
|
---|
56 | * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
|
---|
57 | * \param offset offset in matrix file to the first force component
|
---|
58 | * \param DoConstrainedMD whether a constrained MD shall be done
|
---|
59 | * \param FixedCenterOfMass whether forces and velocities are correct to have fixed center of mass
|
---|
60 | * \return true - file found and parsed, false - no atoms, file not found or imparsable
|
---|
61 | * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
|
---|
62 | */
|
---|
63 | bool operator()(const int NextStep, const size_t offset, const int DoConstrainedMD, const bool FixedCenterOfMass)
|
---|
64 | {
|
---|
65 | Info FunctionInfo(__func__);
|
---|
66 |
|
---|
67 | // make sum of forces equal zero
|
---|
68 | if (FixedCenterOfMass)
|
---|
69 | AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(offset,NextStep-1);
|
---|
70 |
|
---|
71 | // solve a constrained potential if we are meant to
|
---|
72 | // if (DoConstrainedMD)
|
---|
73 | // performConstraintMinimization(DoConstrainedMD,NextStep-1);
|
---|
74 |
|
---|
75 | if (NextStep > 0) {
|
---|
76 | for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
77 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
78 | //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
|
---|
79 | (*iter)->VelocityVerletUpdateU(
|
---|
80 | (*iter)->getId(),
|
---|
81 | NextStep-1,
|
---|
82 | AtomicForceManipulator<T>::Deltat,
|
---|
83 | AtomicForceManipulator<T>::IsAngstroem);
|
---|
84 | }
|
---|
85 |
|
---|
86 | // make sum of velocities equal zero
|
---|
87 | if (FixedCenterOfMass)
|
---|
88 | AtomicForceManipulator<T>::correctVelocitiesForFixedCenterOfMass(NextStep-1);
|
---|
89 |
|
---|
90 | // thermostat
|
---|
91 | AtomicForceManipulator<T>::performThermostatControl(NextStep-1);
|
---|
92 | }
|
---|
93 |
|
---|
94 | //std::cout << "Force before velocity verlet, " << Force << std::endl;
|
---|
95 | // and perform Verlet integration for each atom with position, velocity and force vector
|
---|
96 | // check size of vectors
|
---|
97 | for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
|
---|
98 | iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
|
---|
99 | //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
|
---|
100 | (*iter)->VelocityVerletUpdateX(
|
---|
101 | (*iter)->getId(),
|
---|
102 | NextStep,
|
---|
103 | AtomicForceManipulator<T>::Deltat,
|
---|
104 | AtomicForceManipulator<T>::IsAngstroem);
|
---|
105 | }
|
---|
106 |
|
---|
107 | // exit
|
---|
108 | return true;
|
---|
109 | }
|
---|
110 | };
|
---|
111 |
|
---|
112 | #endif /* VERLETFORCEINTEGRATION_HPP_ */
|
---|