[gmx-users] grompp is using a very large amount of memory on a modestly-sized system

[Sean Marks]

Hi, everyone,

I am running into an issue where grompp is using a tremendous amount of
memory and crashing, even though my system is not especially large (63976
atoms).

I am using GROMACS 2016.3.

My system consists of liquid water (7,930 molecules) next to a block of ice
(8,094 molecules). The ice oxygens are restrained to their lattice position
with a harmonic potential with strength k = 4,000 kJ/mol/nm^2. I am using
the TIP4P/Ice model, which is a rigid 4-site model with a negative partial
charge located on a virtual site rather than the oxygen.

My goal is to systematically reduce the electrostatic interactions between
the water molecules and the position-restrained ice, while leaving
water-water and ice-ice interactions unaffected.

To accomplish this, I am modeling all of the ice molecules using a single
moleculetype so that I can take advantages of GROMACS' FEP features to
selectively scale interactions. I explicitly specify all constraints and
exclusions in the topology file. This moleculetype contains one virtual
site, 3 constraints, and 4 exclusions per "residue" (ice molecule).

When I run grompp, I receive the following error, which I think means that
a huge block of memory (~9 GB) was requested but could not be allocated:

=====
Command line:
  gmx grompp -f npt.mdp -c md.gro -p topol.top -n index.ndx -r
initconf_packmol.gro -o input.tpr -maxwarn 2 -pp processed.top

...

Generated 21 of the 21 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 21 of the 21 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'ICE'
turning H bonds into constraints...
Excluding 3 bonded neighbours molecule type 'SOL'
turning H bonds into constraints...
Coupling 1 copies of molecule type 'ICE'
Setting gen_seed to 1021640799
Velocities were taken from a Maxwell distribution at 273 K
Cleaning up constraints and constant bonded interactions with virtual sites
Removing all charge groups because cutoff-scheme=Verlet

-------------------------------------------------------
Program:     gmx grompp, version 2016.3
Source file: src/gromacs/utility/smalloc.cpp (line 226)

Fatal error:
Not enough memory. Failed to realloc -8589934588 bytes for il->iatoms,
il->iatoms=25e55010
(called from file
/home/semarks/source/gromacs/2016.3/icc/serial/gromacs-2016.3/src/gromacs/
gmxpreprocess/convparm.cpp,
line 565)

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
=======

In the hope that it helps with diagnosing the problem, here is my mdp file:

=======
; Water and ice

; RUN CONTROL PARAMETERS

; Molecular dynamics
integrator               = md       ; Leapfrog
dt                       = 0.002    ; 2 fs - time step [ps]
nsteps                   = 5000000  ; 10 ns - run time [steps]

; Center of mass motion (COMM) removal
comm-mode                = none     ; Linear
nstcomm                  = 10            ; Removal frequency (>=
nstcalcenergy) [steps]
comm-grps                = System   ; Groups for COMM removal (blank =>
whole system)

; Initial velocities
gen_vel                  = yes    ; Generate velocities using Boltzmann
distribution
gen_temp                 = 273    ; Temperature for Boltzmann distribution
[K]
gen_seed                 = -1     ; Seed for RNG is the job ID

; OUTPUT CONTROL OPTIONS
nstcalcenergy            = 1  ; Frequency of energy calculation [steps]
; Output frequency for coords (x), velocities (v) and forces (f)
;     trr
nstxout                  = 0       ; never - print coordinates [steps]
nstvout                  = 0       ; never - print velocities [steps]
nstfout                  = 0       ; never - print forces [steps]
; Log end edr files
nstlog                   = 2500    ; 5 ps - print energies to log file
[steps]
nstenergy                = 2500    ; 5 ps - print energies to energy file
[steps]
;     xtc
nstxout-compressed       = 2500    ; 5 ps - print coordinates to xtc file
[steps]
compressed-x-precision   = 1000    ; Number of zeros = number of places
after decimal point
compressed-x-grps        = System  ; Groups to write to xtc file

; BOUNDARY CONDITIONS
pbc                      = xyz
periodic_molecules       = no      ; Rigid graphene has no intramolecule
interactions

; NEIGHBOR LIST
cutoff-scheme            = Verlet
nstlist                  = 10      ; Neighbor list update frequency [steps]
ns-type                  = grid    ; More efficient than simple search
verlet-buffer-tolerance  = 0.005   ; [kJ/mol/ps] (-1.0 --> use rlist)
rlist                    = 1.0     ; Cutoff distance for short-range
neighbor list [nm]

; VAN DER WAALS
vdwtype                  = Cut-off
vdw-modifier             = Potential-shift-Verlet
rvdw                     = 1.0        ; Radius for vdW cutoff [nm]
DispCorr                 = no

; ELECTROSTATICS
coulombtype              = PME        ; Fast, smooth, particle mesh Ewald
(PME)
coulomb-modifier         = Potential-shift-Verlet
rcoulomb                 = 1.0        ; Short-range Coulomb cutoff [nm]
epsilon-r                = 1      ; Relative dielectric constant

; Ewald/PME/PPPM parameters
fourierspacing           = 0.12   ; FFT grid spacing [nm]
pme_order                = 4      ; Cubic interpolation of PME grid
ewald-rtol               = 1e-05  ;
ewald-geometry           = 3d     ; 3-D charge treatment
epsilon-surface          = 0      ; No dipole correction

; THERMOSTAT
Tcoupl                   = V-rescale     ; stochastic velocity-rescale
nstTcouple               = 1
tc-grps                  = ICE    SOL    ; Groups which are controlled
together
tau_T                    = 0.5    0.5    ; Time constant for coupling [ps]
ref_T                    = 273.0  273.0  ; Temperature setpoint [K]

; BAROSTAT
Pcoupl                   = Parrinello-Rahman      ; (options:
no/Berendsen/Parrinello-Rahman)
Pcoupltype               = anisotropic
nstPcouple               = 1
tau-P                    = 4.0         ; Time constant for coupling [ps]
ref-P                    = 1.0   1.0     1.0     0.0   0.0   0.0  ;
Pressure setpoint [bar]
compressibility          = 0.0   4.8e-5  0.0     0.0   0.0   0.0  ;
Approximate compressibility [bar^-1]
refcoord-scaling         = no

; CONSTRAINTS
constraints              = H-bonds ; Convert bonds with H to constraints
constraint_algorithm     = LINCS   ; Best choice if you don't have angle
constraints
lincs-order              = 4       ; Highest order in expansion of
constraint coupling matrix
lincs-iter               = 1       ; More iterations --> higher accuracy
continuation             = no      ;

; Use GROMACS' FEP features to scale ice-water electrostatics, while leaving
; ice-ice electrostatics unaffected.
; - This setup: lambda = scaling factor for interactions (ranges from 0 to
1)
free-energy     = yes   ; turns on FEP
couple-moltype  = ICE   ; group for coupling
couple-intramol = no    ; don't adjust ice-ice interactions

couple-lambda0 = vdw     ; ice-water electrostatics are off at lambda=0
couple-lambda1 = vdw-q   ; all interactions are on at lambda=1
init-lambda    = 0.5     ; value of lambda (stays fixed with this setup)
delta-lambda   = 0.0     ; don't change lambda

nstdhdl          = 0    ; no extra output needed
dhdl-derivatives = no   ; no extra output needed
=======

Thank you for your time.

Sincerely,
Sean

-- 
Sean M. Marks
PhD Candidate
Dept. of Chemical and Biomolecular Engineering
University of Pennsylvania