[gmx-users] Energy minimization failure

[Justin Lemkul]

On Tue, Mar 19, 2013 at 9:16 AM, Jeff Woodford <
jwoodford at missouriwestern.edu> wrote:

> 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.
>
>
I would suggest a few things:

1. Simplify as much as possible, using small units of structure to diagnose.
2. Look at various energy terms in the .edr file. One will likely stick out
like a sore thumb in terms of being huge and positive, hence the large
forces. My suspicion would be short-range Coulombic forces. Including
nonbonded terms in 1-2 and 1-3 interactions is very weird, since they may
oppose the normal bonded forces.

The "bonds" you see forming are not actually there. Your visualization
program is just guessing what's going on. It's certainly indicative of bad
geometry, but no bonds are breaking or forming, unless you're doing QM and
not telling us ;)

-Justin

-- 

========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540)
231-9080http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================