#include <Sequencer.h>
|
|
|
Definition at line 69 of file Sequencer.C. 00070 {
00071 delete broadcast;
00072 delete reduction;
00073 delete pressureProfileReduction;
00074 delete random;
00075 }
|
|
||||||||||||||||||||
|
Definition at line 1047 of file Sequencer.C. References HomePatch::addForceToMomentum(), and patch. Referenced by integrate(). 01049 {
01050 #if CMK_BLUEGENEL
01051 CmiNetworkProgressAfter (0);
01052 #endif
01053 patch->addForceToMomentum(dt,ftag,useSaved);
01054 }
|
|
|
Definition at line 395 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(). 00395 {
00396 FullAtom *atom = patch->atom.begin();
00397 int numAtoms = patch->numAtoms;
00398 Molecule *molecule = Node::Object()->molecule; // need its methods
00399 const BigReal movDragGlobVel = simParams->movDragGlobVel;
00400 const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
00401 Vector movDragVel, dragIncrement;
00402 for ( int i = 0; i < numAtoms; ++i )
00403 {
00404 // skip if fixed atom or zero drag attribute
00405 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00406 || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00407 molecule->get_movdrag_params(movDragVel, atom[i].id);
00408 dragIncrement = movDragGlobVel * movDragVel * dt;
00409 atom[i].position += dragIncrement;
00410 }
00411 }
|
|
|
Definition at line 414 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(). 00414 {
00415 FullAtom *atom = patch->atom.begin();
00416 int numAtoms = patch->numAtoms;
00417 Molecule *molecule = Node::Object()->molecule; // need its methods
00418 const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00419 const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
00420 BigReal rotDragVel, dAngle;
00421 Vector atomRadius;
00422 Vector rotDragAxis, rotDragPivot, dragIncrement;
00423 for ( int i = 0; i < numAtoms; ++i )
00424 {
00425 // skip if fixed atom or zero drag attribute
00426 if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
00427 || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00428 molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00429 dAngle = rotDragGlobVel * rotDragVel * dt;
00430 rotDragAxis /= rotDragAxis.length();
00431 atomRadius = atom[i].position - rotDragPivot;
00432 dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00433 atom[i].position += dragIncrement;
00434 }
00435 }
|
|
|
Definition at line 1056 of file Sequencer.C. References HomePatch::addVelocityToPosition(), and patch. Referenced by integrate(). 01057 {
01058 #if CMK_BLUEGENEL
01059 CmiNetworkProgressAfter (0);
01060 #endif
01061 patch->addVelocityToPosition(dt);
01062 }
|
|
|
|
Definition at line 28 of file Sequencer.h. Referenced by LdbCoordinator::awakenSequencers(), HomePatch::boxClosed(), HomePatch::depositMigration(), and run(). 00028 { CthAwaken(thread); }
|
|
|
Definition at line 785 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(). 00786 {
00787 if ( simParams->berendsenPressureOn ) {
00788 berendsenPressure_count += 1;
00789 const int freq = simParams->berendsenPressureFreq;
00790 if ( ! (berendsenPressure_count % freq ) ) {
00791 berendsenPressure_count = 0;
00792 FullAtom *a = patch->atom.begin();
00793 int numAtoms = patch->numAtoms;
00794 Tensor factor = broadcast->positionRescaleFactor.get(step);
00795 patch->lattice.rescale(factor);
00796 if ( simParams->useGroupPressure )
00797 {
00798 int hgs;
00799 for ( int i = 0; i < numAtoms; i += hgs ) {
00800 int j;
00801 hgs = a[i].hydrogenGroupSize;
00802 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00803 for ( j = i; j < (i+hgs); ++j ) {
00804 a[j].position = patch->lattice.apply_transform(
00805 a[j].fixedPosition,a[j].transform);
00806 }
00807 continue;
00808 }
00809 BigReal m_cm = 0;
00810 Position x_cm(0,0,0);
00811 for ( j = i; j < (i+hgs); ++j ) {
00812 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00813 m_cm += a[j].mass;
00814 x_cm += a[j].mass * a[j].position;
00815 }
00816 x_cm /= m_cm;
00817 Position new_x_cm = x_cm;
00818 patch->lattice.rescale(new_x_cm,factor);
00819 Position delta_x_cm = new_x_cm - x_cm;
00820 for ( j = i; j < (i+hgs); ++j ) {
00821 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00822 a[j].position = patch->lattice.apply_transform(
00823 a[j].fixedPosition,a[j].transform);
00824 continue;
00825 }
00826 a[j].position += delta_x_cm;
00827 }
00828 }
00829 }
00830 else
00831 {
00832 for ( int i = 0; i < numAtoms; ++i )
00833 {
00834 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00835 a[i].position = patch->lattice.apply_transform(
00836 a[i].fixedPosition,a[i].transform);
00837 continue;
00838 }
00839 patch->lattice.rescale(a[i].position,factor);
00840 }
00841 }
00842 }
00843 } else {
00844 berendsenPressure_count = 0;
00845 }
00846 }
|
|
||||||||||||
|
Definition at line 573 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(). 00573 {
00574
00575 if ( simParams->fixedAtomsOn )
00576 NAMD_die("Cannot zero momentum when fixed atoms are present.");
00577
00578 const Vector dv = broadcast->momentumCorrection.get(step);
00579 const Vector dx = dv * ( drifttime / TIMEFACTOR );
00580
00581 FullAtom *a = patch->atom.begin();
00582 const int numAtoms = patch->numAtoms;
00583
00584 if ( simParams->zeroMomentumAlt ) {
00585 for ( int i = 0; i < numAtoms; ++i ) {
00586 BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00587 a[i].velocity += dv * rmass;
00588 a[i].position += dx * rmass;
00589 }
00590 } else {
00591 for ( int i = 0; i < numAtoms; ++i ) {
00592 a[i].velocity += dv;
00593 a[i].position += dx;
00594 }
00595 }
00596
00597 }
|
|
||||||||||||
|
Definition at line 1695 of file Sequencer.C. References broadcast. Referenced by integrate(). 01695 {
01696 #if USE_BARRIER
01697 if (doBarrier)
01698 broadcast->cycleBarrier.get(step);
01699 #endif
01700 }
|
|
|
Definition at line 139 of file Sequencer.C. References addForceToMomentum(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), berendsenPressure(), BigReal, Bool, SimParameters::colvarsOn, 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(), redistrib_swm4_forces(), redistrib_tip4p_forces(), rescaleVelocities(), SimParameters::rotDragOn, runComputeObjects(), saveForce(), ComputeGlobal::saveTotalForces(), simParams, slowFreq, Flags::step, SimParameters::stepsPerCycle, submitCollections(), submitHalfstep(), submitMomentum(), submitReductions(), SimParameters::tclForcesOn, tcoupleVelocities(), SimParameters::watmodel, and SimParameters::zeroMomentum. Referenced by algorithm(). 00139 {
00140
00141 // patch->write_tip4_props();
00142
00143 int &step = patch->flags.step;
00144 step = simParams->firstTimestep;
00145
00146 // drag switches
00147 const Bool rotDragOn = simParams->rotDragOn;
00148 const Bool movDragOn = simParams->movDragOn;
00149
00150 const int commOnly = simParams->commOnly;
00151
00152 int &maxForceUsed = patch->flags.maxForceUsed;
00153 int &maxForceMerged = patch->flags.maxForceMerged;
00154 maxForceUsed = Results::normal;
00155 maxForceMerged = Results::normal;
00156
00157 const int numberOfSteps = simParams->N;
00158 const int stepsPerCycle = simParams->stepsPerCycle;
00159 const BigReal timestep = simParams->dt;
00160
00161 // what MTS method?
00162 const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00163
00164 const int nonbondedFrequency = simParams->nonbondedFrequency;
00165 slowFreq = nonbondedFrequency;
00166 const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00167 int &doNonbonded = patch->flags.doNonbonded;
00168 doNonbonded = !(step%nonbondedFrequency);
00169 if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00170 if ( doNonbonded ) maxForceUsed = Results::nbond;
00171
00172 // Do we do full electrostatics?
00173 const int dofull = ( simParams->fullDirectOn ||
00174 simParams->FMAOn || simParams->PMEOn );
00175 const int fullElectFrequency = simParams->fullElectFrequency;
00176 if ( dofull ) slowFreq = fullElectFrequency;
00177 const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00178 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00179 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00180 if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00181 maxForceMerged = Results::slow;
00182 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00183 int &doMolly = patch->flags.doMolly;
00184 doMolly = simParams->mollyOn && doFullElectrostatics;
00185
00186 int zeroMomentum = simParams->zeroMomentum;
00187
00188 // Do we need to return forces to TCL script or Colvar module?
00189 int doTcl = simParams->tclForcesOn;
00190 int doColvars = simParams->colvarsOn;
00191 ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00192
00193 // Bother to calculate energies?
00194 int &doEnergy = patch->flags.doEnergy;
00195 int energyFrequency = simParams->outputEnergies;
00196
00197 // printf("Doing initial rattle\n");
00198 rattle1(0.,0); // enforce rigid bond constraints on initial positions
00199
00200 const int reassignFreq = simParams->reassignFreq;
00201 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00202 reassignVelocities(timestep,step);
00203 }
00204
00205 doEnergy = ! ( step % energyFrequency );
00206 runComputeObjects(1,step<numberOfSteps); // must migrate here!
00207
00208 // Redistribute forces, if needed for lonepairs
00209 if (simParams->watmodel == WAT_TIP4) {
00210 redistrib_tip4p_forces(Results::normal, 0);
00211 redistrib_tip4p_forces(Results::nbond, 0);
00212 redistrib_tip4p_forces(Results::slow, 0);
00213 }
00214 else if (simParams->watmodel == WAT_SWM4) {
00215 redistrib_swm4_forces(Results::normal, 0);
00216 redistrib_swm4_forces(Results::nbond, 0);
00217 redistrib_swm4_forces(Results::slow, 0);
00218 }
00219
00220 if ( staleForces || doTcl || doColvars ) {
00221 if ( doNonbonded ) saveForce(Results::nbond);
00222 if ( doFullElectrostatics ) saveForce(Results::slow);
00223 }
00224 if ( ! commOnly ) {
00225 addForceToMomentum(-0.5*timestep);
00226 if (staleForces || doNonbonded)
00227 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00228 if (staleForces || doFullElectrostatics)
00229 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces);
00230 }
00231 minimizationQuenchVelocity();
00232 rattle1(-timestep,0);
00233 submitHalfstep(step);
00234 if ( ! commOnly ) {
00235 addForceToMomentum(timestep);
00236 if (staleForces || doNonbonded)
00237 addForceToMomentum(nbondstep,Results::nbond,staleForces,0);
00238 if (staleForces || doFullElectrostatics)
00239 addForceToMomentum(slowstep,Results::slow,staleForces,0);
00240 }
00241 rattle1(timestep,1);
00242 if (doTcl || doColvars) // include constraint forces
00243 computeGlobal->saveTotalForces(patch);
00244 submitHalfstep(step);
00245 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00246 if ( ! commOnly ) {
00247 addForceToMomentum(-0.5*timestep);
00248 if (staleForces || doNonbonded)
00249 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00250 if (staleForces || doFullElectrostatics)
00251 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00252 }
00253 submitReductions(step);
00254 rebalanceLoad(step);
00255
00256 for ( ++step; step <= numberOfSteps; ++step )
00257 {
00258
00259 rescaleVelocities(step);
00260 tcoupleVelocities(timestep,step);
00261 berendsenPressure(step);
00262
00263 if ( ! commOnly ) {
00264 addForceToMomentum(0.5*timestep);
00265 if (staleForces || doNonbonded)
00266 addForceToMomentum(0.5*nbondstep,Results::nbond,staleForces,1);
00267 if (staleForces || doFullElectrostatics)
00268 addForceToMomentum(0.5*slowstep,Results::slow,staleForces,1);
00269 }
00270
00271 /* reassignment based on half-step velocities
00272 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00273 addVelocityToPosition(0.5*timestep);
00274 reassignVelocities(timestep,step);
00275 addVelocityToPosition(0.5*timestep);
00276 rattle1(0.,0);
00277 rattle1(-timestep,0);
00278 addVelocityToPosition(-1.0*timestep);
00279 rattle1(timestep,0);
00280 } */
00281
00282 maximumMove(timestep);
00283 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00284 langevinPiston(step);
00285 if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00286
00287 minimizationQuenchVelocity();
00288
00289 doNonbonded = !(step%nonbondedFrequency);
00290 doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00291
00292 if ( zeroMomentum && doFullElectrostatics )
00293 correctMomentum(step,slowstep);
00294
00295 submitHalfstep(step);
00296
00297 doMolly = simParams->mollyOn && doFullElectrostatics;
00298
00299 maxForceUsed = Results::normal;
00300 if ( doNonbonded ) maxForceUsed = Results::nbond;
00301 if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00302
00303 // Migrate Atoms on stepsPerCycle
00304 doEnergy = ! ( step % energyFrequency );
00305 runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00306
00307 // Redistribute forces, if needed for lonepairs
00308 if (simParams->watmodel == WAT_TIP4) {
00309 redistrib_tip4p_forces(Results::normal, 1);
00310 if (doNonbonded) redistrib_tip4p_forces(Results::nbond, 1);
00311 if (doFullElectrostatics) redistrib_tip4p_forces(Results::slow, 1);
00312 }
00313 else if (simParams->watmodel == WAT_SWM4) {
00314 redistrib_swm4_forces(Results::normal, 1);
00315 if (doNonbonded) redistrib_swm4_forces(Results::nbond, 1);
00316 if (doFullElectrostatics) redistrib_swm4_forces(Results::slow, 1);
00317 }
00318
00319 if ( staleForces || doTcl || doColvars ) {
00320 if ( doNonbonded ) saveForce(Results::nbond);
00321 if ( doFullElectrostatics ) saveForce(Results::slow);
00322 }
00323
00324 // reassignment based on full-step velocities
00325 if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00326 reassignVelocities(timestep,step);
00327 addForceToMomentum(-0.5*timestep);
00328 if (staleForces || doNonbonded)
00329 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00330 if (staleForces || doFullElectrostatics)
00331 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,0);
00332 rattle1(-timestep,0);
00333 }
00334
00335 if ( ! commOnly ) {
00336 langevinVelocitiesBBK1(timestep);
00337 addForceToMomentum(timestep);
00338 if (staleForces || doNonbonded)
00339 addForceToMomentum(nbondstep,Results::nbond,staleForces,1);
00340 if (staleForces || doFullElectrostatics)
00341 addForceToMomentum(slowstep,Results::slow,staleForces,1);
00342 langevinVelocitiesBBK2(timestep);
00343 }
00344
00345 // add drag to each atom's positions
00346 if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00347 if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00348
00349 rattle1(timestep,1);
00350 if (doTcl || doColvars) // include constraint forces
00351 computeGlobal->saveTotalForces(patch);
00352
00353 submitHalfstep(step);
00354 if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00355
00356 if ( ! commOnly ) {
00357 addForceToMomentum(-0.5*timestep);
00358 if (staleForces || doNonbonded)
00359 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00360 if (staleForces || doFullElectrostatics)
00361 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00362 }
00363
00364 // rattle2(timestep,step);
00365
00366 submitReductions(step);
00367 submitCollections(step);
00368
00369 #if CYCLE_BARRIER
00370 cycleBarrier(!((step+1) % stepsPerCycle), step);
00371 #elif PME_BARRIER
00372 cycleBarrier(doFullElectrostatics, step);
00373 #elif STEP_BARRIER
00374 cycleBarrier(1, step);
00375 #endif
00376
00377 rebalanceLoad(step);
00378
00379 #if PME_BARRIER
00380 // a step before PME
00381 cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00382 #endif
00383
00384 #if USE_HPM
00385 if(step == START_HPM_STEP)
00386 (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
00387
00388 if(step == STOP_HPM_STEP)
00389 (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
00390 #endif
00391 }
00392 }
|
|
|
Definition at line 848 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(). 00849 {
00850 if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
00851 {
00852 FullAtom *a = patch->atom.begin();
00853 int numAtoms = patch->numAtoms;
00854 Tensor factor = broadcast->positionRescaleFactor.get(step);
00855 // JCP FIX THIS!!!
00856 Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
00857 patch->lattice.rescale(factor);
00858 Molecule *mol = Node::Object()->molecule;
00859 if ( simParams->useGroupPressure )
00860 {
00861 int hgs;
00862 for ( int i = 0; i < numAtoms; i += hgs ) {
00863 int j;
00864 hgs = a[i].hydrogenGroupSize;
00865 if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00866 for ( j = i; j < (i+hgs); ++j ) {
00867 a[j].position = patch->lattice.apply_transform(
00868 a[j].fixedPosition,a[j].transform);
00869 }
00870 continue;
00871 }
00872 BigReal m_cm = 0;
00873 Position x_cm(0,0,0);
00874 Velocity v_cm(0,0,0);
00875 for ( j = i; j < (i+hgs); ++j ) {
00876 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00877 m_cm += a[j].mass;
00878 x_cm += a[j].mass * a[j].position;
00879 v_cm += a[j].mass * a[j].velocity;
00880 }
00881 x_cm /= m_cm;
00882 Position new_x_cm = x_cm;
00883 patch->lattice.rescale(new_x_cm,factor);
00884 Position delta_x_cm = new_x_cm - x_cm;
00885 v_cm /= m_cm;
00886 Velocity delta_v_cm;
00887 delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
00888 delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
00889 delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
00890 for ( j = i; j < (i+hgs); ++j ) {
00891 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00892 a[j].position = patch->lattice.apply_transform(
00893 a[j].fixedPosition,a[j].transform);
00894 continue;
00895 }
00896 if ( mol->is_atom_exPressure(a[j].id) ) continue;
00897 a[j].position += delta_x_cm;
00898 a[j].velocity += delta_v_cm;
00899 }
00900 }
00901 }
00902 else
00903 {
00904 for ( int i = 0; i < numAtoms; ++i )
00905 {
00906 if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00907 a[i].position = patch->lattice.apply_transform(
00908 a[i].fixedPosition,a[i].transform);
00909 continue;
00910 }
00911 if ( mol->is_atom_exPressure(a[i].id) ) continue;
00912 patch->lattice.rescale(a[i].position,factor);
00913 a[i].velocity.x *= velFactor.x;
00914 a[i].velocity.y *= velFactor.y;
00915 a[i].velocity.z *= velFactor.z;
00916 }
00917 }
00918 }
00919 }
|
|
|
Definition at line 599 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, Random::gaussian_vector(), CompAtomExt::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. 00600 {
00601 if ( simParams->langevinOn )
00602 {
00603 FullAtom *a = patch->atom.begin();
00604 int numAtoms = patch->numAtoms;
00605 Molecule *molecule = Node::Object()->molecule;
00606 BigReal dt = dt_fs * 0.001; // convert to ps
00607 BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00608
00609 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00610 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00611
00612 for ( int i = 0; i < numAtoms; ++i )
00613 {
00614 int aid = a[i].id;
00615 BigReal f1 = exp( -1. * dt * molecule->langevin_param(aid) );
00616 BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
00617 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00618
00619 a[i].velocity *= f1;
00620 a[i].velocity += f2 * random->gaussian_vector();
00621 }
00622 }
00623 }
|
|
|
Definition at line 625 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, SimParameters::drudeOn, CompAtomExt::id, Molecule::langevin_param(), SimParameters::langevinOn, FullAtom::mass, Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, and FullAtom::velocity. Referenced by integrate(). 00626 {
00627 if ( simParams->langevinOn )
00628 {
00629 FullAtom *a = patch->atom.begin();
00630 int numAtoms = patch->numAtoms;
00631 Molecule *molecule = Node::Object()->molecule;
00632 BigReal dt = dt_fs * 0.001; // convert to ps
00633 int i;
00634
00635 if (simParams->drudeOn) {
00636 for (i = 0; i < numAtoms; i++) {
00637
00638 if (i < numAtoms-1 &&
00639 a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
00640 //printf("*** Found Drude particle %d\n", a[i+1].id);
00641 // i+1 is a Drude particle with parent i
00642
00643 // convert from Cartesian coordinates to (COM,bond) coordinates
00644 BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
00645 Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
00646 Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
00647 BigReal dt_gamma;
00648
00649 // use Langevin damping factor i for v_com
00650 dt_gamma = dt * molecule->langevin_param(a[i].id);
00651 if (dt_gamma != 0.0) {
00652 v_com *= ( 1. - 0.5 * dt_gamma );
00653 }
00654
00655 // use Langevin damping factor i+1 for v_bnd
00656 dt_gamma = dt * molecule->langevin_param(a[i+1].id);
00657 if (dt_gamma != 0.0) {
00658 v_bnd *= ( 1. - 0.5 * dt_gamma );
00659 }
00660
00661 // convert back
00662 a[i].velocity = v_com - m * v_bnd;
00663 a[i+1].velocity = v_bnd + a[i].velocity;
00664
00665 i++; // +1 from loop, we've updated both particles
00666 }
00667 else {
00668 int aid = a[i].id;
00669 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00670 if ( ! dt_gamma ) continue;
00671
00672 a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00673 }
00674
00675 } // end for
00676 } // end if drudeOn
00677 else {
00678
00679 for ( i = 0; i < numAtoms; ++i )
00680 {
00681 int aid = a[i].id;
00682 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00683 if ( ! dt_gamma ) continue;
00684
00685 a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00686 }
00687
00688 } // end else
00689
00690 } // end if langevinOn
00691 }
|
|
|
Definition at line 694 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, SimParameters::drudeOn, SimParameters::drudeTemp, Random::gaussian_vector(), CompAtomExt::id, Molecule::langevin_param(), SimParameters::langevinOn, 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(). 00695 {
00696 if ( simParams->langevinOn )
00697 {
00698 rattle1(dt_fs,1); // conserve momentum if gammas differ
00699
00700 FullAtom *a = patch->atom.begin();
00701 int numAtoms = patch->numAtoms;
00702 Molecule *molecule = Node::Object()->molecule;
00703 BigReal dt = dt_fs * 0.001; // convert to ps
00704 BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00705
00706 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00707 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00708 int i;
00709
00710 if (simParams->drudeOn) {
00711 BigReal kbT_bnd = BOLTZMAN*(simParams->drudeTemp); // drude bond Temp
00712
00713 for (i = 0; i < numAtoms; i++) {
00714
00715 if (i < numAtoms-1 &&
00716 a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
00717 //printf("*** Found Drude particle %d\n", a[i+1].id);
00718 // i+1 is a Drude particle with parent i
00719
00720 // convert from Cartesian coordinates to (COM,bond) coordinates
00721 BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
00722 Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
00723 Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
00724 BigReal dt_gamma;
00725
00726 // use Langevin damping factor i for v_com
00727 dt_gamma = dt * molecule->langevin_param(a[i].id);
00728 if (dt_gamma != 0.0) {
00729 BigReal mass = a[i].mass + a[i+1].mass;
00730 v_com += random->gaussian_vector() *
00731 sqrt( 2 * dt_gamma * kbT *
00732 ( a[i].partition ? tempFactor : 1.0 ) / mass );
00733 v_com /= ( 1. + 0.5 * dt_gamma );
00734 }
00735
00736 // use Langevin damping factor i+1 for v_bnd
00737 dt_gamma = dt * molecule->langevin_param(a[i+1].id);
00738 if (dt_gamma != 0.0) {
00739 BigReal mass = a[i+1].mass * (1. - m);
00740 v_bnd += random->gaussian_vector() *
00741 sqrt( 2 * dt_gamma * kbT_bnd *
00742 ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
00743 v_bnd /= ( 1. + 0.5 * dt_gamma );
00744 }
00745
00746 // convert back
00747 a[i].velocity = v_com - m * v_bnd;
00748 a[i+1].velocity = v_bnd + a[i].velocity;
00749
00750 i++; // +1 from loop, we've updated both particles
00751 }
00752 else {
00753 int aid = a[i].id;
00754 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00755 if ( ! dt_gamma ) continue;
00756
00757 a[i].velocity += random->gaussian_vector() *
00758 sqrt( 2 * dt_gamma * kbT *
00759 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00760 a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00761 }
00762
00763 } // end for
00764 } // end if drudeOn
00765 else {
00766
00767 for ( i = 0; i < numAtoms; ++i )
00768 {
00769 int aid = a[i].id;
00770 BigReal dt_gamma = dt * molecule->langevin_param(aid);
00771 if ( ! dt_gamma ) continue;
00772
00773 a[i].velocity += random->gaussian_vector() *
00774 sqrt( 2 * dt_gamma * kbT *
00775 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00776 a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00777 }
00778
00779 } // end else
00780
00781 } // end if langevinOn
00782 }
|
|
|
Definition at line 1093 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, PDBVELFACTOR, simParams, terminate(), and FullAtom::velocity. Referenced by integrate(). 01094 {
01095 FullAtom *a = patch->atom.begin();
01096 int numAtoms = patch->numAtoms;
01097 if ( simParams->maximumMove ) {
01098 const BigReal dt = timestep / TIMEFACTOR;
01099 const BigReal maxvel = simParams->maximumMove / dt;
01100 const BigReal maxvel2 = maxvel * maxvel;
01101 for ( int i=0; i<numAtoms; ++i ) {
01102 if ( a[i].velocity.length2() > maxvel2 ) {
01103 a[i].velocity *= ( maxvel / a[i].velocity.length() );
01104 }
01105 }
01106 } else {
01107 const BigReal dt = timestep / TIMEFACTOR;
01108 const BigReal maxvel = simParams->cutoff / dt;
01109 const BigReal maxvel2 = maxvel * maxvel;
01110 int killme = 0;
01111 for ( int i=0; i<numAtoms; ++i ) {
01112 killme = killme || ( a[i].velocity.length2() > maxvel2 );
01113 }
01114 if ( killme ) {
01115 for ( int i=0; i<numAtoms; ++i ) {
01116 if ( a[i].velocity.length2() > maxvel2 ) {
01117 iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
01118 << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
01119 << ( PDBVELFACTOR * maxvel ) << ")\n" << endi;
01120 }
01121 }
01122 iout << iERROR <<
01123 "Atoms moving too fast; simulation has become unstable.\n" << endi;
01124 Node::Object()->enableEarlyExit();
01125 terminate();
01126 }
01127 }
01128 }
|
|
|
Definition at line 1130 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::minimizeOn, Patch::numAtoms, patch, simParams, and FullAtom::velocity. Referenced by integrate(). 01131 {
01132 if ( simParams->minimizeOn ) {
01133 FullAtom *a = patch->atom.begin();
01134 int numAtoms = patch->numAtoms;
01135 for ( int i=0; i<numAtoms; ++i ) {
01136 a[i].velocity = 0.;
01137 }
01138 }
01139 }
|
|
|
Definition at line 437 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, minimizeMoveDownhill(), SimParameters::mollyOn, SimParameters::N, newMinimizeDirection(), newMinimizePosition(), patch, SimParameters::PMEOn, quenchVelocities(), rebalanceLoad(), runComputeObjects(), simParams, Flags::step, submitCollections(), and submitMinimizeReductions(). Referenced by algorithm(). 00437 {
00438 const int numberOfSteps = simParams->N;
00439 int &step = patch->flags.step;
00440 step = simParams->firstTimestep;
00441
00442 int &maxForceUsed = patch->flags.maxForceUsed;
00443 int &maxForceMerged = patch->flags.maxForceMerged;
00444 maxForceUsed = Results::normal;
00445 maxForceMerged = Results::normal;
00446 int &doNonbonded = patch->flags.doNonbonded;
00447 doNonbonded = 1;
00448 maxForceUsed = Results::nbond;
00449 maxForceMerged = Results::nbond;
00450 const int dofull = ( simParams->fullDirectOn ||
00451 simParams->FMAOn || simParams->PMEOn );
00452 int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00453 doFullElectrostatics = dofull;
00454 if ( dofull ) {
00455 maxForceMerged = Results::slow;
00456 maxForceUsed = Results::slow;
00457 }
00458 int &doMolly = patch->flags.doMolly;
00459 doMolly = simParams->mollyOn && doFullElectrostatics;
00460 int &doEnergy = patch->flags.doEnergy;
00461 doEnergy = 1;
00462
00463 runComputeObjects(1,0); // must migrate here!
00464
00465 submitMinimizeReductions(step);
00466 rebalanceLoad(step);
00467
00468 int downhill = 1; // start out just fixing bad contacts
00469 int minSeq = 0;
00470 for ( ++step; step <= numberOfSteps; ++step ) {
00471 BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00472 if ( downhill ) {
00473 if ( c ) minimizeMoveDownhill();
00474 else downhill = 0;
00475 }
00476 if ( ! downhill ) {
00477 if ( ! c ) { // new direction
00478 c = broadcast->minimizeCoefficient.get(minSeq++);
00479 newMinimizeDirection(c); // v = c * v + f
00480 c = broadcast->minimizeCoefficient.get(minSeq++);
00481 } // same direction
00482 newMinimizePosition(c); // x = x + c * v
00483 }
00484
00485 runComputeObjects(1,0);
00486 submitMinimizeReductions(step);
00487 submitCollections(step, 1); // write out zeros for velocities
00488 rebalanceLoad(step);
00489 }
00490 quenchVelocities(); // zero out bogus velocity
00491 }
|
|
|
Definition at line 498 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), Patch::f, SimParameters::fixedAtomsOn, Force, Vector::length2(), Patch::numAtoms, patch, CompAtom::position, simParams, and Vector::unit(). Referenced by minimize(). 00498 {
00499
00500 FullAtom *a = patch->atom.begin();
00501 Force *f1 = patch->f[Results::normal].begin();
00502 Force *f2 = patch->f[Results::nbond].begin();
00503 Force *f3 = patch->f[Results::slow].begin();
00504 int numAtoms = patch->numAtoms;
00505
00506 for ( int i = 0; i < numAtoms; ++i ) {
00507 if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
00508 Force f = f1[i] + f2[i] + f3[i];
00509 if ( f.length2() > FMAX2 ) {
00510 a[i].position += ( 0.1 * f.unit() );
00511 }
00512 }
00513 }
|
|
|
Definition at line 516 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), Patch::f, SimParameters::fixedAtomsOn, Force, Patch::numAtoms, patch, simParams, and FullAtom::velocity. Referenced by minimize(). 00516 {
00517 FullAtom *a = patch->atom.begin();
00518 Force *f1 = patch->f[Results::normal].begin();
00519 Force *f2 = patch->f[Results::nbond].begin();
00520 Force *f3 = patch->f[Results::slow].begin();
00521 int numAtoms = patch->numAtoms;
00522
00523 for ( int i = 0; i < numAtoms; ++i ) {
00524 a[i].velocity *= c;
00525 a[i].velocity += f1[i] + f2[i] + f3[i];
00526 if ( simParams->fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00527 }
00528 }
|
|
|
Definition at line 531 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, patch, CompAtom::position, and FullAtom::velocity. Referenced by minimize(). 00531 {
00532 FullAtom *a = patch->atom.begin();
00533 int numAtoms = patch->numAtoms;
00534
00535 for ( int i = 0; i < numAtoms; ++i ) {
00536 a[i].position += c * a[i].velocity;
00537 }
00538 }
|
|
|
Definition at line 541 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity. Referenced by minimize(). 00541 {
00542 FullAtom *a = patch->atom.begin();
00543 int numAtoms = patch->numAtoms;
00544
00545 for ( int i = 0; i < numAtoms; ++i ) {
00546 a[i].velocity = 0;
00547 }
00548 }
|
|
||||||||||||
|
Definition at line 1064 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(). 01065 {
01066 if ( simParams->rigidBonds != RIGID_NONE ) {
01067 Tensor virial;
01068 Tensor *vp = ( pressure ? &virial : 0 );
01069 if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
01070 iout << iERROR <<
01071 "Constraint failure; simulation has become unstable.\n" << endi;
01072 Node::Object()->enableEarlyExit();
01073 terminate();
01074 }
01075 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01076 }
01077 }
|
|
||||||||||||
|
Definition at line 1079 of file Sequencer.C. References ADD_TENSOR_OBJECT, patch, HomePatch::rattle2(), reduction, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_VIRIAL_NORMAL, SimParameters::rigidBonds, and simParams. 01080 {
01081 if ( simParams->rigidBonds != RIGID_NONE ) {
01082 Tensor virial;
01083 patch->rattle2(dt, &virial);
01084 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01085 // we need to add to alt and int virial because not included in forces
01086 #ifdef ALTVIRIAL
01087 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01088 #endif
01089 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,virial);
01090 }
01091 }
|
|
||||||||||||
|
Definition at line 939 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, 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(). 00940 {
00941 const int reassignFreq = simParams->reassignFreq;
00942 if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
00943 FullAtom *a = patch->atom.begin();
00944 int numAtoms = patch->numAtoms;
00945 BigReal newTemp = simParams->reassignTemp;
00946 newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
00947 if ( simParams->reassignIncr > 0.0 ) {
00948 if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
00949 newTemp = simParams->reassignHold;
00950 } else {
00951 if ( newTemp < simParams->reassignHold )
00952 newTemp = simParams->reassignHold;
00953 }
00954 BigReal kbT = BOLTZMAN * newTemp;
00955
00956 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00957 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00958
00959 for ( int i = 0; i < numAtoms; ++i )
00960 {
00961 a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
00962 sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
00963 * random->gaussian_vector() );
00964 }
00965 } else {
00966 NAMD_bug("Sequencer::reassignVelocities called improperly!");
00967 }
00968 }
|
|
|
Definition at line 1684 of file Sequencer.C. References LdbCoordinator::getNumStepsToRun(), Patch::getPatchID(), ldbSteps, LdbCoordinator::Object(), pairlistsAreValid, patch, LdbCoordinator::rebalance(), and HomePatch::submitLoadStats(). Referenced by integrate(), and minimize(). 01684 {
01685 if ( ! ldbSteps ) {
01686 ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
01687 }
01688 if ( ! --ldbSteps ) {
01689 patch->submitLoadStats(timestep);
01690 ldbCoordinator->rebalance(this,patch->getPatchID());
01691 pairlistsAreValid = 0;
01692 }
01693 }
|
|
||||||||||||
|
Definition at line 1040 of file Sequencer.C. References ADD_TENSOR_OBJECT, patch, HomePatch::redistrib_swm4_forces(), reduction, and REDUCTION_VIRIAL_NORMAL. Referenced by integrate(). 01040 {
01041 Tensor virial;
01042 Tensor *vp = ( pressure ? &virial : 0 );
01043 patch->redistrib_swm4_forces(ftag, vp);
01044 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01045 }
|
|
||||||||||||
|
Definition at line 1033 of file Sequencer.C. References ADD_TENSOR_OBJECT, patch, HomePatch::redistrib_tip4p_forces(), reduction, and REDUCTION_VIRIAL_NORMAL. Referenced by integrate(). 01033 {
01034 Tensor virial;
01035 Tensor *vp = ( pressure ? &virial : 0 );
01036 patch->redistrib_tip4p_forces(ftag, vp);
01037 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01038 }
|
|
|
Definition at line 970 of file Sequencer.C. References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, SimParameters::fixedAtomsOn, Random::gaussian_vector(), SimParameters::initialTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, FullAtom::mass, Patch::numAtoms, patch, random, simParams, and FullAtom::velocity. Referenced by algorithm(). 00971 {
00972 FullAtom *a = patch->atom.begin();
00973 int numAtoms = patch->numAtoms;
00974 BigReal newTemp = simParams->initialTemp;
00975 BigReal kbT = BOLTZMAN * newTemp;
00976
00977 int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00978 BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00979
00980 for ( int i = 0; i < numAtoms; ++i )
00981 {
00982 a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
00983 sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
00984 * random->gaussian_vector() );
00985 }
00986 }
|
|
|
Definition at line 998 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(). 00999 {
01000 FullAtom *a = patch->atom.begin();
01001 int numAtoms = patch->numAtoms;
01002 Molecule *molecule = Node::Object()->molecule;
01003 for ( int i = 0; i < numAtoms; ++i )
01004 {
01005 a[i].charge = molecule->atomcharge(a[i].id);
01006 }
01007 }
|
|
|
Definition at line 921 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(). 00922 {
00923 const int rescaleFreq = simParams->rescaleFreq;
00924 if ( rescaleFreq > 0 ) {
00925 FullAtom *a = patch->atom.begin();
00926 int numAtoms = patch->numAtoms;
00927 ++rescaleVelocities_numTemps;
00928 if ( rescaleVelocities_numTemps == rescaleFreq ) {
00929 BigReal factor = broadcast->velocityRescaleFactor.get(step);
00930 for ( int i = 0; i < numAtoms; ++i )
00931 {
00932 a[i].velocity *= factor;
00933 }
00934 rescaleVelocities_numTemps = 0;
00935 }
00936 }
00937 }
|
|
|
Definition at line 988 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity. Referenced by algorithm(). 00989 {
00990 FullAtom *a = patch->atom.begin();
00991 int numAtoms = patch->numAtoms;
00992 for ( int i = 0; i < numAtoms; ++i )
00993 {
00994 a[i].velocity *= factor;
00995 }
00996 }
|
|
|
Definition at line 84 of file Sequencer.C. References awaken(), DebugM, and SEQ_STK_SZ. Referenced by HomePatch::runSequencer(). 00085 {
00086 // create a Thread and invoke it
00087 DebugM(4, "::run() - this = " << this << "\n" );
00088 thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
00089 CthSetStrategyDefault(thread);
00090 awaken();
00091 }
|
|
||||||||||||
|
Definition at line 1631 of file Sequencer.C. References ADD_TENSOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), Flags::doMolly, Flags::doNonbonded, Patch::f, Patch::flags, HomePatch::mollyMollify(), Patch::numAtoms, pairlistsAge, pairlistsAreValid, SimParameters::pairlistsPerCycle, patch, HomePatch::positionsReady(), reduction, REDUCTION_VIRIAL_SLOW, Flags::savePairlists, simParams, SimParameters::stepsPerCycle, suspend(), Flags::usePairlists, and SimParameters::usePairlists. Referenced by integrate(), and minimize(). 01632 {
01633 if ( migration ) pairlistsAreValid = 0;
01634 if ( pairlistsAreValid && ( pairlistsAge > (
01635 (simParams->stepsPerCycle - 1) / simParams->pairlistsPerCycle ) ) ) {
01636 pairlistsAreValid = 0;
01637 }
01638 if ( ! simParams->usePairlists ) pairlists = 0;
01639 patch->flags.usePairlists = pairlists || pairlistsAreValid;
01640 patch->flags.savePairlists =
01641 pairlists && ! pairlistsAreValid;
01642 patch->positionsReady(migration);
01643 suspend(); // until all deposit boxes close
01644 if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
01645 pairlistsAreValid = 1;
01646 pairlistsAge = 0;
01647 }
01648 if ( pairlistsAreValid ) ++pairlistsAge;
01649 if ( patch->flags.doMolly ) {
01650 Tensor virial;
01651 patch->mollyMollify(&virial);
01652 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
01653 }
01654 #ifdef NAMD_CUDA_XXX
01655 int numAtoms = patch->numAtoms;
01656 FullAtom *a = patch->atom.begin();
01657 for ( int i=0; i<numAtoms; ++i ) {
01658 CkPrintf("%d %g %g %g\n", a[i].id,
01659 patch->f[Results::normal][i].x +
01660 patch->f[Results::nbond][i].x +
01661 patch->f[Results::slow][i].x,
01662 patch->f[Results::normal][i].y +
01663 patch->f[Results::nbond][i].y +
01664 patch->f[Results::slow][i].y,
01665 patch->f[Results::normal][i].z +
01666 patch->f[Results::nbond][i].z +
01667 patch->f[Results::slow][i].z);
01668 CkPrintf("%d %g %g %g\n", a[i].id,
01669 patch->f[Results::normal][i].x,
01670 patch->f[Results::nbond][i].x,
01671 patch->f[Results::slow][i].x);
01672 CkPrintf("%d %g %g %g\n", a[i].id,
01673 patch->f[Results::normal][i].y,
01674 patch->f[Results::nbond][i].y,
01675 patch->f[Results::slow][i].y);
01676 CkPrintf("%d %g %g %g\n", a[i].id,
01677 patch->f[Results::normal][i].z,
01678 patch->f[Results::nbond][i].z,
01679 patch->f[Results::slow][i].z);
01680 }
01681 #endif
01682 }
|
|
|
Definition at line 1028 of file Sequencer.C. References patch, and HomePatch::saveForce(). Referenced by integrate().
|
|
||||||||||||
|
Definition at line 1622 of file Sequencer.C. References HomePatch::atom, collection, Output::coordinateNeeded(), Patch::lattice, patch, CollectionMgr::submitPositions(), CollectionMgr::submitVelocities(), and Output::velocityNeeded(). Referenced by algorithm(), integrate(), and minimize(). 01623 {
01624 int prec = Output::coordinateNeeded(step);
01625 if ( prec )
01626 collection->submitPositions(step,patch->atom,patch->lattice,prec);
01627 if ( Output::velocityNeeded(step) )
01628 collection->submitVelocities(step,zeroVel,patch->atom);
01629 }
|
|
|
Definition at line 1141 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(). 01142 {
01143 // velocity-dependent quantities *** ONLY ***
01144 // positions are not at half-step when called
01145 FullAtom *a = patch->atom.begin();
01146 int numAtoms = patch->numAtoms;
01147
01148 #if CMK_BLUEGENEL
01149 CmiNetworkProgressAfter (0);
01150 #endif
01151
01152 {
01153 BigReal kineticEnergy = 0;
01154 Tensor virial;
01155 if ( simParams->pairInteractionOn ) {
01156 if ( simParams->pairInteractionSelf ) {
01157 for ( int i = 0; i < numAtoms; ++i ) {
01158 if ( a[i].partition != 1 ) continue;
01159 kineticEnergy += a[i].mass * a[i].velocity.length2();
01160 virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01161 }
01162 }
01163 } else {
01164 for ( int i = 0; i < numAtoms; ++i ) {
01165 if (a[i].mass < 0.01) continue;
01166 kineticEnergy += a[i].mass * a[i].velocity.length2();
01167 virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01168 }
01169 }
01170
01171 kineticEnergy *= 0.5 * 0.5;
01172 reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
01173 virial *= 0.5;
01174 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01175 #ifdef ALTVIRIAL
01176 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01177 #endif
01178 }
01179
01180 if (pressureProfileReduction) {
01181 int nslabs = simParams->pressureProfileSlabs;
01182 const Lattice &lattice = patch->lattice;
01183 BigReal idz = nslabs/lattice.c().z;
01184 BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01185 int useGroupPressure = simParams->useGroupPressure;
01186
01187 // Compute kinetic energy partition, possibly subtracting off
01188 // internal kinetic energy if group pressure is enabled.
01189 // Since the regular pressure is 1/2 mvv and the internal kinetic
01190 // term that is subtracted off for the group pressure is
01191 // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
01192 // 1/2 m * v * v_cm. The factor of 1/2 is because submitHalfstep
01193 // gets called twice per timestep.
01194 int hgs;
01195 for (int i=0; i<numAtoms; i += hgs) {
01196 int j, ppoffset;
01197 hgs = a[i].hydrogenGroupSize;
01198 int partition = a[i].partition;
01199
01200 BigReal m_cm = 0;
01201 Velocity v_cm(0,0,0);
01202 for (j=i; j< i+hgs; ++j) {
01203 m_cm += a[j].mass;
01204 v_cm += a[j].mass * a[j].velocity;
01205 }
01206 v_cm /= m_cm;
01207 for (j=i; j < i+hgs; ++j) {
01208 BigReal mass = a[j].mass;
01209 if (! (useGroupPressure && j != i)) {
01210 BigReal z = a[j].position.z;
01211 int slab = (int)floor((z-zmin)*idz);
01212 if (slab < 0) slab += nslabs;
01213 else if (slab >= nslabs) slab -= nslabs;
01214 ppoffset = 3*(slab + partition*nslabs);
01215 }
01216 BigReal wxx, wyy, wzz;
01217 if (useGroupPressure) {
01218 wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
01219 wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
01220 wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
01221 } else {
01222 wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
01223 wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
01224 wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
01225 }
01226 pressureProfileReduction->item(ppoffset ) += wxx;
01227 pressureProfileReduction->item(ppoffset+1) += wyy;
01228 pressureProfileReduction->item(ppoffset+2) += wzz;
01229 }
01230 }
01231 }
01232
01233 {
01234 BigReal intKineticEnergy = 0;
01235 Tensor intVirialNormal;
01236
01237 int hgs;
01238 for ( int i = 0; i < numAtoms; i += hgs ) {
01239
01240 #if CMK_BLUEGENEL
01241 CmiNetworkProgress ();
01242 #endif
01243
01244 hgs = a[i].hydrogenGroupSize;
01245 int j;
01246 BigReal m_cm = 0;
01247 Velocity v_cm(0,0,0);
01248 for ( j = i; j < (i+hgs); ++j ) {
01249 m_cm += a[j].mass;
01250 v_cm += a[j].mass * a[j].velocity;
01251 }
01252 v_cm /= m_cm;
01253 if ( simParams->pairInteractionOn ) {
01254 if ( simParams->pairInteractionSelf ) {
01255 for ( j = i; j < (i+hgs); ++j ) {
01256 if ( a[j].partition != 1 ) continue;
01257 BigReal mass = a[j].mass;
01258 Vector v = a[j].velocity;
01259 Vector dv = v - v_cm;
01260 intKineticEnergy += mass * (v * dv);
01261 intVirialNormal.outerAdd (mass, v, dv);
01262 }
01263 }
01264 } else {
01265 for ( j = i; j < (i+hgs); ++j ) {
01266 BigReal mass = a[j].mass;
01267 Vector v = a[j].velocity;
01268 Vector dv = v - v_cm;
01269 intKineticEnergy += mass * (v * dv);
01270 intVirialNormal.outerAdd(mass, v, dv);
01271 }
01272 }
01273 }
01274
01275 intKineticEnergy *= 0.5 * 0.5;
01276 reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
01277 intVirialNormal *= 0.5;
01278 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01279 }
01280
01281 }
|
|
|
Definition at line 1521 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, FMAX2, Force, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, FullAtom::mass, Patch::numAtoms, Tensor::outerAdd(), patch, CompAtom::position, Position, 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, simParams, SubmitReduction::submit(), and FullAtom::velocity. Referenced by minimize(). 01522 {
01523 FullAtom *a = patch->atom.begin();
01524 Force *f1 = patch->f[Results::normal].begin();
01525 Force *f2 = patch->f[Results::nbond].begin();
01526 Force *f3 = patch->f[Results::slow].begin();
01527 int numAtoms = patch->numAtoms;
01528
01529 reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01530
01531 BigReal fdotf = 0;
01532 BigReal fdotv = 0;
01533 BigReal vdotv = 0;
01534 int numHuge = 0;
01535 for ( int i = 0; i < numAtoms; ++i ) {
01536 if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
01537 Force f = f1[i] + f2[i] + f3[i];
01538 BigReal ff = f * f;
01539 if ( ff > FMAX2 ) {
01540 ++numHuge;
01541 // pad scaling so minimizeMoveDownhill() doesn't miss them
01542 BigReal fmult = 1.01 * sqrt(FMAX2/ff);
01543 f *= fmult; ff = f * f;
01544 f1[i] *= fmult;
01545 f2[i] *= fmult;
01546 f3[i] *= fmult;
01547 }
01548 fdotf += ff;
01549 fdotv += f * a[i].velocity;
01550 vdotv += a[i].velocity * a[i].velocity;
01551 }
01552
01553 reduction->item(REDUCTION_MIN_F_DOT_F) += fdotf;
01554 reduction->item(REDUCTION_MIN_F_DOT_V) += fdotv;
01555 reduction->item(REDUCTION_MIN_V_DOT_V) += vdotv;
01556 reduction->item(REDUCTION_MIN_HUGE_COUNT) += numHuge;
01557
01558 {
01559 Tensor intVirialNormal;
01560 Tensor intVirialNbond;
01561 Tensor intVirialSlow;
01562
01563 int hgs;
01564 for ( int i = 0; i < numAtoms; i += hgs ) {
01565 hgs = a[i].hydrogenGroupSize;
01566 int j;
01567 BigReal m_cm = 0;
01568 Position x_cm(0,0,0);
01569 for ( j = i; j < (i+hgs); ++j ) {
01570 m_cm += a[j].mass;
01571 x_cm += a[j].mass * a[j].position;
01572 }
01573 x_cm /= m_cm;
01574 for ( j = i; j < (i+hgs); ++j ) {
01575 BigReal mass = a[j].mass;
01576 // net force treated as zero for fixed atoms
01577 if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01578 Vector dx = a[j].position - x_cm;
01579 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01580 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01581 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01582 }
01583 }
01584
01585 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01586 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
01587 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
01588 }
01589
01590 if ( simParams->fixedAtomsOn ) {
01591 Tensor fixVirialNormal;
01592 Tensor fixVirialNbond;
01593 Tensor fixVirialSlow;
01594 Vector fixForceNormal = 0;
01595 Vector fixForceNbond = 0;
01596 Vector fixForceSlow = 0;
01597
01598 for ( int j = 0; j < numAtoms; j++ ) {
01599 if ( a[j].atomFixed ) {
01600 Vector dx = a[j].fixedPosition;
01601 // all negative because fixed atoms cancels these forces
01602 fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j],dx);
01603 fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j],dx);
01604 fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j],dx);
01605 fixForceNormal -= patch->f[Results::normal][j];
01606 fixForceNbond -= patch->f[Results::nbond][j];
01607 fixForceSlow -= patch->f[Results::slow][j];
01608 }
01609 }
01610
01611 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
01612 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
01613 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
01614 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
01615 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
01616 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
01617 }
01618
01619 reduction->submit();
01620 }
|
|
|
Definition at line 550 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(). 00550 {
00551
00552 FullAtom *a = patch->atom.begin();
00553 const int numAtoms = patch->numAtoms;
00554
00555 Vector momentum = 0;
00556 BigReal mass = 0;
00557 if ( simParams->zeroMomentumAlt ) {
00558 for ( int i = 0; i < numAtoms; ++i ) {
00559 momentum += a[i].mass * a[i].velocity;
00560 mass += 1.;
00561 }
00562 } else {
00563 for ( int i = 0; i < numAtoms; ++i ) {
00564 momentum += a[i].mass * a[i].velocity;
00565 mass += a[i].mass;
00566 }
00567 }
00568
00569 ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
00570 reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
00571 }
|
|
|
Definition at line 1283 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(). 01284 {
01285 FullAtom *a = patch->atom.begin();
01286 int numAtoms = patch->numAtoms;
01287
01288 #if CMK_BLUEGENEL
01289 CmiNetworkProgressAfter(0);
01290 #endif
01291
01292 reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01293 reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
01294
01295 {
01296 BigReal kineticEnergy = 0;
01297 Vector momentum = 0;
01298 Vector angularMomentum = 0;
01299 Vector o = patch->lattice.origin();
01300 int i;
01301 if ( simParams->pairInteractionOn ) {
01302 if ( simParams->pairInteractionSelf ) {
01303 for (i = 0; i < numAtoms; ++i ) {
01304 if ( a[i].partition != 1 ) continue;
01305 kineticEnergy += a[i].mass * a[i].velocity.length2();
01306 momentum += a[i].mass * a[i].velocity;
01307 angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01308 }
01309 }
01310 } else {
01311 for (i = 0; i < numAtoms; ++i ) {
01312 kineticEnergy += a[i].mass * a[i].velocity.length2();
01313 momentum += a[i].mass * a[i].velocity;
01314 angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01315 }
01316 if (simParams->drudeOn) {
01317 BigReal drudeComKE = 0.;
01318 BigReal drudeBondKE = 0.;
01319
01320 for (i = 0; i < numAtoms; i++) {
01321 if (i < numAtoms-1 &&
01322 a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01323 // i+1 is a Drude particle with parent i
01324
01325 // convert from Cartesian coordinates to (COM,bond) coordinates
01326 BigReal m_com = (a[i].mass + a[i+1].mass); // mass of COM
01327 BigReal m = a[i+1].mass / m_com; // mass ratio
01328 BigReal m_bond = a[i+1].mass * (1. - m); // mass of bond
01329 Vector v_bond = a[i+1].velocity - a[i].velocity; // vel of bond
01330 Vector v_com = a[i].velocity + m * v_bond; // vel of COM
01331
01332 drudeComKE += m_com * v_com.length2();
01333 drudeBondKE += m_bond * v_bond.length2();
01334
01335 i++; // +1 from loop, we've updated both particles
01336 }
01337 else {
01338 drudeComKE += a[i].mass * a[i].velocity.length2();
01339 }
01340 } // end for
01341
01342 drudeComKE *= 0.5;
01343 drudeBondKE *= 0.5;
01344 reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY)
01345 += drudeComKE;
01346 reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY)
01347 += drudeBondKE;
01348 } // end drudeOn
01349
01350 } // end else
01351
01352 kineticEnergy *= 0.5;
01353 reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
01354 ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
01355 ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
01356 }
01357
01358 #ifdef ALTVIRIAL
01359 // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
01360 {
01361 Tensor altVirial;
01362 for ( int i = 0; i < numAtoms; ++i ) {
01363 altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
01364 }
01365 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
01366 }
01367 {
01368 Tensor altVirial;
01369 for ( int i = 0; i < numAtoms; ++i ) {
01370 altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
01371 }
01372 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
01373 }
01374 {
01375 Tensor altVirial;
01376 for ( int i = 0; i < numAtoms; ++i ) {
01377 altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
01378 }
01379 ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
01380 }
01381 #endif
01382
01383 {
01384 BigReal intKineticEnergy = 0;
01385 Tensor intVirialNormal;
01386 Tensor intVirialNbond;
01387 Tensor intVirialSlow;
01388
01389 int hgs;
01390 for ( int i = 0; i < numAtoms; i += hgs ) {
01391 #if CMK_BLUEGENEL
01392 CmiNetworkProgress();
01393 #endif
01394 hgs = a[i].hydrogenGroupSize;
01395 int j;
01396 BigReal m_cm = 0;
01397 Position x_cm(0,0,0);
01398 Velocity v_cm(0,0,0);
01399 for ( j = i; j < (i+hgs); ++j ) {
01400 m_cm += a[j].mass;
01401 x_cm += a[j].mass * a[j].position;
01402 v_cm += a[j].mass * a[j].velocity;
01403 }
01404 x_cm /= m_cm;
01405 v_cm /= m_cm;
01406 int fixedAtomsOn = simParams->fixedAtomsOn;
01407 if ( simParams->pairInteractionOn ) {
01408 int pairInteractionSelf = simParams->pairInteractionSelf;
01409 for ( j = i; j < (i+hgs); ++j ) {
01410 if ( a[j].partition != 1 &&
01411 ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
01412 // net force treated as zero for fixed atoms
01413 if ( fixedAtomsOn && a[j].atomFixed ) continue;
01414 BigReal mass = a[j].mass;
01415 Vector v = a[j].velocity;
01416 Vector dv = v - v_cm;
01417 intKineticEnergy += mass * (v * dv);
01418 Vector dx = a[j].position - x_cm;
01419 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01420 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01421 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01422 }
01423 } else {
01424 for ( j = i; j < (i+hgs); ++j ) {
01425 // net force treated as zero for fixed atoms
01426 if ( fixedAtomsOn && a[j].atomFixed ) continue;
01427 BigReal mass = a[j].mass;
01428 Vector v = a[j].velocity;
01429 Vector dv = v - v_cm;
01430 intKineticEnergy += mass * (v * dv);
01431 Vector dx = a[j].position - x_cm;
01432 intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01433 intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01434 intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01435 }
01436 }
01437 }
01438
01439 intKineticEnergy *= 0.5;
01440 reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
01441 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01442 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
01443 ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
01444 }
01445
01446 if (pressureProfileReduction && simParams->useGroupPressure) {
01447 // subtract off internal virial term, calculated as for intVirial.
01448 int nslabs = simParams->pressureProfileSlabs;
01449 const Lattice &lattice = patch->lattice;
01450 BigReal idz = nslabs/lattice.c().z;
01451 BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01452 int useGroupPressure = simParams->useGroupPressure;
01453
01454 int hgs;
01455 for (int i=0; i<numAtoms; i += hgs) {
01456 int j;
01457 hgs = a[i].hydrogenGroupSize;
01458 BigReal m_cm = 0;
01459 Position x_cm(0,0,0);
01460 for (j=i; j< i+hgs; ++j) {
01461 m_cm += a[j].mass;
01462 x_cm += a[j].mass * a[j].position;
01463 }
01464 x_cm /= m_cm;
01465
01466 BigReal z = a[i].position.z;
01467 int slab = (int)floor((z-zmin)*idz);
01468 if (slab < 0) slab += nslabs;
01469 else if (slab >= nslabs) slab -= nslabs;
01470 int partition = a[i].partition;
01471 int ppoffset = 3*(slab + nslabs*partition);
01472 for (j=i; j < i+hgs; ++j) {
01473 BigReal mass = a[j].mass;
01474 Vector dx = a[j].position - x_cm;
01475 const Vector &fnormal = patch->f[Results::normal][j];
01476 const Vector &fnbond = patch->f[Results::nbond][j];
01477 const Vector &fslow = patch->f[Results::slow][j];
01478 BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
01479 BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
01480 BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
01481 pressureProfileReduction->item(ppoffset ) -= wxx;
01482 pressureProfileReduction->item(ppoffset+1) -= wyy;
01483 pressureProfileReduction->item(ppoffset+2) -= wzz;
01484 }
01485 }
01486 }
01487
01488 if ( simParams->fixedAtomsOn ) {
01489 Tensor fixVirialNormal;
01490 Tensor fixVirialNbond;
01491 Tensor fixVirialSlow;
01492 Vector fixForceNormal = 0;
01493 Vector fixForceNbond = 0;
01494 Vector fixForceSlow = 0;
01495
01496 for ( int j = 0; j < numAtoms; j++ ) {
01497 if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01498 Vector dx = a[j].fixedPosition;
01499 // all negative because fixed atoms cancels these forces
01500 fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
01501 fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
01502 fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
01503 fixForceNormal -= patch->f[Results::normal][j];
01504 fixForceNbond -= patch->f[Results::nbond][j];
01505 fixForceSlow -= patch->f[Results::slow][j];
01506 }
01507 }
01508
01509 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
01510 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
01511 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
01512 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
01513 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
01514 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
01515 }
01516
01517 reduction->submit();
01518 if (pressureProfileReduction) pressureProfileReduction->submit();
01519 }
|
|
|
Definition at line 29 of file Sequencer.h. Referenced by HomePatch::doAtomMigration(), LdbCoordinator::rebalance(), and runComputeObjects(). 00029 { CthSuspend(); }
|
|
||||||||||||
|
Definition at line 1009 of file Sequencer.C. References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, broadcast, SimpleBroadcastObject< T >::get(), CompAtomExt::id, Molecule::langevin_param(), Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, and FullAtom::velocity. Referenced by integrate(). 01010 {
01011 if ( simParams->tCoupleOn )
01012 {
01013 FullAtom *a = patch->atom.begin();
01014 int numAtoms = patch->numAtoms;
01015 BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
01016 Molecule *molecule = Node::Object()->molecule;
01017 BigReal dt = dt_fs * 0.001; // convert to ps
01018 coefficient *= dt;
01019 for ( int i = 0; i < numAtoms; ++i )
01020 {
01021 int aid = a[i].id;
01022 BigReal f1 = exp( coefficient * molecule->langevin_param(aid) );
01023 a[i].velocity *= f1;
01024 }
01025 }
01026 }
|
|
|
Definition at line 1702 of file Sequencer.C. Referenced by algorithm(), maximumMove(), and rattle1(). 01702 {
01703 CthFree(thread);
01704 CthSuspend();
01705 }
|
|
|
Definition at line 23 of file Sequencer.h. |
|
|
Definition at line 83 of file Sequencer.h. Referenced by algorithm(), berendsenPressure(), and Sequencer(). |
|
|
Definition at line 101 of file Sequencer.h. Referenced by algorithm(), berendsenPressure(), correctMomentum(), cycleBarrier(), langevinPiston(), minimize(), rescaleVelocities(), Sequencer(), and tcoupleVelocities(). |
|
|
Definition at line 84 of file Sequencer.h. Referenced by algorithm(). |
|
|
Definition at line 100 of file Sequencer.h. Referenced by submitCollections(). |
|
|
Definition at line 103 of file Sequencer.h. Referenced by rebalanceLoad(). |
|
|
Definition at line 39 of file Sequencer.h. Referenced by runComputeObjects(). |
|
|
Definition at line 38 of file Sequencer.h. Referenced by algorithm(), rebalanceLoad(), and runComputeObjects(). |
|
|
|
Definition at line 98 of file Sequencer.h. Referenced by rattle1(), Sequencer(), submitHalfstep(), and submitReductions(). |
|
|
Definition at line 94 of file Sequencer.h. Referenced by langevinVelocities(), langevinVelocitiesBBK2(), reassignVelocities(), reinitVelocities(), and Sequencer(). |
|
|
Definition at line 97 of file Sequencer.h. Referenced by rattle1(), rattle2(), redistrib_swm4_forces(), redistrib_tip4p_forces(), runComputeObjects(), Sequencer(), submitHalfstep(), submitMinimizeReductions(), submitMomentum(), and submitReductions(). |
|
|
Definition at line 77 of file Sequencer.h. Referenced by rescaleVelocities(), and Sequencer(). |
|
|
|
Definition at line 86 of file Sequencer.h. Referenced by integrate(), and langevinPiston(). |
1.3.9.1