1 #include "Report.h" 2 #include "mathutilities.h" 3 #include "pmconstants.h" 4 #include "Vector3DBlock.h" 5 #include "ScalarStructure.h" 6 #include "topologyutilities.h" 7 #include "GenericTopology.h" 8 #include "ModifierUpdateBeta.h" 9 #include "ModifierUpdateBetaAndPush.h" 10 #include "ShadowHMCIntegrator.h" 11 12 #include <iomanip> 13 using namespace std; 14 15 using namespace ProtoMol::Report; 16 using std::vector; 17 using std::string; 18 using std::deque; 19 20 namespace ProtoMol { 21 22 //____________________________________________________________ ShadowHMCIntegrator 23 24 const string ShadowHMCIntegrator::keyword("ShadowHMC"); 25 ShadowHMCIntegrator()26 ShadowHMCIntegrator::ShadowHMCIntegrator() 27 : MCIntegrator(), 28 myOrder(0), 29 myShadowK(0), 30 myShadowKover2(0), 31 myC(0.0), 32 myBeta(0.0), 33 sumTotalEnergy(0.0), 34 myPrevVelocities(NULL), 35 myPrevPositions(NULL), 36 myPrevBeta(NULL), 37 myModifier1(NULL), 38 myModifier2(NULL), 39 myOldForces(NULL){} 40 41 // --------------------------------------------------------------------- // 42 ShadowHMCIntegrator(int cycles,bool randomCycLen,Real initialTemperature,unsigned int shadowOrder,Real c,ForceGroup * overloadedForces,StandardIntegrator * nextIntegrator)43 ShadowHMCIntegrator::ShadowHMCIntegrator(int cycles, 44 bool randomCycLen, 45 Real initialTemperature, 46 unsigned int shadowOrder, 47 Real c, 48 ForceGroup *overloadedForces, 49 StandardIntegrator *nextIntegrator) 50 51 : MCIntegrator(cycles,randomCycLen,initialTemperature,overloadedForces,nextIntegrator), 52 myOrder(shadowOrder), 53 54 // This only works for (order == 2k), where k is even. If/when we 55 // implement k odd, this will need to be changed. 56 myShadowK(myOrder >> 1), 57 myShadowKover2(myOrder >> 2), 58 59 myC(c), 60 myBeta(0.0), 61 sumTotalEnergy(0.0), 62 myPrevVelocities(new deque<Vector3DBlock>()), 63 myPrevPositions(new deque<Vector3DBlock>()), 64 myPrevBeta(new deque<Real>()), 65 myModifier1(NULL), 66 myModifier2(NULL), 67 myOldForces(new Vector3DBlock()) { 68 } 69 70 // --------------------------------------------------------------------- // 71 ~ShadowHMCIntegrator()72 ShadowHMCIntegrator::~ShadowHMCIntegrator(){ 73 if(myPrevVelocities != NULL) 74 delete myPrevVelocities; 75 if(myPrevPositions != NULL) 76 delete myPrevPositions; 77 if(myPrevBeta != NULL) 78 delete myPrevBeta; 79 if(myOldForces != NULL) 80 delete myOldForces; 81 } 82 83 // --------------------------------------------------------------------- // 84 initialize(GenericTopology * topo,Vector3DBlock * positions,Vector3DBlock * velocities,ScalarStructure * energies)85 void ShadowHMCIntegrator::initialize( GenericTopology *topo, 86 Vector3DBlock *positions, 87 Vector3DBlock *velocities, 88 ScalarStructure *energies ) { 89 90 MCIntegrator::initialize(topo,positions,velocities,energies); 91 Report::report << debug(10) <<"[ShadowHMCIntegrator::initialize]"<<Report::endr; 92 93 // Add callbacks 94 myModifier1 = new ModifierUpdateBeta(this); 95 next()->adoptPreDriftOrNextModifier(myModifier1); 96 97 // This one is important because it calls the function to save the shadow 98 // data. 99 myModifier2 = new ModifierUpdateBetaAndPush(this); 100 next()->adoptPostStepModifier(myModifier2); 101 102 // FIXME This needs to be hardcoded eventually. FIXME 103 // Number of samples to estimate myC. 104 // int numSamples = 10; 105 // 106 // myC = -0.1 * ( myEnergies->potentialEnergy() + 107 // kineticEnergy(myTopo,myVelocities) ); 108 // 109 // report << debug(1) << "myC = " << myC << endr; 110 // 111 // run(numSamples); 112 // 113 // myC = 0.90 * ( sumTotalEnergy / numSamples ); 114 // 115 // report << debug(1) << "myC = " << myC << endr; 116 117 } 118 119 // --------------------------------------------------------------------- // 120 run(int numTimesteps)121 void ShadowHMCIntegrator::run(int numTimesteps) { 122 123 Real initTotEnergy = 0.0, 124 initShadowEnergy = 0.0, 125 finShadowEnergy = 0.0, 126 finTotEnergy = 0.0, 127 MDfinEner = 0.0, 128 MDinitEner = 0.0, 129 momentumEnerDiff = 0.0; 130 131 bool acceptMomenta = false; 132 133 // See note in header file. 134 unsigned int localCycleLength = ( myRandomCycLen ? static_cast<int>( 135 myCycleLength * ( 1.3 - 0.6 * randomNumber() ) ) : myCycleLength ); 136 137 // Check that cyclelength is long enough. It is possible to compute the 138 // shadow with a short cyclelength, but it requires a lot of uncessesary 139 // checks. In general, cyclelength should be much larger than 140 // myShadowK. 141 if(localCycleLength < myShadowK) { 142 report << warning << "*** The cyclelength for this step is below the " 143 << "number of steps required for accurate computation of the " 144 << "shadow energy. Adjusting to minimal cycle legnth. ***" << endr; 145 localCycleLength = myShadowK; 146 } 147 148 // Start algorithm. 149 for(int i = 0; i < numTimesteps; i++) { 150 151 preStepModify(); 152 153 // Save current positions, forces, and energies in case we reject. 154 saveValues(); 155 156 // Repeat until new momenta are accepted. 157 acceptMomenta = false; 158 159 while( !acceptMomenta ) { 160 161 // Clear the history queues and reset beta. 162 resetHistory(); 163 164 // Calculate new random velocities. 165 perturbSystem(); 166 167 // Save the initial total energy for the acceptance step. 168 initTotEnergy = myEnergies->potentialEnergy() + kineticEnergy(myTopo,myVelocities); 169 170 // Run steps -k/2 to -1. 171 runPreSteps( myShadowKover2 ); 172 173 // Save data at step 0. 174 pushShadowHistory(); 175 176 // Run steps 1 to k/2. 177 next()->run( myShadowKover2 ); 178 179 // We now have data for steps -k/2 -> k/2, so we can calculate the 180 // shadow at step 0 and save it as the initial shadow. 181 calculateShadow(); 182 183 initShadowEnergy = (*myEnergies)[ScalarStructure::SHADOW]; 184 185 // Accept these momenta with probability proportional to: 186 // min{ 1, exp^(-\beta(H_[2k] - c - H)) }. 187 momentumEnerDiff = initShadowEnergy - initTotEnergy - myC; 188 189 // Decide whether or not to accept new momenta. 190 if( metropolisTest( momentumEnerDiff, 0.0 ) ) { 191 192 acceptMomenta = true; 193 194 report.precision(8); 195 report << debug(1) << "Take momenta (diff = " << momentumEnerDiff << 196 ")" << endr; 197 198 } 199 else { 200 201 // Restore saved values from step 0. 202 restoreValues(); 203 204 report.precision(8); 205 report << debug(1) << "Deny momenta (diff = " << momentumEnerDiff << 206 ")" << endr; 207 208 } 209 210 } 211 212 // Choose max{ H_[2k] - c, H } to be the initial 'energy'. 213 MDinitEner = std::max(initShadowEnergy - myC, initTotEnergy); 214 215 sumTotalEnergy += initShadowEnergy - initTotEnergy; 216 217 // In between the beginning and end of the cycle, we don't need to 218 // save previous values or calculate the shadow so run using LF. Run 219 // from step k/2 + 1 -> n - k/2 - 1. 220 disableBetaUpdate(); 221 next()->run(localCycleLength - myShadowK - 1); 222 enableBetaUpdate(); 223 224 // Finish by running steps n - k/2 -> n. 225 next()->run( myShadowKover2 + 1 ); 226 227 // Run k/2 more steps so we can compute the shadow at step n. 228 runPostSteps( myShadowKover2 ); 229 230 // Make note of the final total and shadow energies. 231 finTotEnergy = myEnergies->potentialEnergy() + kineticEnergy(myTopo,myVelocities); 232 finShadowEnergy = (*myEnergies)[ScalarStructure::SHADOW]; 233 234 235 // FIXME Temporary. Used for figuring out a good value of c. 236 237 // cout.setf( ios::fixed ); 238 // cout << endl << "DEBUG: dSE: " << setw(10) << setprecision(4) << 239 // finShadowEnergy; 240 // cout << " - " << setw(10) << setprecision(4) << initShadowEnergy; 241 // cout << " = " << setw(10) << setprecision(4) << finShadowEnergy - 242 // initShadowEnergy << endl; 243 // 244 // cout << "DEBUG: dTE: " << setw(10) << setprecision(4) << finTotEnergy; 245 // cout << " - " << setw(10) << setprecision(4) << initTotEnergy; 246 // cout << " = " << setw(10) << setprecision(4) << finTotEnergy - 247 // initTotEnergy << endl; 248 // 249 // cout << "DEBUG: dE: " << setw(10) << setprecision(5) << finShadowEnergy 250 // - finTotEnergy; 251 // cout << " - " << setw(10) << setprecision(5) << initShadowEnergy - 252 // initTotEnergy; 253 // cout << " = " << setw(10) << setprecision(5) << finShadowEnergy - 254 // finTotEnergy - initShadowEnergy + initTotEnergy << endl << endl; 255 256 // Choose max{ H_[2k] - c, H } to be the final 'energy'. 257 MDfinEner = std::max(finShadowEnergy - myC, finTotEnergy); 258 259 // Based on the difference in the max energies, decide to accept 260 // the new positions or not. 261 if( metropolisTest( MDfinEner, MDinitEner ) ) { 262 263 // Save the new positions. 264 saveValues(); 265 266 // Print debug information if desired. 267 report << debug(1) << "Accepted (diff = " << MDfinEner - MDinitEner <<")"<< endr; 268 269 } 270 else { 271 272 // If we rejected, the we have a lot of clean up to do. 273 restoreValues(); 274 275 report << debug(1) << "Rejected (diff = " << MDfinEner - MDinitEner <<")" << endr; 276 277 } 278 279 postStepModify(); 280 281 } 282 283 } 284 285 // --------------------------------------------------------------------- // 286 runPreSteps(int numTimesteps)287 void ShadowHMCIntegrator::runPreSteps(int numTimesteps){ 288 289 // Save all values for current timestep (probably time 0). 290 Vector3DBlock positions(*myPositions); 291 Vector3DBlock velocities(*myVelocities); 292 ScalarStructure energies(*myEnergies); 293 Vector3DBlock forces(*(next()->getForces())); 294 295 // Flip clock to run backwards. 296 backward(); 297 298 next()->run( numTimesteps ); 299 300 // Flip clock to run forward. 301 forward(); 302 303 // Restore original values. 304 myPositions->intoAssign(positions); 305 myVelocities->intoAssign(velocities); 306 (*myEnergies) = energies; 307 next()->getForces()->intoAssign(forces); 308 myBeta = 0.0; 309 310 } 311 312 // --------------------------------------------------------------------- // 313 runPostSteps(int numTimesteps)314 void ShadowHMCIntegrator::runPostSteps(int numTimesteps){ 315 316 // We need to save the current state of the system because the extra steps 317 // we are running are only to compute the shadow. All other values must 318 // remain at step n. 319 320 Vector3DBlock positions(*myPositions); 321 Vector3DBlock velocities(*myVelocities); 322 ScalarStructure energies(*myEnergies); 323 Vector3DBlock forces(*(next()->getForces())); 324 325 // Run the extra steps needed to compute the shadow at step n. 326 next()->run( numTimesteps ); 327 calculateShadow(); 328 329 // Restore the energies, forces, positions, velocities and shadow. 330 Real finalShEnergy = (*myEnergies)[ScalarStructure::SHADOW]; 331 myPositions->intoAssign(positions); 332 myVelocities->intoAssign(velocities); 333 (*myEnergies) = energies; 334 next()->getForces()->intoAssign(forces); 335 (*myEnergies)[ScalarStructure::SHADOW] = finalShEnergy; 336 337 } 338 339 // --------------------------------------------------------------------- // 340 updateBeta()341 void ShadowHMCIntegrator::updateBeta() { 342 Report::report << debug(10) <<"[ShadowHMCIntegrator::updateBeta]"<<Report::endr; 343 // Compute the timestep for a 'half' kick. 344 Real h = 0.5 * next()->getTimestep() * Constant::INV_TIMEFACTOR; 345 346 // Compute the dot product of the positions and forces. 347 Real posDotF = 0.0; 348 const Vector3DBlock& forces(*next()->getForces()); 349 for(unsigned int i = 0; i < myPositions->size(); i++) { 350 posDotF += (*myPositions)[i].dot(forces[i]); 351 } 352 353 // Update beta. 354 myBeta -= h * ( posDotF + 2.0 * myEnergies->potentialEnergy() ); 355 356 } 357 358 // --------------------------------------------------------------------- // 359 pushShadowHistory()360 void ShadowHMCIntegrator::pushShadowHistory(){ 361 362 Report::report<< debug(10) <<"[ShadowHMCIntegrator::pushShadowHistory]"<<Report::endr; 363 364 // When simulating forward in time, add new values to end of queue. 365 if(getTimestep() > 0){ 366 367 myPrevBeta->push_back( myBeta ); 368 myPrevVelocities->push_back( *myVelocities ); 369 myPrevPositions->push_back( *myPositions ); 370 371 } 372 // When simulating backward in time, push new values to front of queue. 373 else { 374 375 myPrevBeta->push_front( myBeta ); 376 myPrevVelocities->push_front( *myVelocities ); 377 myPrevPositions->push_front( *myPositions ); 378 379 } 380 381 // For an order (2k) shadow Hamiltonian, we only need to store (k)+1 // 382 // values. Release old values. // 383 if( myPrevBeta->size() > ( myShadowK + 1 ) ) { 384 385 // The most recent values are added to the end of the queue, so we remove 386 // from the front. 387 if(getTimestep() > 0){ 388 389 myPrevBeta->pop_front(); 390 myPrevVelocities->pop_front(); 391 myPrevPositions->pop_front(); 392 393 } 394 else { 395 396 myPrevBeta->pop_back(); 397 myPrevVelocities->pop_back(); 398 myPrevPositions->pop_back(); 399 400 } 401 402 } 403 404 } 405 406 // --------------------------------------------------------------------- // 407 calculateShadow()408 void ShadowHMCIntegrator::calculateShadow() { 409 410 // Make sure that we have enough saved values to compute the shadow. 411 if( myPrevBeta->size() != ( myShadowK + 1 ) ) { 412 report << warning << "Called calculateShadow() with less than " << 413 ( myShadowK + 1 ) << " values." << endr; 414 return; 415 } 416 417 // Return the appropriate shadow. 418 switch( myOrder ) { 419 case(4) : 420 (*myEnergies)[ScalarStructure::SHADOW] = calcShadow4(); 421 break; 422 case(8) : 423 (*myEnergies)[ScalarStructure::SHADOW] = calcShadow8(); 424 break; 425 default : 426 return; 427 } 428 } 429 430 // --------------------------------------------------------------------- // 431 calcShadow4()432 Real ShadowHMCIntegrator::calcShadow4() { 433 434 Real h = next()->getTimestep() * Constant::INV_TIMEFACTOR; 435 Real pos1DotVel0 = 0.0; 436 Real pos1DotVel2 = 0.0; 437 Real vel1DotPos0 = 0.0; 438 Real vel1DotPos2 = 0.0; 439 Real a10 = 0.0; 440 Real a12 = 0.0; 441 const Real c = 1. / (2 * h); 442 Real mass_div_2 = 0.; 443 444 // The shadow Hamiltonian is a symmetric method about a timestep 'n'. The 445 // way that the data is stored and the fact that myShadowK is 0.5 * the 446 // order of the shadow allows me to access the data in a form consistent 447 // with the algorithm. i.e. data[n-2], data[n-1] ... 448 int n = myShadowKover2; 449 450 // Compute the shadow on my range. 451 for( unsigned int i = 0; i < myPositions->size(); i++ ) { 452 453 // Get the mass for this atom. 454 mass_div_2 = 0.5 * myTopo->atoms[i].scaledMass; 455 456 // Compute the dot products. 457 pos1DotVel0 += mass_div_2 * (((*myPrevPositions)[n+1])[i] - 458 ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n])[i]); 459 460 vel1DotPos0 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] - 461 ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n])[i]); 462 463 pos1DotVel2 += mass_div_2 * (((*myPrevPositions)[n+1])[i] - 464 ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n+1])[i] - 465 ((*myPrevVelocities)[n])[i] * 2.0 + 466 ((*myPrevVelocities)[n-1])[i]); 467 468 vel1DotPos2 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] - 469 ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n+1])[i] - 470 ((*myPrevPositions)[n])[i] * 2.0 + 471 ((*myPrevPositions)[n-1])[i]); 472 473 } 474 475 // Compute the a_i_j. 476 a10 = pos1DotVel0 - vel1DotPos0 - 0.5 * ( (*myPrevBeta)[n+1] - (*myPrevBeta)[n-1] ); 477 a12 = pos1DotVel2 - vel1DotPos2; 478 479 report << debug(2) << "calcShadow4() = " << ( c * ( a10 - a12 / 6.0 ) ) << 480 endr; 481 482 // Calculate the 4th order shadow. 483 return( c * ( a10 - a12 / 6.0 ) ); 484 485 } 486 487 // --------------------------------------------------------------------- // 488 calcShadow8()489 Real ShadowHMCIntegrator::calcShadow8() { 490 Real h = next()->getTimestep() * Constant::INV_TIMEFACTOR; 491 Real pos1DotVel0 = 0.0; 492 Real pos1DotVel2 = 0.0; 493 Real pos1DotVel4 = 0.0; 494 Real vel1DotPos0 = 0.0; 495 Real vel1DotPos2 = 0.0; 496 Real vel1DotPos4 = 0.0; 497 Real pos3DotVel0 = 0.0; 498 Real pos3DotVel2 = 0.0; 499 Real pos3DotVel4 = 0.0; 500 Real vel3DotPos0 = 0.0; 501 Real vel3DotPos2 = 0.0; 502 Real vel3DotPos4 = 0.; 503 Real a10 = 0.0; 504 Real a12 = 0.0; 505 Real a14 = 0.0; 506 Real a30 = 0.0; 507 Real a32 = 0.0; 508 Real a34 = 0.; 509 const Real c = 1. / (2 * h); 510 Real mass_div_2 = 0.; 511 512 // The shadow Hamiltonian is a symmetric method about a timestep 'n'. The 513 // way that the data is stored and the fact that myShadowK is 0.5 * the 514 // order of the shadow allows me to access the data in a form consistent 515 // with the algorithm. i.e. data[n-2], data[n-1] ... 516 int n = myShadowKover2; 517 518 // Compute the shadow on my range. 519 for( unsigned int i = 0; i < myPositions->size(); i++ ) { 520 521 mass_div_2 = 0.5 * myTopo->atoms[i].scaledMass; 522 523 // Compute the dot products. 524 pos1DotVel0 += mass_div_2 * (((*myPrevPositions)[n+1])[i] - ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n])[i]); 525 vel1DotPos0 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] - ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n])[i]); 526 pos1DotVel2 += mass_div_2 * (((*myPrevPositions)[n+1])[i] - 527 ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n+1])[i] - 528 ((*myPrevVelocities)[n])[i] * 2.0 + 529 ((*myPrevVelocities)[n-1])[i]); 530 vel1DotPos2 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] - 531 ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n+1])[i] - 532 ((*myPrevPositions)[n])[i] * 2.0 + 533 ((*myPrevPositions)[n-1])[i]); 534 pos1DotVel4 += mass_div_2 * (((*myPrevPositions)[n+1])[i] - 535 ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n+2])[i] - 536 ((*myPrevVelocities)[n+1])[i] * 4.0 + 537 ((*myPrevVelocities)[n])[i] * 6.0 - 538 ((*myPrevVelocities)[n-1])[i] * 4.0 + 539 ((*myPrevVelocities)[n-2])[i]); 540 vel1DotPos4 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] - 541 ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n+2])[i] - 542 ((*myPrevPositions)[n+1])[i] * 4.0 + 543 ((*myPrevPositions)[n])[i] * 6.0 - 544 ((*myPrevPositions)[n-1])[i] * 4.0 + 545 ((*myPrevPositions)[n-2])[i]); 546 pos3DotVel0 += mass_div_2 * (((*myPrevPositions)[n+2])[i] - 547 ((*myPrevPositions)[n+1])[i] * 2.0 + 548 ((*myPrevPositions)[n-1])[i] * 2.0 - 549 ((*myPrevPositions)[n-2])[i]).dot(((*myPrevVelocities)[n])[i]); 550 vel3DotPos0 += mass_div_2 * (((*myPrevVelocities)[n+2])[i] - 551 ((*myPrevVelocities)[n+1])[i] * 2.0 + 552 ((*myPrevVelocities)[n-1])[i] * 2.0 - 553 ((*myPrevVelocities)[n-2])[i]).dot(((*myPrevPositions)[n])[i]); 554 pos3DotVel2 += mass_div_2 * (((*myPrevPositions)[n+2])[i] - 555 ((*myPrevPositions)[n+1])[i] * 2.0 + 556 ((*myPrevPositions)[n-1])[i] * 2.0 - 557 ((*myPrevPositions)[n-2])[i]).dot(((*myPrevVelocities)[n+1])[i] - 558 ((*myPrevVelocities)[n])[i] * 2.0 + 559 ((*myPrevVelocities)[n-1])[i]); 560 vel3DotPos2 += mass_div_2 * (((*myPrevVelocities)[n+2])[i] - 561 ((*myPrevVelocities)[n+1])[i] * 2.0 + 562 ((*myPrevVelocities)[n-1])[i] * 2.0 - 563 ((*myPrevVelocities)[n-2])[i]).dot(((*myPrevPositions)[n+1])[i] - 564 ((*myPrevPositions)[n])[i] * 2.0 + 565 ((*myPrevPositions)[n-1])[i]); 566 pos3DotVel4 += mass_div_2 * (((*myPrevPositions)[n+2])[i] - 567 ((*myPrevPositions)[n+1])[i] * 2.0 + 568 ((*myPrevPositions)[n-1])[i] * 2.0 - 569 ((*myPrevPositions)[n-2])[i]).dot(((*myPrevVelocities)[n+2])[i] - 570 ((*myPrevVelocities)[n+1])[i] * 4.0 + 571 ((*myPrevVelocities)[n])[i] * 6.0 - 572 ((*myPrevVelocities)[n-1])[i] * 4.0 + 573 ((*myPrevVelocities)[n-2])[i]); 574 vel3DotPos4 += mass_div_2 * (((*myPrevVelocities)[n+2])[i] - 575 ((*myPrevVelocities)[n+1])[i] * 2.0 + 576 ((*myPrevVelocities)[n-1])[i] * 2.0 - 577 ((*myPrevVelocities)[n-2])[i]).dot(((*myPrevPositions)[n+2])[i] - 578 ((*myPrevPositions)[n+1])[i] * 4.0 + 579 ((*myPrevPositions)[n])[i] * 6.0 - 580 ((*myPrevPositions)[n-1])[i] * 4.0 + 581 ((*myPrevPositions)[n-2])[i]); 582 } 583 584 // Compute the a_i_j. 585 a10 = pos1DotVel0 - vel1DotPos0 - 0.5 * ( (*myPrevBeta)[n+1] - (*myPrevBeta)[n-1] ); 586 a12 = pos1DotVel2 - vel1DotPos2; 587 a14 = pos1DotVel4 - vel1DotPos4; 588 a30 = pos3DotVel0 - vel3DotPos0 - 0.5 * (*myPrevBeta)[n+2] + (*myPrevBeta)[n+1] - (*myPrevBeta)[n-1] + 0.5 * (*myPrevBeta)[n-2]; 589 a32 = pos3DotVel2 - vel3DotPos2; 590 a34 = pos3DotVel4 - vel3DotPos4; 591 592 report << debug(2) << "calcShadow8() = " << ( (c / 210.) * (210. * a10 + 25. 593 * a30 + 26. * a32 - (60. * a12 + 19. * a14 + 1.5 * a34)) ) << 594 endr; 595 596 // Calculate the 8th order shadow. I multiplied the algorithm by 210 in 597 // order to easily remove the 5 fractions. I now only have to do 1 598 // division instead of 5. I also grouped the additions and subtractions to 599 // avoid as much round off error as possible. 600 601 return ((c / 210.) * (210. * a10 + 25. * a30 + 26. * a32 - (60. * a12 + 19. * a14 + 1.5 * a34))); 602 } 603 604 // --------------------------------------------------------------------- // 605 doMake(string & errMsg,const vector<Value> & values,ForceGroup * fg,StandardIntegrator * nextIntegrator) const606 MTSIntegrator* ShadowHMCIntegrator::doMake(string& errMsg, const vector<Value>& values, 607 ForceGroup* fg, StandardIntegrator *nextIntegrator)const{ 608 609 unsigned int order = values[3]; 610 if(order != 4 && order != 8){ 611 errMsg += " Expecting order 4 or 8, but not "+values[3].getString(); 612 return NULL; 613 } 614 615 return new ShadowHMCIntegrator(values[0],values[1],values[2],values[3],values[4],fg,nextIntegrator); 616 } 617 618 // --------------------------------------------------------------------- // 619 getParameters(vector<Parameter> & parameters) const620 void ShadowHMCIntegrator::getParameters(vector< Parameter> ¶meters) const { 621 MCIntegrator::getParameters(parameters); 622 parameters.push_back(Parameter("order", Value(myOrder,ConstraintValueType::Positive()))); 623 parameters.push_back(Parameter("c", Value(myC))); 624 } 625 626 // --------------------------------------------------------------------- // 627 perturbSystem()628 void ShadowHMCIntegrator::perturbSystem() { 629 630 report << debug(10) << "[ShadowHMCIntegrator::perturbSystem]"<<endr; 631 632 randomVelocity(getInitialTemperature(),myTopo,myVelocities); 633 634 buildMolecularMomentum(myVelocities,myTopo); 635 636 } 637 638 // --------------------------------------------------------------------- // 639 resetHistory()640 void ShadowHMCIntegrator::resetHistory() { 641 // Reset the queues and beta term. 642 myPrevBeta->clear(); 643 myPrevVelocities->clear(); 644 myPrevPositions->clear(); 645 myBeta = 0.; 646 } 647 648 // --------------------------------------------------------------------- // 649 disableBetaUpdate()650 void ShadowHMCIntegrator::disableBetaUpdate(){ 651 if(myModifier1 != NULL) 652 myModifier1->disable(); 653 if(myModifier2 != NULL) 654 myModifier2->disable(); 655 } 656 657 // --------------------------------------------------------------------- // 658 enableBetaUpdate()659 void ShadowHMCIntegrator::enableBetaUpdate(){ 660 if(myModifier1 != NULL) 661 myModifier1->enable(); 662 if(myModifier2 != NULL) 663 myModifier2->enable(); 664 665 } 666 667 // --------------------------------------------------------------------- // 668 saveValues()669 void ShadowHMCIntegrator::saveValues() { 670 (*myOldForces) = (*(next()->getForces())); 671 MCIntegrator::saveValues(); 672 } 673 674 // --------------------------------------------------------------------- // 675 restoreValues()676 void ShadowHMCIntegrator::restoreValues(){ 677 next()->getForces()->intoAssign(*myOldForces); 678 MCIntegrator::restoreValues(); 679 } 680 681 } 682 683