[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 |
|
---|
| 31 | template <class T>
|
---|
| 32 | class VerletForceIntegration
|
---|
| 33 | {
|
---|
| 34 | public:
|
---|
[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 |
|
---|
| 136 | private:
|
---|
| 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 |
|
---|
| 207 | private:
|
---|
| 208 | double Deltat;
|
---|
| 209 | bool IsAngstroem;
|
---|
| 210 | AtomSetMixin<T> atoms;
|
---|
| 211 | };
|
---|
| 212 |
|
---|
| 213 | #endif /* VERLETFORCEINTEGRATION_HPP_ */
|
---|