amber软件中动力学模拟的步骤

Molecular Dynamics simulation 从能量最小化到实际模拟 1 基本流程图 Minimization Equilibrium Production Run 1）概述 前面我们已经得到了 Amber 用来动力学模拟的 prmtop 和 inpcrd 文件，它们分别是参数文件和坐标文件。 我们先从一条命令说起来解释 Amber 是如何做动力学模拟的： sander O i mdin o mdout p prmtop c inpcrd r rst x mdcrd 动力学过程是一个连续地解牛顿运动方程的过程：上一个牛顿方 程结束时，蛋白质中各原子的位置和速度保留给下一个牛顿方程， 惟一改变的是原子的加速度，它会根据各种势能函数重新计算（势 能随原子坐标改变：E=f （ r,） ）。只不过每个牛顿方程的时间很 短，短到 fs（10 -15s）级，Amber 软件提供的 sander 主程序可以用来 自动地做这样的数值计算。 它需要参数文件（prmtop ）、坐标文件（inpcrd）、sander 程序配 置文件（mdin）来启动运行，我们已经有了前两种文件，本节内容 最主要的就是讲解如何配置我们需要的动力学模拟。sander 程序运 行过程中会输出临时文件（rst）保存坐标和速度，还有轨迹文件 （mdcrd）。 2）动力学过程 从基本流程图可以知道，一般的动力学过程也就可以分为三步： 能量最小化（minimization） 、体系平衡（equilibrium） 、实际动力学 模拟。 由于我们进行的初始结构来自晶体结构或同源建模，所以在分子 内部存在着一定的结构张力，能量最小化就是真正的动力学之前释 放这些张力，如果没有这个步骤，在动力学模拟开始之后，整个体 系可能会因此变形、散架。 另外，由于动力学模拟的是真实的生物体环境，因此必须使研究 对象升温升压到临界值，体系达到平衡，才能做实际的动力学模拟。 2 各流程输入文件 要通过 Amber 软件做动力学模拟，需要明白如何去配置上述过 程中的每一步。一般来说就是指定一些键/ 值对。 1） Minimization 第一行是标题，&cntrl 是起始符， ”/”是结束符，中间的键/值对就是参 数配置。 上述的参数配置可以归纳如下： imin=1 Choose a minimization run，指定做能量最小化的动力学。 ntx=1 Read coordinates but not velocities from ASCII formatted inpcrd coordinate file ntx 指定如何获得坐标信息，这里直接从 inpcrd 中读取坐标信息。 maxcyc=2000 Maximum minimization cycles 能量最小化的算法涉及循环迭代，这里指定迭代次数。 ncyc=1000 The steepest descent algorithm for the first 0-ncyc cycles, then switches the conjugate gradient algorithm for ncyc-maxcyc cycles 如上所说，循环迭代的算法不同，此处指定到哪一步第一种迭代算法结 束。 ntpr=100 Print to the Amber mdout output file every ntpr cycles 这个表明采集计算信息的频率，输出到 mdout 文件中。 cut=8.0 Nonbonded cutoff distance in Angstroms 由于计算能量时需要有一个截断距离，这个参数指定截断距离。 其余两个参数都是能量最小化时不可用的参数，下面会具体说明。 2） Heating（Equilibrium step1） 上述的参数配置可以归纳如下： imin=0 Choose a molecular dynamics (MD) run no minimization 开始做动力学，而不是能量最小化。 nstlim=10000 Number of MD steps in run (nstlim * dt = run length in ps) 动力学过程的步长，前面说过动力学的原理，它是连续地解牛顿运动方程。动力 学的时长会等于步长（nstlim）乘以每一步的时间间隔（如下，dt） 。 dt=0.002 Time step in picoseconds (ps). The time length of each MD step 每一步的时间间隔，单位是 ps。 ntf=2 Setting to not calculate force for SHAKE constrained bonds ntc=2 Enable SHAKE to constrain all bonds involving hydrogen 以上两个参数针对氢原子做 shake 限制，主要由于氢原子的振动频率过高，氢原子 的运动还存在量子效应。 tempi=0.0 Initial thermostat temperature in K (see NMROPT section) temp0=300.0 Final thermostat temperature in K (see NMROPT section) 以上两个参数指定初始温度和升温后的温度。 ntwx=1000 Write Amber trajectory file mdcrd every ntwx steps 将坐标信息写入轨迹文件的频率。 ntb=1 Periodic boundaries for constant volume 采用周期性边界的恒容条件，表明这是一个 NVT 系综，周期性边界在附录中说明。 ntp=0 No pressure control 暂时不考虑控制压力。 ntt=3 Temperature control with Langevin thermostat gamma_ln=2.0 Langevin thermostat collision frequency 以上两个参数指定如何控制温度。 nmropt=1 NMR restraints and weight changes read (see NMROPT section) 这个参数原本用于指定 nmr 限制，现在被机智地借来用于控制升温过程。 升温过程的控制： 这三行配置表示，在 10000 步的升温中，在最初的 9000 步中， 温度将从 0K 升到 300K，最后 9001 步到 10000 步，温度会保持在 300K。 3） Equilibrium step2、Production 注意：这里配置的动力学时间只有 60ps，主要是为了演示；一 般来说视体系的大小和研究的需要，模拟的时间一般会更长。 第二步平衡（体系升压）与实际动力学的参数配置是一样的，因 为实际动力学便是在恒温恒压下进行的。 上述的参数配置可以归纳如下： ntx=5 Read coordinates and velocities from unformatted inpcrd coordinate file ntx 等于 5 表示将读取另一种格式的 inpcrd 文件，它往往配合 irest=1（如下）使用。 irest=1 Restart previous MD run This means velocities are expected in the inpcrd file and will be used to provide initial atom velocities irest 等于 1 表明将重启上一次的动力学，也就意味着希望能从输入文件中读取到上 一次动力学结束时对应的（除了坐标信息外的）原子速度信息（储存在 rst 文件中） 。 temp0=300.0 Thermostat temperature. Run at 300K 动力学过程中需要控制保持的温度，一般在 300K 左右 ntb=2 Use periodic boundary conditions with constant pressure 采用周期性边界的恒压条件。一般动力学即采用此 NPT 系综。 ntp=1 Use the Berendsen barostat for constant pressure simulation 设定控压的算法。 3.进入 Amber 实际操作 1） min1：限制蛋白，溶剂部分能量最小化 终端运行命令： $AMBERHOME/bin/sander -O -i min1.in -o min1.out p NAME.prmtop -c NAME.inpcrd r min1.rst -ref NAME.inpcrd $AMBERHOME 指定 Amber 软件安装包的目录信息。 NAME.prmtop、NAME.inpcrd 即表示 prmtop 文件盒 inpcrd 文件。 -ref 表示按照给定的坐标信息来限制蛋白，这里按照蛋白的 inpcrd 内的坐标信息。 2） min2：放松蛋白，整个体系能量最小化 终端运行命令： $AMBERHOME/bin/sander -O -i min2.in -o min2.out p NAME.prmtop -c min1.rst -r min2.rst 输入的坐标信息来源于上一步输出的 rst 文件，其储存了坐标信息。 下面的步骤是同理的，上一步输出的 rst 文件将作为下一步输入的坐 标文件。 3） Equilibration step1：（heating：系统在约束蛋白下升 温） 给系统加热，从 0K 到 300K，运行 20ps 动力学（nstlim 的值可以 适当提高，视体系而定）： 终端运行命令： $AMBERHOME/bin/sander -O -i md1.in -o md1.out p NAME.prmtop -c min2.rst -r md1.rst -x md1.mdcrd -ref min2.rst 这一步开始会收集轨迹文件（-x 指定的输出），用于后续可能的分 析。 4） Equilibration step2：（体系升压） 为了平衡体系，将运行 100ps 的恒温恒压的动力学计算： 终端运行命令： $AMBERHOME/bin/sander -O -i md2.in -o md2.out -p NAME.prmtop c md1.rst -r md2.rst x md2.mdcrd 5） Production Run：实际模拟 获取轨迹，将运行 10ns 的分子动力学计算。 一般来说，实际模拟中的配置文件与系统平衡时的配置文件是一致 的： 这里只是”nstlim”这个配置参数的值改变了。 终端运行命令： $AMBERHOME/bin/sander O i prod.in o prod.out p NAME.prmtop c md2.rst r prod.rst x prod.mdcrd 最终我们将对 prod.out、 prod.mdcrd 进行分析、整合，这个将再下一 节中讲述。 4. 不同的 sander 运行命令：集群限制 本阶段可以实现多种运行方式，如串行、并行、多集群并行和 GPU 计算等等。整理执行命令如下。 1） 串行： $AMBERHOME/bin/sander O i mdin o mdout p prmtop c r rst x mdcrd $AMBERHOME/bin/pmemd O i mdin o mdout p prmtop c r rst x mdcrd Amber 中有两个模块程序支持动力学模拟：sander 和 pmemd。 2） 并行： mpirun np N $AMBERHOME/bin/sander.MPI O i mdin o mdout p prmtop c r rst x mdcrd mpirun np N $AMBERHOME/bin/pmemd.MPI O i mdin o mdout p prmtop c r rst x mdcrd N 值为需要加载的服务器的核数，注意，N 值最大值为运行服务 器的 CPU 个数，例如 8、12 等。 3） 多集群并行： mpirun hostfile hosts np N $AMBERHOME/bin/sander.MPI O i mdin o mdout p prmtop c r rst x mdcrd mpirun hostfile hosts np N $AMBERHOME/bin/pmemd.MPI O i mdin o mdout p prmtop c r rst x mdcrd hosts 用于指定并行计算的节点，这种方法多用于跨节点的并行计算 中，一般会提供专门用于提交任务的脚本。 4） 单块 GPU 计算 $AMBERHOME/bin/pmemd.cuda O i mdin o mdout p prmtop c r rst x mdcrd 只有 pmemd 模块支持 GPU。 5.附录： 1）周期性边界条件简介 一个二维盒子排列，分子 1 将从中心盒子（有阴影的盒子）运动到盒子 B 中，这时为了保 持中心盒子的粒子守恒，盒子镜像 F 中有一个相应的粒子会进入到中心盒子中来，就好像 分子 1 从上缘出去，而从下缘进入。 周期性边界条件使得粒子包裹在无穷多的溶液分子中，这样在使 用相对少量的粒子的情况下模拟了真实的生物体系，使动力学模拟 的成本大大减少。假设一个有粒子在内部运动的立方体盒子从各个 方向被复制，形成了周期性排列的盒子，就像上图的二维盒子展示 的一样。 在二维的情况下，每个盒子周围都有 8 个相邻的盒子，而在三维 的情况下，每个盒子周围则会有 26 个相邻的盒子。对于镜像盒子，