search

PiRSquared17 / sire

Automatically exported from code.google.com/p/sire
0 stars 0 forks source link

Gradients are not zero for null perturbation at lambda 0.00 only using OpenMMFrEnergyST #9

Open GoogleCodeExporter opened 9 years ago

GoogleCodeExporter commented 9 years ago 
Using a MORPH file in explicit solvent with no actual changes in parameters 
(initial=final). 

python debug.py 0.0 1; python debug.py 0.1 1; python debug.py 0.2 1; python 
debug.py 0.3 1; python debug.py 0.4 1; python debug.py 0.5 1; python debug.py 
0.6 1; python debug.py 0.7 1; python debug.py 0.8 1; python debug.py 0.9 1; 
python debug.py 1.0 1; cat gradient*
    1         0.0031697021
    1        -2.1469100520
    1         2.0698155911
    1         0.0698637728
    1        -0.0729800251
    1        -0.0000000000
    1        -0.0000000000
    1        -0.0000000000
    1        -0.0000000000
    1        -0.0000000000
    1        -0.0000000000

The gradients are not zero for lambda 0.0, 0.1, 0.2, 0.3, 0.4 and strictly 0.0 
for lambda >= 0.5
Running with lambda=0.49 -->  
    1        -0.0014132330
With lambda 0.51 --> 
      1        -0.

Original issue reported on code.google.com by julienm...@googlemail.com on 12 Jul 2013 at 10:58

Attachments:

GoogleCodeExporter commented 9 years ago 
Here is an update. 
Everything run using module sire/jm that has merged the sire/gac openmm code 
and is linked with OpenMM5.1 on section6.

The problem seems to be related to numerical errors with the calculation of the 
total energy. Varying delta_lambda has a large effect, with low values creating 
larger errors. Essentially, a low value of delta_lambda magnifies numerical 
errors. All lambda values are affected, but it is more apparent in general for 
the first half of the perturbation.  The selection of an OpenCL or CUDA 
single/mixed/double platform does not change dramatically the behavior. The 
details of the gradients computed using different delta_lambda values are 
pasted below. 

Delta 0.005 works best here, but his is likely system dependent. 

Numerical errors on gradients should be on the order of the second decimal to 
have a negligible effect on computed free energy changes. 

The error may grow for larger systems (protein simulations) because the total 
energy of the system is recomputed from scratch. 
A workaround could be to define a custom hamltonian that only computes the 
energy of the perturbed solute (since solvent-solvent terms cancel out in the 
calculation of the gradients). Alternatively, one that runs through everything 
but reset to 0 any solute-solvent term may be more numerically stable. 

Running with a Reference platform doesn't work (I get Nan energies).

julien@node004:methanol-discharge-solvated$ cat 
delta-0.02-mixed-opencl-repeat.dat 
 0.00000  0.00000
 0.10000 -0.10386
 0.20000 -0.00001
 0.30000  0.00047
 0.40000  0.00010
 0.50000  0.00113
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00000
 0.90000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-solvated$ cat delta-0.01-mixed-opencl.dat 
 0.00000 -0.28384
 0.10000  0.14098
 0.20000  0.00000
 0.30000  0.11158
 0.40000  0.20868
 0.50000  0.00000
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00000
 0.90000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-solvated$ cat delta-0.005-mixed-opencl.dat 
 0.00000  0.00000
 0.10000  0.00033
 0.20000  0.00000
 0.30000 -0.00132
 0.40000  0.00005
 0.50000  0.00000
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00000
 0.90000 -0.00007
 1.00000  0.00000
julien@node004:methanol-discharge-solvated$ cat delta-0.001-mixed-opencl.dat 
 0.00000  0.02284
 0.10000 -2.13764
 0.20000  2.10732
 0.30000  0.07431
 0.40000 -0.06951
 0.50000  0.00000
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00000
 0.90000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-solvated$ cat delta-0.0001-mixed-opencl.dat 
 0.00000 20.79240
 0.10000  0.00114
 0.20000  0.81962
 0.30000  0.00057
 0.40000 -0.02878
 0.50000  0.03818
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00057
 0.90000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-solvated$ cat delta-0.00001-mixed-opencl.dat 
cat: delta-0.00001-mixed-opencl.dat: No such file or directory
julien@node004:methanol-discharge-solvated$ cat delta-0.000001-mixed-opencl.dat 
cat: delta-0.000001-mixed-opencl.dat: No such file or directory
julien@node004:methanol-discharge-solvated$ cat delta-0.0001-mixed-opencl.dat 
 0.00000 20.79240
 0.10000  0.00114
 0.20000  0.81962
 0.30000  0.00057
 0.40000 -0.02878
 0.50000  0.03818
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00057
 0.90000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-solvated$ cat delta-0.00001-mixed-opencl.dat 
cat: delta-0.00001-mixed-opencl.dat: No such file or directory
julien@node004:methanol-discharge-solvated$ cat 
delta-0.000001-mixed-opencl-repeat.dat 
 0.00000 -0.61257
 0.10000 208.14986
 0.20000 -79.88378
 0.30000 -147.14732
 0.40000  1.15320
 0.50000  0.00000
 0.60000  0.00000
 0.70000  0.00000
 0.80000  0.00000
 0.90000  0.00000
 1.00000  0.00000

Original comment by julienm...@googlemail.com on 12 Jul 2013 at 12:12

GoogleCodeExporter commented 9 years ago 
Here is an update. I have run the same null perturbation, but this time in 
vacuum. The major difference is that there aren't solute-solvent and 
solvent-solvent interactions, thus the number of pairwise interactions is much 
smaller. The results for three values of delta-lambda are below (all OpenCL 
mixed precisio)n

julien@node004:methanol-discharge-vacuum-testdeltalambda$ cat 
grads-delta-0.01.dat 
 0.00000  0.00000
 0.05000  0.00000
 0.10000  0.00002
 0.15000  0.00002
 0.20000  0.00000
 0.25000  0.00000
 0.30000  0.00000
 0.35000  0.00002
 0.40000  0.00002
 0.45000  0.00000
 0.50000  0.00000
 0.55000  0.00000
 0.60000  0.00000
 0.65000  0.00000
 0.70000  0.00000
 0.75000  0.00000
 0.80000  0.00000
 0.85000  0.00000
 0.90000  0.00000
 0.95000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-vacuum-testdeltalambda$ cat 
grads-delta-0.001.dat 
 0.00000  0.00000
 0.05000  0.00000
 0.10000 -0.00017
 0.15000 -0.00017
 0.20000 -0.00001
 0.25000  0.00017
 0.30000  0.00000
 0.35000 -0.00001
 0.40000 -0.00017
 0.45000  0.00000
 0.50000  0.00000
 0.55000  0.00000
 0.60000  0.00000
 0.65000  0.00000
 0.70000  0.00000
 0.75000  0.00000
 0.80000  0.00001
 0.85000  0.00000
 0.90000  0.00000
 0.95000  0.00000
 1.00000  0.00000
julien@node004:methanol-discharge-vacuum-testdeltalambda$ cat 
grads-delta-0.0001.dat 
 0.00000  0.00000
 0.05000  0.00000
 0.10000  0.00169
 0.15000  0.00169
 0.20000  0.00168
 0.25000 -0.00177
 0.30000  0.00178
 0.35000  0.00000
 0.40000 -0.00006
 0.45000  0.00000
 0.50000 -0.00176
 0.55000  0.00000
 0.60000  0.00000
 0.65000  0.00000
 0.70000  0.00000
 0.75000  0.00000
 0.80000  0.00173
 0.85000  0.00000
 0.90000  0.00000
 0.95000  0.00000
 1.00000  0.00000

* The error is much smaller than for the solvated runs and essentially 
negligible. It grows as the value of delta-lambda is decreased. 

Recommendation: 

To consider a revised implementation where, to evaluate gradients, the system 
energy at lambda, lambda+delta_lambda and lambda-delta_lambda is computed with 
a custom force field expression that multiplies by 0 any solvent-solvent terms. 
This may remove some round off error due to accumulation of many 
solvent-solvent pairwise interaction energies that end up causing a significant 
round-off error. There will still be potentially round-off error due to the 
solute-solvent terms, but it should be smaller than the current one.

Original comment by julienm...@googlemail.com on 15 Jul 2013 at 10:51