version 1.1229 | version 1.1230 |
---|
| |
/***************************************************************************** | /***************************************************************************** |
* $Source: /home/cvs/namd/cvsroot/namd2/src/Sequencer.C,v $ | * $Source: /home/cvs/namd/cvsroot/namd2/src/Sequencer.C,v $ |
* $Author: jim $ | * $Author: jim $ |
* $Date: 2016/03/02 21:33:06 $ | * $Date: 2016/08/26 19:40:32 $ |
* $Revision: 1.1229 $ | * $Revision: 1.1230 $ |
*****************************************************************************/ | *****************************************************************************/ |
| |
//for gbis debugging; print net force on each atom | //for gbis debugging; print net force on each atom |
| |
if ( doFullElectrostatics ) saveForce(Results::slow); | if ( doFullElectrostatics ) saveForce(Results::slow); |
} | } |
if ( ! commOnly ) { | if ( ! commOnly ) { |
addForceToMomentum(-0.5*timestep); | newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1); |
if (staleForces || doNonbonded) | |
addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0); | |
if (staleForces || doFullElectrostatics) | |
addForceToMomentum(-0.5*slowstep,Results::slow,staleForces); | |
} | } |
minimizationQuenchVelocity(); | minimizationQuenchVelocity(); |
rattle1(-timestep,0); | rattle1(-timestep,0); |
submitHalfstep(step); | submitHalfstep(step); |
if ( ! commOnly ) { | if ( ! commOnly ) { |
addForceToMomentum(timestep); | newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1); |
if (staleForces || doNonbonded) | |
addForceToMomentum(nbondstep,Results::nbond,staleForces,0); | |
if (staleForces || doFullElectrostatics) | |
addForceToMomentum(slowstep,Results::slow,staleForces,0); | |
} | } |
rattle1(timestep,1); | rattle1(timestep,1); |
if (doTcl || doColvars) // include constraint forces | if (doTcl || doColvars) // include constraint forces |
| |
submitHalfstep(step); | submitHalfstep(step); |
if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step); | if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step); |
if ( ! commOnly ) { | if ( ! commOnly ) { |
addForceToMomentum(-0.5*timestep); | newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1); |
if (staleForces || doNonbonded) | |
addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1); | |
if (staleForces || doFullElectrostatics) | |
addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1); | |
} | } |
submitReductions(step); | submitReductions(step); |
if(traceIsOn()){ | if(traceIsOn()){ |
| |
berendsenPressure(step); | berendsenPressure(step); |
| |
if ( ! commOnly ) { | if ( ! commOnly ) { |
if (staleForces || (doNonbonded && doFullElectrostatics)) { | newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); |
addForceToMomentum3(0.5*timestep,Results::normal,0, | |
0.5*nbondstep,Results::nbond,staleForces, | |
0.5*slowstep,Results::slow,staleForces); | |
} else { | |
addForceToMomentum(0.5*timestep); | |
if (staleForces || doNonbonded) | |
addForceToMomentum(0.5*nbondstep,Results::nbond,staleForces,1); | |
if (staleForces || doFullElectrostatics) | |
addForceToMomentum(0.5*slowstep,Results::slow,staleForces,1); | |
} | |
} | } |
| |
// We do RATTLE here if multigrator thermostat was applied in the previous step | // We do RATTLE here if multigrator thermostat was applied in the previous step |
| |
if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) { | if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) { |
if ( ! commOnly ) addVelocityToPosition(0.5*timestep); | if ( ! commOnly ) addVelocityToPosition(0.5*timestep); |
// We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin) | // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin) |
if ( simParams->langevinOn && simParams->langevin_useBAOAB ) langevinVelocities(timestep); | langevinVelocities(timestep); |
langevinPiston(step); | langevinPiston(step); |
if ( ! commOnly ) addVelocityToPosition(0.5*timestep); | if ( ! commOnly ) addVelocityToPosition(0.5*timestep); |
} else { | } else { |
| |
// reassignment based on full-step velocities | // reassignment based on full-step velocities |
if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) { | if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) { |
reassignVelocities(timestep,step); | reassignVelocities(timestep,step); |
if (staleForces || (doNonbonded && doFullElectrostatics)) { | newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); |
addForceToMomentum3(-0.5*timestep,Results::normal,0, | |
-0.5*nbondstep,Results::nbond,staleForces, | |
-0.5*slowstep,Results::slow,staleForces); | |
} else { | |
addForceToMomentum(-0.5*timestep); | |
if (staleForces || doNonbonded) | |
addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0); | |
if (staleForces || doFullElectrostatics) | |
addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,0); | |
} | |
rattle1(-timestep,0); | rattle1(-timestep,0); |
} | } |
| |
if ( ! commOnly ) { | if ( ! commOnly ) { |
langevinVelocitiesBBK1(timestep); | langevinVelocitiesBBK1(timestep); |
if (staleForces || (doNonbonded && doFullElectrostatics)) { | newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); |
addForceToMomentum3(timestep,Results::normal,0, | |
nbondstep,Results::nbond,staleForces, | |
slowstep,Results::slow,staleForces); | |
} else { | |
addForceToMomentum(timestep); | |
if (staleForces || doNonbonded) { | |
addForceToMomentum(nbondstep,Results::nbond,staleForces,1); | |
} | |
if (staleForces || doFullElectrostatics) { | |
addForceToMomentum(slowstep,Results::slow,staleForces,1); | |
} | |
} | |
langevinVelocitiesBBK2(timestep); | langevinVelocitiesBBK2(timestep); |
} | } |
| |
| |
if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step); | if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step); |
| |
if ( ! commOnly ) { | if ( ! commOnly ) { |
if (staleForces || (doNonbonded && doFullElectrostatics)) { | newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); |
addForceToMomentum3(-0.5*timestep,Results::normal,0, | |
-0.5*nbondstep,Results::nbond,staleForces, | |
-0.5*slowstep,Results::slow,staleForces); | |
} else { | |
addForceToMomentum(-0.5*timestep); | |
if (staleForces || doNonbonded) | |
addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1); | |
if (staleForces || doFullElectrostatics) | |
addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1); | |
} | |
} | } |
| |
// rattle2(timestep,step); | // rattle2(timestep,step); |
| |
} | } |
// --------- End Multigrator --------- | // --------- End Multigrator --------- |
| |
| void Sequencer::newtonianVelocities(BigReal stepscale, const BigReal timestep, |
| const BigReal nbondstep, |
| const BigReal slowstep, |
| const int staleForces, |
| const int doNonbonded, |
| const int doFullElectrostatics) |
| { |
| // Deterministic velocity update, account for multigrator |
| if (staleForces || (doNonbonded && doFullElectrostatics)) { |
| addForceToMomentum3(stepscale*timestep, Results::normal, 0, |
| stepscale*nbondstep, Results::nbond, staleForces, |
| stepscale*slowstep, Results::slow, staleForces); |
| } else { |
| addForceToMomentum(stepscale*timestep); |
| if (staleForces || doNonbonded) |
| addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces); |
| if (staleForces || doFullElectrostatics) |
| addForceToMomentum(stepscale*slowstep, Results::slow, staleForces); |
| } |
| } |
| |
void Sequencer::langevinVelocities(BigReal dt_fs) | void Sequencer::langevinVelocities(BigReal dt_fs) |
{ | { |
// This routine is used for the BAOAB integrator, | // This routine is used for the BAOAB integrator, |
| |
// See B. Leimkuhler and C. Matthews, AMRX (2012) | // See B. Leimkuhler and C. Matthews, AMRX (2012) |
// Routine originally written by JPhillips, with fresh errors by CMatthews June2012 | // Routine originally written by JPhillips, with fresh errors by CMatthews June2012 |
| |
if ( simParams->langevinOn ) | if ( simParams->langevinOn && simParams->langevin_useBAOAB ) |
{ | { |
FullAtom *a = patch->atom.begin(); | FullAtom *a = patch->atom.begin(); |
int numAtoms = patch->numAtoms; | int numAtoms = patch->numAtoms; |
| |
| |
for ( int i = 0; i < numAtoms; ++i ) | for ( int i = 0; i < numAtoms; ++i ) |
{ | { |
BigReal f1 = exp( -1. * dt * a[i].langevinParam ); | BigReal dt_gamma = dt * a[i].langevinParam; |
| if ( ! dt_gamma ) continue; |
| |
| BigReal f1 = exp( -dt_gamma ); |
BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT * | BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT * |
( a[i].partition ? tempFactor : 1.0 ) / a[i].mass ); | ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass ); |
| |
a[i].velocity *= f1; | a[i].velocity *= f1; |
a[i].velocity += f2 * random->gaussian_vector(); | a[i].velocity += f2 * random->gaussian_vector(); |
} | } |
| |
patch->saveForce(ftag); | patch->saveForce(ftag); |
} | } |
| |
void Sequencer::addForceToMomentum(BigReal dt, const int ftag, | void Sequencer::addForceToMomentum(BigReal dt, const int ftag, const int useSaved) |
const int useSaved, int pressure) | |
{ | { |
#if CMK_BLUEGENEL | #if CMK_BLUEGENEL |
CmiNetworkProgressAfter (0); | CmiNetworkProgressAfter (0); |