1 #include "NoseNVTLeapfrogIntegrator.h"
2 #include "ScalarStructure.h"
3 #include "Vector3DBlock.h"
4 #include "ForceGroup.h"
5 #include "GenericTopology.h"
6 #include "pmconstants.h"
7 #include "topologyutilities.h"
8 #include "ModifierFriction.h"
9 
10 using std::string;
11 using std::vector;
12 
13 using namespace ProtoMol::Report;
14 
15 namespace ProtoMol {
16 
17   //_________________________________________________________________ NoseNVTLeapfrogIntegrator
18 
19   const string NoseNVTLeapfrogIntegrator::keyword("NoseNVTLeapfrog");
20 
NoseNVTLeapfrogIntegrator()21   NoseNVTLeapfrogIntegrator::NoseNVTLeapfrogIntegrator(): STSIntegrator(),
22 							  myTemperature(0.0),
23 							  myThermalInertia(0.0),
24 							  myBathPosition(0.0){}
25 
NoseNVTLeapfrogIntegrator(Real timestep,Real temperature,Real thermalInertia,Real bathPosition,ForceGroup * overloadedForces)26   NoseNVTLeapfrogIntegrator::NoseNVTLeapfrogIntegrator(Real timestep,
27 						       Real temperature,
28 						       Real thermalInertia,
29 						       Real bathPosition,
30 						       ForceGroup *overloadedForces)
31 
32     : STSIntegrator(timestep, overloadedForces),
33       myTemperature(temperature),
34       myThermalInertia(thermalInertia),
35       myBathPosition(bathPosition){}
36 
friction()37   void NoseNVTLeapfrogIntegrator::friction() {
38 
39     const unsigned int numberOfAtoms   = myTopo->atoms.size();
40     const Real h = getTimestep() * Constant::INV_TIMEFACTOR;
41     Real kineticEnergy = 0.0;
42     for(unsigned int i = 0; i < numberOfAtoms; ++i){
43       Real m = myTopo->atoms[i].scaledMass;
44       kineticEnergy += ((*myVelocities)[i]).normSquared()*m;//+(*myForces)[i]*h/m
45     }
46     myBathPosition += (0.5*kineticEnergy-myTargetKE)*myThermalInertia/(h*numberOfAtoms);
47 
48     for (unsigned int i=0;i<numberOfAtoms;i++) {
49       (*myForces)[i]  -= (*myVelocities)[i]*myBathPosition*myTopo->atoms[i].scaledMass;
50     }
51 
52   }
53 
initialize(GenericTopology * topo,Vector3DBlock * positions,Vector3DBlock * velocities,ScalarStructure * energies)54   void NoseNVTLeapfrogIntegrator::initialize(GenericTopology *topo,
55 					     Vector3DBlock   *positions,
56 					     Vector3DBlock   *velocities,
57 					     ScalarStructure *energies){
58     STSIntegrator::initialize(topo,positions,velocities,energies);
59 
60     myTargetKE = (3.0/2.0*positions->size() * Constant::BOLTZMANN * myTemperature);
61     mySumMass = 0.0;
62     for (unsigned int i=0;i<positions->size();i++)
63       mySumMass += myTopo->atoms[i].scaledMass;
64 
65     initializeForces();
66   }
67 
addModifierAfterInitialize()68   void NoseNVTLeapfrogIntegrator::addModifierAfterInitialize(){
69     adoptPostForceModifier(new ModifierFriction(this));
70     STSIntegrator::addModifierAfterInitialize();
71   }
72 
getParameters(vector<Parameter> & parameters) const73   void NoseNVTLeapfrogIntegrator::getParameters(vector<Parameter>& parameters) const {
74     STSIntegrator::getParameters(parameters);
75     parameters.push_back(Parameter("temperature",Value(myTemperature,ConstraintValueType::NotNegative()),Text("preferred system temperature")));
76     parameters.push_back(Parameter("thermal",Value(myThermalInertia),Text("heat bath coupling: 1.0 very, very strong, 0.0 none")));
77     parameters.push_back(Parameter("bathPos",Value(myBathPosition),0.0,Text("history of the difference of system and heat bath")));
78   }
79 
doMake(string &,const vector<Value> & values,ForceGroup * fg) const80   STSIntegrator* NoseNVTLeapfrogIntegrator::doMake(string&, const vector<Value>& values,ForceGroup* fg)const{
81     return new NoseNVTLeapfrogIntegrator(values[0],values[1],values[2],values[3],fg);
82   }
83 
84 
85 }
86