#include <Sequencer.h>
|
|
|
Definition at line 60 of file Sequencer.C. 00061 {
00062 delete broadcast;
00063 delete reduction;
00064 delete pressureProfileReduction;
00065 delete random;
00066 }
|
|
||||||||||||||||||||
|
Definition at line 844 of file Sequencer.C. References ADD_TENSOR_OBJECT, HomePatch::addForceToMomentum(), patch, reduction, REDUCTION_VIRIAL_NORMAL, simParams, and SimParameters::watmodel. Referenced by integrate(). 00846 {
00847 #if CMK_VERSION_BLUEGENE
00848 CmiNetworkProgressAfter (0);
00849 #endif
00850 Tensor virial;
00851 Tensor *vp = ( pressure ? &virial : 0 );
00852 patch->addForceToMomentum(dt,ftag,useSaved,vp);
00853 if (simParams->watmodel == WAT_TIP4) ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
00854 }
|
|
|
Definition at line 346 of file Sequencer.C. References HomePatch::atom, CompAtom::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(). 00346 {
00347 FullAtom *atom = patch->atom.begin();
00348 int numAtoms = patch->numAtoms;
00349 Molecule *molecule = Node::Object()->molecule; // need its methods
00350 const BigReal movDragGlobVel = simParams->movDragGlobVel;
00351 const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
00352 Vector movDragVel, dragIncrement;
00353 for ( int i = 0; i < numAtoms; ++i )
00354 {
00355 // skip if fixed atom or zero drag attribute
00356 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00357 || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00358 molecule->get_movdrag_params(movDragVel, atom[i].id);
00359 dragIncrement = movDragGlobVel * movDragVel * dt;
00360 atom[i].position += dragIncrement;
00361 }
00362 }
|
|
|
Definition at line 365 of file Sequencer.C. References HomePatch::atom, CompAtom::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(). 00365 {
00366 FullAtom *atom = patch->atom.begin();
00367 int numAtoms = patch->numAtoms;
00368 Molecule *molecule = Node::Object()->molecule; // need its methods
00369 const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00370 const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
00371 BigReal rotDragVel, dAngle;
00372 Vector atomRadius;
00373 Vector rotDragAxis, rotDragPivot, dragIncrement;
00374 for ( int i = 0; i < numAtoms; ++i )
00375 {
00376 // skip if fixed atom or zero drag attribute
00377 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00378 || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00379 molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00380 dAngle = rotDragGlobVel * rotDragVel * dt;
00381 rotDragAxis /= rotDragAxis.length();
00382 atomRadius = atom[i].position - rotDragPivot;
00383 dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00384 atom[i].position += dragIncrement;
00385 }
00386 }
|
|
|
Definition at line 856 of file Sequencer.C. References HomePatch::addVelocityToPosition(), and patch. Referenced by integrate(). 00857 {
00858 #if CMK_VERSION_BLUEGENE
00859 CmiNetworkProgressAfter (0);
00860 #endif
00861 patch->addVelocityToPosition(dt);
00862 }
|
|
|
|
Definition at line 28 of file Sequencer.h. Referenced by LdbCoordinator::awakenSequencers(), HomePatch::boxClosed(), HomePatch::depositMigration(), and run(). 00028 { CthAwaken(thread); }
|
|
|
Definition at line 596 of file Sequencer.C. References Lattice::apply_transform(), HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), berendsenPressure_count, SimParameters::berendsenPressureFreq, SimParameters::berendsenPressureOn, BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), CompAtom::groupFixed, CompAtom::hydrogenGroupSize, Patch::lattice, FullAtom::mass, Patch::numAtoms, patch, Position, CompAtom::position, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, and SimParameters::useGroupPressure. Referenced by integrate(). 00597 {
00598 if ( simParams->berendsenPressureOn ) {
00599 berendsenPressure_count += 1;
00600 const int freq = simParams->berendsenPressureFreq;
00601 if ( ! (berendsenPressure_count % freq ) ) {
00602 berendsenPressure_count = 0;
00603 FullAtom *a = patch->atom.begin();
00604 int numAtoms = patch->numAtoms;
00605 Tensor factor = broadcast->positionRescaleFactor.get(step);
00606 patch->lattice.rescale(factor);
00607 if ( simParams->useGroupPressure )
00608 {
00609 int hgs;
00610 for ( int i = 0; i < numAtoms; i += hgs ) {
00611 int j;
00612 hgs = a[i].hydrogenGroupSize;
00613 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00614 for ( j = i; j < (i+hgs); ++j ) {
00615 a[j].position = patch->lattice.apply_transform(
00616 a[j].fixedPosition,a[j].transform);
00617 }
00618 continue;
00619 }
00620 BigReal m_cm = 0;
00621 Position x_cm(0,0,0);
00622 for ( j = i; j < (i+hgs); ++j ) {
00623 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00624 m_cm += a[j].mass;
00625 x_cm += a[j].mass * a[j].position;
00626 }
00627 x_cm /= m_cm;
00628 Position new_x_cm = x_cm;
00629 patch->lattice.rescale(new_x_cm,factor);
00630 Position delta_x_cm = new_x_cm - x_cm;
00631 for ( j = i; j < (i+hgs); ++j ) {
00632 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00633 a[j].position = patch->lattice.apply_transform(
00634 a[j].fixedPosition,a[j].transform);
00635 continue;
00636 }
00637 a[j].position += delta_x_cm;
00638 }
00639 }
00640 }
00641 else
00642 {
00643 for ( int i = 0; i < numAtoms; ++i )
00644 {
00645 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00646 a[i].position = patch->lattice.apply_transform(
00647 a[i].fixedPosition,a[i].transform);
00648 continue;
00649 }
00650 patch->lattice.rescale(a[i].position,factor);
00651 }
00652 }
00653 }
00654 } else {
00655 berendsenPressure_count = 0;
00656 }
00657 }
|
|
||||||||||||
|
Definition at line 495 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(). 00495 {
00496
00497 if ( simParams->fixedAtomsOn )
00498 NAMD_die("Cannot zero momentum when fixed atoms are present.");
00499
00500 const Vector dv = broadcast->momentumCorrection.get(step);
00501 const Vector dx = dv * ( drifttime / TIMEFACTOR );
00502
00503 FullAtom *a = patch->atom.begin();
00504 const int numAtoms = patch->numAtoms;
00505
00506 if ( simParams->zeroMomentumAlt ) {
00507 for ( int i = 0; i < numAtoms; ++i ) {
00508 BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00509 a[i].velocity += dv * rmass;
00510 a[i].position += dx * rmass;
00511 }
00512 } else {
00513 for ( int i = 0; i < numAtoms; ++i ) {
00514 a[i].velocity += dv;
00515 a[i].position += dx;
00516 }
00517 }
00518
00519 }
|
|
||||||||||||
|
Definition at line 1433 of file Sequencer.C. References broadcast. Referenced by integrate(). 01433 {
01434 #if USE_BARRIER
01435 if (doBarrier)
01436 broadcast->cycleBarrier.get(step);
01437 #endif
01438 }
|
|
|
Definition at line 130 of file Sequencer.C. References addForceToMomentum(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), berendsenPressure(), BigReal, Bool, SimParameters::commOnly, ComputeMgr::computeGlobalObject, Node::computeMgr, correctMomentum(), cycleBarrier(), Flags::doEnergy, Flags::doFullElectrostatics, Flags::doMolly, Flags::doNonbonded, SimParameters::dt, SimParameters::firstTimestep, Patch::flags, SimParameters::FMAOn, SimParameters::fullDirectOn, SimParameters::fullElectFrequency, langevinPiston(), langevinVelocitiesBBK1(), langevinVelocitiesBBK2(), Flags::maxForceMerged, Flags::maxForceUsed, maximumMove(), minimizationQuenchVelocity(), SimParameters::mollyOn, SimParameters::movDragOn, SimParameters::MTSAlgorithm, SimParameters::N, SimParameters::nonbondedFrequency, Node::Object(), SimParameters::outputEnergies, patch, SimParameters::PMEOn, rattle1(), SimParameters::reassignFreq, reassignVelocities(), rebalanceLoad(), rescaleVelocities(), SimParameters::rotDragOn, runComputeObjects(), saveForce(), ComputeGlobal::saveTotalForces(), simParams, slowFreq, Flags::step, SimParameters::stepsPerCycle, submitCollections(), submitHalfstep(), submitMomentum(), submitReductions(), SimParameters::tclForcesOn, tcoupleVelocities(), and SimParameters::zeroMomentum. Referenced by algorithm(). 00130 {
00131
00132 int &step = patch->flags.step;
00133 step = simParams->firstTimestep;
00134
00135 // drag switches
00136 const Bool rotDragOn = simParams->rotDragOn;
00137 const Bool movDragOn = simParams->movDragOn;
00138
00139 const int commOnly = simParams->commOnly;
00140
00141 int &maxForceUsed = patch->flags.maxForceUsed;
00142 int &maxForceMerged = patch->flags.maxForceMerged;
00143 maxForceUsed = Results::normal;
00144 maxForceMerged = Results::normal;
00145
00146 const int numberOfSteps = simParams->N;
00147 const int stepsPerCycle = simParams->stepsPerCycle;
00148 const BigReal timestep = simParams->dt;
00149
00150 // what MTS method?
00151 const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00152
00153 const int nonbondedFrequency = simParams->nonbondedFrequency;
00154 slowFreq = nonbondedFrequency;
00155 const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00156 int &doNonbonded = patch->flags.doNonbonded;
00157 doNonbonded = !(step%nonbondedFrequency);
00158 if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00159 if ( doNonbonded ) maxForceUsed = Results::nbond;
00160
00161 // Do we do full electrostatics?
00162 const int dofull = ( simParams->fullDirectOn ||
00163 simParams->FMAOn || simParams->PMEOn );
00164 const int fullElectFrequency = simParams->fullElectFrequency;
00165 if ( dofull ) slowFreq = fullElectFrequency;
00166 const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00167 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00168 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00169 if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00170 maxForceMerged = Results::slow;
00171 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00172 int &doMolly = patch->flags.doMolly;
00173 doMolly = simParams->mollyOn && doFullElectrostatics;
00174
00175 int zeroMomentum = simParams->zeroMomentum;
00176
00177 // Do we need to return forces to TCL script?
00178 int doTcl = simParams->tclForcesOn;
00179 ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00180
00181 // Bother to calculate energies?
00182 int &doEnergy = patch->flags.doEnergy;
00183 int energyFrequency = simParams->outputEnergies;
00184
00185 rattle1(0.,0); // enforce rigid bond constraints on initial positions
00186
00187 const int reassignFreq = simParams->reassignFreq;
00188 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00189 reassignVelocities(timestep,step);
00190 }
00191
00192 doEnergy = ! ( step % energyFrequency );
00193 runComputeObjects(1,step<numberOfSteps); // must migrate here!
00194 if ( staleForces || doTcl ) {
00195 if ( doNonbonded ) saveForce(Results::nbond);
00196 if ( doFullElectrostatics ) saveForce(Results::slow);
00197 }
00198 if ( ! commOnly ) {
00199 addForceToMomentum(-0.5*timestep);
00200 if (staleForces || doNonbonded)
00201 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00202 if (staleForces || doFullElectrostatics)
00203 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces);
00204 }
00205 minimizationQuenchVelocity();
00206 rattle1(-timestep,0);
00207 submitHalfstep(step);
00208 if ( ! commOnly ) {
00209 addForceToMomentum(timestep);
00210 if (staleForces || doNonbonded)
00211 addForceToMomentum(nbondstep,Results::nbond,staleForces,0);
00212 if (staleForces || doFullElectrostatics)
00213 addForceToMomentum(slowstep,Results::slow,staleForces,0);
00214 }
00215 rattle1(timestep,1);
00216 if (doTcl) // include constraint forces
00217 computeGlobal->saveTotalForces(patch);
00218 submitHalfstep(step);
00219 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00220 if ( ! commOnly ) {
00221 addForceToMomentum(-0.5*timestep);
00222 if (staleForces || doNonbonded)
00223 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00224 if (staleForces || doFullElectrostatics)
00225 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00226 }
00227 submitReductions(step);
00228 rebalanceLoad(step);
00229
00230 for ( ++step; step <= numberOfSteps; ++step )
00231 {
00232
00233 rescaleVelocities(step);
00234 tcoupleVelocities(timestep,step);
00235 berendsenPressure(step);
00236
00237 if ( ! commOnly ) {
00238 addForceToMomentum(0.5*timestep);
00239 if (staleForces || doNonbonded)
00240 addForceToMomentum(0.5*nbondstep,Results::nbond,staleForces,1);
00241 if (staleForces || doFullElectrostatics)
00242 addForceToMomentum(0.5*slowstep,Results::slow,staleForces,1);
00243 }
00244
00245 /* reassignment based on half-step velocities
00246 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00247 addVelocityToPosition(0.5*timestep);
00248 reassignVelocities(timestep,step);
00249 addVelocityToPosition(0.5*timestep);
00250 rattle1(0.,0);
00251 rattle1(-timestep,0);
00252 addVelocityToPosition(-1.0*timestep);
00253 rattle1(timestep,0);
00254 } */
00255
00256 maximumMove(timestep);
00257 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00258 langevinPiston(step);
00259 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00260
00261 minimizationQuenchVelocity();
00262
00263 doNonbonded = !(step%nonbondedFrequency);
00264 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00265
00266 if ( zeroMomentum && doFullElectrostatics )
00267 correctMomentum(step,slowstep);
00268
00269 submitHalfstep(step);
00270
00271 doMolly = simParams->mollyOn && doFullElectrostatics;
00272
00273 maxForceUsed = Results::normal;
00274 if ( doNonbonded ) maxForceUsed = Results::nbond;
00275 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00276
00277 // Migrate Atoms on stepsPerCycle
00278 doEnergy = ! ( step % energyFrequency );
00279 runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00280 if ( staleForces || doTcl ) {
00281 if ( doNonbonded ) saveForce(Results::nbond);
00282 if ( doFullElectrostatics ) saveForce(Results::slow);
00283 }
00284
00285 // reassignment based on full-step velocities
00286 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00287 reassignVelocities(timestep,step);
00288 addForceToMomentum(-0.5*timestep);
00289 if (staleForces || doNonbonded)
00290 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00291 if (staleForces || doFullElectrostatics)
00292 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,0);
00293 rattle1(-timestep,0);
00294 }
00295
00296 if ( ! commOnly ) {
00297 langevinVelocitiesBBK1(timestep);
00298 addForceToMomentum(timestep);
00299 if (staleForces || doNonbonded)
00300 addForceToMomentum(nbondstep,Results::nbond,staleForces,1);
00301 if (staleForces || doFullElectrostatics)
00302 addForceToMomentum(slowstep,Results::slow,staleForces,1);
00303 langevinVelocitiesBBK2(timestep);
00304 }
00305
00306 // add drag to each atom's positions
00307 if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00308 if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00309
00310 rattle1(timestep,1);
00311 if (doTcl) // include constraint forces
00312 computeGlobal->saveTotalForces(patch);
00313
00314 submitHalfstep(step);
00315 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00316
00317 if ( ! commOnly ) {
00318 addForceToMomentum(-0.5*timestep);
00319 if (staleForces || doNonbonded)
00320 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00321 if (staleForces || doFullElectrostatics)
00322 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00323 }
00324
00325 // rattle2(timestep,step);
00326
00327 submitReductions(step);
00328 submitCollections(step);
00329
00330 #if CYCLE_BARRIER
00331 cycleBarrier(!((step+1) % stepsPerCycle), step);
00332 #elif PME_BARRIER
00333 cycleBarrier(doFullElectrostatics, step);
00334 #endif
00335
00336 rebalanceLoad(step);
00337
00338 #if PME_BARRIER
00339 // a step before PME
00340 cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00341 #endif
00342 }
00343 }
|
|
|
Definition at line 659 of file Sequencer.C. References Lattice::apply_transform(), HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), CompAtom::groupFixed, CompAtom::hydrogenGroupSize, Molecule::is_atom_exPressure(), 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(). 00660 {
00661 if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
00662 {
00663 FullAtom *a = patch->atom.begin();
00664 int numAtoms = patch->numAtoms;
00665 Tensor factor = broadcast->positionRescaleFactor.get(step);
00666 // JCP FIX THIS!!!
00667 Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
00668 patch->lattice.rescale(factor);
00669 Molecule *mol = Node::Object()->molecule;
00670 if ( simParams->useGroupPressure )
00671 {
00672 int hgs;
00673 for ( int i = 0; i < numAtoms; i += hgs ) {
00674 int j;
00675 hgs = a[i].hydrogenGroupSize;
00676 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00677 for ( j = i; j < (i+hgs); ++j ) {
00678 a[j].position = patch->lattice.apply_transform(
00679 a[j].fixedPosition,a[j].transform);
00680 }
00681 continue;
00682 }
00683 BigReal m_cm = 0;
00684 Position x_cm(0,0,0);
00685 Velocity v_cm(0,0,0);
00686 for ( j = i; j < (i+hgs); ++j ) {
00687 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00688 m_cm += a[j].mass;
00689 x_cm += a[j].mass * a[j].position;
00690 v_cm += a[j].mass * a[j].velocity;
00691 }
00692 x_cm /= m_cm;
00693 Position new_x_cm = x_cm;
00694 patch->lattice.rescale(new_x_cm,factor);
00695 Position delta_x_cm = new_x_cm - x_cm;
00696 v_cm /= m_cm;
00697 Velocity delta_v_cm;
00698 delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
00699 delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
00700 delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
00701 for ( j = i; j < (i+hgs); ++j ) {
00702 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00703 a[j].position = patch->lattice.apply_transform(
00704 a[j].fixedPosition,a[j].transform);
00705 continue;
00706 }
00707 if ( mol->is_atom_exPressure(a[j].id) ) continue;
00708 a[j].position += delta_x_cm;
00709 a[j].velocity += delta_v_cm;
00710 }
00711 }
00712 }
00713 else
00714 {
00715 for ( int i = 0; i < numAtoms; ++i )
00716 {
00717 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00718 a[i].position = patch->lattice.apply_transform(
00719 a[i].fixedPosition,a[i].transform);
00720 continue;
00721 }
00722 if ( mol->is_atom_exPressure(a[i].id) ) continue;
00723 patch->lattice.rescale(a[i].position,factor);
00724 a[i].velocity.x *= velFactor.x;
00725 a[i].velocity.y *= velFactor.y;
00726 a[i].velocity.z *= velFactor.z;
00727 }
00728 }
00729 }
00730 }
|
|
|
Definition at line 521 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, Random::gaussian_vector(), CompAtom::id, Molecule::langevin_param(), SimParameters::langevinOn, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, Node::molecule, Patch::numAtoms, Node::Object(), patch, random, simParams, and FullAtom::velocity. 00522 {
00523 if ( simParams->langevinOn )
00524 {
00525 FullAtom *a = patch->atom.begin();
00526 int numAtoms = patch->numAtoms;
00527 Molecule *molecule = Node::Object()->molecule;
00528 BigReal dt = dt_fs * 0.001; // convert to ps
00529 BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00530
00531 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00532 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00533
00534 for ( int i = 0; i < numAtoms; ++i )
00535 {
00536 int aid = a[i].id;
00537 BigReal f1 = exp( -1. * dt * molecule->langevin_param(aid) );
00538 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
00539 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00540
00541 a[i].velocity *= f1;
00542 a[i].velocity += f2 * random->gaussian_vector();
00543 }
00544 }
00545 }
|
|
|
Definition at line 547 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, CompAtom::id, Molecule::langevin_param(), SimParameters::langevinOn, Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, and FullAtom::velocity. Referenced by integrate(). 00548 {
00549 if ( simParams->langevinOn )
00550 {
00551 FullAtom *a = patch->atom.begin();
00552 int numAtoms = patch->numAtoms;
00553 Molecule *molecule = Node::Object()->molecule;
00554 BigReal dt = dt_fs * 0.001; // convert to ps
00555 for ( int i = 0; i < numAtoms; ++i )
00556 {
00557 int aid = a[i].id;
00558 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00559 if ( ! dt_gamma ) continue;
00560
00561 a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00562 }
00563 }
00564 }
|
|
|
Definition at line 567 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, Random::gaussian_vector(), CompAtom::id, Molecule::langevin_param(), SimParameters::langevinOn, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, Node::molecule, Patch::numAtoms, Node::Object(), patch, random, rattle1(), simParams, and FullAtom::velocity. Referenced by integrate(). 00568 {
00569 if ( simParams->langevinOn )
00570 {
00571 rattle1(dt_fs,1); // conserve momentum if gammas differ
00572
00573 FullAtom *a = patch->atom.begin();
00574 int numAtoms = patch->numAtoms;
00575 Molecule *molecule = Node::Object()->molecule;
00576 BigReal dt = dt_fs * 0.001; // convert to ps
00577 BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00578
00579 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00580 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00581
00582 for ( int i = 0; i < numAtoms; ++i )
00583 {
00584 int aid = a[i].id;
00585 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00586 if ( ! dt_gamma ) continue;
00587
00588 a[i].velocity += random->gaussian_vector() *
00589 sqrt( 2 * dt_gamma * kbT *
00590 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00591 a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00592 }
00593 }
00594 }
|
|
|
Definition at line 893 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, SimParameters::cutoff, Node::enableEarlyExit(), CompAtom::id, iERROR(), iout, Vector::length(), Vector::length2(), SimParameters::maximumMove, Patch::numAtoms, Node::Object(), patch, PDBVELFACTOR, simParams, terminate(), and FullAtom::velocity. Referenced by integrate(). 00894 {
00895 FullAtom *a = patch->atom.begin();
00896 int numAtoms = patch->numAtoms;
00897 if ( simParams->maximumMove ) {
00898 const BigReal dt = timestep / TIMEFACTOR;
00899 const BigReal maxvel = simParams->maximumMove / dt;
00900 const BigReal maxvel2 = maxvel * maxvel;
00901 for ( int i=0; i<numAtoms; ++i ) {
00902 if ( a[i].velocity.length2() > maxvel2 ) {
00903 a[i].velocity *= ( maxvel / a[i].velocity.length() );
00904 }
00905 }
00906 } else {
00907 const BigReal dt = timestep / TIMEFACTOR;
00908 const BigReal maxvel = simParams->cutoff / dt;
00909 const BigReal maxvel2 = maxvel * maxvel;
00910 int killme = 0;
00911 for ( int i=0; i<numAtoms; ++i ) {
00912 killme = killme || ( a[i].velocity.length2() > maxvel2 );
00913 }
00914 if ( killme ) {
00915 for ( int i=0; i<numAtoms; ++i ) {
00916 if ( a[i].velocity.length2() > maxvel2 ) {
00917 iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
00918 << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
00919 << ( PDBVELFACTOR * maxvel ) << ")\n" << endi;
00920 }
00921 }
00922 iout << iERROR <<
00923 "Atoms moving too fast; simulation has become unstable.\n" << endi;
00924 Node::Object()->enableEarlyExit();
00925 terminate();
00926 }
00927 }
00928 }
|
|
|
Definition at line 930 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::minimizeOn, Patch::numAtoms, patch, simParams, and FullAtom::velocity. Referenced by integrate(). 00931 {
00932 if ( simParams->minimizeOn ) {
00933 FullAtom *a = patch->atom.begin();
00934 int numAtoms = patch->numAtoms;
00935 for ( int i=0; i<numAtoms; ++i ) {
00936 a[i].velocity = 0.;
00937 }
00938 }
00939 }
|
|
|
Definition at line 388 of file Sequencer.C. References BigReal, broadcast, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doMolly, Flags::doNonbonded, SimParameters::firstTimestep, Patch::flags, SimParameters::FMAOn, SimParameters::fullDirectOn, SimpleBroadcastObject< T >::get(), Flags::maxForceMerged, Flags::maxForceUsed, ControllerBroadcasts::minimizeCoefficient, SimParameters::mollyOn, SimParameters::N, newMinimizeDirection(), newMinimizePosition(), patch, SimParameters::PMEOn, quenchVelocities(), rebalanceLoad(), runComputeObjects(), simParams, Flags::step, submitCollections(), and submitMinimizeReductions(). Referenced by algorithm(). 00388 {
00389 const int numberOfSteps = simParams->N;
00390 int &step = patch->flags.step;
00391 step = simParams->firstTimestep;
00392
00393 int &maxForceUsed = patch->flags.maxForceUsed;
00394 int &maxForceMerged = patch->flags.maxForceMerged;
00395 maxForceUsed = Results::normal;
00396 maxForceMerged = Results::normal;
00397 int &doNonbonded = patch->flags.doNonbonded;
00398 doNonbonded = 1;
00399 maxForceUsed = Results::nbond;
00400 maxForceMerged = Results::nbond;
00401 const int dofull = ( simParams->fullDirectOn ||
00402 simParams->FMAOn || simParams->PMEOn );
00403 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00404 doFullElectrostatics = dofull;
00405 if ( dofull ) {
00406 maxForceMerged = Results::slow;
00407 maxForceUsed = Results::slow;
00408 }
00409 int &doMolly = patch->flags.doMolly;
00410 doMolly = simParams->mollyOn && doFullElectrostatics;
00411 int &doEnergy = patch->flags.doEnergy;
00412 doEnergy = 1;
00413
00414 runComputeObjects(1,0); // must migrate here!
00415
00416 submitMinimizeReductions(step);
00417 rebalanceLoad(step);
00418
00419 int minSeq = 0;
00420 for ( ++step; step <= numberOfSteps; ++step ) {
00421 BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00422 if ( ! c ) { // new direction
00423 c = broadcast->minimizeCoefficient.get(minSeq++);
00424 newMinimizeDirection(c); // v = c * v + f
00425 c = broadcast->minimizeCoefficient.get(minSeq++);
00426 } // same direction
00427 newMinimizePosition(c); // x = x + c * v
00428
00429 runComputeObjects(1,0);
00430 submitMinimizeReductions(step);
00431 submitCollections(step, 1); // write out zeros for velocities
00432 rebalanceLoad(step);
00433 }
00434 quenchVelocities(); // zero out bogus velocity
00435 }
|
|
|
Definition at line 438 of file Sequencer.C. References HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), Patch::f, SimParameters::fixedAtomsOn, Force, Patch::numAtoms, |