verbose pdb analysis and energy minimization

[NIKOLAY PLOTNIKOV]

Hi everyone,
does anyone know how to get extra insight from
A.
energy minimization

B.
reading pdb

I'm using a protonated and optimized protein with a script created based on Robert's app http://builder.openmm.org/

the problem is that everything explodes shortly after I start relaxation run following the energy minimization with a very tight threshold.


The output of MD relaxation looks like something (probably hydrogens) are not connected in the topology file and probably are causing nonbonded clashed but pdb.getTopology() gives a very llimited output

Code: Select all

from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout
from simtk.unit import *

pdb = app.PDBFile('MOE_FIX.pdb')
forcefield = app.ForceField('amber10.xml', 'amber10_obc.xml')

system = forcefield.createSystem(pdb.topology,
    nonbondedMethod=app.CutoffNonPeriodic, nonbondedCutoff=2.0*nanometers)
integrator = mm.LangevinIntegrator(300*kelvin, 1.0/picoseconds,
    0.5*femtoseconds)
#integrator.setConstraintTolerance(0.00001)

platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed', 'CudaDeviceIndex': '0'}
simulation = app.Simulation(pdb.topology, system, integrator, platform,
    properties)
simulation.context.setPositions(pdb.positions)

print('Minimizing...')
simulation.minimizeEnergy(tolerance=0.1*kilojoule/mole)
#simulation.reporters.append(app.DCDReporter('equilibrate.dcd', 1000))
#simulation.reporters.append(app.PDBReporter('equilibrate.pdb', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1, step=True,
    potentialEnergy=True, temperature=True,
    totalSteps=10000,separator ='\t' ))
simulation.context.setVelocitiesToTemperature(10*kelvin)
print('Equilibrating...')
simulation.step(10000)

#simulation.reporters.append(app.DCDReporter('trajectory.dcd', 10000))
#simulation.reporters.append(app.PDBReporter('trajectory.pdb', 10000))
#simulation.reporters.append(app.StateDataReporter(stdout, 10000, step=True,
#    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
#    speed=True, totalSteps=100000, separator='\t'))

#print('Running Production...')
#simulation.step(10000)
print('Done!')

the problem is that everything explodes shortly after I start relaxation run following the energy minimization with a very tight threshold.


The output of MD relaxation looks like something (probably hydrogens) are not connected in the topology file and probably are causing nonbonded clashed but pdb.getTopology() gives a very llimited output

Code: Select all

Minimizing...
Equilibrating...
#"Step"	"Potential Energy (kJ/mole)"	"Temperature (K)"
1	4.72611915361e+11	12326086179.3

[Peter Eastman]

You can call getBonds() on the Topology to get a list of all the bonds it contains. If some bonds are missing, you should be able to determine from that.

If that is really the problem, there are two possible causes that immediately occur to me.

1. Does the PDB file contain anything other than standard amino acids and nucleotides? If so, it needs to include CONECT records defining their bonds.

2. Perhaps the atoms are named in a nonstandard way? PDBFile can only fill in the bonds for standard amino acids and nucleotides if it can tell which atoms are which. If they have nonstandard names, it might not be able to figure out.

On the other hand, I don't see how either of those could actually be the problem. The reason is that, if any residues were missing bonds, then createSystem() would fail. It would report that it didn't have any templates matching those residues. So I suspect the problem is actually something else.

Peter