[gmx-users] User defined tabulated potential for a particle with only repulsive interactions, in water

[Braden Kelly]

I would like to simulate a particle that has only a repulsive (no attraction) terms in its potential energy interaction with all other particles. I am simulating this particle in tip3p water.


I am trying to use the tabulated potential option since I know its potential form, and force as well. In python lingo the potential and force are


def cavityPot(r):
    pot = potA * np.exp(-r/potB + Lambda )
    return pot
def cavityForce(r):
    force = - potA / (-potB) * np.exp( -r/potB + Lambda )
    return force

potA  = 1255 # kj/mol

potB = 0.1 # nm

Lambda = -2

r = distance between two atoms, nm


The simulation does run, but is not correct and terminates prematurely, but without raising an error. I say this because the energy minimization converges in around 10 steps which does not happen even for a simulation with only water in it, and the NVT equilibration starts without error, runs, but quits without actually running i.e., the atom coordinates are the same after as before.


Unfortunately there is no error during running, so I must show my files and explain my procedure, and hope a flaw is picked out of these. One note is that gromacs complains that the force it calculates is 52% different from my tabulated force. The equation is analytical, the derivative of potential energy is straightforward. When I test it numerically as well, I still get results saying my force/energy are correct. Some values from a table are shown near the bottom.


If you are still reading, I will go into more depth on my topology,mdp file used, and table generation: here is the top portion of the .top file


---------------------------------------------------------------
[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
  1             2               yes             0.5       0.8333

[ atomtypes ]
; name   mass      charge     ptype  sigma(nm)     epsilon(kJ/mol)
; tip3p atoms
OW        OW    15.99940   -0.8340      A   0.315061     0.6364000
HW        HW    1.008000   0.4170       A   0.0000       0.0000000


;name  bond_type    mass    charge   ptype          sigma      epsilon
 CAV       CAV          0.00000  0.00000   D     0.000000   0.000000 ;

[ nonbond_params ] ; for tabulated potential A and C parameters need to be set to 1
;  i  j    func    V(c6)    W(c12)
  CAV OW    1       1          1
  CAV HW    1       1          1

;=====================================================================
;                              CAVITY - DUMMY
;=====================================================================
; Include CAVITY topology
#include "CAVITY_GMX.itp"

;=====================================================================
;                             WATER tip3p
;=====================================================================
[ moleculetype ]
; molname    nrexcl
SOL    3

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1  OW    1    SOL      OW      1      -0.8340        15.99940
     2  HW    1    SOL      HW1      2       0.4170        1.0080
     3  HW    1    SOL      HW2      3       0.4170        1.0080

(Removed rest of tip3p info to save space... it involves settles etc)
;=====================================================================
[ system ]
CAVITY in Water

[ molecules ]
CAVITY 1
SOL 500
--------------------------------------------------------------------------------------------------------
The corresponding #include "CAVITY_GMX.itp" is
--------------------------------------------------------------------------------------------------------
[ moleculetype ]
; Name            nrexcl
CAVITY  3

[ atoms ] ;
;   nr  type  resi  res  atom  cgnr     charge      mass       ; qtot   bond_type
     1   CAV     1   CAV   CAV    1    0.0000     0.00000 ;
--------------------------------------------------------------------------------------------------------
The mdp parameters I use in addition to a typical minimization/NVT run is:

--------------------------------------------------------------------------------------------------------

cutoff-scheme             = group

; USER DEFINED POTENTIAL TABLES
energygrps                    = CAV OW HW1 HW2
energygrp_table          = CAV OW  CAV HW1 CAV HW2
table-extension            = 1           ; Extension of the potential lookup tables beyond the cut-off
; Electrostatics
coulombtype                = PME-User
; method for vdW (LJ)
vdwtype                        = user
--------------------------------------------------------------------------------------------------------

For the NVT run i use deffnm = NVTEq, so my tables are respectively called

NVTEq.xvg, NVTEq_CAV_CAV.xvg, NVTEq_CAV_OW.xvg,NVTEq_CAV_HW1.xvg,NVTEq_CAV_HW2.xvg


All tables are the same since the cavity does not care what particle it is interacting with and looks like


---------------------------------------------------------------------------------------------------------

0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0345582 0.3455824
0.0020000 0.0000000 0.0000000 0.0000000 0.0000000 0.0338739 0.3387394
0.0040000 0.0000000 0.0000000 0.0000000 0.0000000 0.0332032 0.3320320
0.0060000 0.0000000 0.0000000 0.0000000 0.0000000 0.0325457 0.3254573
0.0080000 0.0000000 0.0000000 0.0000000 0.0000000 0.0319013 0.3190128
0.0100000 0.0000000 0.0000000 0.0000000 0.0000000 0.0312696 0.3126959
0.0120000 0.0000000 0.0000000 0.0000000 0.0000000 0.0306504 0.3065041
(and so on to the r = 2.5)


I make an index file with groups [ system ], OW, HW1,HW2, CAV.


I call gromacs using

gmx grompp -c aqueous.g96 -p run.top -f index.ndx -maxwarn 2 -o NVTEq.tpr
gmx mdrun -deffnm NVTEq -table NVTEq.xvg -c aqueous.g96


To recap, the energy minimization runs but stops way to soon. The NVT runs, but doesn't actually seem to do anything. The coordinates are not updated and it stops right away. No errors are generated.


I appreciate any help,


Braden