Adding ligand parameters for AMOEBA

[Saurabh Belsare]

Hi,

I am trying to simulate a system of protein + ligand using the AMOEBA force field in OpenMM 7. I have the AMOEBA parameters for the ligand, and I have previously simulated the system in TINKER. I am now trying to do the simulation in OpenMM.

I created a new AMOEBA parameter file from the existing amoeba2013.xml called amoeba2013_lig.xml, which has all the existing AMOEBA parameters, and added the ligand parameters. I'm trying to use this new parameter file to simulate my protein-ligand system.

When I tried to minimize my system with this file, it gave me the following error:

Traceback (most recent call last):
File "minimizePdb.py", line 14, in <module>
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
File "/usr/lib64/python2.6/site-packages/simtk/openmm/app/forcefield.py", line 380, in createSystem
raise ValueError('No template found for residue %d (%s). %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
ValueError: No template found for residue 254 (LIG). The set of atoms matches LIG, but the bonds are different.

Tracing back this error, I realized that the bond definitions for the ligand would need to be added separately from the amoeba2013_lig.xml file. Based on the function loadBondDefinitions in wrappers/python/simtk/openmm/app/topology.py, I thought I would need to have a lig_residues.xml, analogous to residues.xml, which I created, and added the bond connectivity of my ligand in the same way as residues.xml has it for the protein residues. However, after I call the loadBondDefinitions function on lig_residues.xml in my minimize.py, the job completes without any error, but the output structure is the same as the input structure, i.e. no minimization has occurred.

What is the correct way to incorporate the ligand parameterization into the xml files to run the simulation using AMOEBA?

Thanks.

Saurabh

[Peter Eastman]

Hi Saurabh,

It's probably actually simpler than you're making it. First, you don't need to modify the existing amoeba2013.xml file (and it's generally better not to). Just put the definitions for your ligand into its own file, then list both files when you create your ForceField.

Also, the simplest way to define bonds is just to include CONECT records in your PDB file. Then the PDB reader will add all the correct bonds automatically. loadBondDefinitions() is a much less usual way of doing it.

As for what's happening with minimization, it's hard to guess. Can you post your files so we can try to reproduce it?

Peter

[Saurabh Belsare]

Hi Dr. Eastman,

I tried using the CONECT keyword and providing the connectivities at the bottom of the PDB file before I tried this convoluted approach, but that still gave me the same ValueError: No template found for residue 254 (LIG). The set of atoms matches LIG, but the bonds are different. error. Maybe I did it wrong somehow. My input files can be found here: https://drive.google.com/open?id=0B9osC ... UIxZ1pBQk0

[Peter Eastman]

The atom indices in your CONECT records are off. The atoms in the ligand are 3971 through 3986:

Code: Select all

ATOM   3971  C7  LIG   254      14.312  42.477  34.152  1.00  0.00           C  
ATOM   3972  H7  LIG   254      13.745  41.695  34.651  1.00  0.00           H  
ATOM   3973  C6  LIG   254      15.672  42.648  34.438  1.00  0.00           C  
ATOM   3974  H6  LIG   254      16.181  41.997  35.144  1.00  0.00           H  
ATOM   3975  C4  LIG   254      16.423  43.557  33.739  1.00  0.00           C  
ATOM   3976  O3  LIG   254      17.705  43.875  33.909  1.00  0.00           O  
ATOM   3977  N2  LIG   254      18.107  44.794  32.941  1.00  0.00           N  
ATOM   3978  C1  LIG   254      16.980  45.116  32.309  1.00  0.00           C  
ATOM   3979  H1  LIG   254      17.058  45.883  31.548  1.00  0.00           H  
ATOM   3980  C5  LIG   254      15.832  44.352  32.824  1.00  0.00           C  
ATOM   3981  C8  LIG   254      14.480  44.049  32.373  1.00  0.00           C  
ATOM   3982  H8  LIG   254      14.023  44.558  31.528  1.00  0.00           H  
ATOM   3983  C9  LIG   254      13.698  43.171  33.102  1.00  0.00           C  
ATOM   3984  N10 LIG   254      12.283  42.926  32.735  1.00  0.00           N  
ATOM   3985  O12 LIG   254      11.646  43.734  32.142  1.00  0.00           O  
ATOM   3986  O11 LIG   254      11.922  41.814  32.916  1.00  0.00           O  

But the indices that appear in your CONECT records are 3972 through 3987:

Code: Select all

CONECT 3972 3973 3974 3984
CONECT 3973 3972
CONECT 3974 3972 3975 3976
CONECT 3975 3974
CONECT 3976 3974 3977 3981
CONECT 3977 3976 3978
CONECT 3978 3977 3979
CONECT 3979 3978 3980 3981
CONECT 3980 3979
CONECT 3981 3976 3979 3982
CONECT 3982 3981 3983 3984
CONECT 3983 3982
CONECT 3984 3972 3982 3985
CONECT 3985 3984 3986 3987
CONECT 3986 3985
CONECT 3987 3985

Peter

[Saurabh Belsare]

Hi Dr. Eastman,

You're right, there was in issue with the numbering. But on fixing that, after the minimization, the RMS between the protein before and after minimization is only 0.002. Also, assuming that that minimization has occurred correctly, when I now take my output of minimization as the input for nvt equilibration, it gives me the following error:

Traceback (most recent call last):
File "equilibrate_nvt.py", line 26, in <module>
simulation.step(1)
File "/usr/lib64/python2.6/site-packages/simtk/openmm/app/simulation.py", line 95, in step
self._simulate(endStep=self.currentStep+steps)
File "/usr/lib64/python2.6/site-packages/simtk/openmm/app/simulation.py", line 172, in _simulate
reporter.report(self, state)
File "/usr/lib64/python2.6/site-packages/simtk/openmm/app/pdbreporter.py", line 78, in report
PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
File "/usr/lib64/python2.6/site-packages/simtk/openmm/app/pdbfile.py", line 285, in writeModel
raise ValueError('Particle position is NaN')
ValueError: Particle position is NaN

I've attached my equilibration python script here:

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout

if __name__=="__main__":

writeFreq = 1
nsteps = 100

platform = Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'}

forcefield = ForceField('amoeba2013.xml', 'lig.xml')
pdb = PDBFile('input_minimized.pdb') # This is the output of the minimization from the previous minimization files.
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.001*picoseconds)

simulation = Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.reporters.append(PDBReporter('input_equilibrated.pdb',writeFreq))
simulation.reporters.append(StateDataReporter(stdout, writeFreq, step=True, potentialEnergy=True, temperature=True))

for i in range(nsteps):
simulation.step(1)

[Peter Eastman]

Could you post the fixed PDB file? Then I can try to reproduce it.

Peter

[Saurabh Belsare]

The input pdb file is here:

https://drive.google.com/open?id=0B9osC ... kotMUp0Wlk The lig.xml which has the parameters for the ligand is in https://drive.google.com/open?id=0B9osC ... UIxZ1pBQk0

Thanks.

[Peter Eastman]

I think you have some incorrect units in your force field definition. For example, you have

Code: Select all

  <Vdw class="301" sigma="3.8000" epsilon="0.1010" reduction="1.0" /> 

That specifies 3.8 nm for sigma. I think you probably meant that to be 3.8 A (0.38 nm) instead? The result is a huge energy. The minimizer tries valiantly to bring it down, but there isn't much it can do.

Peter

[Peter Eastman]

Your polarizabilities are also clearly in the wrong units. They're about three orders of magnitude larger than the ones in amoeba2013.xml.

Peter

[Saurabh Belsare]

Hi Dr. Eastman,

That's solved the problem. I thought I had converted all the quantities from the TINKER units to OpenMM units, but clearly had missed these two. That was a stupid mistake.

Thanks a lot!