Difference for src/Sequencer.C from version 1.1229 to 1.1230

version 1.1229version 1.1230
Line 7
Line 7
 /***************************************************************************** /*****************************************************************************
  * $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
Line 292
Line 292
       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
Line 314
Line 306
     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()){
Line 340
Line 328
       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
Line 371
Line 349
       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 {
Line 431
Line 409
       // 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);
       }       }
  
Line 475
Line 431
       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);
Line 1099
Line 1045
 } }
 // --------- 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,
Line 1106
Line 1073
 // 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;
Line 1123
Line 1090
  
     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();
     }     }
Line 1600
Line 1569
   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);

Legend:
Removed in v.1.1229 
changed lines
 Added in v.1.1230


Made by using version 1.53 of cvs2html