[gmx-users] Possible bug with -rerun option

[Alexandre Suman de Araujo]

Hi GMXers

I observed a strange behavior in a simulation of an equilibrated box
with 1409 tip3p water molecules using OPLS-AA, PME and GROMACS 3.3.1 in
double precision.
My .mdp file is as fallow:

;
;       Input file
;
title               =  Simulacao de calculo   
cpp                 =  /lib/cpp               
dt                  =  0.001                 
integrator          =  md                    
nsteps              =  10000                  
nstcomm             =  1                       
nstxout             =  100                     
nstvout             =  100                    
nstlog              =  100                    
nstenergy           =  100                 
nstlist             =  10                      
ns_type             =  grid                  
pbc                 =  xyz                   
rlist               =  0.9               
coulombtype         =  PME                   
rcoulomb            =  0.9                   
vdwtype             =  cut-off             
rvdw                =  0.9                   
tcoupl              =  berendsen        
ref_t               =  300                  
tc-grps             =  System               
tau_t               =  0.1                   
pcoupl              =  berendsen            
pcoupltype          =  isotropic             
tau_p               =  1                    
compressibility     =  4.5e-5
ref_p               =  1.0

Finishing the simulation I run g_energy in the .edr, obtaining the
follow results:

Statistics over 10001 steps [ 0.0000 thru 10.0000 ps ], 3 data sets

Energy                      Average       RMSD     Fluct.      Drift 
Tot-Drift
-------------------------------------------------------------------------------
Potential                  *-55852.6*    136.692        136  
-4.75933    -47.598
Total-Energy               *-45314.7*    78.2564    77.8627   -2.71554  
-27.1581
Temperature                 *299.942*    3.15335    3.14887   0.058173  
0.581788
Heat Capacity Cv:      12.4738 J/mol K (factor = 0.000110527)

Until here everything is ok. So I used the same .tpr file and -rerun
option in mdrun_d in the trajectory from this first simulation. So I ran
g_energy in the new .edr file and obtained:

Statistics over 10001 steps [ 0.0000 thru 10.0000 ps ], 3 data sets

Energy                      Average       RMSD     Fluct.      Drift 
Tot-Drift
-------------------------------------------------------------------------------
Potential                  *-563.871*    5526.06    5526.05  
-4.75933    -47.598
Total-Energy               *-457.605*    4484.56    4484.55   
-2.9575    -29.578
Temperature                 *3.02467*    29.6455    29.6451  0.0512859  
0.512911
Heat Capacity Cv:   -0.0871562 J/mol K (factor = 96.0644)

It can be seen that there are differences of approximately a factor 100
between Average values.


In the .log files I saw the follow results:

First Simulation:

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    8.90375e+03   -5.73880e+04   -7.13930e+03   *-5.56236e+04 *  
1.04571e+04
   Total Energy    Temperature Pressure (bar)
   *-4.51665e+04*    *2.97643e+02*    4.12898e+02

           Step           Time         Lambda
          10000       10.00000        0.00000

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    8.36435e+03   -5.68959e+04   -7.10656e+03   *-5.56381e+04*   
1.05288e+04
   Total Energy    Temperature Pressure (bar)
   *-4.51094e+04*    *2.99683e+02*   -2.69511e+02

        <======  ###############  ==>
        <====  A V E R A G E S  ====>
        <==  ###############  ======>

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    8.68022e+03   -5.74019e+04   -7.13085e+03   *-5.58526e+04*   
1.05378e+04
   Total Energy    Temperature Pressure (bar)
   *-4.53147e+04*    *2.99942e+02*    1.07997e+01





Second Simulation (rerun):

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    8.90375e+03   -5.73880e+04   -7.13930e+03   *-5.56236e+04*   
1.04594e+04
   Total Energy    Temperature Pressure (bar)
   *-4.51641e+04*    *2.97709e+02*    4.13037e+02

           Step           Time         Lambda
          10000       10.00000        0.00000

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    8.36435e+03   -5.68959e+04   -7.10656e+03   *-5.56381e+04*   
1.05293e+04
   Total Energy    Temperature Pressure (bar)
   *-4.51088e+04*    *2.99699e+02*   -2.70138e+02

        <======  ###############  ==>
        <====  A V E R A G E S  ====>
        <==  ###############  ======>

   Energies (kJ/mol)
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    8.79288e+01   -5.79793e+02   -7.20063e+01   *-5.63871e+02*   
1.06266e+02
   Total Energy    Temperature Pressure (bar)
   *-4.57605e+02*    *3.02467e+00*    4.68805e-01

Here we observed again that in the first simulation everything is ok,
but in the "reran" simulation the Average values are approximately 100
times lesser than the values obtained in every step (the values are in
bold)!!!
I made this test with GROMACS 2.2.1 in single precision and this
behavior happened too.

This could be a bug??? Or am I doing something wrong????

Best regards.

-- 
Alexandre Suman de Araujo
asaraujo at if.sc.usp.br
UIN: 6194055
IFSC - USP - São Carlos - Brasil