OpenMM——教程学习(1)

1. 如何从零开始做一个蛋白小分子动力学模拟

AmberTools将被用来生成输入文件，OpenMM 将被用来运行模拟，模拟平台为在线百度AI Stuio, 并使用GPU加速。
First thing’s first, 到PDB 蛋白数据库下载一需要模拟的靶点晶体，备用。
1. H++ web server ，输入蛋白编码，预处理蛋白，下载 后缀为 crd和top 两个文件备用 (Note: If you use it regularly, please register)
（如果第一步生成失败了，请阅读一下原因，如果是因为氨基酸序列缺失，backbone boken，or missing residue，调到第5步，使用openmm-setup修复下载，再回来走完整合流程。)
2. 在线（AI Studio）使用conda 安装一个独立的python3.9 环境，取名可为AmberTools (可能需要重启命令终端) ,如有询问y/n，输入y
(请自行用手机号提前注册一个百度AI studio平台的账号，新建一个"项目"，点击"启动环境"，这样，你在云端就有了一台性能不错的虚拟机)

conda create --name AmberTools python=3.9

3. 激活python环境, 并安装AmberTools21,如有询问y/n，输入y (！注意，如果你是用本地计算机需要用conda activate 但是在百度的平台我们需要用sourcre activate)

source activate AmberTools
conda install -c conda-forge ambertools=21 compilers 

运行AmberTools 所包含的一些列工具处理小分子，合并小分子和蛋白，生成新的crd和top文件（相对于第1步）

4. Windows本地电脑安装Anaconda，然后启动Anaconda命令行，输入(与第3步类似)，并安装一个用于准备OpenMM动力学输入文件的图形界面工具,如有询问y/n，输入y

conda create --name openmm python=3.9 # 这一句命令后，可能需要关掉再打开
conda activate openmm
（openmm）conda install -c conda-forge openmm-setup

5.启动openmm-setup图形界面，你的默认浏览器将会自动弹出成为一个动力学准备软件界面, 输入openmm-setup，回车即可

（openmm）openmm-setup

6. 回到百度AI Stuido, 在公开项目中搜索”OpenMM-蛋白-小分子动力学模拟“项目，喜欢，点赞后，fork到你自己的账号中，点击，”运行一下“。上传第5步生成的输入文件，开始模拟。
7. 下载模拟生成的output.pdb或者output.dcd到本地计算机，使用VMD进行分析，分析的教程请自行搜索，再次不在赘述。

步骤3：（在继续之前，请验证你已经成功完成了ambertools的安装）
conda activate ambertools
注意，第1步下载的两个文件的名称分别为 0.15_80_10_pH6.5_1ODX.top，和0.15_80_10_pH6.5_1ODX.crd, 数字代表了生成这些文件的H++web server条件。这些名字是什么不重要，你可以随意修改。下面一条命令是利用这两个文件来重新生成一个蛋白pdb,我们后面的模拟就用这个新的蛋白。新蛋白的名字为0.15_80_10_pH6.5_1ODX.pdb

ambpdb -p 0.15_80_10_pH6.5_1ODX.top -c 0.15_80_10_pH6.5_1ODX.crd > 0.15_80_10_pH6.5_1ODX.pdb

接下来，我们需要用到最最开始我们下载的那个PDB结构了，虽然这个结构中的蛋白部分我们不再需要（因为我们刚刚新生成了一个），但是小分子还是要从这个结构中提取。并对提取的这个小分子进行加氢，或者还原，并对其格式进行正规化矫正。从这里开始，我们下面将一共生成6个与小分子有关的结构，即ligand1到ligand6，注意命名方式。

awk '$4=="0E8"' 1ODX.pdb > ligand1.pdb # 该命令意思为提取pdb文件的第四列，小分子在这里
reduce ligand1.pdb > ligand2.pdb # 还原，也就是向小分子加氢
#如果加氢失败，请使用open babel 加氢。
#conda install -c conda-forge openbabel
#obabel -ipdb ligand1.pdb -opdb -O ligand2.pdb -h

pdb4amber -i ligand2.pdb -o ligand3.pdb # 格式美化，或者格式矫正

（选择性步骤：用txt编辑器打开确认你的clean版本的小分子pdb文件已经加氢，并且其格式不同于格式标准化之前的版本）
接下来我们需要生成小分子的mol2格式文件
(注意！如果小分子含有氧原子如酮，需要特别注意加氢是否正确，特别是如果使用pymol代替reduce加氢的情况下)

antechamber -fi pdb -i ligand3.pdb -fo mol2 -o ligand4.mol2 -c bcc -pf y
#This step may take a little bit of time, upon finish, check "sqm.out", the last line should be
# -------------calculation complete---------------------
antechamber -i ligand4.mol2 -fi mol2 -o ligand5.prepi -fo prepi -pf y
parmchk2 -f prepi -i ligand5.prepi -o ligand6.frcmod

(ATTENTION: You DON’T have to generate prepi file, you could move on with ONLY mol2 file, but then you need modify “http://tleap.in” file accordingly.)
……………………………………………………………………………………………………………………………………
合并小分子和蛋白，并对复合物pdb格式进行标准化

cat 0.15_80_10_pH6.5_1ODX.pdb 0E8_clean_H.pdb > 1ODX_H.pdb
pdb4amber -i 1ODX_H.pdb -o 1ODX_clean_H.pdb

新建一个文件，名称与格式为 tleap.in, 写入下面的内容，保存

source leaprc.protein.ff14SB #Source leaprc file for ff14SB protein force field
source leaprc.gaff #Source leaprc file for gaff
source leaprc.water.tip3p #Source leaprc file for TIP3P water model
loadamberprep ligand5.prepi #Load the prepi file for the ligand
loadamberparams ligand6.frcmod #Load the additional frcmod file for ligand
mol = loadpdb 1ODX_clean_H.pdb #Load PDB file for protein-ligand complex
solvatebox mol TIP3PBOX 8 #Solvate the complex with a cubic water box
addions mol Cl- 0 #Add Cl- ions to neutralize the system
saveamberparm mol 1ODX.prmtop 1ODX.inpcrd #Save AMBER topology and coordinate files
quit #Quit tleap program

(As we have mentioned earlier, if you do not have a “prepi” file, that is fine, what you need do is delete the “loadamberprep” line and add a new line for loading the mol2 file directly.)

在命令终端使用AmberTools的tleap 运行这个文件

tleap -s -f 1ODX_tleap.in > 1ODX_tleap.out

这样我们就得到了新的两个文件1ODX.prmtop 1ODX.inpcrd。这两个文件将是第5步的输入文件。
步骤3是最容易出错也是最重要的一个步骤，请参考
Tutorial for the LEaP Program

选择 Amber 文件，然后上传上面得到的两个文件，OpenMM 会要求用户输入模拟的步数，即时长，模拟的溶剂类型，以及很多其他参数，正常情况下是都有一个默认值已经填好的，依次点击下一步，平台 platform 选择 CUDA，如果你可用的平台没有GPU则更改CUDA为CPU即可。最后保存， ”save all files" ，得到一个压缩文件。里面包含了我们openmm 模拟需要的文件和脚本。如果你不知道如何修改界面上的参数，请全部保持默认即可！

AmberInpcrdFile('SYS_gaff2.crd')

你需要把上一步生成的你的prmtop文件和inpcrd文件(或者crd)替换到以上两个括号里。

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

prmtop = AmberPrmtopFile('SYS_gaff2.prmtop')
inpcrd = AmberInpcrdFile('SYS_gaff2.crd')
system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
        constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)
if inpcrd.boxVectors is not None:
    simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('amber_output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,
        potentialEnergy=True, temperature=True))
simulation.step(10000)

你可以使用上面这段代码，也可以使用下面这一段，作用是一样的，模拟的参数文件稍有区别。 如果使用下面这一段，我们需要替换的两个文件是 “complex.prmtop" 和 “complex.inpcrd”

from simtk.openmm import *
from simtk.openmm.app import *
from simtk.unit import *

# Input Files

prmtop = AmberPrmtopFile('complex.prmtop')
inpcrd = AmberInpcrdFile('complex.inpcrd')

# 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.004*picoseconds
temperature = 300*kelvin
friction = 1.0/picosecond
pressure = 1.0*atmospheres
barostatInterval = 25

# Simulation Options

steps = 1000000
equilibrationSteps = 1000
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'single'}
dcdReporter = DCDReporter('trajectory.dcd', 10000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=steps,
    step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator='\t')
checkpointReporter = CheckpointReporter('checkpoint.chk', 10000)

# Prepare the Simulation

print('Building system...')
topology = prmtop.topology
positions = inpcrd.positions
system = prmtop.createSystem(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)
if inpcrd.boxVectors is not None:
    simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)

# 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)

最理想的情况是，你自己根据自己的需要或者理解用Openmm生成适合自己蛋白体系的代码区块然后替换上面的例子。如果你会使用本地Jupyter notebook, 理论上也是可以打开以上脚本运行的，但墙裂不建议这么做，因为那样速度将非常非常慢，还容易出错。自己有Linux系统的，可以考虑把步骤2和3挪到本地，这样会更容易一些。如:

https://github.com/quantaosun/Ambertools-OpenMM-MD