Adding ligand parameters for AMOEBA
- Saurabh Belsare
- Posts: 32
- Joined: Sat Aug 14, 2010 8:43 am
Adding ligand parameters for AMOEBA
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
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Adding ligand parameters for AMOEBA
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
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
- Posts: 32
- Joined: Sat Aug 14, 2010 8:43 am
Re: Adding ligand parameters for AMOEBA
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
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Adding ligand parameters for AMOEBA
The atom indices in your CONECT records are off. The atoms in the ligand are 3971 through 3986:
But the indices that appear in your CONECT records are 3972 through 3987:
Peter
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
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
- Saurabh Belsare
- Posts: 32
- Joined: Sat Aug 14, 2010 8:43 am
Re: Adding ligand parameters for AMOEBA
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)
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Adding ligand parameters for AMOEBA
Could you post the fixed PDB file? Then I can try to reproduce it.
Peter
Peter
- Saurabh Belsare
- Posts: 32
- Joined: Sat Aug 14, 2010 8:43 am
Re: Adding ligand parameters for AMOEBA
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.
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
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Adding ligand parameters for AMOEBA
I think you have some incorrect units in your force field definition. For example, you have
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
Code: Select all
<Vdw class="301" sigma="3.8000" epsilon="0.1010" reduction="1.0" />
Peter
- Peter Eastman
- Posts: 2564
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Adding ligand parameters for AMOEBA
Your polarizabilities are also clearly in the wrong units. They're about three orders of magnitude larger than the ones in amoeba2013.xml.
Peter
Peter
- Saurabh Belsare
- Posts: 32
- Joined: Sat Aug 14, 2010 8:43 am
Re: Adding ligand parameters for AMOEBA
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!
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!