Error in implementing LJPME in gromacs-2020.1

User discussions
1

GROMACS version: 2020.1
GROMACS modification: No
Dear Gromacs Users,

I am just following up on my previous email which unfortunately didn’t get any response. Now everything is moved to forum. I would highly appreciate if any of you can comment on this if you have implemented long-range vdw interaction in any of the version of gromacs.
I am trying to implement PME for vdw interaction for an umbrella sampling simulation containing nanoparticles. There are two nano-particles; one is frozen in its place and the other is allowed to move along z-direction. Below is my pulling .mdp file (at the end of the message). When I run the same system with vdw=cut-off ; I see no issue running in any version of the gromacs.

grompp warning:
unlike previous versions(2019 and before), grompp doesn’t process the topology file unless i exclude the frozen(all-directions) nanoparticle from the comm-grps. Here is the warning after processing the topology.

WARNING 1 [file pull.mdp]:
There are 249 atoms that are frozen along less then 3 dimensions and part
of COMM removal group(s), due to limitations in the code these still
contribute to the mass of the COM along frozen dimensions and therefore
the COMM correction will be too small.

mdrun:
The energy minimization seems to run OK with LJPME which I had issues with previous version of gromacs
For the production run, after running few hundred steps the simulation crashes with following error: Decreasing the step size to 0.5femto second doesn’t lead to anything different.

500000 steps, 500.0 ps.
step 400, will finish Tue May 12 04:14:34 2020
step 482: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates

Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.

mpirun noticed that process rank 0 with PID 0 on node scc-k04 exited on signal 11 (Segmentation fault).

pull.mdp file

integrator = md
tinit = 0
dt = 0.001
nsteps = 500000
comm-mode = Linear
nstcomm = 100
comm-grps = !NP1
emtol = 10
emstep = 0.01
nstlog = 10000
nstcalcenergy = 1000
nstenergy = 1000
nstxout-compressed = 50000
compressed-x-precision = 50000
; This selects the subset of atoms for the compressed
; trajectory file. You can select multiple groups. By
; default, all atoms will be written.
compressed-x-grps =
; Selection of energy groups
energygrps =

; NEIGHBORSEARCHING PARAMETERS
cutoff-scheme = Verlet
nstlist = 10
pbc = xyz
periodic-molecules = no
verlet-buffer-tolerance = 0.005
rlist = 1.2
coulombtype = PME
coulomb-modifier = Potential-shift-Verlet
rcoulomb-switch = 0
rcoulomb = 1.2
epsilon-r = 1
epsilon-rf = 0
vdw-type = PME
vdw-modifier = Potential-shift-Verlet
rvdw-switch = 1
rvdw = 1.2
DispCorr = No
table-extension = 1
energygrp-table =
fourierspacing = 0.12
pme-order = 4
ewald-rtol = 1e-05
ewald-rtol-lj = 0.001
lj-pme-comb-rule = Geometric
ewald-geometry = 3d
epsilon-surface = 0

tcoupl = Nose-hoover; Berendsen
tc-grps = SOL non-water
ref-t = 300 300
tau-t = 1 1
pcoupl = Parrinello-Rahman; Berendsen
pcoupltype = Isotropic
nstpcouple = -1
tau-p = 4
compressibility = 4.5e-05
ref-p = 1
refcoord-scaling = No
gen-vel = no;
gen-temp = 300
gen-seed = -1
constraints = none;
constraint-algorithm = Lincs
continuation = no
Shake-SOR = no
shake-tol = 0.0001
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
morse = no
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp-excl =
; COM PULLING
pull = yes
; Cylinder radius for dynamic reaction force groups (nm)
pull-cylinder-r = 1.5
pull-constr-tol = 1e-06
pull-print-com = yes
pull-print-ref-value = yes
pull-print-components = yes
pull-nstxout = 50
pull-nstfout = 50
; Number of pull groups
pull_ngroups = 2
; Number of pull coordinates
pull_ncoords = 1
; Group and coordinate parameters
pull_group1_name = NP1
pull-group1-weights =
pull-group1-pbcatom = 0
pull_group2_name = NP2
pull-group2-weights =
pull-group2-pbcatom = 0
pull_coord1_type = umbrella
pull-coord1-potential-provider =
pull_coord1_geometry = distance
pull_coord1_groups = 1 2
pull_coord1_dim = N N Y
pull-coord1-origin = 0.0 0.0 0.0
pull-coord1-vec = 0.0 0.0 0.0
pull_coord1_start = yes
pull-coord1-init = 0
pull_coord1_rate = 0
pull_coord1_k = 1000
pull-coord1-kB = 1000

; Non-equilibrium MD stuff
acc-grps =
accelerate =
freezegrps = NP1 NP2
freezedim = Y Y Y Y Y N
cos-acceleration = 0

Thank you,
Udaya Dahal

2

I added the warning about the COMM removal to grompp in version 2020.
I don’t know what you want to do exactly, but you need to think what should happen with the COMM removal of your nano-particles. I guess both should not be part of COMM removal groups, which would solve your issues. You can also use your old settings, ignore the warning and get the 2019 behavior.

3

Note that your topic is incorrect. Your error has nothing to do with LJPME, but rather with COMM removal (combined with pulling).

4

Thank you for the comment. The reason I put LJPME in title is because when I switch to vdw=cut-off; I don’t see the error for pulling simulation. I actually tried removing both particle from COMM removal groups but I am unable to run the equilibration for pulling. It crashes almost immediately. As mentioned this behavior is not observed when using vdw = cut-off for the same setup. That’s why I thought the error is stemming out from the LJPME.

If I run the two particles without ions and water – everything runs fine with LJPME
If I run the water and ions without nanoparticle – everything runs fine with LJPME
When I mix nanoparticle with ions or water – the system crashes almost immediately after few steps only when using LJPME but works fine with cut-off.

I am not sure why that’s the case.

5

Indeed that is what you wrote. I almost can’t believe that the non-bonded treatment would affect a COMM option. Can you attach all input files needed for running grompp, so I can check what is going on here?

6

Please find the attached file. https://drive.google.com/file/d/1nJi3T6-x8S-r15c9j2WWmdciwUnq_u-r/view?usp=sharing

I missed the index file please use this new link.

7

When I change vdwtype from PME to cut-off, I don’t see any difference in warnings, as I expected. So this is not an LJ-PME issue.

I don’t understand why you freeze NP1 along x,y and z and NP2 only along x and y. Then atoms within NP2 can move with respect to each other along z. Can you not simply remove all freezing (and have the whole system as COMM group)?

8

Thank you for the information. I am not worried about the initial warning but the md simulation. I actually checked with non-frozen nanoparticles but the final outcome is the same. If you try “mdrun” ; may be for short period of time <100ps ; you will quickly know that implementing LJPME will crash the simulation immediately. But if you switch to cut-off; it just works fine.

If you think it is not due to LJPME; what do you suspect for the source of error?

I need to freeze the atoms eventually because I am running umbrella sampling simulations. The reason I froze along XY is to let the particle drift along Z-direction. I wouldn’t worry about individual atoms motion along z-direction because the intra-NP bonds are stiff and the NP2 will drift eventually as a whole nanoparticle. I didn’t have any issue calculating PMF using LJ cut-off for the same system.

Thanks,

9

Apologies for the late answer.

I didn’t realize the issue is not related to the COMM removal at all.
My guess is that the issue is that the LJ parameter for your AU atomtype are very far off from the combination rules, whereas LJ-PME PME assumes combination rules, so the long-range attraction between AU and other atoms is an order of magnitude too large. Anyhow, you do not want to use LJ-PME when you interactions are so far off from combination rules.

10

Thank you for the reply. I will look into if changing the LJ parameter actually fixes the error. Then we will know for sure , it’s the AU parameter that’s what causing the error.

11

Dear Berk,

I tested the simulation by turning off all the non-bonded parameters defined in pair types so that all the parameters follow the combination rules (using “gen pairs = yes”) but still I got the similar error. Does it mean that I won’t be able to use LJ-PME at all for any system that has similar LJ potential like AU? If so, what is the source of the error?

Thanks,