Fixing Langevin piston to allow off-diagonal stress (and why does the source say "FIX THIS"?)

From: FX (fxcoudert_at_gmail.com)
Date: Wed Mar 06 2013 - 04:43:49 CST

Hi all,

I am currently trying to modify the NAMD sources related to the Langevin piston barostat method, in order to allow off-diagonal stress fluctuations (shearing). I think I have gotten the parts in Controller.C okay (my "first attempt" patch is attached, but I have not yet tested it), but I do not understand well what is done in Sequencer::langevinPiston:

> void Sequencer::langevinPiston(int step)
> {
> if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
> {
> FullAtom *a = patch->atom.begin();
> int numAtoms = patch->numAtoms;
> Tensor factor = broadcast->positionRescaleFactor.get(step);
> // JCP FIX THIS!!!
> Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
> patch->lattice.rescale(factor);

It seems to me that the position rescaling code will work fine, even with my modifications to the calculation of the rescale factor. However, I have the following questions:

  1. Why is the velocity scaled by the inverse of the position scaling factor? I cannot get that from the Langevin piston equations…
  2. Who is JCP? (James C. Phillips?) Does he read this list?
  3. How much should I be scared of the "FIX THIS" comment?

Thanks in advance for any help you can give me.

François-Xavier Coudert

diff -rpu NAMD_2.9_Source-orig/src/Controller.C NAMD_2.9_Source/src/Controller.C
--- NAMD_2.9_Source-orig/src/Controller.C 2012-02-24 02:25:01.000000000 +0100
+++ NAMD_2.9_Source/src/Controller.C 2013-03-06 11:31:07.000000000 +0100
@@ -704,14 +704,14 @@ void Controller::langevinPiston1(int ste
       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
       strainRate *= f1;
       if ( simParams->useFlexibleCell ) {
- // We only use on-diagonal terms (for now)
         if ( simParams->useConstantRatio ) {
           BigReal r = f2 * random->gaussian();
           strainRate.xx += r;
           strainRate.yy += r;
           strainRate.zz += f2 * random->gaussian();
         } else {
- strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
+ // FXC
+ strainRate += f2 * Tensor::symmetric(random->gaussian_vector(), random->gaussian_vector());
         }
       } else {
         strainRate += f2 * Tensor::identity(random->gaussian());
@@ -746,15 +746,22 @@ void Controller::langevinPiston1(int ste
 
     if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
     {
- // We only use on-diagonal terms (for now)
       Tensor factor;
- if ( !simParams->useConstantArea ) {
+ if ( simParams->useConstantArea ) {
+ factor.xx = factor.yy = 1;
+ factor.zz = exp( dt_long * strainRate.zz );
+ } else {
+ // FXC
         factor.xx = exp( dt_long * strainRate.xx );
+ factor.xy = exp( dt_long * strainRate.xy );
+ factor.xz = exp( dt_long * strainRate.xz );
+ factor.yx = exp( dt_long * strainRate.yx );
         factor.yy = exp( dt_long * strainRate.yy );
- } else {
- factor.xx = factor.yy = 1;
+ factor.yz = exp( dt_long * strainRate.yz );
+ factor.zx = exp( dt_long * strainRate.zx );
+ factor.zy = exp( dt_long * strainRate.zy );
+ factor.zz = exp( dt_long * strainRate.zz );
       }
- factor.zz = exp( dt_long * strainRate.zz );
       broadcast->positionRescaleFactor.publish(step,factor);
       state->lattice.rescale(factor);
 #ifdef DEBUG_PRESSURE
@@ -854,14 +861,14 @@ void Controller::langevinPiston2(int ste
       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
       strainRate *= f1;
       if ( simParams->useFlexibleCell ) {
- // We only use on-diagonal terms (for now)
         if ( simParams->useConstantRatio ) {
           BigReal r = f2 * random->gaussian();
           strainRate.xx += r;
           strainRate.yy += r;
           strainRate.zz += f2 * random->gaussian();
         } else {
- strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
+ // FXC
+ strainRate += f2 * Tensor::symmetric(random->gaussian_vector(), random->gaussian_vector());
         }
       } else {
         strainRate += f2 * Tensor::identity(random->gaussian());

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:00 CST