question for PDBFixer

[yhgon]

Hi,

I want to make simulation pipeline from RCSB PDB file to OpenMM input with pdbfixer.

so I get get simple PDB file.

wget -N -O short.pdb -q https://files.rcsb.org/download/1CQU.pdb

I check PDB file and it has only one chain and 56 residues.

pdb2fasta(filename='short.pdb')
0 A 56 MKVIFLKDVKGKGKKGEIKNVADGYANNFLFKQGLAIEATPANLKALEAQKQKEQR

when I try to get openMM system, I got error in PDB file with OpenMMException: NonbondedForce: The cutoff distance cannot be greater than half the periodic box size.

pdb = PDBFile('short.pdb')
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, 
                                 nonbondedCutoff=1.0*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)

---------------------------------------------------------------------------
OpenMMException                           Traceback (most recent call last)
<ipython-input-75-7c91fa975bcf> in <module>
----> 1 simulation = Simulation(pdb.topology, system, integrator)

/opt/conda/lib/python3.8/site-packages/simtk/openmm/app/simulation.py in __init__(self, topology, system, integrator, platform, platformProperties, state)
     99         if platform is None:
    100             ## The Context containing the current state of the simulation
--> 101             self.context = mm.Context(self.system, self.integrator)
    102         elif platformProperties is None:
    103             self.context = mm.Context(self.system, self.integrator, platform)

/opt/conda/lib/python3.8/site-packages/simtk/openmm/openmm.py in __init__(self, *args)
  13230             a set of values for platform-specific properties. Keys are the property names.
  13231         """
> 13232         _openmm.Context_swiginit(self, _openmm.new_Context(*args))
  13233 
  13234         self._system = args[0]

OpenMMException: NonbondedForce: The cutoff distance cannot be greater than half the periodic box size.

so I try to fix PDB file with PDBFixer same as manual however I still have problem. how to fix this problem?

from pdbfixer import PDBFixer
fixer = PDBFixer(filename='short.pdb')
fixer.findMissingResidues()
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()
fixer.removeHeterogens(True)
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
fixer.addSolvent(fixer.topology.getUnitCellDimensions())
PDBFile.writeFile(fixer.topology, fixer.positions, open('short_fixed.pdb', 'w'))

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-71-f6f97a79181e> in <module>
      7 fixer.addMissingAtoms()
      8 fixer.addMissingHydrogens(7.0)
----> 9 fixer.addSolvent(fixer.topology.getUnitCellDimensions())
     10 PDBFile.writeFile(fixer.topology, fixer.positions, open('short_fixed.pdb', 'w'))

/opt/conda/lib/python3.8/site-packages/pdbfixer/pdbfixer.py in addSolvent(self, boxSize, padding, boxVectors, positiveIon, negativeIon, ionicStrength)
   1081         modeller = app.Modeller(self.topology, self.positions)
   1082         forcefield = self._createForceField(self.topology, True)
-> 1083         modeller.addSolvent(forcefield, padding=padding, boxSize=boxSize, boxVectors=boxVectors, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength)
   1084         chains = list(modeller.topology.chains())
   1085         if len(chains) == 1:

/opt/conda/lib/python3.8/site-packages/simtk/openmm/app/modeller.py in addSolvent(self, forcefield, model, boxSize, boxVectors, padding, numAdded, positiveIon, negativeIon, ionicStrength, neutralize)
    630 
    631         # Add ions to neutralize the system.
--> 632         self._addIons(forcefield, numTotalWaters, waterPos, positiveIon=positiveIon, negativeIon=negativeIon, ionicStrength=ionicStrength, neutralize=neutralize)
    633 
    634     class _ResidueData:

/opt/conda/lib/python3.8/site-packages/simtk/openmm/app/modeller.py in _addIons(self, forcefield, numWaters, replaceableMols, ionCutoff, positiveIon, negativeIon, ionicStrength, neutralize)
    323         if neutralize:
    324             if abs(totalCharge) > numReplaceableMols:
--> 325                 raise Exception('Cannot neutralize the system because the charge is greater than the number of available positions for ions')
    326             if totalCharge > 0:
    327                 numNegative += totalCharge

Exception: Cannot neutralize the system because the charge is greater than the number of available positions for ions

openMM and PDBFixer from conda-forge

openmm             conda-forge/linux-64::openmm-7.5.1-py38hafe6fa4_1
pdbfixer           conda-forge/noarch::pdbfixer-1.7-pyhd3deb0d_0

[peastman]

The error is exactly what it says. If you look in the PDB file you'll find the line

CRYST1    1.000    1.000    1.000  90.00  90.00  90.00 P 1           1          

That means this model is a 1 A cube. That means the cutoff distance can't be more than 0.5 A. You tell it to use 1 nm (that is, 10 A).

I'm pretty sure that CRYST1 record is incorrect. The coordinates in the file span more than 1 A. But that's what the file says, so that's what it uses.

The second error is caused by the same problem:

fixer.addSolvent(fixer.topology.getUnitCellDimensions())

You tell it to add solvent, keeping the total size of the periodic box to the value loaded from the PDB file, that is, a 1 A cube. That doesn't give it room to actually add anything (including any neutralizing counterions). I would change that line to

fixer.addSolvent(padding=1*nanometer)

That will tell it to make sure there is at least 1 nm of water surrounding the solute on all sides.

[yhgon]

thanks for your kind help. it works well.