free_energy=yes affects energies at lambda=0
I'm using mdrun -rerun to reprocess a (trr) trajectory with several different free energy settings for
purposes of reweighting using my own WHAM installation. However, I found something rather alarming:
The total potential energies I get, even with lambda=0, depend on whether free_energy=yes in the tpr
files I use for reprocessing. For example, from the same trajectory, here are some energy values from
With free_energy=yes, at lambda=0
< 0.000000 -278065.687500
< 100.000008 -278276.562500
< 200.000015 -278689.375000
< 300.000000 -278253.187500
< 400.000031 -278992.375000
< 500.000031 -278813.750000
< 600.000000 -278452.562500
As you can see, the values differ by 0.2 kJ/mol or so -- but as far as I can tell from the documentation,
at lambda=0, it should be irrelevant whether free energy is turned on or off.
This turns out to be a huge problem for what I'm doing, as it's absolutely essential that lambda=0 be
equivalent to having no perturbation on the system, so I really need some help here.
I am providing a link to a tarball to reproduce the problem. Files are as follows:
- L7.trr: An example (short) trajectory.
- orig.mdp, orig.tpr, orig.ene, orig_pot.xvg: "Original" mdp file, input for mdrun
rerun, and output noFE.mdp, noFE.tpr, noFE.ene, noFE_pot.xvg: Same mdp file but with free_energy=no; associated tpr
energies from mdrun -rerun; in these, free_energy=yes
files and output from mdrun
rerun L7.top, L7.gro: Top and gro files used to generate the tpr files.
I am using single precision 3.3 with various bugfixes applied. I assume the same thing still happens in
#1 Updated by David van der Spoel over 13 years ago
as always, examples to reproduce the problem will be needed. I also wonder
wheter you have tried it in double precision? It could be that another inner
loop is used which although mathematically equivalent is nt numerically
equivalent. DP should fix that. Also, have you investigated which term in the
energy causes the difference?
#3 Updated by David Mobley over 13 years ago
Yes, I tried in double precision (after posting this) and the difference is gone, as well.
I like Berk's suggestion (on list) of making stuff internally be double precision -- I do still find it alarming
that I could get such substantial (for me) energy differences just because stuff is in single precision.
#4 Updated by David van der Spoel over 13 years ago
what you call significant is in fact the last bit in the float representation
which is 7 decimals, i.e. we're looking at roundoff errors. Maybe you could look
into defining your problem in a different manner such that you compute
differences. Of course, this is why people usually shy away from computing
absolute free energies. Either way if you know which term of the energy causes
this it might still be cured by modifying the code for the free energy
calculation although it will be very difficult to make it water tight.