Molecular Simulation with OpenMM

OpenMM is a toolkit for molecular simulation. It is a physic based libraries that is useful for refining the structure and exploring functional interactions with other molecules. It provides a combination of extreme flexibility (through custom forces and integrators), openness, and high performance (especially on recent GPUs) that make it truly unique among simulation codes.

Protein data

We use a processed 2DRI dataset that represents the ribose binding protein in bacterial transport and chemotaxis. The source organism is the Escherichia coli bacteria. You can find more details on this protein at the related RCSB Protein Data Bank page.

In [ ]:
from IPython.display import display
from PIL import Image

display(Image.open('2dri-image.png'))

Protein data can be stored in a .pdb file, this is a human readable format. It provides for description and annotation of protein and nucleic acid structures including atomic coordinates, secondary structure assignments, as well as atomic connectivity. See more information about PDB format here.

We are printing the first 10 lines of the file. The output contains a number of ATOM records. These describe the coordinates of the atoms that are part of the protein.

In [ ]:
%%bash
head ./2dri-processed.pdb
REMARK   1 CREATED WITH OPENMM 7.6, 2022-07-12
CRYST1   81.309   81.309   81.309  90.00  90.00  90.00 P 1           1 
ATOM      1  N   LYS A   1      64.731   9.461  59.430  1.00  0.00           N  
ATOM      2  CA  LYS A   1      63.588  10.286  58.927  1.00  0.00           C  
ATOM      3  HA  LYS A   1      62.707   9.486  59.038  1.00  0.00           H  
ATOM      4  C   LYS A   1      63.790  10.671  57.468  1.00  0.00           C  
ATOM      5  O   LYS A   1      64.887  11.089  57.078  1.00  0.00           O  
ATOM      6  CB  LYS A   1      63.458  11.567  59.749  1.00  0.00           C  
ATOM      7  HB2 LYS A   1      63.333  12.366  58.879  1.00  0.00           H  
ATOM      8  HB3 LYS A   1      64.435  11.867  60.372  1.00  0.00           H  

Write the script

To run the script above all we need is a Python environment with the OpenMM library installed.

In [ ]:
%%bash
conda install -c conda-forge openmm
In [ ]:
%%bash
python -m openmm.testInstallation
OpenMM Version: 8.0
Git Revision: a7800059645f4471f4b91c21e742fe5aa4513cda

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.30471e-06
Reference vs. OpenCL: 6.74827e-06
CPU vs. OpenCL: 7.07625e-07

All differences are within tolerance.
In [ ]:
import os
from openmm import *
from openmm.app import *
from openmm.unit import *

# Input Files
input_path = './2dri-processed.pdb'
os.path.exists(input_path) # check if input file exists
pdb = PDBFile(input_path)
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

# Output
output_path = './final_state.pdbx'
if not os.path.exists(os.path.dirname(output_path)): # check if ouput dir exists
    os.makedirs(os.path.dirname(output_path))

# System Configuration

nonbondedMethod = PME
nonbondedCutoff = 1.0*nanometers
ewaldErrorTolerance = 0.0005
constraints = HBonds
rigidWater = True
constraintTolerance = 0.000001
hydrogenMass = 1.5*amu

# Integration Options

dt = 0.002*picoseconds
temperature = 310*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25

# Simulation Options

steps = 10
equilibrationSteps = 0
#platform = Platform.getPlatformByName('CUDA')
platform = Platform.getPlatformByName('CPU')
#platformProperties = {'Precision': 'single'}
platformProperties = {}
dcdReporter = DCDReporter('trajectory.dcd', 1000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=steps,
    step=True, time=True, speed=True, progress=True, elapsedTime=True, remainingTime=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True, volume=True, density=True, separator='\t')
checkpointReporter = CheckpointReporter('checkpoint.chk', 1000)

# Prepare the Simulation

print('Building system...')
topology = pdb.topology
positions = pdb.positions
system = forcefield.createSystem(topology, nonbondedMethod=nonbondedMethod, nonbondedCutoff=nonbondedCutoff,
    constraints=constraints, rigidWater=rigidWater, ewaldErrorTolerance=ewaldErrorTolerance, hydrogenMass=hydrogenMass)
system.addForce(MonteCarloBarostat(pressure, temperature, barostatInterval))
integrator = LangevinMiddleIntegrator(temperature, friction, dt)
integrator.setConstraintTolerance(constraintTolerance)
simulation = Simulation(topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(positions)

# Minimize and Equilibrate

print('Performing energy minimization...')
simulation.minimizeEnergy()
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(equilibrationSteps)

# Simulate

print('Simulating...')
simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0
simulation.step(steps)

# Write file with final simulation state

state = simulation.context.getState(getPositions=True, enforcePeriodicBox=system.usesPeriodicBoundaryConditions())
with open(output_path, mode='w+') as file:
    PDBxFile.writeFile(simulation.topology, state.getPositions(), file)
print('Simulation complete, file written to disk at: {}'.format(output_path))
Building system...
Performing energy minimization...
Equilibrating...
Simulating...
Simulation complete, file written to disk at: ./final_state.pdbx
In [ ]:
%%bash
head ./final_state.pdbx
data_cell
# Created with OpenMM 8.0, 2023-03-12
#
_cell.length_a        81.3090
_cell.length_b        81.3090
_cell.length_c        81.3090
_cell.angle_alpha     90.0000
_cell.angle_beta      90.0000
_cell.angle_gamma     90.0000
#


什么是 OpenMM?

OpenMM 是一个高性能开源软件库,用于分子动力学 (MD) 模拟。它主要用于模拟和研究生物分子系统,如蛋白质、核酸和小分子。在这些模拟中,OpenMM 计算原子之间的相互作用力,并使用这些力来预测分子随时间的运动。

OpenMM 的主要特点

  1. 高性能计算

    • 支持 GPU 加速,显著提高模拟速度,适用于大规模系统和长时间尺度的模拟。
  2. 可扩展性

    • 允许用户定义自定义力场和模拟协议,适应不同的研究需求。
  3. 易于使用

    • 提供友好的 Python API,使得脚本化和自动化模拟变得容易。
  4. 跨平台

    • 支持多种操作系统,包括 Windows、macOS 和 Linux。
  5. 社区支持和文档

    • 提供丰富的文档和社区支持,帮助用户快速上手并解决问题。

OpenMM 使用量子化学吗?

OpenMM 主要用于经典分子动力学模拟,不直接使用量子化学方法。它依赖经典力场来描述分子间的相互作用,如 AMBER、CHARMM、OPLS 和 GROMOS 等。这些力场基于经验参数和半经验公式,用于计算原子间的相互作用力。

然而,OpenMM 可以与量子化学计算相结合,用于特定场景下的多尺度模拟。例如,在 QM/MM(量子力学/分子力学)方法中,系统的一部分(如反应中心)使用量子化学方法计算,而其余部分使用经典力场。这种结合可以通过与其他量子化学软件(如 Gaussian、NWChem 或 Psi4)配合使用来实现。

使用 OpenMM 进行分子动力学模拟的示例

以下是一个使用 OpenMM 进行简单分子动力学模拟的 Python 示例。假设你有一个 PDB 文件(input.pdb),描述了你的分子系统。

安装 OpenMM

首先,确保你已经安装了 OpenMM。你可以使用 conda 安装:

In [ ]:
%%bash
conda install -c conda-forge openmm

Python 脚本

下面是一个 Python 脚本,演示如何使用 OpenMM 读取 PDB 文件并进行分子动力学模拟:

In [ ]:
from openmm import *
from openmm.app import *
from openmm.unit import *
import sys

# 读取 PDB 文件
pdb = PDBFile('input.pdb')

# 使用 AMBER FF14SB 力场
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

# 创建系统对象
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, 
                                 nonbondedCutoff=1*nanometer, constraints=HBonds)

# 创建模拟对象
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)

# 设置初始位置
simulation.context.setPositions(pdb.positions)

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

# 运行分子动力学模拟
simulation.reporters.append(StateDataReporter(sys.stdout, 1000, step=True, 
                                              potentialEnergy=True, temperature=True))
simulation.reporters.append(PDBReporter('output.pdb', 1000))
print('Running simulation...')
simulation.step(10000)

print('Simulation completed')

解释代码

  1. 导入库

    • 导入 OpenMM 和相关库。
  2. 读取 PDB 文件

    • 使用 PDBFile 类读取输入的 PDB 文件。
  3. 选择力场

    • 使用 AMBER FF14SB 力场来创建系统。
  4. 创建系统对象

    • 使用力场创建系统,指定非键合相互作用和约束。
  5. 创建模拟对象

    • 使用 Langevin 积分器创建模拟对象,设置温度、摩擦系数和时间步长。
  6. 设置初始位置

    • 将初始原子位置设置为 PDB 文件中的位置。
  7. 能量最小化

    • 执行能量最小化以消除初始结构中的高能量冲突。
  8. 运行分子动力学模拟

    • 设置状态报告器和 PDB 报告器,进行分子动力学模拟并输出结果。

总结

OpenMM 是一个功能强大的分子动力学模拟工具,适用于经典分子力学模拟。尽管它不直接使用量子化学方法,但可以与量子化学软件结合,进行多尺度模拟。通过易于使用的 Python API,研究人员可以高效地进行分子动力学模拟和分析。


获取 .pdb 文件通常有几种途径,具体取决于你需要的分子结构和研究目的。以下是一些常见的获取途径:

  1. 蛋白质数据银行 (Protein Data Bank, PDB)

    • PDB 是一个公开的数据库,包含了大量的蛋白质、核酸和其他生物大分子的三维结构。你可以在 PDB 的网站上搜索并下载你感兴趣的分子的 .pdb 文件。
  2. 分子建模软件

    • 一些分子建模软件,如PyMOL、VMD 和 ChimeraX,可以从 PDB 数据库中直接获取分子结构并保存为 .pdb 文件。
  3. 分子模拟软件

    • 使用分子模拟软件进行分子动力学模拟时,通常会生成包含模拟结果的 .pdb 文件。你可以从模拟结果中提取所需的分子结构并保存为 .pdb 文件。
  4. 实验数据

    • 如果你在实验室进行生物化学或结构生物学研究,可能会获得实验测定的分子结构数据。这些数据通常可以转换为 .pdb 格式。
  5. 在线数据库和资源

    • 除了 PDB,还有一些其他数据库和资源提供了分子结构数据,例如 PubChem、ChemSpider 和 DrugBank。你可以在这些数据库中搜索并下载 .pdb 文件。
  6. 文献和研究文章

    • 有时候,科学研究论文中会提供分子结构数据的补充材料,包括 .pdb 文件。你可以通过阅读相关的文献或联系作者来获取这些数据。

通过以上途径,你应该能够获取到你需要的 .pdb 文件,从而进行后续的分子建模、分子动力学模拟或其他分子结构分析工作。


The MMTF Python API is available from pip:

In [ ]:
%%bash
pip3 install mmtf-python
  1. Get the data for a PDB structure and print the number of chains:
In [ ]:
from mmtf import fetch
# Get the data for 4CUP
decoded_data = fetch('4CUP')
In [ ]:
print(f'PDB Code: {decoded_data.structure_id} has {decoded_data.num_chains} chains')
PDB Code: 4CUP has 6 chains
  1. Show the charge information for the first group:
In [ ]:
atomic_charges = ','.join([str(x) for x in decoded_data.group_list[0]['formalChargeList']])
In [ ]:
print(f'Group name: {decoded_data.group_list[0]["groupName"]} has the following atomic charges: {atomic_charges}')
Group name: GLU has the following atomic charges: 0,0,0,0,0,0,0,0,0,0,0,0,0,0
  1. Show how many bioassemblies it has:
In [ ]:
print(f'PDB Code: {decoded_data.structure_id} has {len(decoded_data.bio_assembly)} bioassemblies')
PDB Code: 4CUP has 1 bioassemblies

In [ ]:
%%bash
pip install markdown
In [ ]:
import markdown
In [ ]:
markdown.markdown(
    '[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/bacalhau-project/examples/blob/main/molecular-dynamics/openmm/index.ipynb)'
)
Out[ ]:
'<p><a href="https://colab.research.google.com/github/bacalhau-project/examples/blob/main/molecular-dynamics/openmm/index.ipynb"><img alt="Open In Colab" src="https://colab.research.google.com/assets/colab-badge.svg" /></a></p>'
In [ ]:
markdown.markdown(
    '[![Open In Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/bacalhau-project/examples/HEAD?labpath=molecular-dynamics/openmm/index.ipynb)'
)
Out[ ]:
'<p><a href="https://mybinder.org/v2/gh/bacalhau-project/examples/HEAD?labpath=molecular-dynamics/openmm/index.ipynb"><img alt="Open In Binder" src="https://mybinder.org/badge.svg" /></a></p>'

=> How to Fold Proteins with Decentralized Computing?

Comments

2023-03-12