#include <Sequencer.h>
|
|
|
Definition at line 72 of file Sequencer.C. 00073 {
00074 delete broadcast;
00075 delete reduction;
00076 delete pressureProfileReduction;
00077 delete random;
00078 }
|
|
|
Definition at line 1099 of file Sequencer.C. References ControllerBroadcasts::adaptTemperature, SimParameters::adaptTempFirstStep, SimParameters::adaptTempFreq, SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, adaptTempT, broadcast, SimParameters::firstTimestep, SimpleBroadcastObject< T >::get(), SimParameters::langevinOn, SimParameters::langevinTemp, and simParams. Referenced by integrate(). 01100 {
01101 //check if adaptive tempering is enabled and in the right timestep range
01102 if (!simParams->adaptTempOn) return;
01103 if ( (step < simParams->adaptTempFirstStep ) ||
01104 ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) {
01105 if (simParams->langevinOn) // restore langevin temperature
01106 adaptTempT = simParams->langevinTemp;
01107 return;
01108 }
01109 // Get Updated Temperature
01110 if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
01111 adaptTempT = broadcast->adaptTemperature.get(step);
01112 }
|
|
||||||||||||||||||||
|
Definition at line 1207 of file Sequencer.C. References HomePatch::addForceToMomentum(), and patch. Referenced by integrate(). 01209 {
01210 #if CMK_BLUEGENEL
01211 CmiNetworkProgressAfter (0);
01212 #endif
01213 patch->addForceToMomentum(dt,ftag,useSaved);
01214 }
|
|
|
Definition at line 468 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, SimParameters::fixedAtomsOn, Molecule::get_movdrag_params(), Molecule::is_atom_movdragged(), Node::molecule, SimParameters::movDragGlobVel, Patch::numAtoms, Node::Object(), patch, CompAtom::position, and simParams. Referenced by integrate(). 00468 {
00469 FullAtom *atom = patch->atom.begin();
00470 int numAtoms = patch->numAtoms;
00471 Molecule *molecule = Node::Object()->molecule; // need its methods
00472 const BigReal movDragGlobVel = simParams->movDragGlobVel;
00473 const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
00474 Vector movDragVel, dragIncrement;
00475 for ( int i = 0; i < numAtoms; ++i )
00476 {
00477 // skip if fixed atom or zero drag attribute
00478 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00479 || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00480 molecule->get_movdrag_params(movDragVel, atom[i].id);
00481 dragIncrement = movDragGlobVel * movDragVel * dt;
00482 atom[i].position += dragIncrement;
00483 }
00484 }
|
|
|
Definition at line 487 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, SimParameters::fixedAtomsOn, Molecule::get_rotdrag_params(), Molecule::is_atom_rotdragged(), Vector::length(), Node::molecule, Patch::numAtoms, Node::Object(), patch, CompAtom::position, SimParameters::rotDragGlobVel, and simParams. Referenced by integrate(). 00487 {
00488 FullAtom *atom = patch->atom.begin();
00489 int numAtoms = patch->numAtoms;
00490 Molecule *molecule = Node::Object()->molecule; // need its methods
00491 const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00492 const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
00493 BigReal rotDragVel, dAngle;
00494 Vector atomRadius;
00495 Vector rotDragAxis, rotDragPivot, dragIncrement;
00496 for ( int i = 0; i < numAtoms; ++i )
00497 {
00498 // skip if fixed atom or zero drag attribute
00499 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00500 || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00501 molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00502 dAngle = rotDragGlobVel * rotDragVel * dt;
00503 rotDragAxis /= rotDragAxis.length();
00504 atomRadius = atom[i].position - rotDragPivot;
00505 dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00506 atom[i].position += dragIncrement;
00507 }
00508 }
|
|
|
Definition at line 1216 of file Sequencer.C. References HomePatch::addVelocityToPosition(), and patch. Referenced by integrate(). 01217 {
01218 #if CMK_BLUEGENEL
01219 CmiNetworkProgressAfter (0);
01220 #endif
01221 patch->addVelocityToPosition(dt);
01222 }
|
|
|
|
Definition at line 29 of file Sequencer.h. References PRIORITY_SIZE. Referenced by LdbCoordinator::awakenSequencers(), HomePatch::boxClosed(), HomePatch::depositMigration(), HomePatch::receiveResult(), and run(). 00029 {
00030 CthAwakenPrio(thread, CK_QUEUEING_IFIFO, PRIORITY_SIZE, &priority);
00031 }
|
|
|
Definition at line 903 of file Sequencer.C. References Lattice::apply_transform(), HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), berendsenPressure_count, SimParameters::berendsenPressureFreq, SimParameters::berendsenPressureOn, BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), CompAtomExt::groupFixed, CompAtom::hydrogenGroupSize, j, Patch::lattice, FullAtom::mass, Patch::numAtoms, patch, Position, CompAtom::position, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, and SimParameters::useGroupPressure. Referenced by integrate(). 00904 {
00905 if ( simParams->berendsenPressureOn ) {
00906 berendsenPressure_count += 1;
00907 const int freq = simParams->berendsenPressureFreq;
00908 if ( ! (berendsenPressure_count % freq ) ) {
00909 berendsenPressure_count = 0;
00910 FullAtom *a = patch->atom.begin();
00911 int numAtoms = patch->numAtoms;
00912 Tensor factor = broadcast->positionRescaleFactor.get(step);
00913 patch->lattice.rescale(factor);
00914 if ( simParams->useGroupPressure )
00915 {
00916 int hgs;
00917 for ( int i = 0; i < numAtoms; i += hgs ) {
00918 int j;
00919 hgs = a[i].hydrogenGroupSize;
00920 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00921 for ( j = i; j < (i+hgs); ++j ) {
00922 a[j].position = patch->lattice.apply_transform(
00923 a[j].fixedPosition,a[j].transform);
00924 }
00925 continue;
00926 }
00927 BigReal m_cm = 0;
00928 Position x_cm(0,0,0);
00929 for ( j = i; j < (i+hgs); ++j ) {
00930 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00931 m_cm += a[j].mass;
00932 x_cm += a[j].mass * a[j].position;
00933 }
00934 x_cm /= m_cm;
00935 Position new_x_cm = x_cm;
00936 patch->lattice.rescale(new_x_cm,factor);
00937 Position delta_x_cm = new_x_cm - x_cm;
00938 for ( j = i; j < (i+hgs); ++j ) {
00939 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00940 a[j].position = patch->lattice.apply_transform(
00941 a[j].fixedPosition,a[j].transform);
00942 continue;
00943 }
00944 a[j].position += delta_x_cm;
00945 }
00946 }
00947 }
00948 else
00949 {
00950 for ( int i = 0; i < numAtoms; ++i )
00951 {
00952 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00953 a[i].position = patch->lattice.apply_transform(
00954 a[i].fixedPosition,a[i].transform);
00955 continue;
00956 }
00957 patch->lattice.rescale(a[i].position,factor);
00958 }
00959 }
00960 }
00961 } else {
00962 berendsenPressure_count = 0;
00963 }
00964 }
|
|
||||||||||||
|
Definition at line 689 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), FullAtom::mass, ControllerBroadcasts::momentumCorrection, NAMD_die(), Patch::numAtoms, patch, CompAtom::position, simParams, FullAtom::velocity, and SimParameters::zeroMomentumAlt. Referenced by integrate(). 00689 {
00690
00691 if ( simParams->fixedAtomsOn )
00692 NAMD_die("Cannot zero momentum when fixed atoms are present.");
00693
00694 const Vector dv = broadcast->momentumCorrection.get(step);
00695 const Vector dx = dv * ( drifttime / TIMEFACTOR );
00696
00697 FullAtom *a = patch->atom.begin();
00698 const int numAtoms = patch->numAtoms;
00699
00700 if ( simParams->zeroMomentumAlt ) {
00701 for ( int i = 0; i < numAtoms; ++i ) {
00702 BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00703 a[i].velocity += dv * rmass;
00704 a[i].position += dx * rmass;
00705 }
00706 } else {
00707 for ( int i = 0; i < numAtoms; ++i ) {
00708 a[i].velocity += dv;
00709 a[i].position += dx;
00710 }
00711 }
00712
00713 }
|
|
||||||||||||
|
Definition at line 1986 of file Sequencer.C. References broadcast. Referenced by integrate(). 01986 {
01987 #if USE_BARRIER
01988 if (doBarrier)
01989 broadcast->cycleBarrier.get(step);
01990 #endif
01991 }
|
|
|
Definition at line 148 of file Sequencer.C. References SimParameters::accelMDdihe, SimParameters::accelMDdual, SimParameters::accelMDOn, SimParameters::adaptTempOn, adaptTempT, adaptTempUpdate(), addForceToMomentum(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), HomePatch::atom, Patch::atomMapper, ResizeArray< Elem >::begin(), berendsenPressure(), BigReal, Bool, SimParameters::colvarsOn, SimParameters::commOnly, ComputeMgr::computeGlobalObject, Node::computeMgr, correctMomentum(), cycleBarrier(), Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, Flags::doNonbonded, SimParameters::dt, ResizeArray< Elem >::end(), eventEndOfTimeStep, SimParameters::firstTimestep, Patch::flags, SimParameters::fullElectFrequency, SimParameters::GBISOn, SimParameters::initialTemp, SimParameters::langevinOn, langevinPiston(), SimParameters::langevinTemp, langevinVelocitiesBBK1(), langevinVelocitiesBBK2(), SimParameters::LCPOOn, HomePatch::ldObjHandle, SimParameters::lonepairs, SimParameters::loweAndersenOn, Flags::maxForceMerged, Flags::maxForceUsed, maximumMove(), minimizationQuenchVelocity(), SimParameters::mollyOn, SimParameters::movDragOn, SimParameters::MTSAlgorithm, SimParameters::N, SimParameters::nonbondedFrequency, SimParameters::numTraceSteps, LdbCoordinator::Object(), Node::Object(), SimParameters::outputEnergies, patch, Patch::patchID, LdbCoordinator::pauseWork(), rattle1(), SimParameters::reassignFreq, reassignVelocities(), rebalanceLoad(), AtomMapper::registerIDsFullAtom(), rescaleaccelMD(), SimParameters::rescaleFreq, SimParameters::rescaleTemp, rescaleVelocities(), SimParameters::rotDragOn, runComputeObjects(), saveForce(), ComputeGlobal::saveTotalForces(), simParams, slowFreq, Node::specialTracing, LdbCoordinator::startWork(), SimParameters::statsOn, Flags::step, SimParameters::stepsPerCycle, submitCollections(), submitHalfstep(), submitMomentum(), submitReductions(), SimParameters::tclForcesOn, tcoupleVelocities(), traceBarrier(), SimParameters::traceStartStep, and SimParameters::zeroMomentum. Referenced by algorithm(). 00148 {
00149 char traceNote[24];
00150 char tracePrefix[20];
00151 sprintf(tracePrefix, "p:%d,s:",patch->patchID);
00152 // patch->write_tip4_props();
00153
00154 int &step = patch->flags.step;
00155 step = simParams->firstTimestep;
00156
00157 // drag switches
00158 const Bool rotDragOn = simParams->rotDragOn;
00159 const Bool movDragOn = simParams->movDragOn;
00160
00161 const int commOnly = simParams->commOnly;
00162
00163 int &maxForceUsed = patch->flags.maxForceUsed;
00164 int &maxForceMerged = patch->flags.maxForceMerged;
00165 maxForceUsed = Results::normal;
00166 maxForceMerged = Results::normal;
00167
00168 const int numberOfSteps = simParams->N;
00169 const int stepsPerCycle = simParams->stepsPerCycle;
00170 const BigReal timestep = simParams->dt;
00171
00172 // what MTS method?
00173 const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00174
00175 const int nonbondedFrequency = simParams->nonbondedFrequency;
00176 slowFreq = nonbondedFrequency;
00177 const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00178 int &doNonbonded = patch->flags.doNonbonded;
00179 doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
00180 if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00181 if ( doNonbonded ) maxForceUsed = Results::nbond;
00182
00183 // Do we do full electrostatics?
00184 const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00185 const int fullElectFrequency = simParams->fullElectFrequency;
00186 if ( dofull ) slowFreq = fullElectFrequency;
00187 const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00188 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00189 doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
00190 if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00191 maxForceMerged = Results::slow;
00192 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00193
00194 const Bool accelMDOn = simParams->accelMDOn;
00195 const Bool accelMDdihe = simParams->accelMDdihe;
00196 const Bool accelMDdual = simParams->accelMDdual;
00197 if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
00198
00199 // Is adaptive tempering on?
00200 const Bool adaptTempOn = simParams->adaptTempOn;
00201 adaptTempT = simParams->initialTemp;
00202 if (simParams->langevinOn)
00203 adaptTempT = simParams->langevinTemp;
00204 else if (simParams->rescaleFreq > 0)
00205 adaptTempT = simParams->rescaleTemp;
00206
00207
00208 int &doMolly = patch->flags.doMolly;
00209 doMolly = simParams->mollyOn && doFullElectrostatics;
00210 // BEGIN LA
00211 int &doLoweAndersen = patch->flags.doLoweAndersen;
00212 doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00213 // END LA
00214
00215 int &doGBIS = patch->flags.doGBIS;
00216 doGBIS = simParams->GBISOn;
00217
00218 int &doLCPO = patch->flags.doLCPO;
00219 doLCPO = simParams->LCPOOn;
00220
00221 int zeroMomentum = simParams->zeroMomentum;
00222
00223 // Do we need to return forces to TCL script or Colvar module?
00224 int doTcl = simParams->tclForcesOn;
00225 int doColvars = simParams->colvarsOn;
00226 ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00227
00228 // Bother to calculate energies?
00229 int &doEnergy = patch->flags.doEnergy;
00230 int energyFrequency = simParams->outputEnergies;
00231
00232 // printf("Doing initial rattle\n");
00233 rattle1(0.,0); // enforce rigid bond constraints on initial positions
00234
00235 if (simParams->lonepairs) {
00236 patch->atomMapper->registerIDsFullAtom(
00237 patch->atom.begin(),patch->atom.end());
00238 }
00239
00240 const int reassignFreq = simParams->reassignFreq;
00241 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00242 reassignVelocities(timestep,step);
00243 }
00244
00245 doEnergy = ! ( step % energyFrequency );
00246 if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00247 //Update energy every timestep for adaptive tempering
00248 if ( adaptTempOn ) doEnergy=1;
00249 runComputeObjects(1,step<numberOfSteps); // must migrate here!
00250 rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
00251 adaptTempUpdate(step); // update adaptive tempering temperature
00252
00253 if ( staleForces || doTcl || doColvars ) {
00254 if ( doNonbonded ) saveForce(Results::nbond);
00255 if ( doFullElectrostatics ) saveForce(Results::slow);
00256 }
00257 if ( ! commOnly ) {
00258 addForceToMomentum(-0.5*timestep);
00259 if (staleForces || doNonbonded)
00260 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00261 if (staleForces || doFullElectrostatics)
00262 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces);
00263 }
00264 minimizationQuenchVelocity();
00265 rattle1(-timestep,0);
00266 submitHalfstep(step);
00267 if ( ! commOnly ) {
00268 addForceToMomentum(timestep);
00269 if (staleForces || doNonbonded)
00270 addForceToMomentum(nbondstep,Results::nbond,staleForces,0);
00271 if (staleForces || doFullElectrostatics)
00272 addForceToMomentum(slowstep,Results::slow,staleForces,0);
00273 }
00274 rattle1(timestep,1);
00275 if (doTcl || doColvars) // include constraint forces
00276 computeGlobal->saveTotalForces(patch);
00277 submitHalfstep(step);
00278 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00279 if ( ! commOnly ) {
00280 addForceToMomentum(-0.5*timestep);
00281 if (staleForces || doNonbonded)
00282 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00283 if (staleForces || doFullElectrostatics)
00284 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00285 }
00286 submitReductions(step);
00287 if(traceIsOn()){
00288 traceUserEvent(eventEndOfTimeStep);
00289 sprintf(traceNote, "%s%d",tracePrefix,step);
00290 traceUserSuppliedNote(traceNote);
00291 }
00292 rebalanceLoad(step);
00293
00294 for ( ++step; step <= numberOfSteps; ++step )
00295 {
00296 #ifndef NAMD_CUDA
00297 LdbCoordinator::Object()->startWork(patch->ldObjHandle);
00298 #endif
00299 rescaleVelocities(step);
00300 tcoupleVelocities(timestep,step);
00301 berendsenPressure(step);
00302
00303 if ( ! commOnly ) {
00304 addForceToMomentum(0.5*timestep);
00305 if (staleForces || doNonbonded)
00306 addForceToMomentum(0.5*nbondstep,Results::nbond,staleForces,1);
00307 if (staleForces || doFullElectrostatics)
00308 addForceToMomentum(0.5*slowstep,Results::slow,staleForces,1);
00309 }
00310
00311 /* reassignment based on half-step velocities
00312 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00313 addVelocityToPosition(0.5*timestep);
00314 reassignVelocities(timestep,step);
00315 addVelocityToPosition(0.5*timestep);
00316 rattle1(0.,0);
00317 rattle1(-timestep,0);
00318 addVelocityToPosition(-1.0*timestep);
00319 rattle1(timestep,0);
00320 } */
00321
00322 maximumMove(timestep);
00323 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00324 langevinPiston(step);
00325 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00326
00327 minimizationQuenchVelocity();
00328
00329 doNonbonded = !(step%nonbondedFrequency);
00330 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00331
00332 if ( zeroMomentum && doFullElectrostatics )
00333 correctMomentum(step,slowstep);
00334
00335 submitHalfstep(step);
00336
00337 doMolly = simParams->mollyOn && doFullElectrostatics;
00338 // BEGIN LA
00339 doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00340 // END LA
00341
00342 maxForceUsed = Results::normal;
00343 if ( doNonbonded ) maxForceUsed = Results::nbond;
00344 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00345 if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
00346
00347 // Migrate Atoms on stepsPerCycle
00348 doEnergy = ! ( step % energyFrequency );
00349 if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00350 if ( adaptTempOn ) doEnergy=1;
00351 #ifndef NAMD_CUDA
00352 LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
00353 #endif
00354 runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00355 #ifndef NAMD_CUDA
00356 LdbCoordinator::Object()->startWork(patch->ldObjHandle);
00357 #endif
00358
00359 rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
00360
00361 if ( staleForces || doTcl || doColvars ) {
00362 if ( doNonbonded ) saveForce(Results::nbond);
00363 if ( doFullElectrostatics ) saveForce(Results::slow);
00364 }
00365
00366 // reassignment based on full-step velocities
00367 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00368 reassignVelocities(timestep,step);
00369 addForceToMomentum(-0.5*timestep);
00370 if (staleForces || doNonbonded)
00371 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00372 if (staleForces || doFullElectrostatics)
00373 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,0);
00374 rattle1(-timestep,0);
00375 }
00376
00377 if ( ! commOnly ) {
00378 langevinVelocitiesBBK1(timestep);
00379 addForceToMomentum(timestep);
00380 if (staleForces || doNonbonded) {
00381 addForceToMomentum(nbondstep,Results::nbond,staleForces,1);
00382 }
00383 if (staleForces || doFullElectrostatics) {
00384 addForceToMomentum(slowstep,Results::slow,staleForces,1);
00385 }
00386 langevinVelocitiesBBK2(timestep);
00387 }
00388
00389 // add drag to each atom's positions
00390 if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00391 if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00392
00393 rattle1(timestep,1);
00394 if (doTcl || doColvars) // include constraint forces
00395 computeGlobal->saveTotalForces(patch);
00396
00397 submitHalfstep(step);
00398 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00399
00400 if ( ! commOnly ) {
00401 addForceToMomentum(-0.5*timestep);
00402 if (staleForces || doNonbonded)
00403 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00404 if (staleForces || doFullElectrostatics)
00405 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00406 }
00407
00408 // rattle2(timestep,step);
00409
00410 submitReductions(step);
00411 submitCollections(step);
00412 //Update adaptive tempering temperature
00413 adaptTempUpdate(step);
00414
00415 #ifndef NAMD_CUDA
00416 LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
00417 #endif
00418
00419 #if CYCLE_BARRIER
00420 cycleBarrier(!((step+1) % stepsPerCycle), step);
00421 #elif PME_BARRIER
00422 cycleBarrier(doFullElectrostatics, step);
00423 #elif STEP_BARRIER
00424 cycleBarrier(1, step);
00425 #endif
00426
00427 if(Node::Object()->specialTracing || simParams->statsOn){
00428 int bstep = simParams->traceStartStep;
00429 int estep = bstep + simParams->numTraceSteps;
00430 if(step == bstep || step == estep){
00431 traceBarrier(step);
00432 }
00433 }
00434
00435 #ifdef MEASURE_NAMD_WITH_PAPI
00436 if(simParams->papiMeasure) {
00437 int bstep = simParams->papiMeasureStartStep;
00438 int estep = bstep + simParams->numPapiMeasureSteps;
00439 if(step == bstep || step==estep) {
00440 papiMeasureBarrier(step);
00441 }
00442 }
00443 #endif
00444
00445 if(traceIsOn()){
00446 traceUserEvent(eventEndOfTimeStep);
00447 sprintf(traceNote, "%s%d",tracePrefix,step);
00448 traceUserSuppliedNote(traceNote);
00449 }
00450 rebalanceLoad(step);
00451
00452 #if PME_BARRIER
00453 // a step before PME
00454 cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00455 #endif
00456
00457 #if USE_HPM
00458 if(step == START_HPM_STEP)
00459 (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
00460
00461 if(step == STOP_HPM_STEP)
00462 (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
00463 #endif
00464 }
00465 }
|
|
|
Definition at line 966 of file Sequencer.C. References Lattice::apply_transform(), HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), CompAtomExt::groupFixed, CompAtom::hydrogenGroupSize, Molecule::is_atom_exPressure(), j, SimParameters::langevinPistonOn, Patch::lattice, FullAtom::mass, Node::molecule, Patch::numAtoms, Node::Object(), patch, Position, CompAtom::position, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, slowFreq, SimParameters::useGroupPressure, FullAtom::velocity, Velocity, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz. Referenced by integrate(). 00967 {
00968 if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
00969 {
00970 FullAtom *a = patch->atom.begin();
00971 int numAtoms = patch->numAtoms;
00972 Tensor factor = broadcast->positionRescaleFactor.get(step);
00973 // JCP FIX THIS!!!
00974 Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
00975 patch->lattice.rescale(factor);
00976 Molecule *mol = Node::Object()->molecule;
00977 if ( simParams->useGroupPressure )
00978 {
00979 int hgs;
00980 for ( int i = 0; i < numAtoms; i += hgs ) {
00981 int j;
00982 hgs = a[i].hydrogenGroupSize;
00983 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00984 for ( j = i; j < (i+hgs); ++j ) {
00985 a[j].position = patch->lattice.apply_transform(
00986 a[j].fixedPosition,a[j].transform);
00987 }
00988 continue;
00989 }
00990 BigReal m_cm = 0;
00991 Position x_cm(0,0,0);
00992 Velocity v_cm(0,0,0);
00993 for ( j = i; j < (i+hgs); ++j ) {
00994 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00995 m_cm += a[j].mass;
00996 x_cm += a[j].mass * a[j].position;
00997 v_cm += a[j].mass * a[j].velocity;
00998 }
00999 x_cm /= m_cm;
01000 Position new_x_cm = x_cm;
01001 patch->lattice.rescale(new_x_cm,factor);
01002 Position delta_x_cm = new_x_cm - x_cm;
01003 v_cm /= m_cm;
01004 Velocity delta_v_cm;
01005 delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
01006 delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
01007 delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
01008 for ( j = i; j < (i+hgs); ++j ) {
01009 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01010 a[j].position = patch->lattice.apply_transform(
01011 a[j].fixedPosition,a[j].transform);
01012 continue;
01013 }
01014 if ( mol->is_atom_exPressure(a[j].id) ) continue;
01015 a[j].position += delta_x_cm;
01016 a[j].velocity += delta_v_cm;
01017 }
01018 }
01019 }
01020 else
01021 {
01022 for ( int i = 0; i < numAtoms; ++i )
01023 {
01024 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01025 a[i].position = patch->lattice.apply_transform(
01026 a[i].fixedPosition,a[i].transform);
01027 continue;
01028 }
01029 if ( mol->is_atom_exPressure(a[i].id) ) continue;
01030 patch->lattice.rescale(a[i].position,factor);
01031 a[i].velocity.x *= velFactor.x;
01032 a[i].velocity.y *= velFactor.y;
01033 a[i].velocity.z *= velFactor.z;
01034 }
01035 }
01036 }
01037 }
|
|
|
Definition at line 715 of file Sequencer.C. References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMANN, Random::gaussian_vector(), SimParameters::langevinOn, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, Node::molecule, Patch::numAtoms, Node::Object(), patch, random, simParams, and FullAtom::velocity. 00716 {
00717 if ( simParams->langevinOn )
00718 {
00719 FullAtom *a = patch->atom.begin();
00720 int numAtoms = patch->numAtoms;
00721 Molecule *molecule = Node::Object()->molecule;
00722 BigReal dt = dt_fs * 0.001; // convert to ps
00723 BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
00724 if (simParams->adaptTempOn && simParams->adaptTempLangevin)
00725 {
00726 kbT = BOLTZMANN*adaptTempT;
00727 }
00728
00729 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00730 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00731
00732 for ( int i = 0; i < numAtoms; ++i )
00733 {
00734 BigReal f1 = exp( -1. * dt * a[i].langevinParam );
00735 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
00736 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00737
00738 a[i].velocity *= f1;
00739 a[i].velocity += f2 * random->gaussian_vector();
00740 }
00741 }
00742 }
|
|
|
Definition at line 744 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, SimParameters::drudeOn, SimParameters::langevinOn, FullAtom::langevinParam, FullAtom::mass, Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, and FullAtom::velocity. Referenced by integrate(). 00745 {
00746 if ( simParams->langevinOn )
00747 {
00748 FullAtom *a = patch->atom.begin();
00749 int numAtoms = patch->numAtoms;
00750 Molecule *molecule = Node::Object()->molecule;
00751 BigReal dt = dt_fs * 0.001; // convert to ps
00752 int i;
00753
00754 if (simParams->drudeOn) {
00755 for (i = 0; i < numAtoms; i++) {
00756
00757 if (i < numAtoms-1 &&
00758 a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
00759 //printf("*** Found Drude particle %d\n", a[i+1].id);
00760 // i+1 is a Drude particle with parent i
00761
00762 // convert from Cartesian coordinates to (COM,bond) coordinates
00763 BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
00764 Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
00765 Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
00766 BigReal dt_gamma;
00767
00768 // use Langevin damping factor i for v_com
00769 dt_gamma = dt * a[i].langevinParam;
00770 if (dt_gamma != 0.0) {
00771 v_com *= ( 1. - 0.5 * dt_gamma );
00772 }
00773
00774 // use Langevin damping factor i+1 for v_bnd
00775 dt_gamma = dt * a[i+1].langevinParam;
00776 if (dt_gamma != 0.0) {
00777 v_bnd *= ( 1. - 0.5 * dt_gamma );
00778 }
00779
00780 // convert back
00781 a[i].velocity = v_com - m * v_bnd;
00782 a[i+1].velocity = v_bnd + a[i].velocity;
00783
00784 i++; // +1 from loop, we've updated both particles
00785 }
00786 else {
00787 BigReal dt_gamma = dt * a[i].langevinParam;
00788 if ( ! dt_gamma ) continue;
00789
00790 a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00791 }
00792
00793 } // end for
00794 } // end if drudeOn
00795 else {
00796
00797 for ( i = 0; i < numAtoms; ++i )
00798 {
00799 BigReal dt_gamma = dt * a[i].langevinParam;
00800 if ( ! dt_gamma ) continue;
00801
00802 a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00803 }
00804
00805 } // end else
00806
00807 } // end if langevinOn
00808 }
|
|
|
Definition at line 811 of file Sequencer.C. References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMANN, SimParameters::drudeOn, SimParameters::drudeTemp, Random::gaussian_vector(), SimParameters::langevinOn, FullAtom::langevinParam, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, FullAtom::mass, Node::molecule, Patch::numAtoms, Node::Object(), patch, random, rattle1(), simParams, and FullAtom::velocity. Referenced by integrate(). 00812 {
00813 if ( simParams->langevinOn )
00814 {
00815 rattle1(dt_fs,1); // conserve momentum if gammas differ
00816
00817 FullAtom *a = patch->atom.begin();
00818 int numAtoms = patch->numAtoms;
00819 Molecule *molecule = Node::Object()->molecule;
00820 BigReal dt = dt_fs * 0.001; // convert to ps
00821 BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
00822 if (simParams->adaptTempOn && simParams->adaptTempLangevin)
00823 {
00824 kbT = BOLTZMANN*adaptTempT;
00825 }
00826 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00827 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00828 int i;
00829
00830 if (simParams->drudeOn) {
00831 BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp); // drude bond Temp
00832
00833 for (i = 0; i < numAtoms; i++) {
00834
00835 if (i < numAtoms-1 &&
00836 a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
00837 //printf("*** Found Drude particle %d\n", a[i+1].id);
00838 // i+1 is a Drude particle with parent i
00839
00840 // convert from Cartesian coordinates to (COM,bond) coordinates
00841 BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
00842 Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
00843 Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
00844 BigReal dt_gamma;
00845
00846 // use Langevin damping factor i for v_com
00847 dt_gamma = dt * a[i].langevinParam;
00848 if (dt_gamma != 0.0) {
00849 BigReal mass = a[i].mass + a[i+1].mass;
00850 v_com += random->gaussian_vector() *
00851 sqrt( 2 * dt_gamma * kbT *
00852 ( a[i].partition ? tempFactor : 1.0 ) / mass );
00853 v_com /= ( 1. + 0.5 * dt_gamma );
00854 }
00855
00856 // use Langevin damping factor i+1 for v_bnd
00857 dt_gamma = dt * a[i+1].langevinParam;
00858 if (dt_gamma != 0.0) {
00859 BigReal mass = a[i+1].mass * (1. - m);
00860 v_bnd += random->gaussian_vector() *
00861 sqrt( 2 * dt_gamma * kbT_bnd *
00862 ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
00863 v_bnd /= ( 1. + 0.5 * dt_gamma );
00864 }
00865
00866 // convert back
00867 a[i].velocity = v_com - m * v_bnd;
00868 a[i+1].velocity = v_bnd + a[i].velocity;
00869
00870 i++; // +1 from loop, we've updated both particles
00871 }
00872 else {
00873 BigReal dt_gamma = dt * a[i].langevinParam;
00874 if ( ! dt_gamma ) continue;
00875
00876 a[i].velocity += random->gaussian_vector() *
00877 sqrt( 2 * dt_gamma * kbT *
00878 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00879 a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00880 }
00881
00882 } // end for
00883 } // end if drudeOn
00884 else {
00885
00886 for ( i = 0; i < numAtoms; ++i )
00887 {
00888 BigReal dt_gamma = dt * a[i].langevinParam;
00889 if ( ! dt_gamma ) continue;
00890
00891 a[i].velocity += random->gaussian_vector() *
00892 sqrt( 2 * dt_gamma * kbT *
00893 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00894 a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00895 }
00896
00897 } // end else
00898
00899 } // end if langevinOn
00900 }
|
|
|
Definition at line 1253 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, SimParameters::cutoff, Node::enableEarlyExit(), CompAtomExt::id, iERROR(), iout, Vector::length(), Vector::length2(), SimParameters::maximumMove, Patch::numAtoms, Node::Object(), patch, Patch::patchID, PDBVELFACTOR, simParams, terminate(), and FullAtom::velocity. Referenced by integrate(). 01254 {
01255 FullAtom *a = patch->atom.begin();
01256 int numAtoms = patch->numAtoms;
01257 if ( simParams->maximumMove ) {
01258 const BigReal dt = timestep / TIMEFACTOR;
01259 const BigReal maxvel = simParams->maximumMove / dt;
01260 const BigReal maxvel2 = maxvel * maxvel;
01261 for ( int i=0; i<numAtoms; ++i ) {
01262 if ( a[i].velocity.length2() > maxvel2 ) {
01263 a[i].velocity *= ( maxvel / a[i].velocity.length() );
01264 }
01265 }
01266 } else {
01267 const BigReal dt = timestep / TIMEFACTOR;
01268 const BigReal maxvel = simParams->cutoff / dt;
01269 const BigReal maxvel2 = maxvel * maxvel;
01270 int killme = 0;
01271 for ( int i=0; i<numAtoms; ++i ) {
01272 killme = killme || ( a[i].velocity.length2() > maxvel2 );
01273 }
01274 if ( killme ) {
01275 killme = 0;
01276 for ( int i=0; i<numAtoms; ++i ) {
01277 if ( a[i].velocity.length2() > maxvel2 ) {
01278 ++killme;
01279 iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
01280 << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
01281 << ( PDBVELFACTOR * maxvel ) << ", atom "
01282 << i << " of " << numAtoms << " on patch "
01283 << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
01284 }
01285 }
01286 iout << iERROR <<
01287 "Atoms moving too fast; simulation has become unstable ("
01288 << killme << " atoms on patch " << patch->patchID
01289 << " pe " << CkMyPe() << ").\n" << endi;
01290 Node::Object()->enableEarlyExit();
01291 terminate();
01292 }
01293 }
01294 }
|
|
|
Definition at line 1296 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::minimizeOn, Patch::numAtoms, patch, simParams, and FullAtom::velocity. Referenced by integrate(). 01297 {
01298 if ( simParams->minimizeOn ) {
01299 FullAtom *a = patch->atom.begin();
01300 int numAtoms = patch->numAtoms;
01301 for ( int i=0; i<numAtoms; ++i ) {
01302 a[i].velocity = 0.;
01303 }
01304 }
01305 }
|
|
|
Definition at line 510 of file Sequencer.C. References HomePatch::atom, Patch::atomMapper, ResizeArray< Elem >::begin(), BigReal, broadcast, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, Flags::doNonbonded, ResizeArray< Elem >::end(), SimParameters::firstTimestep, Patch::flags, SimParameters::fullElectFrequency, SimParameters::GBISOn, SimpleBroadcastObject< T >::get(), SimParameters::LCPOOn, SimParameters::lonepairs, Flags::maxForceMerged, Flags::maxForceUsed, ControllerBroadcasts::minimizeCoefficient, minimizeMoveDownhill(), SimParameters::mollyOn, SimParameters::N, newMinimizeDirection(), newMinimizePosition(), patch, quenchVelocities(), rebalanceLoad(), AtomMapper::registerIDsFullAtom(), runComputeObjects(), simParams, Flags::step, SimParameters::stepsPerCycle, submitCollections(), submitMinimizeReductions(), and TIMEFACTOR. Referenced by algorithm(). 00510 {
00511 const int numberOfSteps = simParams->N;
00512 const int stepsPerCycle = simParams->stepsPerCycle;
00513 int &step = patch->flags.step;
00514 step = simParams->firstTimestep;
00515
00516 int &maxForceUsed = patch->flags.maxForceUsed;
00517 int &maxForceMerged = patch->flags.maxForceMerged;
00518 maxForceUsed = Results::normal;
00519 maxForceMerged = Results::normal;
00520 int &doNonbonded = patch->flags.doNonbonded;
00521 doNonbonded = 1;
00522 maxForceUsed = Results::nbond;
00523 maxForceMerged = Results::nbond;
00524 const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00525 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00526 doFullElectrostatics = dofull;
00527 if ( dofull ) {
00528 maxForceMerged = Results::slow;
00529 maxForceUsed = Results::slow;
00530 }
00531 int &doMolly = patch->flags.doMolly;
00532 doMolly = simParams->mollyOn && doFullElectrostatics;
00533 // BEGIN LA
00534 int &doLoweAndersen = patch->flags.doLoweAndersen;
00535 doLoweAndersen = 0;
00536 // END LA
00537
00538 int &doGBIS = patch->flags.doGBIS;
00539 doGBIS = simParams->GBISOn;
00540
00541 int &doLCPO = patch->flags.doLCPO;
00542 doLCPO = simParams->LCPOOn;
00543
00544 int &doEnergy = patch->flags.doEnergy;
00545 doEnergy = 1;
00546
00547 if (simParams->lonepairs) {
00548 patch->atomMapper->registerIDsFullAtom(
00549 patch->atom.begin(),patch->atom.end());
00550 }
00551
00552 runComputeObjects(1,step<numberOfSteps); // must migrate here!
00553
00554 BigReal fmax2 = TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00555
00556 submitMinimizeReductions(step,fmax2);
00557 rebalanceLoad(step);
00558
00559 int downhill = 1; // start out just fixing bad contacts
00560 int minSeq = 0;
00561 for ( ++step; step <= numberOfSteps; ++step ) {
00562 BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00563 if ( downhill ) {
00564 if ( c ) minimizeMoveDownhill(fmax2);
00565 else {
00566 downhill = 0;
00567 fmax2 *= 10000.;
00568 }
00569 }
00570 if ( ! downhill ) {
00571 if ( ! c ) { // new direction
00572 c = broadcast->minimizeCoefficient.get(minSeq++);
00573 newMinimizeDirection(c); // v = c * v + f
00574 c = broadcast->minimizeCoefficient.get(minSeq++);
00575 } // same direction
00576 newMinimizePosition(c); // x = x + c * v
00577 }
00578
00579 runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00580 submitMinimizeReductions(step,fmax2);
00581 submitCollections(step, 1); // write out zeros for velocities
00582 rebalanceLoad(step);
00583 }
00584 quenchVelocities(); // zero out bogus velocity
00585 }
|
|
|
Definition at line 588 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), Patch::f, SimParameters::fixedAtomsOn, Force, CompAtom::hydrogenGroupSize, j, Vector::length2(), Patch::numAtoms, patch, CompAtom::position, simParams, and Vector::unit(). Referenced by minimize(). 00588 {
00589
00590 FullAtom *a = patch->atom.begin();
00591 Force *f1 = patch->f[Results::normal].begin();
00592 Force *f2 = patch->f[Results::nbond].begin();
00593 Force *f3 = patch->f[Results::slow].begin();
00594 int numAtoms = patch->numAtoms;
00595
00596 for ( int i = 0; i < numAtoms; ++i ) {
00597 if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
00598 Force f = f1[i] + f2[i] + f3[i];
00599 if ( f.length2() > fmax2 ) {
00600 a[i].position += ( 0.1 * f.unit() );
00601 int hgs = a[i].hydrogenGroupSize; // 0 if not parent
00602 for ( int j=1; j<hgs; ++j ) {
00603 a[++i].position += ( 0.1 * f.unit() );
00604 }
00605 }
00606 }
00607 }
|
|
|
Definition at line 610 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, Patch::f, SimParameters::fixedAtomsOn, Force, CompAtom::hydrogenGroupSize, j, Vector::length2(), Patch::numAtoms, patch, simParams, TIMEFACTOR, and FullAtom::velocity. Referenced by minimize(). 00610 {
00611 FullAtom *a = patch->atom.begin();
00612 Force *f1 = patch->f[Results::normal].begin();
00613 Force *f2 = patch->f[Results::nbond].begin();
00614 Force *f3 = patch->f[Results::slow].begin();
00615 int numAtoms = patch->numAtoms;
00616
00617 for ( int i = 0; i < numAtoms; ++i ) {
00618 a[i].velocity *= c;
00619 a[i].velocity += f1[i] + f2[i] + f3[i];
00620 if ( simParams->fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00621 }
00622
00623 // prevent hydrogens from being left behind
00624 BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00625 // int adjustCount = 0;
00626 int hgs;
00627 for ( int i = 0; i < numAtoms; i += hgs ) {
00628 hgs = a[i].hydrogenGroupSize;
00629 BigReal minChildVel = a[i].velocity.length2();
00630 if ( minChildVel < fmax2 ) continue;
00631 int adjustChildren = 1;
00632 for ( int j = i+1; j < (i+hgs); ++j ) {
00633 if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
00634 }
00635 if ( adjustChildren ) {
00636 // if ( hgs > 1 ) ++adjustCount;
00637 for ( int j = i+1; j < (i+hgs); ++j ) {
00638 a[j].velocity = a[i].velocity;
00639 }
00640 }
00641 }
00642 // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
00643
00644 }
|
|
|
Definition at line 647 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, patch, CompAtom::position, and FullAtom::velocity. Referenced by minimize(). 00647 {
00648 FullAtom *a = patch->atom.begin();
00649 int numAtoms = patch->numAtoms;
00650
00651 for ( int i = 0; i < numAtoms; ++i ) {
00652 a[i].position += c * a[i].velocity;
00653 }
00654 }
|
|
|
Definition at line 657 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity. Referenced by minimize(). 00657 {
00658 FullAtom *a = patch->atom.begin();
00659 int numAtoms = patch->numAtoms;
00660
00661 for ( int i = 0; i < numAtoms; ++i ) {
00662 a[i].velocity = 0;
00663 }
00664 }
|
|
||||||||||||
|
Definition at line 1224 of file Sequencer.C. References ADD_TENSOR_OBJECT, Node::enableEarlyExit(), iERROR(), iout, Node::Object(), patch, pressureProfileReduction, HomePatch::rattle1(), reduction, REDUCTION_VIRIAL_NORMAL, SimParameters::rigidBonds, simParams, and terminate(). Referenced by integrate(), and langevinVelocitiesBBK2(). 01225 {
01226 if ( simParams->rigidBonds != RIGID_NONE ) {
01227 Tensor virial;
01228 Tensor *vp = ( pressure ? &virial : 0 );
01229 if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
01230 iout << iERROR <<
01231 "Constraint failure; simulation has become unstable.\n" << endi;
01232 Node::Object()->enableEarlyExit();
01233 terminate();
01234 }
01235 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01236 }
01237 }
|
|
||||||||||||
|
Definition at line 1239 of file Sequencer.C. References ADD_TENSOR_OBJECT, patch, HomePatch::rattle2(), reduction, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_VIRIAL_NORMAL, SimParameters::rigidBonds, and simParams. 01240 {
01241 if ( simParams->rigidBonds != RIGID_NONE ) {
01242 Tensor virial;
01243 patch->rattle2(dt, &virial);
01244 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01245 // we need to add to alt and int virial because not included in forces
01246 #ifdef ALTVIRIAL
01247 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01248 #endif
01249 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,virial);
01250 }
01251 }
|
|
||||||||||||
|
Definition at line 1114 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, BOLTZMANN, SimParameters::fixedAtomsOn, Random::gaussian_vector(), SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, FullAtom::mass, NAMD_bug(), Patch::numAtoms, patch, random, SimParameters::reassignFreq, SimParameters::reassignHold, SimParameters::reassignIncr, SimParameters::reassignTemp, simParams, and FullAtom::velocity. Referenced by integrate(). 01115 {
01116 const int reassignFreq = simParams->reassignFreq;
01117 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01118 FullAtom *a = patch->atom.begin();
01119 int numAtoms = patch->numAtoms;
01120 BigReal newTemp = simParams->reassignTemp;
01121 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01122 if ( simParams->reassignIncr > 0.0 ) {
01123 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01124 newTemp = simParams->reassignHold;
01125 } else {
01126 if ( newTemp < simParams->reassignHold )
01127 newTemp = simParams->reassignHold;
01128 }
01129 BigReal kbT = BOLTZMANN * newTemp;
01130
01131 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01132 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01133
01134 for ( int i = 0; i < numAtoms; ++i )
01135 {
01136 a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
01137 sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
01138 * random->gaussian_vector() );
01139 }
01140 } else {
01141 NAMD_bug("Sequencer::reassignVelocities called improperly!");
01142 }
01143 }
|
|
|
Definition at line 1975 of file Sequencer.C. References LdbCoordinator::getNumStepsToRun(), Patch::getPatchID(), ldbSteps, LdbCoordinator::Object(), pairlistsAreValid, patch, LdbCoordinator::rebalance(), and HomePatch::submitLoadStats(). Referenced by integrate(), and minimize(). 01975 {
01976 if ( ! ldbSteps ) {
01977 ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
01978 }
01979 if ( ! --ldbSteps ) {
01980 patch->submitLoadStats(timestep);
01981 ldbCoordinator->rebalance(this,patch->getPatchID());
01982 pairlistsAreValid = 0;
01983 }
01984 }
|
|
|
|
Definition at line 1173 of file Sequencer.C. References HomePatch::atom, Molecule::atomcharge(), ResizeArray< Elem >::begin(), CompAtom::charge, Node::molecule, Patch::numAtoms, Node::Object(), and patch. Referenced by algorithm(). 01174 {
01175 FullAtom *a = patch->atom.begin();
01176 int numAtoms = patch->numAtoms;
01177 Molecule *molecule = Node::Object()->molecule;
01178 for ( int i = 0; i < numAtoms; ++i )
01179 {
01180 a[i].charge = molecule->atomcharge(a[i].id);
01181 }
01182 }
|
|
||||||||||||||||
|
Definition at line 1057 of file Sequencer.C. References SimParameters::accelMDdihe, SimParameters::accelMDdual, SimParameters::accelMDFirstStep, SimParameters::accelMDLastStep, SimParameters::accelMDOn, ControllerBroadcasts::accelMDRescaleFactor, BigReal, broadcast, Patch::f, SimpleBroadcastObject< T >::get(), NAMD_die(), Patch::numAtoms, patch, and simParams. Referenced by integrate(). 01058 {
01059 if (!simParams->accelMDOn) return;
01060 if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
01061
01062 Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
01063 const BigReal factor_dihe = accelMDfactor[0];
01064 const BigReal factor_tot = accelMDfactor[1];
01065 const int numAtoms = patch->numAtoms;
01066
01067 if (simParams->accelMDdihe && factor_tot <1 )
01068 NAMD_die("accelMD broadcasting error!\n");
01069 if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
01070 NAMD_die("accelMD broadcasting error!\n");
01071
01072 if (simParams->accelMDdihe && factor_dihe < 1) {
01073 for (int i = 0; i < numAtoms; ++i)
01074 if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01075 patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);
01076 }
01077
01078 if ( !simParams->accelMDdihe && factor_tot < 1) {
01079 for (int i = 0; i < numAtoms; ++i)
01080 patch->f[Results::normal][i] *= factor_tot;
01081 if (doNonbonded) {
01082 for (int i = 0; i < numAtoms; ++i)
01083 patch->f[Results::nbond][i] *= factor_tot;
01084 }
01085 if (doFullElectrostatics) {
01086 for (int i = 0; i < numAtoms; ++i)
01087 patch->f[Results::slow][i] *= factor_tot;
01088 }
01089 }
01090
01091 if (simParams->accelMDdual && factor_dihe < 1) {
01092 for (int i = 0; i < numAtoms; ++i)
01093 if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01094 patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
01095 }
01096
01097 }
|
|
|
Definition at line 1039 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, broadcast, SimpleBroadcastObject< T >::get(), Patch::numAtoms, patch, SimParameters::rescaleFreq, rescaleVelocities_numTemps, simParams, FullAtom::velocity, and ControllerBroadcasts::velocityRescaleFactor. Referenced by integrate(). 01040 {
01041 const int rescaleFreq = simParams->rescaleFreq;
01042 if ( rescaleFreq > 0 ) {
01043 FullAtom *a = patch->atom.begin();
01044 int numAtoms = patch->numAtoms;
01045 ++rescaleVelocities_numTemps;
01046 if ( rescaleVelocities_numTemps == rescaleFreq ) {
01047 BigReal factor = broadcast->velocityRescaleFactor.get(step);
01048 for ( int i = 0; i < numAtoms; ++i )
01049 {
01050 a[i].velocity *= factor;
01051 }
01052 rescaleVelocities_numTemps = 0;
01053 }
01054 }
01055 }
|
|
|
Definition at line 1163 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity. Referenced by algorithm(). 01164 {
01165 FullAtom *a = patch->atom.begin();
01166 int numAtoms = patch->numAtoms;
01167 for ( int i = 0; i < numAtoms; ++i )
01168 {
01169 a[i].velocity *= factor;
01170 }
01171 }
|
|
|
Definition at line 87 of file Sequencer.C. References awaken(), DebugM, Patch::getPatchID(), patch, PATCH_PRIORITY, and SEQ_STK_SZ. Referenced by HomePatch::runSequencer(). 00088 {
00089 // create a Thread and invoke it
00090 DebugM(4, "::run() - this = " << this << "\n" );
00091 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
00092 CthSetStrategyDefault(thread);
00093 priority = PATCH_PRIORITY(patch->getPatchID());
00094 awaken();
00095 }
|
|
||||||||||||
|
Definition at line 1805 of file Sequencer.C. References ADD_TENSOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), Flags::doFullElectrostatics, Flags::doGBIS, Flags::doLoweAndersen, Flags::doMolly, Flags::doNonbonded, Patch::f, Patch::flags, SimParameters::fullElectFrequency, HomePatch::gbisComputeAfterP1(), HomePatch::gbisComputeAfterP2(), Patch::getPatchID(), SimParameters::lonepairs, HomePatch::loweAndersenFinish(), HomePatch::mollyMollify(), Patch::numAtoms, Patch::p, pairlistsAge, pairlistsAreValid, SimParameters::pairlistsPerCycle, patch, PATCH_PRIORITY, Patch::pExt, HomePatch::positionsReady(), HomePatch::redistrib_lonepair_forces(), HomePatch::redistrib_swm4_forces(), HomePatch::redistrib_tip4p_forces(), reduction, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, HomePatch::reposition_all_lonepairs(), Flags::savePairlists, Flags::sequence, simParams, SimParameters::stepsPerCycle, suspend(), Flags::usePairlists, SimParameters::usePairlists, and SimParameters::watmodel. Referenced by integrate(), and minimize(). 01806 {
01807 if ( migration ) pairlistsAreValid = 0;
01808 #ifdef NAMD_CUDA
01809 if ( pairlistsAreValid &&
01810 ( patch->flags.doFullElectrostatics || ! simParams->fullElectFrequency )
01811 && ( pairlistsAge > (
01812 #else
01813 if ( pairlistsAreValid && ( pairlistsAge > (
01814 #endif
01815 (simParams->stepsPerCycle - 1) / simParams->pairlistsPerCycle ) ) ) {
01816 pairlistsAreValid = 0;
01817 }
01818 if ( ! simParams->usePairlists ) pairlists = 0;
01819 patch->flags.usePairlists = pairlists || pairlistsAreValid;
01820 patch->flags.savePairlists =
01821 pairlists && ! pairlistsAreValid;
01822
01823 if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
01824
01825 patch->positionsReady(migration); // updates flags.sequence
01826 int seq = patch->flags.sequence;
01827 int basePriority = ( (seq & 0xffff) << 15 )
01828 + PATCH_PRIORITY(patch->getPatchID());
01829 if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
01830 priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
01831 suspend(); // until all deposit boxes close
01832 patch->gbisComputeAfterP1();
01833 priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
01834 suspend();
01835 patch->gbisComputeAfterP2();
01836 priority = basePriority + COMPUTE_HOME_PRIORITY;
01837 suspend();
01838 } else {
01839 priority = basePriority + COMPUTE_HOME_PRIORITY;
01840 suspend(); // until all deposit boxes close
01841 }
01842
01843 if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
01844 pairlistsAreValid = 1;
01845 pairlistsAge = 0;
01846 }
01847 if ( pairlistsAreValid ) ++pairlistsAge;
01848
01849 if (simParams->lonepairs) {
01850 {
01851 Tensor virial;
01852 patch->redistrib_lonepair_forces(Results::normal, &virial);
01853 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
01854 }
01855 if (patch->flags.doNonbonded) {
01856 Tensor virial;
01857 patch->redistrib_lonepair_forces(Results::nbond, &virial);
01858 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
01859 }
01860 if (patch->flags.doFullElectrostatics) {
01861 Tensor virial;
01862 patch->redistrib_lonepair_forces(Results::slow, &virial);
01863 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
01864 }
01865 } else if (simParams->watmodel == WAT_TIP4) {
01866 {
01867 Tensor virial;
01868 patch->redistrib_tip4p_forces(Results::normal, &virial);
01869 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
01870 }
01871 if (patch->flags.doNonbonded) {
01872 Tensor virial;
01873 patch->redistrib_tip4p_forces(Results::nbond, &virial);
01874 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
01875 }
01876 if (patch->flags.doFullElectrostatics) {
01877 Tensor virial;
01878 patch->redistrib_tip4p_forces(Results::slow, &virial);
01879 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
01880 }
01881 } else if (simParams->watmodel == WAT_SWM4) {
01882 {
01883 Tensor virial;
01884 patch->redistrib_swm4_forces(Results::normal, &virial);
01885 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
01886 }
01887 if (patch->flags.doNonbonded) {
01888 Tensor virial;
01889 patch->redistrib_swm4_forces(Results::nbond, &virial);
01890 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
01891 }
01892 if (patch->flags.doFullElectrostatics) {
01893 Tensor virial;
01894 patch->redistrib_swm4_forces(Results::slow, &virial);
01895 ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
01896 }
01897 }
01898
01899 if ( patch->flags.doMolly ) {
01900 Tensor virial;
01901 patch->mollyMollify(&virial);
01902 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
01903 }
01904
01905
01906 // BEGIN LA
01907 if (patch->flags.doLoweAndersen) {
01908 patch->loweAndersenFinish();
01909 }
01910 // END LA
01911 #ifdef NAMD_CUDA_XXX
01912 int numAtoms = patch->numAtoms;
01913 FullAtom *a = patch->atom.begin();
01914 for ( int i=0; i<numAtoms; ++i ) {
01915 CkPrintf("%d %g %g %g\n", a[i].id,
01916 patch->f[Results::normal][i].x +
01917 patch->f[Results::nbond][i].x +
01918 patch->f[Results::slow][i].x,
01919 patch->f[Results::normal][i].y +
01920 patch->f[Results::nbond][i].y +
01921 patch->f[Results::slow][i].y,
01922 patch->f[Results::normal][i].z +
01923 patch->f[Results::nbond][i].z +
01924 patch->f[Results::slow][i].z);
01925 CkPrintf("%d %g %g %g\n", a[i].id,
01926 patch->f[Results::normal][i].x,
01927 patch->f[Results::nbond][i].x,
01928 patch->f[Results::slow][i].x);
01929 CkPrintf("%d %g %g %g\n", a[i].id,
01930 patch->f[Results::normal][i].y,
01931 patch->f[Results::nbond][i].y,
01932 patch->f[Results::slow][i].y);
01933 CkPrintf("%d %g %g %g\n", a[i].id,
01934 patch->f[Results::normal][i].z,
01935 patch->f[Results::nbond][i].z,
01936 patch->f[Results::slow][i].z);
01937 }
01938 #endif
01939
01940 #if PRINT_FORCES
01941 int numAtoms = patch->numAtoms;
01942 FullAtom *a = patch->atom.begin();
01943 for ( int i=0; i<numAtoms; ++i ) {
01944 float fxNo = patch->f[Results::normal][i].x;
01945 float fxNb = patch->f[Results::nbond][i].x;
01946 float fxSl = patch->f[Results::slow][i].x;
01947 float fyNo = patch->f[Results::normal][i].y;
01948 float fyNb = patch->f[Results::nbond][i].y;
01949 float fySl = patch->f[Results::slow][i].y;
01950 float fzNo = patch->f[Results::normal][i].z;
01951 float fzNb = patch->f[Results::nbond][i].z;
01952 float fzSl = patch->f[Results::slow][i].z;
01953 float fx = fxNo+fxNb+fxSl;
01954 float fy = fyNo+fyNb+fySl;
01955 float fz = fzNo+fzNb+fzSl;
01956
01957 float f = sqrt(fx*fx+fy*fy+fz*fz);
01958 int id = patch->pExt[i].id;
01959 int seq = patch->flags.sequence;
01960 float x = patch->p[i].position.x;
01961 float y = patch->p[i].position.y;
01962 float z = patch->p[i].position.z;
01963 //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
01964 CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
01965 //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
01966 //fxNo,fyNo,fzNo,
01967 //fxNb,fyNb,fzNb,
01968 //fxSl,fySl,fzSl,
01969 fx,fy,fz
01970 );
01971 }
01972 #endif
01973 }
|
|
|
Definition at line 1202 of file Sequencer.C. References patch, and HomePatch::saveForce(). Referenced by integrate().
|
|
||||||||||||
|
Definition at line 1791 of file Sequencer.C. References HomePatch::atom, collection, Output::coordinateNeeded(), Patch::f, Patch::flags, Output::forceNeeded(), Patch::lattice, Flags::maxForceUsed, patch, CollectionMgr::submitForces(), CollectionMgr::submitPositions(), CollectionMgr::submitVelocities(), and Output::velocityNeeded(). Referenced by algorithm(), integrate(), and minimize(). 01792 {
01793 int prec = Output::coordinateNeeded(step);
01794 if ( prec )
01795 collection->submitPositions(step,patch->atom,patch->lattice,prec);
01796 if ( Output::velocityNeeded(step) )
01797 collection->submitVelocities(step,zeroVel,patch->atom);
01798 if ( Output::forceNeeded(step) ) {
01799 int maxForceUsed = patch->flags.maxForceUsed;
01800 if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
01801 collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
01802 }
01803 }
|
|
|
Definition at line 1307 of file Sequencer.C. References ADD_TENSOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, Lattice::c(), CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, Patch::lattice, Vector::length2(), FullAtom::mass, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, CompAtom::partition, patch, CompAtom::position, pressureProfileReduction, SimParameters::pressureProfileSlabs, reduction, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_VIRIAL_NORMAL, simParams, SimParameters::useGroupPressure, Velocity, FullAtom::velocity, Vector::x, Vector::y, and Vector::z. Referenced by integrate(). 01308 {
01309 // velocity-dependent quantities *** ONLY ***
01310 // positions are not at half-step when called
01311 FullAtom *a = patch->atom.begin();
01312 int numAtoms = patch->numAtoms;
01313
01314 #if CMK_BLUEGENEL
01315 CmiNetworkProgressAfter (0);
01316 #endif
01317
01318 {
01319 BigReal kineticEnergy = 0;
01320 Tensor virial;
01321 if ( simParams->pairInteractionOn ) {
01322 if ( simParams->pairInteractionSelf ) {
01323 for ( int i = 0; i < numAtoms; ++i ) {
01324 if ( a[i].partition != 1 ) continue;
01325 kineticEnergy += a[i].mass * a[i].velocity.length2();
01326 virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01327 }
01328 }
01329 } else {
01330 for ( int i = 0; i < numAtoms; ++i ) {
01331 if (a[i].mass < 0.01) continue;
01332 kineticEnergy += a[i].mass * a[i].velocity.length2();
01333 virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01334 }
01335 }
01336
01337 kineticEnergy *= 0.5 * 0.5;
01338 reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
01339 virial *= 0.5;
01340 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01341 #ifdef ALTVIRIAL
01342 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01343 #endif
01344 }
01345
01346 if (pressureProfileReduction) {
01347 int nslabs = simParams->pressureProfileSlabs;
01348 const Lattice &lattice = patch->lattice;
01349 BigReal idz = nslabs/lattice.c().z;
01350 BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01351 int useGroupPressure = simParams->useGroupPressure;
01352
01353 // Compute kinetic energy partition, possibly subtracting off
01354 // internal kinetic energy if group pressure is enabled.
01355 // Since the regular pressure is 1/2 mvv and the internal kinetic
01356 // term that is subtracted off for the group pressure is
01357 // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
01358 // 1/2 m * v * v_cm. The factor of 1/2 is because submitHalfstep
01359 // gets called twice per timestep.
01360 int hgs;
01361 for (int i=0; i<numAtoms; i += hgs) {
01362 int j, ppoffset;
01363 hgs = a[i].hydrogenGroupSize;
01364 int partition = a[i].partition;
01365
01366 BigReal m_cm = 0;
01367 Velocity v_cm(0,0,0);
01368 for (j=i; j< i+hgs; ++j) {
01369 m_cm += a[j].mass;
01370 v_cm += a[j].mass * a[j].velocity;
01371 }
01372 v_cm /= m_cm;
01373 for (j=i; j < i+hgs; ++j) {
01374 BigReal mass = a[j].mass;
01375 if (! (useGroupPressure && j != i)) {
01376 BigReal z = a[j].position.z;
01377 int slab = (int)floor((z-zmin)*idz);
01378 if (slab < 0) slab += nslabs;
01379 else if (slab >= nslabs) slab -= nslabs;
01380 ppoffset = 3*(slab + partition*nslabs);
01381 }
01382 BigReal wxx, wyy, wzz;
01383 if (useGroupPressure) {
01384 wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
01385 wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
01386 wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
01387 } else {
01388 wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
01389 wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
01390 wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
01391 }
01392 pressureProfileReduction->item(ppoffset ) += wxx;
01393 pressureProfileReduction->item(ppoffset+1) += wyy;
01394 pressureProfileReduction->item(ppoffset+2) += wzz;
01395 }
01396 }
01397 }
01398
01399 {
01400 BigReal intKineticEnergy = 0;
01401 Tensor intVirialNormal;
01402
01403 int hgs;
01404 for ( int i = 0; i < numAtoms; i += hgs ) {
01405
01406 #if CMK_BLUEGENEL
01407 CmiNetworkProgress ();
01408 #endif
01409
01410 hgs = a[i].hydrogenGroupSize;
01411 int j;
01412 BigReal m_cm = 0;
01413 Velocity v_cm(0,0,0);
01414 for ( j = i; j < (i+hgs); ++j ) {
01415 m_cm += a[j].mass;
01416 v_cm += a[j].mass * a[j].velocity;
01417 }
01418 v_cm /= m_cm;
01419 if ( simParams->pairInteractionOn ) {
01420 if ( simParams->pairInteractionSelf ) {
01421 for ( j = i; j < (i+hgs); ++j ) {
01422 if ( a[j].partition != 1 ) continue;
01423 BigReal mass = a[j].mass;
01424 Vector v = a[j].velocity;
01425 Vector dv = v - v_cm;
01426 intKineticEnergy += mass * (v * dv);
01427 intVirialNormal.outerAdd (mass, v, dv);
01428 }
01429 }
01430 } else {
01431 for ( j = i; j < (i+hgs); ++j ) {
01432 BigReal mass = a[j].mass;
01433 Vector v = a[j].velocity;
01434 Vector dv = v - v_cm;
01435 intKineticEnergy += mass * (v * dv);
01436 intVirialNormal.outerAdd(mass, v, dv);
01437 }
01438 }
01439 }
01440
01441 intKineticEnergy *= 0.5 * 0.5;
01442 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
01443 intVirialNormal *= 0.5;
01444 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01445 }
01446
01447 }
|
|
||||||||||||
|
Definition at line 1687 of file Sequencer.C. References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, Patch::f, SimParameters::fixedAtomsOn, FullAtom::fixedPosition, Patch::flags, Force, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, FullAtom::mass, Patch::numAtoms, Tensor::outerAdd(), patch, Patch::pExt, CompAtom::position, Position, SimParameters::printBadContacts, reduction, REDUCTION_ATOM_CHECKSUM, REDUCTION_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_MIN_F_DOT_F, REDUCTION_MIN_F_DOT_V, REDUCTION_MIN_HUGE_COUNT, REDUCTION_MIN_V_DOT_V, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, Flags::sequence, simParams, SubmitReduction::submit(), and FullAtom::velocity. Referenced by minimize(). 01688 {
01689 FullAtom *a = patch->atom.begin();
01690 Force *f1 = patch->f[Results::normal].begin();
01691 Force *f2 = patch->f[Results::nbond].begin();
01692 Force *f3 = patch->f[Results::slow].begin();
01693 int numAtoms = patch->numAtoms;
01694
01695 reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01696
01697 BigReal fdotf = 0;
01698 BigReal fdotv = 0;
01699 BigReal vdotv = 0;
01700 int numHuge = 0;
01701 for ( int i = 0; i < numAtoms; ++i ) {
01702 if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
01703 Force f = f1[i] + f2[i] + f3[i];
01704 BigReal ff = f * f;
01705 if ( ff > fmax2 ) {
01706 if (simParams->printBadContacts) {
01707 CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
01708 }
01709 ++numHuge;
01710 // pad scaling so minimizeMoveDownhill() doesn't miss them
01711 BigReal fmult = 1.01 * sqrt(fmax2/ff);
01712 f *= fmult; ff = f * f;
01713 f1[i] *= fmult;
01714 f2[i] *= fmult;
01715 f3[i] *= fmult;
01716 }
01717 fdotf += ff;
01718 fdotv += f * a[i].velocity;
01719 vdotv += a[i].velocity * a[i].velocity;
01720 }
01721
01722 reduction->item(REDUCTION_MIN_F_DOT_F) += fdotf;
01723 reduction->item(REDUCTION_MIN_F_DOT_V) += fdotv;
01724 reduction->item(REDUCTION_MIN_V_DOT_V) += vdotv;
01725 reduction->item(REDUCTION_MIN_HUGE_COUNT) += numHuge;
01726
01727 {
01728 Tensor intVirialNormal;
01729 Tensor intVirialNbond;
01730 Tensor intVirialSlow;
01731
01732 int hgs;
01733 for ( int i = 0; i < numAtoms; i += hgs ) {
01734 hgs = a[i].hydrogenGroupSize;
01735 int j;
01736 BigReal m_cm = 0;
01737 Position x_cm(0,0,0);
01738 for ( j = i; j < (i+hgs); ++j ) {
01739 m_cm += a[j].mass;
01740 x_cm += a[j].mass * a[j].position;
01741 }
01742 x_cm /= m_cm;
01743 for ( j = i; j < (i+hgs); ++j ) {
01744 BigReal mass = a[j].mass;
01745 // net force treated as zero for fixed atoms
01746 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01747 Vector dx = a[j].position - x_cm;
01748 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01749 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01750 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01751 }
01752 }
01753
01754 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01755 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
01756 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
01757 }
01758
01759 if ( simParams->fixedAtomsOn ) {
01760 Tensor fixVirialNormal;
01761 Tensor fixVirialNbond;
01762 Tensor fixVirialSlow;
01763 Vector fixForceNormal = 0;
01764 Vector fixForceNbond = 0;
01765 Vector fixForceSlow = 0;
01766
01767 for ( int j = 0; j < numAtoms; j++ ) {
01768 if ( a[j].atomFixed ) {
01769 Vector dx = a[j].fixedPosition;
01770 // all negative because fixed atoms cancels these forces
01771 fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j],dx);
01772 fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j],dx);
01773 fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j],dx);
01774 fixForceNormal -= patch->f[Results::normal][j];
01775 fixForceNbond -= patch->f[Results::nbond][j];
01776 fixForceSlow -= patch->f[Results::slow][j];
01777 }
01778 }
01779
01780 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
01781 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
01782 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
01783 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
01784 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
01785 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
01786 }
01787
01788 reduction->submit();
01789 }
|
|
|
Definition at line 666 of file Sequencer.C. References ADD_VECTOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, SubmitReduction::item(), FullAtom::mass, Patch::numAtoms, patch, reduction, REDUCTION_HALFSTEP_MOMENTUM, REDUCTION_MOMENTUM_MASS, simParams, FullAtom::velocity, and SimParameters::zeroMomentumAlt. Referenced by integrate(). 00666 {
00667
00668 FullAtom *a = patch->atom.begin();
00669 const int numAtoms = patch->numAtoms;
00670
00671 Vector momentum = 0;
00672 BigReal mass = 0;
00673 if ( simParams->zeroMomentumAlt ) {
00674 for ( int i = 0; i < numAtoms; ++i ) {
00675 momentum += a[i].mass * a[i].velocity;
00676 mass += 1.;
00677 }
00678 } else {
00679 for ( int i = 0; i < numAtoms; ++i ) {
00680 momentum += a[i].mass * a[i].velocity;
00681 mass += a[i].mass;
00682 }
00683 }
00684
00685 ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
00686 reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
00687 }
|
|
|
Definition at line 1449 of file Sequencer.C. References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, Lattice::c(), SimParameters::drudeOn, Patch::f, SimParameters::fixedAtomsOn, FullAtom::fixedPosition, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, Patch::lattice, Vector::length2(), HomePatch::marginViolations, FullAtom::mass, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, CompAtom::partition, patch, Position, CompAtom::position, pressureProfileReduction, SimParameters::pressureProfileSlabs, reduction, REDUCTION_ANGULAR_MOMENTUM, REDUCTION_ATOM_CHECKSUM, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY, REDUCTION_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_MARGIN_VIOLATIONS, REDUCTION_MOMENTUM, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, simParams, SubmitReduction::submit(), SimParameters::useGroupPressure, Velocity, FullAtom::velocity, Vector::x, Vector::y, and Vector::z. Referenced by integrate(). 01450 {
01451 FullAtom *a = patch->atom.begin();
01452 int numAtoms = patch->numAtoms;
01453
01454 #if CMK_BLUEGENEL
01455 CmiNetworkProgressAfter(0);
01456 #endif
01457
01458 reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01459 reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
01460
01461 {
01462 BigReal kineticEnergy = 0;
01463 Vector momentum = 0;
01464 Vector angularMomentum = 0;
01465 Vector o = patch->lattice.origin();
01466 int i;
01467 if ( simParams->pairInteractionOn ) {
01468 if ( simParams->pairInteractionSelf ) {
01469 for (i = 0; i < numAtoms; ++i ) {
01470 if ( a[i].partition != 1 ) continue;
01471 kineticEnergy += a[i].mass * a[i].velocity.length2();
01472 momentum += a[i].mass * a[i].velocity;
01473 angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01474 }
01475 }
01476 } else {
01477 for (i = 0; i < numAtoms; ++i ) {
01478 kineticEnergy += a[i].mass * a[i].velocity.length2();
01479 momentum += a[i].mass * a[i].velocity;
01480 angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01481 }
01482 if (simParams->drudeOn) {
01483 BigReal drudeComKE = 0.;
01484 BigReal drudeBondKE = 0.;
01485
01486 for (i = 0; i < numAtoms; i++) {
01487 if (i < numAtoms-1 &&
01488 a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01489 // i+1 is a Drude particle with parent i
01490
01491 // convert from Cartesian coordinates to (COM,bond) coordinates
01492 BigReal m_com = (a[i].mass + a[i+1].mass); // mass of COM
01493 BigReal m = a[i+1].mass / m_com; // mass ratio
01494 BigReal m_bond = a[i+1].mass * (1. - m); // mass of bond
01495 Vector v_bond = a[i+1].velocity - a[i].velocity; // vel of bond
01496 Vector v_com = a[i].velocity + m * v_bond; // vel of COM
01497
01498 drudeComKE += m_com * v_com.length2();
01499 drudeBondKE += m_bond * v_bond.length2();
01500
01501 i++; // +1 from loop, we've updated both particles
01502 }
01503 else {
01504 drudeComKE += a[i].mass * a[i].velocity.length2();
01505 }
01506 } // end for
01507
01508 drudeComKE *= 0.5;
01509 drudeBondKE *= 0.5;
01510 reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY)
01511 += drudeComKE;
01512 reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY)
01513 += drudeBondKE;
01514 } // end drudeOn
01515
01516 } // end else
01517
01518 kineticEnergy *= 0.5;
01519 reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
01520 ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
01521 ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
01522 }
01523
01524 #ifdef ALTVIRIAL
01525 // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
01526 {
01527 Tensor altVirial;
01528 for ( int i = 0; i < numAtoms; ++i ) {
01529 altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
01530 }
01531 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
01532 }
01533 {
01534 Tensor altVirial;
01535 for ( int i = 0; i < numAtoms; ++i ) {
01536 altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
01537 }
01538 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
01539 }
01540 {
01541 Tensor altVirial;
01542 for ( int i = 0; i < numAtoms; ++i ) {
01543 altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
01544 }
01545 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
01546 }
01547 #endif
01548
01549 {
01550 BigReal intKineticEnergy = 0;
01551 Tensor intVirialNormal;
01552 Tensor intVirialNbond;
01553 Tensor intVirialSlow;
01554
01555 int hgs;
01556 for ( int i = 0; i < numAtoms; i += hgs ) {
01557 #if CMK_BLUEGENEL
01558 CmiNetworkProgress();
01559 #endif
01560 hgs = a[i].hydrogenGroupSize;
01561 int j;
01562 BigReal m_cm = 0;
01563 Position x_cm(0,0,0);
01564 Velocity v_cm(0,0,0);
01565 for ( j = i; j < (i+hgs); ++j ) {
01566 m_cm += a[j].mass;
01567 x_cm += a[j].mass * a[j].position;
01568 v_cm += a[j].mass * a[j].velocity;
01569 }
01570 x_cm /= m_cm;
01571 v_cm /= m_cm;
01572 int fixedAtomsOn = simParams->fixedAtomsOn;
01573 if ( simParams->pairInteractionOn ) {
01574 int pairInteractionSelf = simParams->pairInteractionSelf;
01575 for ( j = i; j < (i+hgs); ++j ) {
01576 if ( a[j].partition != 1 &&
01577 ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
01578 // net force treated as zero for fixed atoms
01579 if ( fixedAtomsOn && a[j].atomFixed ) continue;
01580 BigReal mass = a[j].mass;
01581 Vector v = a[j].velocity;
01582 Vector dv = v - v_cm;
01583 intKineticEnergy += mass * (v * dv);
01584 Vector dx = a[j].position - x_cm;
01585 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01586 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01587 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01588 }
01589 } else {
01590 for ( j = i; j < (i+hgs); ++j ) {
01591 // net force treated as zero for fixed atoms
01592 if ( fixedAtomsOn && a[j].atomFixed ) continue;
01593 BigReal mass = a[j].mass;
01594 Vector v = a[j].velocity;
01595 Vector dv = v - v_cm;
01596 intKineticEnergy += mass * (v * dv);
01597 Vector dx = a[j].position - x_cm;
01598 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01599 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01600 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01601 }
01602 }
01603 }
01604
01605 intKineticEnergy *= 0.5;
01606 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
01607 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01608 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
01609 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
01610 }
01611
01612 if (pressureProfileReduction && simParams->useGroupPressure) {
01613 // subtract off internal virial term, calculated as for intVirial.
01614 int nslabs = simParams->pressureProfileSlabs;
01615 const Lattice &lattice = patch->lattice;
01616 BigReal idz = nslabs/lattice.c().z;
01617 BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01618 int useGroupPressure = simParams->useGroupPressure;
01619
01620 int hgs;
01621 for (int i=0; i<numAtoms; i += hgs) {
01622 int j;
01623 hgs = a[i].hydrogenGroupSize;
01624 BigReal m_cm = 0;
01625 Position x_cm(0,0,0);
01626 for (j=i; j< i+hgs; ++j) {
01627 m_cm += a[j].mass;
01628 x_cm += a[j].mass * a[j].position;
01629 }
01630 x_cm /= m_cm;
01631
01632 BigReal z = a[i].position.z;
01633 int slab = (int)floor((z-zmin)*idz);
01634 if (slab < 0) slab += nslabs;
01635 else if (slab >= nslabs) slab -= nslabs;
01636 int partition = a[i].partition;
01637 int ppoffset = 3*(slab + nslabs*partition);
01638 for (j=i; j < i+hgs; ++j) {
01639 BigReal mass = a[j].mass;
01640 Vector dx = a[j].position - x_cm;
01641 const Vector &fnormal = patch->f[Results::normal][j];
01642 const Vector &fnbond = patch->f[Results::nbond][j];
01643 const Vector &fslow = patch->f[Results::slow][j];
01644 BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
01645 BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
01646 BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
01647 pressureProfileReduction->item(ppoffset ) -= wxx;
01648 pressureProfileReduction->item(ppoffset+1) -= wyy;
01649 pressureProfileReduction->item(ppoffset+2) -= wzz;
01650 }
01651 }
01652 }
01653
01654 if ( simParams->fixedAtomsOn ) {
01655 Tensor fixVirialNormal;
01656 Tensor fixVirialNbond;
01657 Tensor fixVirialSlow;
01658 Vector fixForceNormal = 0;
01659 Vector fixForceNbond = 0;
01660 Vector fixForceSlow = 0;
01661
01662 for ( int j = 0; j < numAtoms; j++ ) {
01663 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01664 Vector dx = a[j].fixedPosition;
01665 // all negative because fixed atoms cancels these forces
01666 fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
01667 fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
01668 fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
01669 fixForceNormal -= patch->f[Results::normal][j];
01670 fixForceNbond -= patch->f[Results::nbond][j];
01671 fixForceSlow -= patch->f[Results::slow][j];
01672 }
01673 }
01674
01675 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
01676 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
01677 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
01678 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
01679 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
01680 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
01681 }
01682
01683 reduction->submit();
01684 if (pressureProfileReduction) pressureProfileReduction->submit();
01685 }
|
|
|
Definition at line 32 of file Sequencer.h. Referenced by HomePatch::doAtomMigration(), LdbCoordinator::rebalance(), and runComputeObjects(). 00032 { CthSuspend(); }
|
|
||||||||||||
|
Definition at line 1184 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, broadcast, SimpleBroadcastObject< T >::get(), Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, and FullAtom::velocity. Referenced by integrate(). 01185 {
01186 if ( simParams->tCoupleOn )
01187 {
01188 FullAtom *a = patch->atom.begin();
01189 int numAtoms = patch->numAtoms;
01190 BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
01191 Molecule *molecule = Node::Object()->molecule;
01192 BigReal dt = dt_fs * 0.001; // convert to ps
01193 coefficient *= dt;
01194 for ( int i = 0; i < numAtoms; ++i )
01195 {
01196 BigReal f1 = exp( coefficient * a[i].langevinParam );
01197 a[i].velocity *= f1;
01198 }
01199 }
01200 }
|
|
|
Definition at line 2003 of file Sequencer.C. Referenced by algorithm(), maximumMove(), and rattle1(). 02003 {
02004 CthFree(thread);
02005 CthSuspend();
02006 }
|
|
|
Definition at line 1993 of file Sequencer.C. References broadcast, SimpleBroadcastObject< T >::get(), and ControllerBroadcasts::traceBarrier. Referenced by integrate(). 01993 {
01994 broadcast->traceBarrier.get(step);
01995 }
|
|
|
Definition at line 24 of file Sequencer.h. |
|
|
Definition at line 73 of file Sequencer.h. Referenced by adaptTempUpdate(), and integrate(). |
|
|
Definition at line 84 of file Sequencer.h. Referenced by algorithm(), berendsenPressure(), and Sequencer(). |
|
|
Definition at line 106 of file Sequencer.h. Referenced by adaptTempUpdate(), algorithm(), berendsenPressure(), correctMomentum(), cycleBarrier(), langevinPiston(), minimize(), rescaleaccelMD(), rescaleVelocities(), Sequencer(), tcoupleVelocities(), and traceBarrier(). |
|
|
Definition at line 85 of file Sequencer.h. Referenced by algorithm(). |
|
|
Definition at line 105 of file Sequencer.h. Referenced by submitCollections(). |
|
|
Definition at line 108 of file Sequencer.h. Referenced by rebalanceLoad(). |
|
|
Definition at line 42 of file Sequencer.h. Referenced by runComputeObjects(). |
|
|
Definition at line 41 of file Sequencer.h. Referenced by algorithm(), rebalanceLoad(), and runComputeObjects(). |
|
|
|
Definition at line 103 of file Sequencer.h. Referenced by rattle1(), Sequencer(), submitHalfstep(), and submitReductions(). |
|
|
Definition at line 99 of file Sequencer.h. Referenced by langevinVelocities(), langevinVelocitiesBBK2(), reassignVelocities(), reinitVelocities(), and Sequencer(). |
|
|
Definition at line 102 of file Sequencer.h. Referenced by rattle1(), rattle2(), runComputeObjects(), Sequencer(), submitHalfstep(), submitMinimizeReductions(), submitMomentum(), and submitReductions(). |
|
|
Definition at line 78 of file Sequencer.h. Referenced by rescaleVelocities(), and Sequencer(). |
|
|
|
Definition at line 87 of file Sequencer.h. Referenced by integrate(), and langevinPiston(). |
1.3.9.1