Geometry optimization with AMOEBA force field

[Pushp Bajaj]

Hello,

I was recently using OpenMM to do a energy minimization with AMOEBA for a simple H2O-Cl- dimer system. I printed out the energy and the forces for the optimized geometry and I got non-zero forces on the atoms in the optimized geometry. I used the following commands to perform the optimization and to get the energy and forces:

"simulation.minimizeEnergy(tolerance=0.0001, maxIterations=2000)"
or
"mm.openmm.LocalEnergyMinimizer.minimize(simulation.context, 0.0001)"

state=simulation.context.getState(getEnergy=True,getForces=True,getPositions=True)
energy=state.getPotentialEnergy()
forces=state.getForces()

I also did an optimization of the same system with AMOEBA using TINKER. To compare with OpenMM I used the optimized geometry from TINKER and calculated the energy and forces with OpenMM; I got the energies to be the same to the order of 10^-2 but there were large forces on the atoms. Does anyone know what might be causing these differences and why there are non-zero forces on the atoms in the optimized structure?

Thanks a lot!

[Jason Swails]

Because you are not at a true local minimum. Note that it is very difficult to get to a true local minimum (where the force on every atom is zero). The minimization criteria for OpenMM is (I believe) satisfied when the difference in energy between successive steps is less than the tolerance. Forces are a lot more sensitive than energies.

The OpenMM energy minimized is intended for use in preparing for MD. Maybe try another minimizer (like one from scipy based on second derivatives) to get closer to the true local minimum.

[Pushp Bajaj]

Thank you for your reply!
Yes, I was to able get a more strict optimization done with TINKER where the rms of the final gradients is of the order of 10^-5. I thought it was just strange that after using the TINKER optimized geometry in OpenMM I still got large forces on the atoms.

[Peter Eastman]

When you say "large forces", how large do you mean? The OpenMM minimizer should be able to get the forces very small. It just depends on the tolerance and maximum number of iterations. You may also need to specify double precision, since there's a limit to how precisely it can minimize in single precision.

[Pushp Bajaj]

Using the optimized geometry from TINKER in OpenMM I got forces of the order of 10^2 kJ/mol/nm. After optimizing with OpenMM the final geometry had forces of the order of 10 kJ/mol/nm.
If I'm not mistaken the tolerance is for the rms of the gradients, I played around with that changing it from 10^-3 to 10^-8 and the results did not change.

[Peter Eastman]

Is that in single or double precision? If single, try switching to double. It should be able to minimize it a lot further than that.