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