[gmx-users] Energy minimization failure

[Jeff Woodford]

Greetings,
I am relatively new to MD, and I am attempting to simulate a metal-organic framework but I don't seem to be even able to get past the energy minimization phase. I could appreciate any insight into what I might be doing wrong. Here is a brief summary of what I have done:
 
- constructed a PDB file from the crystal structure coordinates of a model compound, capping the organic linker with phenyl and capping the metal atoms with water and hydroxide ions to yield a full coordination sphere and to yield a net charge of zero
- constructed a custom force field using published parameters
- run pdb2gmx to construct the .gro and .top files
- run editconf to construct a box around the MOF, 5 nm in each dimension
- run genbox to fill the box with water
- run grompp to prepare the box for addition of ions
- run genion to replace 20 solvent molecules with 10 Na+ and 10 Cl- ions
- manually edit the topology file (per the suggestion here on the listserv of a few days ago) to include, in the [ pairs ] directive, all pair interactions for purposes of electrostatic computation, but with parameters of "0 0" for the van der Waals parameters, so as to have the van der Waals interactions remain the usual case of neglected for 1-2 and 1-3 interactions
- create an index file containing a group that has the 60 atoms constituting the terminal water and hydroxide ion atoms, so as to keep their coordinates frozen during the minimization
- run grompp again to prepare the box for energy minimization (below is the minim.mdp file used for this)
- run mdrun to perform the energy minimization
 
The energy minimization does not yield the maximum force falling below the desired 1000, but stays as high as 10^6. Visualization of the output .gro file yields that many atoms within the core of the MOF have moved into "strange" positions, i.e., the formation of seeming O-O bonds between adjacent carboxylates when instead both O atoms should be coordinated to the metal atom.
 
The calculations yielded a great deal of data, and so as not to clog everybody's inboxes, I only included the one that seemed most relevant (minim.mdp). If there is something else required, I'll happily share it.
 
I am not sure what's going on and I could appreciate any insight into what might be causing this and how to fix it.
 
Thanks in advance ,
Jeff Woodford
Assistant Professor of Chemistry
Missouri Western State University
 
minim.mdp:
 
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep  ; Algorithm (steep = steepest descent minimization)
emtol  = 1000.0   ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep          = 0.01          ; Energy step size
nsteps  = 50000    ; Maximum number of (minimization) steps to perform
 
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist  = 1  ; Frequency to update the neighbor list and long range forces
ns_type  = grid  ; Method to determine neighbor list (simple, grid)
rlist  = 1.5  ; Cut-off for making neighbor list (short range forces)
coulombtype = PME  ; Treatment of long range electrostatic interactions
rcoulomb = 1.5  ; Short-range electrostatic cut-off
rvdw  = 1.5  ; Short-range Van der Waals cut-off
pbc  = xyz   ; Periodic Boundary Conditions (yes/no)
 
freezegrps = Terminal_atoms
freezedim = Y Y Y