[gmx-users] How Coulomb energy terms are calculated?

[Mark Abraham]

On 21/06/2012 9:11 PM, vincenzo moreschini wrote:
> Dear users,
>
> I'm trying to understand the way gmx calculates non-bonded interactions
> and I set up a toy system composed of two water molecules with initial
> coordinates:
>
>      1SOL     OW    1   0.3758265   0.5507358   0.8927467
>      1SOL    HW1    2   0.3234106   0.4695494   0.8670300
>      1SOL    HW2    3   0.3936553   0.6061094   0.8114085
>      2SOL     OW    4   0.0000000   0.0000000   0.0000000
>      2SOL    HW1    5   0.1000000   0.0000000   0.0000000
>      2SOL    HW2    6  -0.0333313   0.0942816   0.0000000
>
> and zero velocities. Modifying spce.itp I put HW charges=0 so that only
> Oxygens matter. Calculating distance between
> atom 1 and 4 and then the Coulomb term as f*q*q/r I get 89.5803168839 kj/mol
> with f taken from the manual (138.935), q=-0.8476 and r=1.1142495905.
>
> The I set nstlist=rcoulomb=rvdw=rlist=0.0 pbc=no and nstype simple. Doing 1
> step of sim I get:
>
>     Energies (kJ/mol)
>             Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>      3.65337e-10    5.85205e-05   -1.36690e-03    8.95803e+01    8.95790e+01
>      Kinetic En.   Total Energy    Temperature Pressure (bar)
>      5.31575e-11    8.95790e+01    8.52446e-10    0.00000e+00
>
> Playing with constraints does changes nothing, as I think it should.

That depends whether constraints are already present, satisfied, and/or 
whether unconstrained_start = no. The way to just compute a no-nonsense 
single point is with mdrun -s yourmethod -rerun yourconfiguration.

> However, when I insert the correct charges (0.4238) for hydrogen atoms,
> using the same set up I get
>   
>   Energies (kJ/mol)
>             Bond          Angle        LJ (SR)   Coulomb (SR)      Potential
>      3.65337e-10    5.85205e-05   -1.36690e-03    2.76599e-01    2.75290e-01
>      Kinetic En.   Total Energy    Temperature Pressure (bar)
>      3.13252e-12    2.75290e-01    5.02337e-11    0.00000e+00
>
> But summing the values of the 9 Coulomb interactions present in the systems
> (6 atoms -> 6*(6-1)/2 -bonded terms)
> by hand I get about twice i.e. 0.479870592 kjmol-1.
> Given that I made no errors (I shouldn't since I pasted formulas in a xls
> sheet and checked but only death and taxes are sure) what am I missing?
>

Not sure. You can check the components of your calculation by breaking 
each atom into its own charge group in the .top and defining a uniquely 
named energy group for each charge group in the .mdp. Then the .edr 
output file will have 9 Coulomb (SR) interaction terms.

Mark