00001
00007
00008
00009
00010
00011
00012
00013
00014 #include "InfoStream.h"
00015 #include "Node.h"
00016 #include "SimParameters.h"
00017 #include "Sequencer.h"
00018 #include "HomePatch.h"
00019 #include "ReductionMgr.h"
00020 #include "CollectionMgr.h"
00021 #include "BroadcastObject.h"
00022 #include "Output.h"
00023 #include "Controller.h"
00024 #include "Broadcasts.h"
00025 #include "Molecule.h"
00026 #include "NamdOneTools.h"
00027 #include "LdbCoordinator.h"
00028 #include "Thread.h"
00029 #include "Random.h"
00030 #include "PatchMap.inl"
00031 #include "ComputeMgr.h"
00032 #include "ComputeGlobal.h"
00033
00034 #define MIN_DEBUG_LEVEL 4
00035
00036 #include "Debug.h"
00037
00038 Sequencer::Sequencer(HomePatch *p) :
00039 simParams(Node::Object()->simParameters),
00040 patch(p),
00041 collection(CollectionMgr::Object()),
00042 ldbSteps(0)
00043 {
00044 broadcast = new ControllerBroadcasts;
00045 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00046 if (simParams->pressureProfileOn) {
00047 int ntypes = simParams->pressureProfileAtomTypes;
00048 int nslabs = simParams->pressureProfileSlabs;
00049 pressureProfileReduction =
00050 ReductionMgr::Object()->willSubmit(
00051 REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
00052 } else {
00053 pressureProfileReduction = NULL;
00054 }
00055 ldbCoordinator = (LdbCoordinator::Object());
00056 random = new Random(simParams->randomSeed);
00057 random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
00058
00059 rescaleVelocities_numTemps = 0;
00060 berendsenPressure_count = 0;
00061 }
00062
00063 Sequencer::~Sequencer(void)
00064 {
00065 delete broadcast;
00066 delete reduction;
00067 delete pressureProfileReduction;
00068 delete random;
00069 }
00070
00071
00072 void Sequencer::threadRun(Sequencer* arg)
00073 {
00074 arg->algorithm();
00075 }
00076
00077
00078 void Sequencer::run(void)
00079 {
00080
00081 DebugM(4, "::run() - this = " << this << "\n" );
00082 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
00083 CthSetStrategyDefault(thread);
00084 awaken();
00085 }
00086
00087
00088
00089
00090 void Sequencer::algorithm(void)
00091 {
00092 int scriptTask;
00093 int scriptSeq = 0;
00094 while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
00095 switch ( scriptTask ) {
00096 case SCRIPT_OUTPUT:
00097 submitCollections(FILE_OUTPUT);
00098 break;
00099 case SCRIPT_MEASURE:
00100 submitCollections(EVAL_MEASURE);
00101 break;
00102 case SCRIPT_REINITVELS:
00103 reinitVelocities();
00104 break;
00105 case SCRIPT_RESCALEVELS:
00106 rescaleVelocitiesByFactor(simParams->scriptArg1);
00107 break;
00108 case SCRIPT_RELOADCHARGES:
00109 reloadCharges();
00110 break;
00111 case SCRIPT_CHECKPOINT:
00112 patch->checkpoint();
00113 checkpoint_berendsenPressure_count = berendsenPressure_count;
00114 break;
00115 case SCRIPT_REVERT:
00116 patch->revert();
00117 berendsenPressure_count = checkpoint_berendsenPressure_count;
00118 pairlistsAreValid = 0;
00119 break;
00120 case SCRIPT_MINIMIZE:
00121 minimize();
00122 break;
00123 case SCRIPT_RUN:
00124 integrate();
00125 break;
00126 }
00127 }
00128 submitCollections(END_OF_RUN);
00129 terminate();
00130 }
00131
00132
00133 void Sequencer::integrate() {
00134
00135 int &step = patch->flags.step;
00136 step = simParams->firstTimestep;
00137
00138
00139 const Bool rotDragOn = simParams->rotDragOn;
00140 const Bool movDragOn = simParams->movDragOn;
00141
00142 const int commOnly = simParams->commOnly;
00143
00144 int &maxForceUsed = patch->flags.maxForceUsed;
00145 int &maxForceMerged = patch->flags.maxForceMerged;
00146 maxForceUsed = Results::normal;
00147 maxForceMerged = Results::normal;
00148
00149 const int numberOfSteps = simParams->N;
00150 const int stepsPerCycle = simParams->stepsPerCycle;
00151 const BigReal timestep = simParams->dt;
00152
00153
00154 const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00155
00156 const int nonbondedFrequency = simParams->nonbondedFrequency;
00157 slowFreq = nonbondedFrequency;
00158 const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00159 int &doNonbonded = patch->flags.doNonbonded;
00160 doNonbonded = !(step%nonbondedFrequency);
00161 if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00162 if ( doNonbonded ) maxForceUsed = Results::nbond;
00163
00164
00165 const int dofull = ( simParams->fullDirectOn ||
00166 simParams->FMAOn || simParams->PMEOn );
00167 const int fullElectFrequency = simParams->fullElectFrequency;
00168 if ( dofull ) slowFreq = fullElectFrequency;
00169 const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00170 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00171 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00172 if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00173 maxForceMerged = Results::slow;
00174 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00175 int &doMolly = patch->flags.doMolly;
00176 doMolly = simParams->mollyOn && doFullElectrostatics;
00177
00178 int zeroMomentum = simParams->zeroMomentum;
00179
00180
00181 int doTcl = simParams->tclForcesOn;
00182 ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00183
00184
00185 int &doEnergy = patch->flags.doEnergy;
00186 int energyFrequency = simParams->outputEnergies;
00187
00188 rattle1(0.,0);
00189
00190 const int reassignFreq = simParams->reassignFreq;
00191 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00192 reassignVelocities(timestep,step);
00193 }
00194
00195 doEnergy = ! ( step % energyFrequency );
00196 runComputeObjects(1,step<numberOfSteps);
00197
00198
00199 if (simParams->watmodel == WAT_TIP4) {
00200 redistrib_tip4p_forces(Results::normal, 0);
00201 redistrib_tip4p_forces(Results::nbond, 0);
00202 redistrib_tip4p_forces(Results::slow, 0);
00203 }
00204
00205 if ( staleForces || doTcl ) {
00206 if ( doNonbonded ) saveForce(Results::nbond);
00207 if ( doFullElectrostatics ) saveForce(Results::slow);
00208 }
00209 if ( ! commOnly ) {
00210 addForceToMomentum(-0.5*timestep);
00211 if (staleForces || doNonbonded)
00212 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00213 if (staleForces || doFullElectrostatics)
00214 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces);
00215 }
00216 minimizationQuenchVelocity();
00217 rattle1(-timestep,0);
00218 submitHalfstep(step);
00219 if ( ! commOnly ) {
00220 addForceToMomentum(timestep);
00221 if (staleForces || doNonbonded)
00222 addForceToMomentum(nbondstep,Results::nbond,staleForces,0);
00223 if (staleForces || doFullElectrostatics)
00224 addForceToMomentum(slowstep,Results::slow,staleForces,0);
00225 }
00226 rattle1(timestep,1);
00227 if (doTcl)
00228 computeGlobal->saveTotalForces(patch);
00229 submitHalfstep(step);
00230 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00231 if ( ! commOnly ) {
00232 addForceToMomentum(-0.5*timestep);
00233 if (staleForces || doNonbonded)
00234 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00235 if (staleForces || doFullElectrostatics)
00236 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00237 }
00238 submitReductions(step);
00239 rebalanceLoad(step);
00240
00241 for ( ++step; step <= numberOfSteps; ++step )
00242 {
00243
00244 rescaleVelocities(step);
00245 tcoupleVelocities(timestep,step);
00246 berendsenPressure(step);
00247
00248 if ( ! commOnly ) {
00249 addForceToMomentum(0.5*timestep);
00250 if (staleForces || doNonbonded)
00251 addForceToMomentum(0.5*nbondstep,Results::nbond,staleForces,1);
00252 if (staleForces || doFullElectrostatics)
00253 addForceToMomentum(0.5*slowstep,Results::slow,staleForces,1);
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 maximumMove(timestep);
00268 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00269 langevinPiston(step);
00270 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00271
00272 minimizationQuenchVelocity();
00273
00274 doNonbonded = !(step%nonbondedFrequency);
00275 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00276
00277 if ( zeroMomentum && doFullElectrostatics )
00278 correctMomentum(step,slowstep);
00279
00280 submitHalfstep(step);
00281
00282 doMolly = simParams->mollyOn && doFullElectrostatics;
00283
00284 maxForceUsed = Results::normal;
00285 if ( doNonbonded ) maxForceUsed = Results::nbond;
00286 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00287
00288
00289 doEnergy = ! ( step % energyFrequency );
00290 runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00291
00292
00293 if (simParams->watmodel == WAT_TIP4) {
00294 redistrib_tip4p_forces(Results::normal, 1);
00295 if (doNonbonded) redistrib_tip4p_forces(Results::nbond, 1);
00296 if (doFullElectrostatics) redistrib_tip4p_forces(Results::slow, 1);
00297 }
00298
00299 if ( staleForces || doTcl ) {
00300 if ( doNonbonded ) saveForce(Results::nbond);
00301 if ( doFullElectrostatics ) saveForce(Results::slow);
00302 }
00303
00304
00305 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00306 reassignVelocities(timestep,step);
00307 addForceToMomentum(-0.5*timestep);
00308 if (staleForces || doNonbonded)
00309 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00310 if (staleForces || doFullElectrostatics)
00311 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,0);
00312 rattle1(-timestep,0);
00313 }
00314
00315 if ( ! commOnly ) {
00316 langevinVelocitiesBBK1(timestep);
00317 addForceToMomentum(timestep);
00318 if (staleForces || doNonbonded)
00319 addForceToMomentum(nbondstep,Results::nbond,staleForces,1);
00320 if (staleForces || doFullElectrostatics)
00321 addForceToMomentum(slowstep,Results::slow,staleForces,1);
00322 langevinVelocitiesBBK2(timestep);
00323 }
00324
00325
00326 if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00327 if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00328
00329 rattle1(timestep,1);
00330 if (doTcl)
00331 computeGlobal->saveTotalForces(patch);
00332
00333 submitHalfstep(step);
00334 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00335
00336 if ( ! commOnly ) {
00337 addForceToMomentum(-0.5*timestep);
00338 if (staleForces || doNonbonded)
00339 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00340 if (staleForces || doFullElectrostatics)
00341 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00342 }
00343
00344
00345
00346 submitReductions(step);
00347 submitCollections(step);
00348
00349 #if CYCLE_BARRIER
00350 cycleBarrier(!((step+1) % stepsPerCycle), step);
00351 #elif PME_BARRIER
00352 cycleBarrier(doFullElectrostatics, step);
00353 #elif STEP_BARRIER
00354 cycleBarrier(1, step);
00355 #endif
00356
00357 rebalanceLoad(step);
00358
00359 #if PME_BARRIER
00360
00361 cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00362 #endif
00363 }
00364 }
00365
00366
00367 void Sequencer::addMovDragToPosition(BigReal timestep) {
00368 FullAtom *atom = patch->atom.begin();
00369 int numAtoms = patch->numAtoms;
00370 Molecule *molecule = Node::Object()->molecule;
00371 const BigReal movDragGlobVel = simParams->movDragGlobVel;
00372 const BigReal dt = timestep / TIMEFACTOR;
00373 Vector movDragVel, dragIncrement;
00374 for ( int i = 0; i < numAtoms; ++i )
00375 {
00376
00377 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00378 || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00379 molecule->get_movdrag_params(movDragVel, atom[i].id);
00380 dragIncrement = movDragGlobVel * movDragVel * dt;
00381 atom[i].position += dragIncrement;
00382 }
00383 }
00384
00385
00386 void Sequencer::addRotDragToPosition(BigReal timestep) {
00387 FullAtom *atom = patch->atom.begin();
00388 int numAtoms = patch->numAtoms;
00389 Molecule *molecule = Node::Object()->molecule;
00390 const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00391 const BigReal dt = timestep / TIMEFACTOR;
00392 BigReal rotDragVel, dAngle;
00393 Vector atomRadius;
00394 Vector rotDragAxis, rotDragPivot, dragIncrement;
00395 for ( int i = 0; i < numAtoms; ++i )
00396 {
00397
00398 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00399 || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00400 molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00401 dAngle = rotDragGlobVel * rotDragVel * dt;
00402 rotDragAxis /= rotDragAxis.length();
00403 atomRadius = atom[i].position - rotDragPivot;
00404 dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00405 atom[i].position += dragIncrement;
00406 }
00407 }
00408
00409 void Sequencer::minimize() {
00410 const int numberOfSteps = simParams->N;
00411 int &step = patch->flags.step;
00412 step = simParams->firstTimestep;
00413
00414 int &maxForceUsed = patch->flags.maxForceUsed;
00415 int &maxForceMerged = patch->flags.maxForceMerged;
00416 maxForceUsed = Results::normal;
00417 maxForceMerged = Results::normal;
00418 int &doNonbonded = patch->flags.doNonbonded;
00419 doNonbonded = 1;
00420 maxForceUsed = Results::nbond;
00421 maxForceMerged = Results::nbond;
00422 const int dofull = ( simParams->fullDirectOn ||
00423 simParams->FMAOn || simParams->PMEOn );
00424 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00425 doFullElectrostatics = dofull;
00426 if ( dofull ) {
00427 maxForceMerged = Results::slow;
00428 maxForceUsed = Results::slow;
00429 }
00430 int &doMolly = patch->flags.doMolly;
00431 doMolly = simParams->mollyOn && doFullElectrostatics;
00432 int &doEnergy = patch->flags.doEnergy;
00433 doEnergy = 1;
00434
00435 runComputeObjects(1,0);
00436
00437 submitMinimizeReductions(step);
00438 rebalanceLoad(step);
00439
00440 int minSeq = 0;
00441 for ( ++step; step <= numberOfSteps; ++step ) {
00442 BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00443 if ( ! c ) {
00444 c = broadcast->minimizeCoefficient.get(minSeq++);
00445 newMinimizeDirection(c);
00446 c = broadcast->minimizeCoefficient.get(minSeq++);
00447 }
00448 newMinimizePosition(c);
00449
00450 runComputeObjects(1,0);
00451 submitMinimizeReductions(step);
00452 submitCollections(step, 1);
00453 rebalanceLoad(step);
00454 }
00455 quenchVelocities();
00456 }
00457
00458
00459 void Sequencer::newMinimizeDirection(BigReal c) {
00460 FullAtom *a = patch->atom.begin();
00461 Force *f1 = patch->f[Results::normal].begin();
00462 Force *f2 = patch->f[Results::nbond].begin();
00463 Force *f3 = patch->f[Results::slow].begin();
00464 int numAtoms = patch->numAtoms;
00465
00466 for ( int i = 0; i < numAtoms; ++i ) {
00467 a[i].velocity *= c;
00468 a[i].velocity += f1[i] + f2[i] + f3[i];
00469 if ( simParams->fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00470 }
00471 }
00472
00473
00474 void Sequencer::newMinimizePosition(BigReal c) {
00475 FullAtom *a = patch->atom.begin();
00476 int numAtoms = patch->numAtoms;
00477
00478 for ( int i = 0; i < numAtoms; ++i ) {
00479 a[i].position += c * a[i].velocity;
00480 }
00481 }
00482
00483
00484 void Sequencer::quenchVelocities() {
00485 FullAtom *a = patch->atom.begin();
00486 int numAtoms = patch->numAtoms;
00487
00488 for ( int i = 0; i < numAtoms; ++i ) {
00489 a[i].velocity = 0;
00490 }
00491 }
00492
00493 void Sequencer::submitMomentum(int step) {
00494
00495 FullAtom *a = patch->atom.begin();
00496 const int numAtoms = patch->numAtoms;
00497
00498 Vector momentum = 0;
00499 BigReal mass = 0;
00500 if ( simParams->zeroMomentumAlt ) {
00501 for ( int i = 0; i < numAtoms; ++i ) {
00502 momentum += a[i].mass * a[i].velocity;
00503 mass += 1.;
00504 }
00505 } else {
00506 for ( int i = 0; i < numAtoms; ++i ) {
00507 momentum += a[i].mass * a[i].velocity;
00508 mass += a[i].mass;
00509 }
00510 }
00511
00512 ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
00513 reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
00514 }
00515
00516 void Sequencer::correctMomentum(int step, BigReal drifttime) {
00517
00518 if ( simParams->fixedAtomsOn )
00519 NAMD_die("Cannot zero momentum when fixed atoms are present.");
00520
00521 const Vector dv = broadcast->momentumCorrection.get(step);
00522 const Vector dx = dv * ( drifttime / TIMEFACTOR );
00523
00524 FullAtom *a = patch->atom.begin();
00525 const int numAtoms = patch->numAtoms;
00526
00527 if ( simParams->zeroMomentumAlt ) {
00528 for ( int i = 0; i < numAtoms; ++i ) {
00529 BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00530 a[i].velocity += dv * rmass;
00531 a[i].position += dx * rmass;
00532 }
00533 } else {
00534 for ( int i = 0; i < numAtoms; ++i ) {
00535 a[i].velocity += dv;
00536 a[i].position += dx;
00537 }
00538 }
00539
00540 }
00541
00542 void Sequencer::langevinVelocities(BigReal dt_fs)
00543 {
00544 if ( simParams->langevinOn )
00545 {
00546 FullAtom *a = patch->atom.begin();
00547 int numAtoms = patch->numAtoms;
00548 Molecule *molecule = Node::Object()->molecule;
00549 BigReal dt = dt_fs * 0.001;
00550 BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00551
00552 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00553 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00554
00555 for ( int i = 0; i < numAtoms; ++i )
00556 {
00557 int aid = a[i].id;
00558 BigReal f1 = exp( -1. * dt * molecule->langevin_param(aid) );
00559 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
00560 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00561
00562 a[i].velocity *= f1;
00563 a[i].velocity += f2 * random->gaussian_vector();
00564 }
00565 }
00566 }
00567
00568 void Sequencer::langevinVelocitiesBBK1(BigReal dt_fs)
00569 {
00570 if ( simParams->langevinOn )
00571 {
00572 FullAtom *a = patch->atom.begin();
00573 int numAtoms = patch->numAtoms;
00574 Molecule *molecule = Node::Object()->molecule;
00575 BigReal dt = dt_fs * 0.001;
00576 for ( int i = 0; i < numAtoms; ++i )
00577 {
00578 int aid = a[i].id;
00579 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00580 if ( ! dt_gamma ) continue;
00581
00582 a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00583 }
00584 }
00585 }
00586
00587
00588 void Sequencer::langevinVelocitiesBBK2(BigReal dt_fs)
00589 {
00590 if ( simParams->langevinOn )
00591 {
00592 rattle1(dt_fs,1);
00593
00594 FullAtom *a = patch->atom.begin();
00595 int numAtoms = patch->numAtoms;
00596 Molecule *molecule = Node::Object()->molecule;
00597 BigReal dt = dt_fs * 0.001;
00598 BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00599
00600 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00601 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00602
00603 for ( int i = 0; i < numAtoms; ++i )
00604 {
00605 int aid = a[i].id;
00606 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00607 if ( ! dt_gamma ) continue;
00608
00609 a[i].velocity += random->gaussian_vector() *
00610 sqrt( 2 * dt_gamma * kbT *
00611 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00612 a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00613 }
00614 }
00615 }
00616
00617 void Sequencer::berendsenPressure(int step)
00618 {
00619 if ( simParams->berendsenPressureOn ) {
00620 berendsenPressure_count += 1;
00621 const int freq = simParams->berendsenPressureFreq;
00622 if ( ! (berendsenPressure_count % freq ) ) {
00623 berendsenPressure_count = 0;
00624 FullAtom *a = patch->atom.begin();
00625 int numAtoms = patch->numAtoms;
00626 Tensor factor = broadcast->positionRescaleFactor.get(step);
00627 patch->lattice.rescale(factor);
00628 if ( simParams->useGroupPressure )
00629 {
00630 int hgs;
00631 for ( int i = 0; i < numAtoms; i += hgs ) {
00632 int j;
00633 hgs = a[i].hydrogenGroupSize;
00634 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00635 for ( j = i; j < (i+hgs); ++j ) {
00636 a[j].position = patch->lattice.apply_transform(
00637 a[j].fixedPosition,a[j].transform);
00638 }
00639 continue;
00640 }
00641 BigReal m_cm = 0;
00642 Position x_cm(0,0,0);
00643 for ( j = i; j < (i+hgs); ++j ) {
00644 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00645 m_cm += a[j].mass;
00646 x_cm += a[j].mass * a[j].position;
00647 }
00648 x_cm /= m_cm;
00649 Position new_x_cm = x_cm;
00650 patch->lattice.rescale(new_x_cm,factor);
00651 Position delta_x_cm = new_x_cm - x_cm;
00652 for ( j = i; j < (i+hgs); ++j ) {
00653 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00654 a[j].position = patch->lattice.apply_transform(
00655 a[j].fixedPosition,a[j].transform);
00656 continue;
00657 }
00658 a[j].position += delta_x_cm;
00659 }
00660 }
00661 }
00662 else
00663 {
00664 for ( int i = 0; i < numAtoms; ++i )
00665 {
00666 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00667 a[i].position = patch->lattice.apply_transform(
00668 a[i].fixedPosition,a[i].transform);
00669 continue;
00670 }
00671 patch->lattice.rescale(a[i].position,factor);
00672 }
00673 }
00674 }
00675 } else {
00676 berendsenPressure_count = 0;
00677 }
00678 }
00679
00680 void Sequencer::langevinPiston(int step)
00681 {
00682 if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
00683 {
00684 FullAtom *a = patch->atom.begin();
00685 int numAtoms = patch->numAtoms;
00686 Tensor factor = broadcast->positionRescaleFactor.get(step);
00687
00688 Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
00689 patch->lattice.rescale(factor);
00690 Molecule *mol = Node::Object()->molecule;
00691 if ( simParams->useGroupPressure )
00692 {
00693 int hgs;
00694 for ( int i = 0; i < numAtoms; i += hgs ) {
00695 int j;
00696 hgs = a[i].hydrogenGroupSize;
00697 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00698 for ( j = i; j < (i+hgs); ++j ) {
00699 a[j].position = patch->lattice.apply_transform(
00700 a[j].fixedPosition,a[j].transform);
00701 }
00702 continue;
00703 }
00704 BigReal m_cm = 0;
00705 Position x_cm(0,0,0);
00706 Velocity v_cm(0,0,0);
00707 for ( j = i; j < (i+hgs); ++j ) {
00708 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00709 m_cm += a[j].mass;
00710 x_cm += a[j].mass * a[j].position;
00711 v_cm += a[j].mass * a[j].velocity;
00712 }
00713 x_cm /= m_cm;
00714 Position new_x_cm = x_cm;
00715 patch->lattice.rescale(new_x_cm,factor);
00716 Position delta_x_cm = new_x_cm - x_cm;
00717 v_cm /= m_cm;
00718 Velocity delta_v_cm;
00719 delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
00720 delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
00721 delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
00722 for ( j = i; j < (i+hgs); ++j ) {
00723 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00724 a[j].position = patch->lattice.apply_transform(
00725 a[j].fixedPosition,a[j].transform);
00726 continue;
00727 }
00728 if ( mol->is_atom_exPressure(a[j].id) ) continue;
00729 a[j].position += delta_x_cm;
00730 a[j].velocity += delta_v_cm;
00731 }
00732 }
00733 }
00734 else
00735 {
00736 for ( int i = 0; i < numAtoms; ++i )
00737 {
00738 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00739 a[i].position = patch->lattice.apply_transform(
00740 a[i].fixedPosition,a[i].transform);
00741 continue;
00742 }
00743 if ( mol->is_atom_exPressure(a[i].id) ) continue;
00744 patch->lattice.rescale(a[i].position,factor);
00745 a[i].velocity.x *= velFactor.x;
00746 a[i].velocity.y *= velFactor.y;
00747 a[i].velocity.z *= velFactor.z;
00748 }
00749 }
00750 }
00751 }
00752
00753 void Sequencer::rescaleVelocities(int step)
00754 {
00755 const int rescaleFreq = simParams->rescaleFreq;
00756 if ( rescaleFreq > 0 ) {
00757 FullAtom *a = patch->atom.begin();
00758 int numAtoms = patch->numAtoms;
00759 ++rescaleVelocities_numTemps;
00760 if ( rescaleVelocities_numTemps == rescaleFreq ) {
00761 BigReal factor = broadcast->velocityRescaleFactor.get(step);
00762 for ( int i = 0; i < numAtoms; ++i )
00763 {
00764 a[i].velocity *= factor;
00765 }
00766 rescaleVelocities_numTemps = 0;
00767 }
00768 }
00769 }
00770
00771 void Sequencer::reassignVelocities(BigReal timestep, int step)
00772 {
00773 const int reassignFreq = simParams->reassignFreq;
00774 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
00775 FullAtom *a = patch->atom.begin();
00776 int numAtoms = patch->numAtoms;
00777 BigReal newTemp = simParams->reassignTemp;
00778 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
00779 if ( simParams->reassignIncr > 0.0 ) {
00780 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
00781 newTemp = simParams->reassignHold;
00782 } else {
00783 if ( newTemp < simParams->reassignHold )
00784 newTemp = simParams->reassignHold;
00785 }
00786 BigReal kbT = BOLTZMAN * newTemp;
00787
00788 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00789 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00790
00791 for ( int i = 0; i < numAtoms; ++i )
00792 {
00793 a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
00794 sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
00795 * random->gaussian_vector() );
00796 }
00797 } else {
00798 NAMD_bug("Sequencer::reassignVelocities called improperly!");
00799 }
00800 }
00801
00802 void Sequencer::reinitVelocities(void)
00803 {
00804 FullAtom *a = patch->atom.begin();
00805 int numAtoms = patch->numAtoms;
00806 BigReal newTemp = simParams->initialTemp;
00807 BigReal kbT = BOLTZMAN * newTemp;
00808
00809 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00810 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00811
00812 for ( int i = 0; i < numAtoms; ++i )
00813 {
00814 a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
00815 sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
00816 * random->gaussian_vector() );
00817 }
00818 }
00819
00820 void Sequencer::rescaleVelocitiesByFactor(BigReal factor)
00821 {
00822 FullAtom *a = patch->atom.begin();
00823 int numAtoms = patch->numAtoms;
00824 for ( int i = 0; i < numAtoms; ++i )
00825 {
00826 a[i].velocity *= factor;
00827 }
00828 }
00829
00830 void Sequencer::reloadCharges()
00831 {
00832 FullAtom *a = patch->atom.begin();
00833 int numAtoms = patch->numAtoms;
00834 Molecule *molecule = Node::Object()->molecule;
00835 for ( int i = 0; i < numAtoms; ++i )
00836 {
00837 a[i].charge = molecule->atomcharge(a[i].id);
00838 }
00839 }
00840
00841 void Sequencer::tcoupleVelocities(BigReal dt_fs, int step)
00842 {
00843 if ( simParams->tCoupleOn )
00844 {
00845 FullAtom *a = patch->atom.begin();
00846 int numAtoms = patch->numAtoms;
00847 BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
00848 Molecule *molecule = Node::Object()->molecule;
00849 BigReal dt = dt_fs * 0.001;
00850 coefficient *= dt;
00851 for ( int i = 0; i < numAtoms; ++i )
00852 {
00853 int aid = a[i].id;
00854 BigReal f1 = exp( coefficient * molecule->langevin_param(aid) );
00855 a[i].velocity *= f1;
00856 }
00857 }
00858 }
00859
00860 void Sequencer::saveForce(const int ftag)
00861 {
00862 patch->saveForce(ftag);
00863 }
00864
00865 void Sequencer::redistrib_tip4p_forces(const int ftag, const int pressure) {
00866 Tensor virial;
00867 Tensor *vp = ( pressure ? &virial : 0 );
00868 patch->redistrib_tip4p_forces(ftag, vp);
00869 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
00870 }
00871
00872 void Sequencer::addForceToMomentum(BigReal dt, const int ftag,
00873 const int useSaved, int pressure)
00874 {
00875 #if CMK_BLUEGENEL
00876 CmiNetworkProgressAfter (0);
00877 #endif
00878 patch->addForceToMomentum(dt,ftag,useSaved);
00879 }
00880
00881 void Sequencer::addVelocityToPosition(BigReal dt)
00882 {
00883 #if CMK_BLUEGENEL
00884 CmiNetworkProgressAfter (0);
00885 #endif
00886 patch->addVelocityToPosition(dt);
00887 }
00888
00889 void Sequencer::rattle1(BigReal dt, int pressure)
00890 {
00891 if ( simParams->rigidBonds != RIGID_NONE ) {
00892 Tensor virial;
00893 Tensor *vp = ( pressure ? &virial : 0 );
00894 if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
00895 iout << iERROR <<
00896 "Constraint failure; simulation has become unstable.\n" << endi;
00897 Node::Object()->enableEarlyExit();
00898 terminate();
00899 }
00900 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
00901 }
00902 }
00903
00904 void Sequencer::rattle2(BigReal dt, int step)
00905 {
00906 if ( simParams->rigidBonds != RIGID_NONE ) {
00907 Tensor virial;
00908 patch->rattle2(dt, &virial);
00909 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
00910
00911 #ifdef ALTVIRIAL
00912 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
00913 #endif
00914 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,virial);
00915 }
00916 }
00917
00918 void Sequencer::maximumMove(BigReal timestep)
00919 {
00920 FullAtom *a = patch->atom.begin();
00921 int numAtoms = patch->numAtoms;
00922 if ( simParams->maximumMove ) {
00923 const BigReal dt = timestep / TIMEFACTOR;
00924 const BigReal maxvel = simParams->maximumMove / dt;
00925 const BigReal maxvel2 = maxvel * maxvel;
00926 for ( int i=0; i<numAtoms; ++i ) {
00927 if ( a[i].velocity.length2() > maxvel2 ) {
00928 a[i].velocity *= ( maxvel / a[i].velocity.length() );
00929 }
00930 }
00931 } else {
00932 const BigReal dt = timestep / TIMEFACTOR;
00933 const BigReal maxvel = simParams->cutoff / dt;
00934 const BigReal maxvel2 = maxvel * maxvel;
00935 int killme = 0;
00936 for ( int i=0; i<numAtoms; ++i ) {
00937 killme = killme || ( a[i].velocity.length2() > maxvel2 );
00938 }
00939 if ( killme ) {
00940 for ( int i=0; i<numAtoms; ++i ) {
00941 if ( a[i].velocity.length2() > maxvel2 ) {
00942 iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
00943 << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
00944 << ( PDBVELFACTOR * maxvel ) << ")\n" << endi;
00945 }
00946 }
00947 iout << iERROR <<
00948 "Atoms moving too fast; simulation has become unstable.\n" << endi;
00949 Node::Object()->enableEarlyExit();
00950 terminate();
00951 }
00952 }
00953 }
00954
00955 void Sequencer::minimizationQuenchVelocity(void)
00956 {
00957 if ( simParams->minimizeOn ) {
00958 FullAtom *a = patch->atom.begin();
00959 int numAtoms = patch->numAtoms;
00960 for ( int i=0; i<numAtoms; ++i ) {
00961 a[i].velocity = 0.;
00962 }
00963 }
00964 }
00965
00966 void Sequencer::submitHalfstep(int step)
00967 {
00968
00969
00970 FullAtom *a = patch->atom.begin();
00971 int numAtoms = patch->numAtoms;
00972
00973 #if CMK_BLUEGENEL
00974 CmiNetworkProgressAfter (0);
00975 #endif
00976
00977 {
00978 BigReal kineticEnergy = 0;
00979 Tensor virial;
00980 if ( simParams->pairInteractionOn ) {
00981 if ( simParams->pairInteractionSelf ) {
00982 for ( int i = 0; i < numAtoms; ++i ) {
00983 if ( a[i].partition != 1 ) continue;
00984 kineticEnergy += a[i].mass * a[i].velocity.length2();
00985 virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
00986 }
00987 }
00988 } else {
00989 for ( int i = 0; i < numAtoms; ++i ) {
00990 if (a[i].mass < 0.01) continue;
00991 kineticEnergy += a[i].mass * a[i].velocity.length2();
00992 virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
00993 }
00994 }
00995
00996 kineticEnergy *= 0.5 * 0.5;
00997 reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
00998 virial *= 0.5;
00999 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01000 #ifdef ALTVIRIAL
01001 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01002 #endif
01003 }
01004
01005 if (pressureProfileReduction) {
01006 int nslabs = simParams->pressureProfileSlabs;
01007 const Lattice &lattice = patch->lattice;
01008 BigReal idz = nslabs/lattice.c().z;
01009 BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01010 int useGroupPressure = simParams->useGroupPressure;
01011
01012
01013
01014
01015
01016
01017
01018
01019 int hgs;
01020 for (int i=0; i<numAtoms; i += hgs) {
01021 int j, ppoffset;
01022 hgs = a[i].hydrogenGroupSize;
01023 int partition = a[i].partition;
01024
01025 BigReal m_cm = 0;
01026 Velocity v_cm(0,0,0);
01027 for (j=i; j< i+hgs; ++j) {
01028 m_cm += a[j].mass;
01029 v_cm += a[j].mass * a[j].velocity;
01030 }
01031 v_cm /= m_cm;
01032 for (j=i; j < i+hgs; ++j) {
01033 BigReal mass = a[j].mass;
01034 if (! (useGroupPressure && j != i)) {
01035 BigReal z = a[j].position.z;
01036 int slab = (int)floor((z-zmin)*idz);
01037 if (slab < 0) slab += nslabs;
01038 else if (slab >= nslabs) slab -= nslabs;
01039 ppoffset = 3*(slab + partition*nslabs);
01040 }
01041 BigReal wxx, wyy, wzz;
01042 if (useGroupPressure) {
01043 wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
01044 wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
01045 wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
01046 } else {
01047 wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
01048 wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
01049 wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
01050 }
01051 pressureProfileReduction->item(ppoffset ) += wxx;
01052 pressureProfileReduction->item(ppoffset+1) += wyy;
01053 pressureProfileReduction->item(ppoffset+2) += wzz;
01054 }
01055 }
01056 }
01057
01058 {
01059 BigReal intKineticEnergy = 0;
01060 Tensor intVirialNormal;
01061
01062 int hgs;
01063 for ( int i = 0; i < numAtoms; i += hgs ) {
01064
01065 #if CMK_BLUEGENEL
01066 CmiNetworkProgress ();
01067 #endif
01068
01069 hgs = a[i].hydrogenGroupSize;
01070 int j;
01071 BigReal m_cm = 0;
01072 Velocity v_cm(0,0,0);
01073 for ( j = i; j < (i+hgs); ++j ) {
01074 m_cm += a[j].mass;
01075 v_cm += a[j].mass * a[j].velocity;
01076 }
01077 v_cm /= m_cm;
01078 if ( simParams->pairInteractionOn ) {
01079 if ( simParams->pairInteractionSelf ) {
01080 for ( j = i; j < (i+hgs); ++j ) {
01081 if ( a[j].partition != 1 ) continue;
01082 BigReal mass = a[j].mass;
01083 Vector v = a[j].velocity;
01084 Vector dv = v - v_cm;
01085 intKineticEnergy += mass * (v * dv);
01086 intVirialNormal.outerAdd (mass, v, dv);
01087 }
01088 }
01089 } else {
01090 for ( j = i; j < (i+hgs); ++j ) {
01091 BigReal mass = a[j].mass;
01092 Vector v = a[j].velocity;
01093 Vector dv = v - v_cm;
01094 intKineticEnergy += mass * (v * dv);
01095 intVirialNormal.outerAdd(mass, v, dv);
01096 }
01097 }
01098 }
01099
01100 intKineticEnergy *= 0.5 * 0.5;
01101 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
01102 intVirialNormal *= 0.5;
01103 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01104 }
01105
01106 }
01107
01108 void Sequencer::submitReductions(int step)
01109 {
01110 FullAtom *a = patch->atom.begin();
01111 int numAtoms = patch->numAtoms;
01112
01113 #if CMK_BLUEGENEL
01114 CmiNetworkProgressAfter(0);
01115 #endif
01116
01117 reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01118 reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
01119
01120 {
01121 BigReal kineticEnergy = 0;
01122 Vector momentum = 0;
01123 Vector angularMomentum = 0;
01124 Vector o = patch->lattice.origin();
01125 if ( simParams->pairInteractionOn ) {
01126 if ( simParams->pairInteractionSelf ) {
01127 for ( int i = 0; i < numAtoms; ++i ) {
01128 if ( a[i].partition != 1 ) continue;
01129 kineticEnergy += a[i].mass * a[i].velocity.length2();
01130 momentum += a[i].mass * a[i].velocity;
01131 angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01132 }
01133 }
01134 } else {
01135 for ( int i = 0; i < numAtoms; ++i ) {
01136 kineticEnergy += a[i].mass * a[i].velocity.length2();
01137 momentum += a[i].mass * a[i].velocity;
01138 angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01139 }
01140 }
01141
01142 kineticEnergy *= 0.5;
01143 reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
01144 ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
01145 ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
01146 }
01147
01148 #ifdef ALTVIRIAL
01149
01150 {
01151 Tensor altVirial;
01152 for ( int i = 0; i < numAtoms; ++i ) {
01153 altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
01154 }
01155 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
01156 }
01157 {
01158 Tensor altVirial;
01159 for ( int i = 0; i < numAtoms; ++i ) {
01160 altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
01161 }
01162 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
01163 }
01164 {
01165 Tensor altVirial;
01166 for ( int i = 0; i < numAtoms; ++i ) {
01167 altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
01168 }
01169 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
01170 }
01171 #endif
01172
01173 {
01174 BigReal intKineticEnergy = 0;
01175 Tensor intVirialNormal;
01176 Tensor intVirialNbond;
01177 Tensor intVirialSlow;
01178
01179 int hgs;
01180 for ( int i = 0; i < numAtoms; i += hgs ) {
01181 #if CMK_BLUEGENEL
01182 CmiNetworkProgress();
01183 #endif
01184 hgs = a[i].hydrogenGroupSize;
01185 int j;
01186 BigReal m_cm = 0;
01187 Position x_cm(0,0,0);
01188 Velocity v_cm(0,0,0);
01189 for ( j = i; j < (i+hgs); ++j ) {
01190 m_cm += a[j].mass;
01191 x_cm += a[j].mass * a[j].position;
01192 v_cm += a[j].mass * a[j].velocity;
01193 }
01194 x_cm /= m_cm;
01195 v_cm /= m_cm;
01196 int fixedAtomsOn = simParams->fixedAtomsOn;
01197 if ( simParams->pairInteractionOn ) {
01198 int pairInteractionSelf = simParams->pairInteractionSelf;
01199 for ( j = i; j < (i+hgs); ++j ) {
01200 if ( a[j].partition != 1 &&
01201 ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
01202
01203 if ( fixedAtomsOn && a[j].atomFixed ) continue;
01204 BigReal mass = a[j].mass;
01205 Vector v = a[j].velocity;
01206 Vector dv = v - v_cm;
01207 intKineticEnergy += mass * (v * dv);
01208 Vector dx = a[j].position - x_cm;
01209 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01210 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01211 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01212 }
01213 } else {
01214 for ( j = i; j < (i+hgs); ++j ) {
01215
01216 if ( fixedAtomsOn && a[j].atomFixed ) continue;
01217 BigReal mass = a[j].mass;
01218 Vector v = a[j].velocity;
01219 Vector dv = v - v_cm;
01220 intKineticEnergy += mass * (v * dv);
01221 Vector dx = a[j].position - x_cm;
01222 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01223 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01224 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01225 }
01226 }
01227 }
01228
01229 intKineticEnergy *= 0.5;
01230 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
01231 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01232 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
01233 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
01234 }
01235
01236 if (pressureProfileReduction && simParams->useGroupPressure) {
01237
01238 int nslabs = simParams->pressureProfileSlabs;
01239 const Lattice &lattice = patch->lattice;
01240 BigReal idz = nslabs/lattice.c().z;
01241 BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01242 int useGroupPressure = simParams->useGroupPressure;
01243
01244 int hgs;
01245 for (int i=0; i<numAtoms; i += hgs) {
01246 int j;
01247 hgs = a[i].hydrogenGroupSize;
01248 BigReal m_cm = 0;
01249 Position x_cm(0,0,0);
01250 for (j=i; j< i+hgs; ++j) {
01251 m_cm += a[j].mass;
01252 x_cm += a[j].mass * a[j].position;
01253 }
01254 x_cm /= m_cm;
01255
01256 BigReal z = a[i].