[gmx-users] Force evaluation in umbrella pulling with nonzero pull_rate

[Lin, Dejun]

Hi,

I'm trying to use umbrella pulling to do a targeted MD and I have the following pull code:
pull                    = umbrella
pull_geometry   = distance     
pull_dim            = N N Y
pull_start           = no
pull_ngroups     = 1
pull_group0      = REF
pull_group1      = PULL
pull_init1           = 5
pull_rate1         = -0.0005      ; nm/ps
pull_k1              = 1000          ; kJ mol^-1 nm^-2

which gives the results:
(position.xvg)
Time      Z(nm)    dZ(nm)
0.0000  12.1868 5.04599
0.0100  12.1868 5.04576
0.0200  12.1869 5.04552
0.0300  12.1869 5.0453
0.0400  12.1869 5.04508
0.0500  12.1869 5.04486
0.0600  12.187  5.04465
0.0700  12.1869 5.04447
0.0800  12.1869 5.04423
0.0900  12.1869 5.044
0.1000  12.1869 5.04376
0.1100  12.1869 5.04351
0.1200  12.1868 5.04333
0.1300  12.1868 5.04304
0.1400  12.1868 5.04274
0.1500  12.1868 5.04242
0.1600  12.1868 5.04208
0.1700  12.1866 5.04183
0.1800  12.1866 5.04145
0.1900  12.1866 5.04105
0.2000  12.1866 5.04063
0.2100  12.1866 5.0402
0.2200  12.1864 5.03985
0.2300  12.1864 5.03939
0.2400  12.1864 5.03891
0.2500  12.1864 5.03843
0.2600  12.1864 5.03795
0.2700  12.186  5.03762
0.2800  12.186  5.03715
0.2900  12.186  5.03669
0.3000  12.186  5.03626
0.3100  12.186  5.03586

(force.xvg):
Time      Force
0.0000  -45.9932
0.0100  -45.7626
0.0200  -45.5344
0.0300  -45.3134
0.0400  -45.0967
0.0500  -44.8854
0.0600  -44.6768
0.0700  -44.5008
0.0800  -44.2741
0.0900  -44.045
0.1000  -43.8105
0.1100  -43.5675
0.1200  -43.3877
0.1300  -43.107
0.1400  -42.8104
0.1500  -42.4967
0.1600  -42.1648
0.1700  -41.915
0.1800  -41.538
0.1900  -41.1422
0.2000  -40.7282
0.2100  -40.3026
0.2200  -39.9637
0.2300  -39.5006
0.2400  -39.0298
0.2500  -38.5553
0.2600  -38.0832
0.2700  -37.7541
0.2800  -37.2857
0.2900  -36.8358
0.3000  -36.4099
0.3100  -36.0136

Can anyone tell me If the equilibrium position of the applied harmonic potential is the same as it is in pull_geometry = position (which is P_eq = P_REF + pull_init1 + time*pull_rate1*pull_vec1, where P_REF is the position of the reference group in my case ) and if the change of the equilibrium position of the applied harmonic potential is evaluated every single step (in my case every 0.01 ps )or every ps? 

What's more, is the force in the output evaluated according to the following equation so that only the z-component of the force ( which is a vector ) is reported in force.xvg?
 f = pull_k1*( P_eq - P_PULL ) 
   = pull_k1*( P_REF - P_PULL + pull_init1 + time*pull_rate1*pull_vec1)
   = pull_k1*( pull_init1 - dZ - time*pull_rate1*pull_vec1) , where P_PULL is the position of the pulled group and dZ is  the values in the 3rd column of position.xvg output.

Finally, is umbrella pulling a good choice to do targeted MD ( the 'fast-growth' method proposed by C.Jarzynski )? Is there better way to do that in Gromacs? 

Thanks,
Dejun Lin