GROMACS的使用(二、自由能计算、加电场模拟) GROMACS的使用(二、自由能计算、加电场模拟) 日期:2018-06-25 标签: 阅读: 简介:在上一篇中,我们介绍了一个比较基础的应用,如果你已经掌握了GROMACS的基础命令和步骤,在此基础上来进一步学习一下GROMACS在计算自由能和加电场的模拟。

            在上一篇中,我们介绍了一个比较基础的应用,如果你已经掌握了GROMACS的基础命令和步骤,在此基础上来进一步学习一下GROMACS在计算自由能和加电场的模拟。


1. 加电场模拟

对于最基础的长方体模拟盒子,在某一个方向上(x或y或z)施加一个恒定强度的电场,只需要在mdp文件中添加下面一行:

 

E-z  =  1  1.8  0


        其中E-z表示在z方向上加电场;等号右边第一个数是余弦的数目,因为只实现了单个余弦项(频率为0),因此为1,不能写成1.0;第二个数为电场大小和方向,正数表示从z=0到z=max,负数则相反,所以在如下图中的体系中施加电场时,应当添加负的电场,或者一开始就把多肽放在生物膜的下方,电场的单位是V/nm;第三个数字是余弦的相位,因为频率为0的余弦没有相位,所以可以填写任意数字。

 

1.jpg

图1. 模拟体系

​另外一种常用电场为高斯脉冲电场,但是在gromacs的calc_f_el函数的代码中并没有高斯脉冲电场的数据类型,需要自己添加:

 

typedefstruct { int n; /* Number of terms*/ real *a; /* Coeffients (V / nm )                     */ real *phi; /* Phase angles*/ } t_cosines; typedefstruct { real E0; /* Field strength (V/nm)                        */ real omega; /* Frequency (1/ps)                             */ real t0; /* Centre of the Gaussian pulse (ps)            */ real sigma; /* Width of the Gaussian pulse (FWHM) (ps)      */ } t_efield;


在mdp文件中加入下面两行:

 

E-z   =  1  5  0

E-zt  =  1  50  0 


在真正模拟之前,最好先使用NVE系综进行测试电场对于体系的影响,即在不加电场时在关心的时间尺度内能否保证能量守恒,也可以根据能量漂移的速率来决定,大致小于相应温度下的平动能就能接受。能量守恒测试通过后,再加入电场,看其对体系总能量的影响。真正的模拟体系仍然为nvt或npt体系,相应命令见上一篇GROMACS应用一。

 

注意事项:由于施加的电场一直存在,故在高电场下,磷脂膜等受到破坏之后,电势并没有降低,因此在此情况下,磷脂膜不会自愈。如果模拟一种物质跨膜,应注意是真的跨膜,还是通过周期边界条件而移动到膜的下方。

2.jpg

图2. 加电场进行md后模拟体系。图中只显示了水孔处某一薄层切片。


2. 自由能计算

我们仍然以多肽与膜的体系为例,在进行完npt平衡模拟之后,给多肽施加一定的作用力,使其匀速运动,产生一系列多肽位于不同位置的构型。以这些构型为初始构型,对质心加以限制,使其发生振动和转动,最终达到一个平衡。

(1)产生构型 

首先,在mdp文件中添加下列牵引代码

; Pull code

pull                 = yes   ;激活牵引代码

pull_ngroups         = 2     ;受伞形偏离势的作用的组

pull_ncoords         = 1     ;反应坐标数目

pull_group1_name     = Chain_A   ;索引文件中的第一个组

pull_group2_name     = Chain_B   ;索引文件中的第二个组

pull_coord1_type      = umbrella   ;第一个反应坐标用简谐势牵引

pull_coord1_geometry  = distance   ; 简单的距离增加 simple distance increase

pull_coord1_groups    = 1 2        ;由1组和2组定义反应坐标

pull_coord1_dim      = N N Y      ;只在z方向上牵引

pull_coord1_rate      = 0.01       ; 牵引组虚拟粒子的牵引速度,单位是 nm·ps-1 

pull_coord1_k        = 1000       ; 牵引力常数,单位是kJ·mol-1nm-2

pull_start            = yes        ; 初始的质心距离是第一帧的参考距离

然后定义分组:

 

make_ndx   -f   npt.tpr   -o   pull.ndx


如果拉取整个多肽,可以选择多肽为一组,磷脂为另一组。实际上这样的分组可以不用重新生成pull.ndx,因为在前面的ndx文件中本来就存在。如果是别的情况,应当选择施加牵引的原子或者基团分为一组,将参考相分为一组。

以npt平衡后的构型文件为初始构型文件,进行md模拟。模拟结束后抽取一些窗口进行下一步的模拟,每隔0.1nm左右抽取一个窗口。间隔不易太小或太大,如果太小则模拟耗时较多且相邻两个窗口之间变化不大;如果太大则后面可能两个窗口位移不相交。为了方便查看每一个构型中两组的相对位置,gromacs教程中给出了一个脚本,只需要执行这个脚本,就可以生成一个距离文件distances.dat:

 

perl   distances.pl


根据distances.dat中的信息,可以选取n个构型。

 

gmx  trjconv  -s  pull.tpr  -f  traj.xtc  -o  conf.gro  –sep


这一步将会产生所有的输出构型文件,根据上面选取出第i个文件。

  

3.jpg

图3. PMF模拟前后多肽位置的变化。

 

(2)伞形抽样模拟

将选择的n个构型文件分别作为初始构型进行伞形抽样模拟。在模拟之前,先进行短暂的npt平衡,再进行md模拟。在这两个过程中mdp文件与上面的类似,将牵引速度改为0即可。这意味着仅仅施加了一个弹簧力,使牵引组分尽量维持在原位置,但是该组分可以发生上下振动,自身的转动调整。这样经过长时间模拟,体系趋于稳定。但与之前的模拟不一样的是在脚本里面改为:

 

-deffnm  umbrella0  -pf  pull-umbrella0.xvg  -px  pullx-umbrella0.xvg

 

(3)数据分析

建立两个文件:tpr-files.dat和pullf-files.dat,在tpr-files.dat中写入所有的伞形抽样模拟的tpr文件名,在pullf-files.dat中写入所有的pullf.dat文件名,如图所示.

 

4.jpg

5.jpg


gmx  wham  -it  tpr-files.dat  -if  pullf-files.dat  -o  -hist  -unit  kCal


其中单位为kCal·mol-1,也可以换成kJ·mol-1或者kBT

下图为输出结果样图(图片来自教程),左图为平均力势的曲线变化图,右图为伞形直方图。根据伞形直方图来看,每两个相邻窗口都必须相交,否则平均力势的曲线是有突变的,因此也是不正确的。

6.jpg

 

7.jpg