Deformed PDB files produced with PDBFixer

[jchodera]

pdbfixer seems to produce a highly deformed structure for PDB id code 1AO7.

You can reproduce this yourself by running the HTML-based pdbfixer on 1AO7 to simply build in the hydrogens. For some reason, they end up with all sorts of crazily acute angles (see images).

Any idea why this is? I haven't seen this happen before...

Resulting PDB structure can be found here: pdbfixer-output.pdb

[jchodera]

Actually, I'm seeing this behavior with a lot of PDB files of this size...

[peastman]

I'd guess something's going wrong with the energy minimization. It first adds hydrogens at fairly arbitrary positions, then runs an energy minimization with everything except hydrogens immobilized to fix them up. But if the minimization fails for some reason, it will just return the original positions.

[jchodera]

Let's check to make sure this isn't still a problem before we cut a release.

[peastman]

It looks like the minimization wasn't failing. The force just wasn't strong enough. If you don't specify a ForceField when you call addHydrogens(), it makes up a potential designed to give vaguely plausible positions. All the hydrogens repel each other (but only over very short distances), and there are also angle terms for a few special cases like water. The force on the hydrogens wasn't large enough, and so the minimizer was exiting before they had moved enough.

The fix is trivial, but unfortunately it's in Modeller, not in PdbFixer.

[jchodera]

So we would need an OpenMM 6.3.1?

[peastman]

Yes, if we decide it's important enough. Remember, as soon as you do a minimization with any real force field the problem will get fixed.

[jchodera]

This would seem to be useful for getting PDBFixer to work well as a standalone tool.

Could the hydrogen force constants in PDBFixer be increased instead as a workaround?

[peastman]

There are no hydrogen force constants in PDBFixer. It just calls addHydrogens().

[ajasja]

@jchodera @peastman This is still present when using OpenMM 7.1.1

Couldn't PDBFixer just accept a forcefield? (and use Amber by default)

Input:

Output:

Script:

pH = 7.5
fixer = pdbfixer.PDBFixer('NB26_APH.pdb')
forcefield = app.ForceField('amber99sbildn.xml', 'amber99_obc.xml')
fixer.findMissingResidues()
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()
fixer.removeHeterogens(False)
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(pH)
with open("NB26_APH.fix.pdb", 'w') as outfile:
    app.PDBFile.writeFile(fixer.topology, fixer.positions, outfile)

[jchodera]

Could you upload your NB26_APH.pdb file so we can take a look? You'll probably have to ZIP it up.

[ajasja]

Sure, here it is: NB26_APH.zip It contains a lot of lysins:)

[ajasja]

Minimization/simulation then fails.

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

pdb = app.PDBFile('struct/NB26_APH.fix.pdb')
#forcefield = app.ForceField('amber99sbildn.xml',  'amber99_obc.xml')
#forcefield = app.ForceField('amber99sbildn.xml')
forcefield = app.ForceField('amber10.xml', 'amber10_obc.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffPeriodic, 
    nonbondedCutoff=1.6*unit.nanometers, constraints=None, rigidWater=True)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 5.0/unit.picoseconds, 
    1.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

platform = mm.Platform.getPlatformByName('OpenCL')
properties = {'OpenCLPrecision': 'mixed'}
simulation = app.Simulation(pdb.topology, system, integrator, platform, 
    properties)

simulation.context.setPositions(pdb.positions)

print('Minimizing...')
simulation.minimizeEnergy(maxIterations=100000)

positions = simulation.context.getState(getPositions=True).getPositions()
with open("struct/NB26_APH.min.pdb", 'w') as outfile:
    app.PDBFile.writeFile(simulation.topology, positions, outfile)

simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(1000)

simulation.reporters.append(app.DCDReporter('out/trajectory.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, 
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True, 
    speed=True, totalSteps=10000, separator='\t'))

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

[jchodera]

I just ran your script with the OpenMM 7.2.0 beta and the corresponding pdbfixer, and things look fine: Can you try again after upgrading to these versions?

If you're using conda:

conda install --yes -c omnia/label/beta openmm==7.2.0 pdbfixer

Let us know if that works for you. If you still have trouble, I wonder if this might be a platform-specific issue.

[jchodera]

P.S. Your renderings look beautiful! Did you use VMD to generate those?

[ajasja]

Hi!

Thanks, that's done using Chimera (which is very beautiful by default:)

The problem remains using OpenMM 7.2 and pdbfixer 1.4-py35_0.

Fixed:

Minimized: PS: why is minimization not working?

Script: Nb26-APH.zip also here.

[jchodera]

The latest pdbfixer is 1.5. Try forcing upgrade with:

conda install --yes -c omnia/label/beta openmm==7.2.0 pdbfixer==1.5

If this still fails, it could be that it's trying to use the GPU and the minimization isn't making progress due to restricted precision.

@peastman: Is there a way to force the CPU platform to be used in pdbfixer?

[ajasja]

Updated pdbfixer to 1.5-py35heda29b4_0 omnia/label/beta

still the same result:

[jchodera]

It doesn't look like there's a way to specify the Platform to use in pdbfixer.

We'll need @peastman to help debug further.

[jchodera]

In the meantime, here's your fixed file!

NB26_APH.fix.pdb.zip

[ajasja]

Thanks! Although the simulation still crashes afterwards, even starting from the fixed file. The installation seems to be OK (can I run some more extensive tests?):

(ommprotocol) D:\data\MD\Nanobodies\Nb26-APH>python -m simtk.testInstallation
There are 3 Platforms available:

1 Reference - Successfully computed forces
2 CPU - Successfully computed forces
3 OpenCL - Successfully computed forces

Median difference in forces between platforms:

Reference vs. CPU: 1.98179e-05
Reference vs. OpenCL: 2.15121e-05
CPU vs. OpenCL: 1.55429e-05

[jchodera]

Can you post your simulation script too?

[ajasja]

Sure:

from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout

pdb = app.PDBFile('struct/NB26_APH.jchodera.fix.pdb')
#forcefield = app.ForceField('amber99sbildn.xml',  'amber99_obc.xml')
#forcefield = app.ForceField('amber99sbildn.xml')
forcefield = app.ForceField('amber10.xml', 'amber10_obc.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffPeriodic, 
    nonbondedCutoff=1.6*unit.nanometers, constraints=None, rigidWater=True)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 5.0/unit.picoseconds, 
    1.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

platform = mm.Platform.getPlatformByName('OpenCL')
properties = {'OpenCLPrecision': 'mixed'}
simulation = app.Simulation(pdb.topology, system, integrator, platform, 
    properties)

simulation.context.setPositions(pdb.positions)

print('Minimizing...')
simulation.minimizeEnergy(maxIterations=100000)

positions = simulation.context.getState(getPositions=True).getPositions()
with open("struct/NB26_APH.min.pdb", 'w') as outfile:
    app.PDBFile.writeFile(simulation.topology, positions, outfile)

simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(1000)

simulation.reporters.append(app.DCDReporter('out/trajectory.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, 
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True, 
    speed=True, totalSteps=10000, separator='\t'))

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

[jchodera]

This explodes for me on the CPU platform with and without HBond constraints. Minimization seems to make progress, so I'm not quite sure what's going on here.

[ajasja]

Should I move this to a separate issue on the OpenMM github tracker?

On 20 December 2017 at 19:09, John Chodera notifications@github.com wrote:

This explodes for me on the CPU platform with and without HBond constraints. Minimization seems to make progress, so I'm not quite sure what's going on here.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/pandegroup/pdbfixer/issues/92#issuecomment-353139992, or mute the thread https://github.com/notifications/unsubscribe-auth/AABN1n8o0Z2i2syXSaXlSZwqORwGaZknks5tCU1wgaJpZM4FHmoc .

[jchodera]

I think we should keep it here and figure out what's going on---I do suspect it might be something we can improve with pdbfixer.

[jaimergp]

I got it to work by changing the forcefield to forcefield = app.ForceField('amber99sbildn.xml', 'amber99_obc.xml') (first commented out line).

from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout

pdb = app.PDBFile('struct/NB26_APH.jchodera.fix.pdb')
forcefield = app.ForceField('amber99sbildn.xml',  'amber99_obc.xml')
#forcefield = app.ForceField('amber99sbildn.xml')
#forcefield = app.ForceField('amber10.xml', 'amber10_obc.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffPeriodic,
    nonbondedCutoff=1.6*unit.nanometers, constraints=None, rigidWater=True)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 5.0/unit.picoseconds,
    1.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

platform = mm.Platform.getPlatformByName('OpenCL')
properties = {'OpenCLPrecision': 'mixed'}
simulation = app.Simulation(pdb.topology, system, integrator, platform,
    properties)

simulation.context.setPositions(pdb.positions)

print('Minimizing...')
simulation.minimizeEnergy(maxIterations=100000)

positions = simulation.context.getState(getPositions=True).getPositions()
with open("struct/NB26_APH.min.pdb", 'w') as outfile:
    app.PDBFile.writeFile(simulation.topology, positions, outfile)

simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(1000)

simulation.reporters.append(app.DCDReporter('out/trajectory.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
    speed=True, totalSteps=10000, separator='\t'))

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

Result (converted to single PDB trajectory with MDTraj): traj.pdb.zip

[jchodera]

Interesting! I wonder why we're having trouble with amber10.

@peastman : Any ideas? Hydrogens without LJ centers?

[ajasja]

Jup, that works for me as well (on windows).

On 21 December 2017 at 16:13, John Chodera notifications@github.com wrote:

Interesting! I wonder why we're having trouble with amber10.

@peastman https://github.com/peastman : Any ideas? Hydrogens without LJ centers?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/pandegroup/pdbfixer/issues/92#issuecomment-353375272, or mute the thread https://github.com/notifications/unsubscribe-auth/AABN1md_iCPj1pSRNGTxrhpRnUQugu0mks5tCnWCgaJpZM4FHmoc .

[peastman]

No idea. It may just be random. Any tiny perturbation to the minimization process can change where it ends up.

[jchodera]

@peastman: Minimization worked fine on my machine with CPU, but Dynamics NaNd with any timesteps I tried.

[peastman]

Minimization worked fine on my machine with CPU

What does "worked fine" mean? It appears to have left the system in a conformation that couldn't be simulated.

If it's blowing up immediately, you can generally tell where the problem is. Just check which atoms have large forces on them.

[jchodera]

What does "worked fine" mean?

Minimization produced a structure with a large negative energy that was free of obvious structural defects by visual inspection.

If it's blowing up immediately, you can generally tell where the problem is. Just check which atoms have large forces on them.

Will do.