[gmx-users] free_energy = yes & init_lambda = 0.00 does not yield identical energies to free_energy = no.

[chris.neale at utoronto.ca]

Hi Berk,

the difference remains when utilizing double precision:
free energy code on with lambda=0.00: Coulomb (SR) = -3.23064e+05
free energy code off: Coulomb (SR) = -3.22653e+05

Here are the full double precision values:

; Free energy control stuff
free_energy              = no

    Energies (kJ/mol)
           Angle    Proper Dih. Ryckaert-Bell.          LJ-14     Coulomb-14
     2.91166e+03    5.90834e+02    1.07820e+03    7.70226e+02    6.22283e+03
         LJ (SR)        LJ (LR)  Disper. corr.   Coulomb (SR)   Coul. recip.
     4.94952e+04   -1.84944e+03   -6.38472e+02   -3.22653e+05   -7.12639e+04
       Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
    -3.35336e+05    5.67661e+04   -2.78570e+05    3.00586e+02    6.48848e-01
   Cons. rmsd () Cons.2 rmsd ()
     0.00000e+00    0.00000e+00


; Free energy control stuff
free_energy              = yes
init_lambda              = 0.00
delta_lambda             = 0
sc_alpha                 =0.0
sc-power                 =1.0
sc-sigma                 = 0.3

           Angle    Proper Dih. Ryckaert-Bell.          LJ-14     Coulomb-14
     2.91166e+03    5.90834e+02    1.07820e+03    7.70226e+02    6.22283e+03
         LJ (SR)        LJ (LR)  Disper. corr.   Coulomb (SR)   Coul. recip.
     4.94952e+04   -1.84944e+03   -6.38472e+02   -3.23064e+05   -7.12639e+04
       Potential    Kinetic En.   Total Energy    Temperature Pressure (bar)
    -3.35747e+05    5.67544e+04   -2.78993e+05    3.00524e+02   -4.14918e+01
   dVpot/dlambda  dEkin/dlambda  dG/dl constr.  Cons. rmsd () Cons.2 rmsd ()
     1.29815e+03    0.00000e+00    0.00000e+00    0.00000e+00    0.00000e+00

I'll try to develop a simplet system to see if this is reproducible  
between systems or is system specific.

I am willing to post my topologies, but what would be most useful at  
this point is some assistance confirming if anybody else encounters  
this difference.

Thanks,
Chris.

--- original message ---

Hi,

I would also expect the same results.
I have no clue what kind of system you have.
But check that it is not an accuracy issue.
Could you recalculate the energies in double precision.

Berk

> Date: Tue, 27 Jan 2009 20:00:23 -0500
> From: chris.neale at utoronto.ca
> To: gmx-users at gromacs.org
> Subject: [gmx-users] free_energy = yes & init_lambda = 0.00 does not  
> yield identical energies to free_energy = no.
>
> Hello,
>
> when applying the free energy code, I expected that the  
> instantaneous energy of a given conformation when free_energy = yes  
> & init_lambda = 0.00 would be identical to that when free_energy =  
> no. However, I have found this not to be the case in zero-step  
> mdrun. This is the same result that I get from gmx v 3.3.1, 3.3.3,  
> and 4.0.3. Is this expected?
>
> With free_energy = yes & init_lambda = 0.00:
> Coulomb (SR)    -3.23064e+05
>
> With free_energy = no:
> Coulomb (SR)    -3.22654e+05
>
> and all other energy components are identical. The Coulomb (SR)  
> difference seems significant to me. Note that the free_energy = no  
> values are identical to those that I get when I use a topology that  
> has not been modified to add a B-state so I don't think that my  
> A-state values are incorrect.
>
> ## free_energy related .mdp options:
> free_energy              = yes
> init_lambda              = 0.00
> delta_lambda             = 0
> sc_alpha                 =0.0
> sc-power                 =1.0
> sc-sigma                 = 0.3
>
> ## Other .mdp options:
> nsteps              = 0
> tinit               =  0
> dt                  =  0.004
> integrator          =  sd
> comm_mode           =  linear
> nstcomm             =  1
> comm_grps           =  System
> nstlog              =  2500
> nstlist             =  5
> ns_type             =  grid
> pbc                 =  xyz
> coulombtype         =  PME
> rcoulomb            =  0.9
> fourierspacing      =  0.12
> pme_order           =  4
> vdwtype             =  cut-off
> rvdw_switch         =  0
> rvdw                =  1.4
> rlist               =  0.9
> DispCorr            =  EnerPres
> Pcoupl              =  Berendsen
> pcoupltype          =  isotropic
> compressibility     =  4.5e-5
> ref_p               =  1.
> tau_p               =  4.0
> tcoupl              =  Berendsen
> tc_grps             =  System
> tau_t               =  0.1       ref_t               =  300.       
> annealing           =  no
> gen_vel             =  no
> unconstrained-start =  yes
> gen_temp            =  300.
> gen_seed            =  9896
> constraints         =  all-bonds
> constraint_algorithm=  lincs
> lincs-iter          =  1
> lincs-order         =  6
>
>
> Thank you,
> Chris.