[gmx-users] Periodic Molecule's Free Energy Calculation Error

[Jason Zhu]

Hello Gromacs Community,

I am trying to calculate the solvation free energy of a hBN sheet following
Justin Lemkul and Alchemistry's tutorials.
Since the hBN sheet is infinitely large, I turned the periodic molecules
flag on.
This runs all fine on one core, but when I try to run NVT in parallel (e.g.
4 ranks), the job would throw the following error:

A list of missing interactions:
        LJC Pairs NB of 890400 missing 338688

-------------------------------------------------------
Program gmx mdrun, VERSION 5.1.4
Source code file:
/gpfs/runtime/opt/gromacs/5.1.4/src/gromacs-5.1.4/src/gromacs/domdec/domdec_topology.cpp,
line: 242

Software inconsistency error:
Some interactions seem to be assigned multiple times
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

Halting parallel program gmx mdrun on rank 1 out of 4
In: PMI_Abort(1, application called MPI_Abort(MPI_COMM_WORLD, 1) - process
1)

-------------------------------------------------------
Program gmx mdrun, VERSION 5.1.4
Source code file:
/gpfs/runtime/opt/gromacs/5.1.4/src/gromacs-5.1.4/src/gromacs/domdec/domdec_topology.cpp,
line: 242

Software inconsistency error:
Some interactions seem to be assigned multiple times
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

It seems like a domain decomposition error. My first thought was that the
system "explode".
However, when I check my topology and pbc condition carefully, there is no
sign of anything wrong.
I also tried NPT & PROD MD. The same error when I ran on multiple MPI
threads.

My question is: why the system could run fine on one MPI, but not if I
increased the number of MPI threads?

Any help on this issue will be really appreciated.


Here below is my .mdp file:

; RUN CONTROL—NVT
;——————————————————————————
define       = -DPOSRES_HBN
integrator   = sd            ; stochastic leap-frog integrator
nsteps       = 5000          ; 2 * 5,000 fs = 10 ps
dt           = 0.002         ; 2 fs
comm-mode    = Linear        ; remove center of mass translation
nstcomm      = 100           ; frequency for center of mass motion removal

;——————————————————————————
; OUTPUT CONTROL
;——————————————————————————
nstxout                = 0          ; don't save coordinates to .trr
nstvout                = 0          ; don't save velocities to .trr
nstfout                = 0          ; don't save forces to .trr
nstxout-compressed     = 5000       ; xtc compressed trajectory output
every 5000 steps
compressed-x-precision = 1000       ; precision with which to write to the
compressed trajectory file
nstlog                 = 5000       ; update log file every 10 ps
nstenergy              = 5000       ; save energies every 10 ps
nstcalcenergy          = 100        ; calculate energies every 100 steps

;——————————————————————————
; BONDS
;——————————————————————————
constraint_algorithm   = lincs      ; holonomic constraints
constraints            = all-bonds  ; hydrogens only are constrained
lincs_iter             = 1          ; accuracy of LINCS (1 is default)
lincs_order            = 4          ; also related to accuracy (4 is
default)
lincs-warnangle        = 30         ; maximum angle that a bond can rotate
before LINCS will complain (30 is default)
continuation           = no         ; formerly known as
'unconstrained-start' - useful for exact continuations and reruns

;——————————————————————————
; NEIGHBOR SEARCHING
;——————————————————————————
cutoff-scheme   = Verlet
ns-type         = grid   ; search neighboring grid cells
nstlist         = 10     ; 20 fs (default is 10)
rlist           = 1.2    ; short-range neighborlist cutoff (in nm)
pbc             = xyz    ; 3D PBC
; PBC: grp is infinite
periodic-molecules = yes

;——————————————————————————
; ELECTROSTATICS
;——————————————————————————
coulombtype      = PME      ; Particle Mesh Ewald for long-range
electrostatics
rcoulomb         = 1.2      ; short-range electrostatic cutoff (in nm)
ewald_geometry   = 3d       ; Ewald sum is performed in all three dimensions
pme-order        = 4        ; interpolation order for PME (default is 4)
fourierspacing   = 0.16     ; grid spacing for FFT
ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct
potential at rcoulomb

;——————————————————————————
; VDW
;——————————————————————————
vdw-type                = PME
rvdw                    = 1.2
vdw-modifier            = Potential-Shift
ewald-rtol-lj           = 1e-3
lj-pme-comb-rule        = Geometric
DispCorr                = EnerPres

;——————————————————————————
; TEMPERATURE & PRESSURE COUPL
;——————————————————————————
tc_grps    =  System
tau_t      =  0.1
ref_t      =  300
pcoupl     =  no

;——————————————————————————
; VELOCITY GENERATION
;——————————————————————————
gen_vel      = yes       ; Velocity generation is on (if gen_vel is 'yes',
continuation should be 'no')
gen_seed     = -1       ; Use random seed
gen_temp     = 300

;——————————————————————————
; FREE ENERGY CALCULATIONS
;——————————————————————————
free-energy              = yes
couple-moltype           = BNT
couple-lambda0           = vdw-q
couple-lambda1           = none
couple-intramol          = no
separate-dhdl-file       = yes
sc-alpha                 = 0.5
sc-power                 = 1
sc-sigma = 0.3
init-lambda-state        = 0
coul-lambdas             = 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0
1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas              = 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5
0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
nstdhdl                  = 100
calc-lambda-neighbors    = -1


Best,
Xuliang