[gmx-developers] tpic energies - COUL_SR

[Rui Fartaria]

Hi all,

I am a new user/developer of GROMACS. I am using GROMACS 4.6.3
My goal is to extend the package GromPy, created by René Pool but I am 
still new to the code.

In the meanwhile I wanted to calculate the excess chemical potential 
using TPI/TPIC. TPIC is preferable as my systems tend to be highly 
anisotropic and therefore TPIC allows me to choose where to sample.

I noticed that some energies for the insertion of molecules were 
extremely low. To check if everything was working ok, I prepared the 
following test:

  * I have a configuration without the inserted molecule and the same
    configuration with the molecule inserted.
  * The molecule in question is monoethanolamine (MEA).
  * MEA topology was modified with nrexcl=5 so that intra molecular
    interactions were not an issue for comparing insertion energies.
  * I run 1 time step in each and then run "echo '1 2 3 4 5 6 0' |
    g_energy" to get the different terms for potential energy.
  * When calculating the differences between the "inserted" and "non
    inserted" configurations I corrected for the potential energy due to
    angles and dihedrals.
  * I hacked the tpi.c source code and forced it to insert the same
    configuration as in the tests described above.

I prepared at tgz with the tests and the code for the hacked tpi.c. I 
can be downloaded from:
https://www.dropbox.com/s/0lnjancduzi2frf/example_to_send.tgz

As I have tested it, if you compile GROMACS 4.6.3 with the my modified 
version of tpi.c and run "bash run_examples.bash" in the example folder 
you should get the file result.txt as shown below (really hope it works 
for you too):


You can see from the results below that the energies between the 
standard calculation for potential energy and the TPIC potential energy 
do not match. I searched the archives and the only reference to a 
problem like this was due to problems in the reciprocal sum in PME. In 
my case the problem is in the Coulomb real space sum. This is affected 
by a cutoff and I could not locate the exact place in the code where 
these calculations are done.

The LJ sums match, apart from long range corrections. I also find it 
strange that the reciprocal part of the coulomb term is 0.0 in TPIC.

I would be really thankful if someone could help me understand where the 
problem is. Hopefully you already know.

Kind regards,

Rui Fartaria



   Angle      RB            LJ    Coul(SR)   Coul(RECIP)  Potential
  15611.243  -9512.769  -9837.916 -13497.961   -174.225 -17411.627
  15649.481  -9521.968  -9842.127 -13493.568   -174.517 -17382.699
     38.238     -9.199     -4.211      4.393     -0.292     28.928
Potential - Angle - RB =     -0.111


******************************************
RESULTS FROM TPI
******************************************
Using GMX_CAVITY_CENTER cavity




WARNING: USING CUSTOM CONF FOR INSERTION


epot = -65.895645
embU = 21686035981.671570
F_COUL_RECIP = 0.000000
F_COUL_SR = -61.993370
F_LJ = -3.902278
beta = 0.361176
BEGIN INSERTED CONFIGURATION
Frame time: 0.000000
14.00 2.058448 1.753919 2.382128
  1.00 2.116501 1.707533 2.312707
12.00 1.961425 1.845027 2.319632
  1.00 2.037138 1.699207 2.464866
  1.00 2.021572 1.899345 2.245694
12.00 1.895999 1.942494 2.418418
  1.00 1.896307 1.780891 2.259156
  1.00 1.853303 1.890639 2.505562
  1.00 1.961367 2.024441 2.449268
16.00 1.780418 1.997807 2.356455
  1.00 1.751046 2.084993 2.389755
END INSERTED CONFIGURATION

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20140221/49b31e21/attachment.html>