OpenMM-简介与案例

DrugSci

发布于 2021-03-04 10:16:26

1.5K0

发布于 2021-03-04 10:16:26

文章被收录于专栏：FindKeyFindKey

简介

openmm是用于分子模拟的高性能工具包。该代码是开源，并在MIT和LGPL license下。是 Omnia(预测生物分子模拟的工具套件)的一部分。

web: http://openmm.org

Install

conda install -c conda-forge openmm
# 如果特殊的cuda版本号
conda install -c conda-forge openmm cudatoolkit==10.0
#ps：支持gpu
#如果要在GPU上运行OpenMM，请确保已安装驱动程序。
#如果是Nvidia GPU，请从https://www.nvidia.com/cn/download/index.aspx下载最新的驱动程序。
#如果是AMD GPU，并且正在使用Linux或Windows，请从https://support.amd.com下载最新版本的驱动程序。
#检查openmm是否已经安装
python -m simtk.testInstallation

# out
OpenMM Version: 7.4.2
Git Revision: dc9d188939ad630d240e89806b185061f7cd661a

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: 6.30679e-06
Reference vs. OpenCL: 6.73408e-06
CPU vs. OpenCL: 8.89223e-07

All differences are within tolerance.

Run Examples

加载一个PDB文件
使用Amber14力场和TIP3P-FB水模型对其进行参数化
将其能量最小化
使用Langevin积分器对其进行10,000步模拟，40ps
并将快照保存到每1000个时间步的PDB文件称为output.pdb。

#所用蛋白为1stp
from pdbfixer import PDBFixer
from simtk.openmm.app.pdbfile import PDBFile
import os
def fix_pdb(pdb_id):
    path = os.getcwd()
    if len(pdb_id) != 4:
        print("Creating PDBFixer...")
        fixer = PDBFixer(pdb_id)
        print("Finding missing residues...")
        fixer.findMissingResidues()

        chains = list(fixer.topology.chains())
        keys = fixer.missingResidues.keys()
        for key in list(keys):
            chain = chains[key[0]]
            if key[1] == 0 or key[1] == len(list(chain.residues())):
                print("ok")
                del fixer.missingResidues[key]

        print("Finding nonstandard residues...")
        fixer.findNonstandardResidues()
        print("Replacing nonstandard residues...")
        fixer.replaceNonstandardResidues()
        print("Removing heterogens...")
        fixer.removeHeterogens(keepWater=True)

        print("Finding missing atoms...")
        fixer.findMissingAtoms()
        print("Adding missing atoms...")
        fixer.addMissingAtoms()
        print("Adding missing hydrogens...")
        fixer.addMissingHydrogens(7)
        print("Writing PDB file...")

        PDBFile.writeFile(
            fixer.topology,
            fixer.positions,
            open(os.path.join(path, "%s_fixed_pH_%s.pdb" % (pdb_id.split('.')[0], 7)),
                 "w"),
            keepIds=True)
        return "%s_fixed_pH_%s.pdb" % (pdb_id.split('.')[0], 7)
fix_pdb('1stp.pdb')
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout

#倒入pdb文件，同样openmm也支持PDBx/mmCIF格式
pdb = app.PDBFile('1stp_fixed_pH_7.pdb')
#确定使用的力场，Amber14：amber14-all.xml；TIP3P-FB water model：amber14/tip3pfb.xml
#使用XML files 来定义标准力场
#可用力场从：http://docs.openmm.org/latest/userguide/application.html#force-fields 进行查看
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
#将力场与分子拓扑文件相组合来创建一个对此模拟系统的完整的数学描述
#一些options：
# 将粒子网格Ewald用于长距离静电相互作用（nonbondedMethod = PME）
# 设置1nm截止值用于直接空间相互作用（nonbondedCutoff = 1 *纳米）
# 并约束所有涉及氢原子的键的长度（约束条件 = HBonds）。
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,nonbondedCutoff=1*nanometer, constraints=HBonds)
# 朗之万动力学 ( Langevin Dynamics ) ，模拟温度（300 K），摩擦系数（1 ps-1）和步长（0.004 ps）
# 还可以使用许多其他集成方法。例如，如果要在恒定能量而不是恒定温度下模拟系统，则可以使用VerletIntegrator。
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
# 将topology，system，integrator进行组合，创建一个Simulation对象，并将其分配给一个名为Simulation的变量。
#模拟对象管理运行模拟所涉及的所有过程
simulation = Simulation(pdb.topology, system, integrator)
#指定了用于模拟的初始原子位置。
simulation.context.setPositions(pdb.positions)
#进行局部能量最小化，建议在模拟开始前进行，减少应力
simulation.minimizeEnergy()
#创建一个“报告程序”以在模拟过程中输出结果， PDBReporter将结构写入PDB文件。
#指定输出文件为output.pdb，并应每1000步写入一个结构。
simulation.reporters.append(PDBReporter('output.pdb', 1000))
#在模拟运行时获取常规状态报告，可以监视其进度
#添加另一个报告程序，以每1000个时间步长打印一些基本信息：当前步长索引，系统的势能和温度。
#将stdout（而不用引号）指定为输出文件，这意味着将结果在终端显示。
#也可以给文件名（用引号引起来），将信息写到文件中。
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,potentialEnergy=True, temperature=True))
#开始模拟
simulation.step(10000)
#result

#"Step","Potential Energy (kJ/mole)","Temperature (K)"
1000,-10439.484774531746,299.56325970834155
2000,-10663.160555781746,298.4476284317964
3000,-11088.430087031746,296.83560819072136
4000,-10801.258212031746,296.276886377841
5000,-10968.391024531746,302.23009618402386
6000,-11008.289462031746,295.76176920559607
7000,-10883.461337031746,302.16333057574894
8000,-11009.863680781746,312.57213252166844
9000,-11113.148837031746,307.5580071086222
10000,-11322.344149531746,301.048477900256
#运行完成之后，可以使用任何可视化的程序对结果进行分析(VMD, PyMol, AmberTools, etc.)
#同样oenmm也支持AMBER，Gromacs，CHARMM文件等等，具体可见：http://docs.openmm.org/latest/userguide/application.html#installing-openmm

注意：

不要将蛋白使用其余工具加氢后再使用openmm进行处理，具体原理理由后续会进行阐述。

本文参与?腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2021-02-06，如有侵权请联系?cloudcommunity@tencent.com 删除

output