To create an input file for FMO calculations, we need both the structural information of the protein with hydrogen atoms included and the charge information for each fragment. However, since PDB files do not contain hydrogen atoms by default, we must add them using an external tool. Additionally, because the PDB format does not provide a field for storing partial charges, we repurposed the temperature factor field for this purpose.
Since the mmCIF format allows partial charge attributes to be added without unnecessary hacks, I wrote the code in openMM and gemmi.
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout
import gemmi
pdb = PDBFile("1AKI.pdb")
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
modeller = Modeller(pdb.topology, pdb.positions)
residues=modeller.addHydrogens(forcefield)
system = forcefield.createSystem(
modeller.topology,
nonbondedMethod=PME,
nonbondedCutoff=1.0*nanometer,
constraints=HBonds)
partial_charges = []
for force in system.getForces():
if isinstance(force, openmm.NonbondedForce):
for atom_index in range(force.getNumParticles()):
charge, sigma, epsilon = force.getParticleParameters(atom_index)
partial_charges.append(str(charge)[:-2])
with open('tmp.cif', 'w') as f:
PDBxFile.writeFile(modeller.topology, modeller.positions, f)
doc = gemmi.cif.read('tmp.cif')
block = doc.sole_block()
loop = block.find_loop('_atom_site.type_symbol').get_loop()
loop.add_columns(['_atom_site.partial_charge'], value='?')
pcs = block.find_values('_atom_site.partial_charge')
for n, x in enumerate(pcs):
pcs[n] = partial_charges[n]
doc.write_file('output.cif')
I’m getting tired of writing FMO input generators every time I change jobs, so I’ve decided to upload them to GitHub as open-source software.
