OpenMM-组蛋白甲基转移酶分子动力学模拟-TIP3P

from pdbfixer import PDBFixer
from simtk.openmm.app import PDBFile

fixer = PDBFixer('1o9s.pdb')
# 看下fixer里面几条链
print(list(fixer.topology.chains()))
# 查看fixer里面链1有几个氨基酸，并转化为列表
residue_list = list(fixer.topology._chains[1].residues())
# 除去不相关的氨基酸
fixer.removeChains([1,2,3,4,5,7,9])
# 查看fixer的topo
print(fixer.topology)

# 寻找缺失残基
fixer.findMissingResidues()

# 在链的中间补充缺失缺失氨基酸，而不是末端
chains = list(fixer.topology.chains())
keys = fixer.missingResidues.keys()
missingResidues = dict()
for key in keys:
    chain = chains[key[0]]
    if not (key[1] == 0 or key[1] == len(list(chain.residues()))):
        missingResidues[key] = fixer.missingResidues[key]
fixer.missingResidues = missingResidues

fixer.findMissingAtoms()
fixer.addMissingAtoms()

# 输出
PDBFile.writeFile(fixer.topology, fixer.positions, open('1o9s_fixed.pdb', 'w'))

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

# 导入准备好的蛋白文件
# 力场文件通过xml导入
pdb = app.PDBFile('1o9s_fixed.pdb')
forcefield = app.ForceField('amber99sbildn.xml', 'tip3p.xml')
# 使用app.Modeller来进行加氢以及溶剂化
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield)
# 溶剂化，默认阳离子 'Na+' ，阴离子 'Cl-'
# 使用1ns的填充距离
modeller.addSolvent(forcefield, model='tip3p', padding=1.0*unit.nanometers)
# 输出文件
app.PDBFile.writeFile(modeller.topology, modeller.positions, open('1o9s_modeller_tip3p.pdb', 'w'))

# 准备系统以及integrator
system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,
    nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True,
    ewaldErrorTolerance=0.0005)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,
    2.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

# 准备模拟
# CPU,混合精度
platform = mm.Platform.getPlatformByName('CPU')
properties = {'CpuThreads': 'mixed'}
simulation = app.Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)

# 最小化
print('Minimizing...')
simulation.minimizeEnergy()

# 平衡 100 steps 
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(100)

# 增加 reporters
simulation.reporters.append(app.DCDReporter('trajectory_tip3p.dcd', 100, enforcePeriodicBox=False))
simulation.reporters.append(app.pdbreporter.PDBReporter('trajectory_tip3p.pdb', 100,enforcePeriodicBox=False))
simulation.reporters.append(app.StateDataReporter('trajectory_tip3p.csv', 100, step=True,
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
    speed=True,totalSteps=5000, separator=','))


# 10ps的模拟
print('Running Production...')
simulation.step(5000)
print('Done!')

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

# load in input PDB file and force field XML files
pdb = app.PDBFile('1o9s_fixed.pdb')
forcefield = app.ForceField('amber99sbildn.xml', 'tip4pew.xml')

# use app.Modeller to add hydrogens, extra particles, and solvent
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens()
modeller.addExtraParticles(forcefield)
modeller.addSolvent(forcefield, model='tip4pew', padding=1.0*unit.nanometers)
app.PDBFile.writeFile(modeller.topology, modeller.positions, open('1o9s_modeller_tip4pew.pdb', 'w'))

# prepare system and integrator
system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,
    nonbondedCutoff=1.0*unit.nanometers, constraints=app.HBonds, rigidWater=True,
    ewaldErrorTolerance=0.0005)
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,
    2.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)

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

# minimize
print('Minimizing...')
simulation.minimizeEnergy()

# equilibrate for 100 steps
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)
print('Equilibrating...')
simulation.step(100)

# append reporters
simulation.reporters.append(app.DCDReporter('trajectory_tip4pew.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True,
    potentialEnergy=True, temperature=True, progress=True, remainingTime=True,
    speed=True, totalSteps=25000000, separator='\t'))

# run 50 ns of production simulation
print('Running Production...')
simulation.step(25000000)
print('Done!')