AMOEBA03 water structure reproduction

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
POST REPLY
User avatar
Kyeong-jun Jeong
Posts: 8
Joined: Mon Aug 07, 2017 1:45 pm

AMOEBA03 water structure reproduction

Post by Kyeong-jun Jeong » Thu Sep 12, 2019 1:13 pm

Hi,

I cannot reproduce the AMOEBA water structure in literature (J. Phys. Chem. B 2003, 107, 5933 and J. Phys. Chem. B 2015, 119, 9423), when I use the amoeba2013.xml force field file basically given by OpenMM installation and perform NPT simulation in room temperature / 1 bar.

Density is underestimated seriously (I'm almost getting 0.92g/cm^3), with shifting O-O g(r) peak to the right (larger distance) by 0.2 angstroms. Even though the vdw and multipole parameters look converted correctly in the given .xml file.

Is there any known issue about this?

User avatar
Peter Eastman
Posts: 1905
Joined: Thu Aug 09, 2007 1:25 pm

Re: AMOEBA03 water structure reproduction

Post by Peter Eastman » Thu Sep 12, 2019 3:11 pm

It should exactly reproduce the literature value. Can you post the code you're using?

Did you remember to set rigidWater=False? AMOEBA water is supposed to be flexible.

User avatar
Kyeong-jun Jeong
Posts: 8
Joined: Mon Aug 07, 2017 1:45 pm

Re: AMOEBA03 water structure reproduction

Post by Kyeong-jun Jeong » Thu Sep 12, 2019 7:33 pm

I'm posting the script and initial structure file used for the test.

The essential part of the code is that:
(By the way, is it insufficient to dismiss constraint argument? As the default is set to none, I thought it's okay to leave without mentioning in the system definition)
-->Feedback : I'll first try with rigidWater=False. Now I can see that the default of rigidWater argument is not turned off, as other constraints are.
-->Still it's not fixed with confirming that water is set flexible.

====

from __future__ import print_function
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
from time import gmtime, strftime
from datetime import datetime

temperature=300.0*kelvin
pressure = 1.0*atmosphere
elstcutoff=1.0*nanometer
vdwcut=1.2*nanometer

barofreq = 100
pdb = PDBFile('amoebareform_water_m1000.pdb')

integ_md = LangevinIntegrator(temperature, 1/picosecond, 0.0005*picoseconds)

pdb.topology.loadBondDefinitions('amoeba_residues_reline.xml')
pdb.topology.createStandardBonds();

modeller = Modeller(pdb.topology, pdb.positions)
forcefield = ForceField('amoeba2013_urea_vdwopt.xml') #force field for water molecule hasn't been adjusted from original file.

system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, nonbondedCutoff=elstcutoff,
vdwCutoff=vdwcut, polarization='extrapolated') #AMOEBA : no constraints!

barostat = MonteCarloBarostat(pressure,temperature,barofreq)
system.addForce(barostat)

simmd = Simulation(modeller.topology, system, integ_md, platform, properties)
simmd.context.setPositions(modeller.positions)
simmd.minimizeEnergy(maxIterations=10000)

...
then I ran 10 ns of simulation, 2 x 10^7 MD steps in total.
Attachments
amoeba_water_test.tar.gz
(101.66 KiB) Downloaded 3 times

User avatar
Peter Eastman
Posts: 1905
Joined: Thu Aug 09, 2007 1:25 pm

Re: AMOEBA03 water structure reproduction

Post by Peter Eastman » Fri Sep 13, 2019 10:19 am

You also need the argument rigidWater=False in the call to createSystem(). If rigidWater=True (the default), water molecules are rigid regardless of the value of the constraints argument. That's because most standard water models are defined to be rigid, so you need that even if you don't want any constraints in the protein. AMOEBA water is unusual in being flexible. See http://docs.openmm.org/latest/api-pytho ... eateSystem for details.

User avatar
Kyeong-jun Jeong
Posts: 8
Joined: Mon Aug 07, 2017 1:45 pm

Re: AMOEBA03 water structure reproduction

Post by Kyeong-jun Jeong » Fri Sep 13, 2019 12:22 pm

Thank you. Actually I've tested with rigidWater=False in system creation, but still the volume expands. "iamoeba" was working fine, but until now I couldn't find reason why the simulation using pristine AMOEBA03 parameters make deviation from publications.

User avatar
Peter Eastman
Posts: 1905
Joined: Thu Aug 09, 2007 1:25 pm

Re: AMOEBA03 water structure reproduction

Post by Peter Eastman » Fri Sep 13, 2019 12:56 pm

I see another problem in your script: you shouldn't call createStandardBonds(). There's almost never a reason to do that yourself. PDBFile does it automatically when a file is loaded. By calling it again, you created a second set of bonds for every molecule, which effectively means all the force constants are doubled.

If you need to load custom bond definitions, just put the call to loadBondDefinitions() at the top, before you create the PDBFile. It will then use them automatically.

Your amoeba_residues_reline.xml file is incorrect. You're using the syntax for force field templates, which are different. Bond definition files don't include atom types, and the atom names need to match standard PDB names (unlike force field templates where the names can be arbitrary). See https://github.com/openmm/openmm/blob/m ... sidues.xml for an example of a bond definition file. You also don't need to include water in your file. It's already built in.

User avatar
Kyeong-jun Jeong
Posts: 8
Joined: Mon Aug 07, 2017 1:45 pm

Re: AMOEBA03 water structure reproduction

Post by Kyeong-jun Jeong » Fri Sep 13, 2019 6:16 pm

Thank you for detailed inspection! The bond creation was the problem.

POST REPLY