1 #include "PaulTrapIntegrator.h" 2 #include "ScalarStructure.h" 3 #include "Vector3DBlock.h" 4 #include "ForceGroup.h" 5 #include "GenericTopology.h" 6 #include "topologyutilities.h" 7 8 using std::string; 9 using std::vector; 10 using std::pair; 11 12 using namespace ProtoMol::Report; 13 14 namespace ProtoMol { 15 16 //________________________________________________________ ThermostatType 17 18 const string ThermostatEnum::str[static_cast<int>(LAST)-static_cast<int>(FIRST)] = { 19 // Order is essential, must be in relation to Enum 20 string("undefined"), // Returned when no enum matches 21 string("NVT"), 22 string("NVT_zero"), 23 string("NVT_ind"), 24 string("NVT_shell"), 25 string("NVT_global"), 26 string("berendsen"), 27 string("berendsen_zero"), 28 string("berendsen_ind"), 29 string("berendsen_shell"), 30 string("berendsen_global") 31 }; 32 33 //_________________________________________________________________ PaulTrapIntegrator 34 35 const string PaulTrapIntegrator::keyword("PaulTrap"); 36 PaulTrapIntegrator()37 PaulTrapIntegrator::PaulTrapIntegrator(): STSIntegrator(), 38 myCached(false), 39 myTemperature(0.0), 40 myThermalInertia(0.0), 41 myBathPosition(0.0), 42 myBathVelocity(1.0){} 43 PaulTrapIntegrator(Real timestep,Real temperature,Real thermalInertia,Real bathPosition,Real bathVelocity,ThermostatType nvttype,Real part,const vector<Real> & time,const vector<Real> & t,ForceGroup * overloadedForces)44 PaulTrapIntegrator::PaulTrapIntegrator(Real timestep, 45 Real temperature, 46 Real thermalInertia, 47 Real bathPosition, 48 Real bathVelocity, 49 ThermostatType nvttype, 50 Real part, 51 const vector<Real>& time, 52 const vector<Real>& t, 53 ForceGroup *overloadedForces) 54 55 : STSIntegrator(timestep, overloadedForces), 56 myCached(false), 57 myTemperature(temperature), 58 myThermalInertia(thermalInertia), 59 myBathPosition(bathPosition), 60 myBathVelocity(bathVelocity), 61 myThermostatType(nvttype), 62 myPart(std::max(std::min(part,1.0),0.0)), 63 myPartReal(std::max(std::min(part,1.0),0.0)), 64 myCount(0){ 65 66 // Test 67 const unsigned int count = time.size(); 68 if(count != t.size()) 69 report << error <<"[PaulTrapIntegrator::PaulTrapIntegrator] size of time("<<count<<") and temperature("<<t.size()<<") vectors differ!"<<endr; 70 71 // Sort 72 vector<pair<Real,Real> > tmp(count); 73 for(unsigned int i=0;i<count;++i){ 74 tmp[i].first = time[i]; 75 tmp[i].second = t[i]; 76 } 77 std::sort(tmp.begin(),tmp.end()); 78 myT.resize(count); 79 myTime.resize(count); 80 for(unsigned int i=0;i<count;++i){ 81 myTime[i] = tmp[i].first; 82 myT[i] = tmp[i].second; 83 } 84 } 85 doDrift(void)86 void PaulTrapIntegrator::doDrift(void) { 87 if(!myCached) 88 init(); 89 90 const unsigned int numberOfAtoms = myTopo->atoms.size(); 91 Real deltaT = getTimestep()/Constant::TIMEFACTOR; 92 93 Real theTemperature = myTemperature; 94 for(int i=0;i<static_cast<int>(myT.size())-1;++i){ 95 if(myTopo->time >= myTime[i] && myTopo->time < myTime[i+1]) { 96 theTemperature = myT[i]+(myTopo->time-myTime[i])/(myTime[i+1]-myTime[i])*(myT[i+1]-myT[i]); 97 break; 98 } 99 } 100 101 // Scaling of velocities 102 unsigned int countZero = 0; 103 switch(myThermostatType){ 104 case ThermostatType::NVT_SHELL: 105 case ThermostatType::BERENDSEN_SHELL: 106 { 107 vector<Real> vavg(myLayer.size(),0.0); 108 vector<int> count(myLayer.size(),0); 109 for (unsigned int i=0;i<myLayer.size();i++) { 110 for (unsigned int j=0;j<myLayer[i].size();j++) { 111 unsigned int k = myLayer[i][j]; 112 Real rx = (*myPositions)[k].norm(); 113 if(rx > 1e-10){ 114 count[i]++; 115 vavg[i] += (*myVelocities)[k].dot((*myPositions)[k])/rx; 116 } 117 } 118 } 119 for (unsigned int i=0;i<vavg.size();i++) { 120 Real v = 0.0; 121 if(count[i] > 0) 122 v = vavg[i]/static_cast<Real>(count[i]); 123 124 //report << hint << "Vavg["<<i<<"] ("<<myLayer[i].size()<<") ="<<v<<"."<<endr; 125 for (unsigned int j=0;j<myLayer[i].size();j++) { 126 unsigned int k = myLayer[i][j]; 127 Real rx = (*myPositions)[k].norm(); 128 if(rx > 1e-10){ 129 (*myVelocities)[k] = (*myPositions)[k]*(v/rx); 130 } 131 else { 132 (*myVelocities)[k] = Vector3D(0.0,0.0,0.0); 133 countZero++; 134 } 135 } 136 } 137 } 138 break; 139 case ThermostatType::NVT_GLOBAL: 140 case ThermostatType::BERENDSEN_GLOBAL: 141 { 142 Real vavg=0.0; 143 int count = 0; 144 for (unsigned int i=0;i<numberOfAtoms;i++) { 145 if(myKeep[i] != 0){ 146 Real rx = (*myPositions)[i].norm(); 147 if(rx > 1e-10){ 148 count++; 149 vavg += (*myVelocities)[i].dot((*myPositions)[i])/rx; 150 } 151 } 152 153 } 154 if(count > 0){ 155 vavg /=static_cast<Real>(count); 156 //report << hint << "Vavg="<<vavg<<"."<<endr; 157 for (unsigned int i=0;i<numberOfAtoms;i++) { 158 if(myKeep[i] != 0){ 159 Real rx = (*myPositions)[i].norm(); 160 if(rx > 1e-10){ 161 (*myVelocities)[i] = (*myPositions)[i]*(vavg/rx); 162 } 163 else { 164 (*myVelocities)[i] = Vector3D(0.0,0.0,0.0); 165 countZero++; 166 } 167 } 168 } 169 } 170 } 171 break; 172 case ThermostatType::BERENDSEN_IND: 173 case ThermostatType::NVT_IND: 174 { 175 for (unsigned int i=0;i<numberOfAtoms;i++) { 176 if(myKeep[i] != 0){ 177 Real rx = (*myPositions)[i].normSquared(); 178 if(rx > 1e-10){ 179 (*myVelocities)[i] = (*myPositions)[i]*((*myVelocities)[i].dot((*myPositions)[i])/rx); 180 } 181 else { 182 (*myVelocities)[i] = Vector3D(0.0,0.0,0.0); 183 countZero++; 184 } 185 } 186 } 187 } 188 break; 189 case ThermostatType::BERENDSEN_ZERO: 190 case ThermostatType::NVT_ZERO: 191 { 192 for (unsigned int i=0;i<numberOfAtoms;i++) { 193 if(myKeep[i] != 0){ 194 (*myVelocities)[i] = Vector3D(0.0,0.0,0.0); 195 countZero++; 196 } 197 } 198 } 199 break; 200 default: 201 break; 202 } 203 204 205 // Integration 206 switch(myThermostatType){ 207 case ThermostatType::BERENDSEN: 208 case ThermostatType::BERENDSEN_IND: 209 case ThermostatType::BERENDSEN_ZERO: 210 case ThermostatType::BERENDSEN_GLOBAL: 211 case ThermostatType::BERENDSEN_SHELL: 212 { 213 //report.precision(13); 214 //report << hint << myBathVelocity; 215 Real t = temperature(myTopo,myVelocities); 216 if(t > 0.0 && countZero < numberOfAtoms) 217 t *= static_cast<Real>(numberOfAtoms)/static_cast<Real>(numberOfAtoms-countZero); 218 myBathVelocity = t > 0.0?(myBathVelocity*(1-myThermalInertia) + theTemperature/t*myThermalInertia):1.0; 219 Real a = sqrt(myBathVelocity); 220 //if(t < theTemperature ){ 221 // a = 1.0; 222 // myBathVelocity = 1.0; 223 //} 224 //report << "a="<<a-1<<", T="<<t<<endr; 225 for (unsigned int i=0;i<numberOfAtoms;i++) { 226 (*myVelocities)[i] *= a; 227 (*myPositions)[i] += (*myVelocities)[i]*deltaT; 228 } 229 } 230 break; 231 default: 232 { 233 Real kineticEnergy = 0.0; 234 //report.precision(13); 235 //report << hint << "a="<<-myBathPosition<<", T="<<temperature(myTopo,myVelocities)<<endr; 236 for (unsigned int i=0;i<numberOfAtoms;i++) { 237 Real m = myTopo->atoms[i].scaledMass; 238 (*myVelocities)[i] -= (*myVelocities)[i]*myBathPosition; 239 (*myPositions)[i] += (*myVelocities)[i]*deltaT; 240 kineticEnergy += ((*myVelocities)[i]).normSquared()*m; 241 } 242 kineticEnergy *=0.5; 243 Real e = (3.0/2.0*(Real)numberOfAtoms * Constant::BOLTZMANN * theTemperature); 244 245 myBathPosition += (kineticEnergy-e)*myThermalInertia/(Real)numberOfAtoms; 246 } 247 break; 248 } 249 buildMolecularCenterOfMass(myPositions,myTopo); 250 buildMolecularMomentum(myVelocities,myTopo); 251 } 252 initialize(GenericTopology * topo,Vector3DBlock * positions,Vector3DBlock * velocities,ScalarStructure * energies)253 void PaulTrapIntegrator::initialize(GenericTopology *topo, 254 Vector3DBlock *positions, 255 Vector3DBlock *velocities, 256 ScalarStructure *energies){ 257 STSIntegrator::initialize(topo,positions,velocities,energies); 258 initializeForces(); 259 myCached = false; 260 init(); 261 262 report.precision(13); 263 report << plain 264 << "Paul Trap I : thermostate = " << myThermostatType <<"\n" 265 << " : part = "<<myPart<<"\n"; 266 for(unsigned int i=0;i<myT.size();++i) 267 report << " : t"<<i<<" = ("<<myTime[i]<<","<<myT[i]<<")\n"; 268 report << " : layers = "<<myLayer.size()<<"\n" 269 << " : holding = "<<myCount<<" ("<<100.0*myPartReal<<"%, r="<<pow(myPartReal,1.0/3.0)<<")"; 270 report << endr; 271 } 272 273 getParameters(vector<Parameter> & parameters) const274 void PaulTrapIntegrator::getParameters(vector<Parameter>& parameters) const { 275 STSIntegrator::getParameters(parameters); 276 parameters.push_back(Parameter("temperature",Value(myTemperature,ConstraintValueType::NotNegative()))); 277 parameters.push_back(Parameter("thermal",Value(myThermalInertia))); 278 parameters.push_back(Parameter("bathPos",Value(myBathPosition),0.0)); 279 parameters.push_back(Parameter("bathVel",Value(myBathVelocity),1.0)); 280 parameters.push_back(Parameter("scheme",Value(myThermostatType.getString(),ConstraintValueType::NotEmpty()),std::string("NVT"),Text(std::string("thermostat scheme (")+ThermostatType::getPossibleValues()+std::string(")")))); 281 parameters.push_back(Parameter("part",Value(myPart,ConstraintValueType::NotNegative()),0.0)); 282 parameters.push_back(Parameter("time",Value(myTime),vector<Real>())); 283 parameters.push_back(Parameter("t",Value(myT,ConstraintValueType::NotNegative()),vector<Real>())); 284 } 285 doMake(string & errMsg,const vector<Value> & values,ForceGroup * fg) const286 STSIntegrator* PaulTrapIntegrator::doMake(string& errMsg, const vector<Value>& values,ForceGroup* fg)const{ 287 ThermostatType nvttype = values[5].getString(); 288 if(!nvttype.valid()){ 289 errMsg += " ThermostatType \'"+values[5].getString()+ 290 "\' not recognized, possible values are: "+ThermostatType::getPossibleValues(",")+"."; 291 return NULL; 292 } 293 if(values[7].size() != values[8].size()){ 294 errMsg += " size of time("+toString(values[7].size())+") and temperature("+toString(values[8].size())+") vectors differ."; 295 } 296 return new PaulTrapIntegrator(values[0],values[1],values[2],values[3],values[4],nvttype,values[6],values[7],values[8],fg); 297 } 298 init()299 void PaulTrapIntegrator::init() { 300 const unsigned int numberOfAtoms = myTopo->atoms.size(); 301 302 myKeep.clear(); 303 myKeep.resize(numberOfAtoms,0); 304 myLayer.clear(); 305 306 myCount = 0; 307 if(myThermostatType != "nvt" && myThermostatType != "berendsen" && myPart > 0.0){ 308 vector<pair<Real,unsigned int> > order(numberOfAtoms); 309 for (unsigned int i=0;i<numberOfAtoms;i++) { 310 order[i].first = (*myPositions)[i].norm(); 311 order[i].second = i; 312 } 313 sort(order.begin(),order.end()); 314 Real maxr = -1.0; 315 Real last = -1.0; 316 for (unsigned int i=0;i<numberOfAtoms;i++) { 317 if(maxr >= 0.0 && fabs(order[i].first-maxr) > 1e-2) 318 break; 319 if(fabs(order[i].first-last) > 1e-2){ 320 myLayer.resize(myLayer.size()+1); 321 last = order[i].first; 322 } 323 myLayer[myLayer.size()-1].push_back(i); 324 if(static_cast<unsigned int>(numberOfAtoms*myPart) <= i && maxr < 0.0){ 325 maxr=order[i].first; 326 } 327 myKeep[order[i].second] = 1; 328 myCount++; 329 } 330 } 331 332 myPartReal = (numberOfAtoms>0?(Real)myCount/(Real)numberOfAtoms:0.0); 333 myCached = true; 334 } 335 doUncache()336 void PaulTrapIntegrator::doUncache(){ 337 myCached = false; 338 } 339 340 } 341