[gmx-users] Metal organic framework energy minimization

[Jordan Chapman]

GROMACS users,

I am attempting pre-processing of a metal organic framework (ZIF-8) in
water, but I am running into trouble with energy minimization. I've added
both bonded and nonbonded forcefield parameters for ZIF-8 into the AMBER03
forcefield library and solvated the system with the TIP3P water model using
"gmx solvate". The steepest descent algorithm converges to the options I
have set in my minimization mdp file, and an analysis of potential energy
over the simulation time using "gmx energy" follows an expected trend.
However the simulation returns a warning that at several time steps in the
minimization one or more water molecules could not be settled.

Depending on the size of the system I defined using "gmx editconf" I find
that the ZIF-8 molecule by the end of the simulation is either crumpled
when the box size is too large or stretched to the dimensions of the box
when the box size is too small. For reference, the size of a unit cell of
ZIF-8 is approximately 1.7 nm x 1.7 nm x 1.7 nm.

I've analyzed the range of box dimensions that return either crumpled or
stretched ZIF-8 conformations. At a diameter of 3.09 nm, energy
minimization returns a stretched ZIF-8 structure. At a diameter of 3.0925
nm, energy minimization returns a crumpled ZIF-8 structure. See below for
the mdp options I have set for energy minimization. Energy minimization of
ZIF-8 at the same mdp options was successful for ZIF-8 in vacuum.

; Energy minimization of ZIF-8 unit cell in TIP3P
include                  =
define                   = -DPOSRES

; RUN CONTROL PARAMETERS
integrator               = steep
dt                       = 0.001
nsteps                   = 100000
comm-mode                = Linear
nstcomm                  = 20
comm-grps                =


; ENERGY MINIMIZATION OPTIONS
emtol                    = 10
emstep                   = 0.01


; OUTPUT CONTROL OPTIONS
nstxout                  = 200
nstvout                  = 200
nstfout                  = 200
nstlog                   = 200
nstcalcenergy            = 200
nstenergy                = 200
nstxout-compressed       = 0
energygrps               =


; NEIGHBORSEARCHING PARAMETERS
cutoff-scheme            = Verlet
nstlist                  = 40
ns_type                  = grid
pbc                      = xyz
periodic-molecules       = yes
verlet-buffer-tolerance  = 0.0001
rlist                    = 1.4


; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype              = PME
coulomb-modifier         =
rcoulomb                 = 1.4
vdwtype                  = Cut-off
vdw-modifier             =
rvdw                     = 1.4
DispCorr                 = EnerPres
; EWALD/PME/PPPM parameters
fourierspacing           = 0.12
pme_order                = 4
ewald_rtol               = 1e-05


; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel                  = no
gen-temp                 = 300
gen-seed                 = 3532


; OPTIONS FOR BONDS
constraints              = none

Any advice or insight into what could be causing this or a recommendation
for appropriate pre-processing is greatly appreciated!

Best,
Jordan Chapman
Graduate Research Assistant
Department of Chemical and Biomedical Enginering
West Virginia Univesity