由于专业原因,执行分子动力学模拟一般不会仅仅局限于标准的氨基酸、脂质、核酸之类的分子。更多的时候会去模拟分子力场所没有定义的小分子,执行蛋白-配体复合物的动力学模拟。

NAMD基础

一般建议从NAMD的官方基础教程开始,从UIUC处下载NAMD Tutorial,非常建议同时下载tutorial files用于动手实际操作一下。本节讨论一下该教程不会论述的东西——把NAMD安装上,当然也会提一下其他需要下载的资源。

VMD、NAMD工具链

VMD和NAMD是联系非常紧密的两个工具,VMD是很重要的一个可视化工具,同时也是构建NAMD计算任务的一个必要工具。从UIUC下载即可。VMD 1.9.3版本之后就可以使用QwikMD,这是一个非常好用的动力学模拟任务生成插件,遗憾的是在处理CHARMM36力场且含配体小分子的任务中,不太好用。另外关于NAMD,有且仅有Nvidia图形卡加速的选项,亦即编译带CUDA的。需要注意的是,NAMD要求CUDA的Compute Capability等级不小于3.x。

Hardware Generation Compute Capability
Turing 7.5
Volta 7.x
Pascal 6.x
Maxwell 5.x
Kepler 3.x
Fermi 2.x

因此需要确认显卡是否为开普勒及以上架构的,费米架构虽然能装上CUDA,但会由于等级不足,无法进行计算。

获得CHARMM36力场文件

虽然讲道理VMD 1.9.3之后QwikMD会自带CHARMM36力场文件,不过姑且提一下,对于常见的氨基酸、脂质以及核酸等的力场文件,可以通过MacKerell Lab的网站获得。而对于小分子,通常是不会出现在CHARMM36力场文件中的,因此需要额外获得,通常文献里会提到 SwissParam 这个网站,但总感觉得到的文件和CHARMM36格格不入。因此建议使用CGenFF或者CHARMM-GUI

构建计算任务的注意点

要构建一个能跑起来的计算任务倒也不是很难的问题,照着tutorial files一通抄就行了,只不过有那么几个小问题需要说明一下。

小分子力场文件微调

从CGenFF得到的小分子力场文件是str做后缀的Stream文件,本质上是拓扑文件和参数文件的混合体,因此在使用psfgen以及写NAMD计算文件的过程中,拓扑文件和参数文件都指向该Stream文件即可。psfgen并不会说Stream文件最后结尾处的Return是语法错误,但NAMD会,并且会把它当成一个Error,而不是一个Warning。因此需要把Stream文件最后一行 Return 删掉。

引入所有参数文件

由于CHARMM36的全新设计,在NAMD中使用Stream文件,需要引入所有参数文件(详见说明)。因此,无论是否体系中包含氨基酸、脂质或者核酸,都需要在NAMD计算任务文件中,引入Stream文件之前,引入其他所有参数文件,即

...
structure system-xplor.psf
paraTypeCharmm on
parameters par_all36_prot.prm
mergeCrossterms yes
parameters par_all36_lipid.prm
parameters par_all36_carb.prm
parameters par_all36_cgenff.prm
parameters molecule.str
parameters toppar_water_ions_namd.str
...

QwikMD报错

这里讨论的报错需要有这么一个前件:小分子力场参数由CGenFF生成。在这种情况下QwikMD提示子程序非正常退出,查看NAMD日志发现力场参数缺失。一般情况下,这是由于QwikMD自带的CHARMM36力场参数文件过于陈旧,重新下载CHARMM36,将力场参数文件覆盖进工作目录即可。

面向高效计算的Tech Tips

多条短轨迹的简易做法

通过QwikMD生成过一个计算任务目录后,用NAMD简单做个验证,看看计算任务是否能正常跑起来,如果能正常跑起来那就可以开始后续的工作了。一般来说分子动力学模拟主要分两个流派:1.单条长轨迹(通常是μs级,长的可以到达ms级);2.多条短轨迹(通常每条轨迹只有数十ns)。个人觉得两种做法各有特色,适用的理论框架可能都不太一样,但无论怎么样,涉及到化学键断裂与生成的酶促反应过程那必然选择多条短轨迹——毕竟常规的分子动力学模拟根本不能描述共价过程,长轨迹毫无意义。

于是问题就只剩怎么方便地得到多条短轨迹。如果没有Minimization这一步的话,亦即QwikMD中Easy Run生成的任务,那就直接从退火开始重复执行就行了。参考shell命令如下

#Suppose AmazingComplex as job directory
#Generate 9 replicas, so 10 replicas in total

export namd=$(which namd2)
for i in {1...9};do cp -r AmazingComplex AmazingComplex$i;done
for i in $(ls);do \
namd $i/job/run/qwikmd_equilibration_0.conf > $i/job/run/qwikmd_equilibration_0.log && \
namd $i/job/run/qwikmd_production_1.conf > $d/job/run/qwikmd_production_1.log;done

当然如果有Minimization这一步,也就是通常得到的分子对接构象没有经过后优化(常见于使用开源分子对接工具链),一般建议做完Minimization后再执行上述的命令。因为Minimization步通常都是确定性算法,理论上得到结果是一致的。