GROMACS version: 2018.1
GROMACS modification: No
Greetings,
I am having some issues implementing user-defined potentials in my GROMACS version. I should preface that I am using the group neighbor searching scheme (even though it is deprecated) due to lack of user-potential compatibility with the Verlet scheme.
I am trying to implement a custom VdW interaction and standard coulomb interaction. The coulomb interactions are fairly long-range so I would like to use the PME (or PME-user) coulomb-type interactions; however, the results coming out of gmx energy suggest that something is wrong with the simulation or that the properties are being miscalculated (for both NVT and NPT equilibrations).
The RDF and diffusion coefficients very nearly match the expected shapes and values, but the temperature and density are poor. The reference temperature is set to 1200 K with the velocity rescaling thermostat but the system equilibrates to about - 240,000 K. The density equilibrates to about twice the expected value (~10 g/mL), but my own rough calculation of the density based on the box vectors of the last frame give the expected density (~4 g/mL). Additionally, visual inspection of the NPT trajectory with VMD shows that the box vectors do not change significantly and that the density should be about 4 g/mL.
Perhaps the issue with temperature is similar to Bug #2286 in that the energies for the energy summation are wrong because outputting the temperature with gmx traj, which (I think) calculates temperature solely on atom velocity, shows that the temperature is equilibrated at the 1200 K reference temperature. NOTE: I get this issue for both the PME and PME-user coulomb-types.
As for density, I have no leads on what might produce this problem. Much like the temperature, using gmx traj to output the box vectors shows that the density is stable near the expected 4 g/mL value.
I would appreciate help on figuring out why the temperature and densities coming out of gmx energy seem wrong. Additionally, are the incorrect temperature values shown by gmx energy affecting the thermostat during the simulation, or are the temperature values only incorrectly calculated after the simulation runs?
I cannot attach my .mdp file (new user) so I’ll put some of the relevant lines below:
;
Nonbonded Settings
cutoff-scheme = group
ns_type = grid
nstlist = 10 ; every 0.01 ps
rlist = 0.80 ; [nm];
Van der Waals
vdw-type = user
rvdw = 0.8;
Electrostatics
coulombtype = pme
rcoulomb = 0.8;
Temperature coupling is ON
tcoupl = v-rescale
tc_grps = System
tau_t = 0.1
ref_t = 1200;
Pressure coupling is ON
pcoupl = parrinello-rahman
pcoupltype = isotropic
tau_p = 1.0
ref_p = 1.0 ; reference pressure [bar]
compressibility = 3.704e-6 ; isothermal compressibility [bar^-1] of UCl3 ~ 1200 K
If there is another way I can help to recreate the issue, I would be glad to help.
Thank you for the help.