半经典分子动力学程序PROGDYN使用教程 - 计算化学公社

（典型的bifurcation势能面。这种情况下无法通过TST预测产物比例，只能进行分子动力学模拟）

PROGDYN目前是基于bash和awk的程序。虽然是公开的，但移植到不同机器上需要修改源代码，并不十分方便。

在介绍如何使用之前，首先讲解一下PROGDYN的工作流程。该程序需要两个输入文件，freqinHP和progdyn.conf。动力学轨迹的初始结构用Gaussian做频率计算（同时注意写上freq=hpmode保留力常数的更多有效数字），得到的输出文件重命名为freqinHP；progdyn.conf则配置动力学的相关信息。对每一条轨迹，程序从freqinHP中读取初始构型，通过不同的方式初始化后做半经典分子动力学模拟（对每一步求QM的力常数用于下一步），直到轨迹达到终止条件。

(, 下载次数 Times of downloads: 203)

程序的主体位于binall400文件夹内（可以随意重命名）；bin里包含很多文件，但progdyn只需要其中的randgen就够了。n400则是与存放输入和输出文件的场所。

(, 下载次数 Times of downloads: 212)

binall400中包含一系列文件。其中freqinHP是输入文件，progdynstarterHP是程序的入口；progdyn.conf是混入的奇怪的东西不需要这个。移植这一程序时，首先需要确保机器上安装了openbabel。接下来修改progdynstarterHP的以下几处：

#export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:'/home/ymma/mopac':'~/bin/lib'
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:'/home/accelerator/progdyn/bin/lib'
#export BABEL_LIBDIR='/scratch/user/ymma/babel/build/lib'
#export BABEL_DATADIR='/scratch/user/ymma/babel/openbabel-2.4.1/data'
export LC_ALL=C
echo $1
scratchdir=$1
export g09root=/home/accelerator
. $g09root/g09/bsd/g09.profile
origdir=`pwd`
cd $origdir
logfile=ymmalog
randdir=/home/accelerator/progdyn/bin
proggramdir=/home/accelerator/progdyn/binall400
freqfile=/home/accelerator/progdyn/binall400/freqinHP
echo
echo ORIGDIR at the beginning of run:
echo $origdir
ls $origdir
echo
echo SCRATCHDIR at the beginning of run:
echo $scratchdir
ls $scratchdir
echo
echo PROGGRAMDIR at the beginning of run::
echo $proggramdir
ls $proggramdir
rm -f nogo # assume that if someone is starting a job, they want it to go.
rm -f diagnostics goingwell tempdone # diagnostics contains extra info from previous runs, other two files are from older versions of progdyn
if (test -s g09.com) then
sed -i '/guess=tcheck/d' g09.com # no chk file on first point

复制代码

由于PROGDYN调用Gaussian进行计算，自然要配置相应的环境变量。randdir即为之前提到的随机数生成器randgen的位置；proggramdir为binall400文件所在地，freqfile指向freqinHP文件。如果需要同时运行很多PROGDYN任务的话一定要注意这两个变量是否正确。

如果在集群上运行的话还需要配置一下babel的相关信息，这里不需要所以注释掉了。

n400目录只需要一个文件即progdyn.conf，指定了轨迹的计算方法，自旋多重度和电荷数，并行数和内存占用，初始化方式，apply force等许多重要信息。其中关于每一条设置的含义都有详细的注释。画风是这样的：

#***method, charge, multiplicity, memory, processors, title
#*** method --The following word is copied exactly to the gaussian input file.
method MP2/cc-pvdz
#To do a nonstandard route, make nonstandard 1. For normal calcs, use nonstandard 0 or else leave it out.
#Then make a file called "nonstandard" containing the nonstandard route with no extra lines.
nonstandard 0
# NMRoptions As is NMRtype=1 will add a section for an NMR calc at every NMRevery intervals. If you want to combine the two use nonstandard
#NMRtype 1
#NMRmethod2 B97D/6-31G*
#NMRmethod LC-wPBE/6-31G*
#NMRmethod3 B3LYP/cc-pvtz
#NMRevery 4
#NMRrand 1
#NMRcc 1
#loadlimit 10.0
#geometry linear
rotationmode 1
#*** method2 --The options here are restricted, unrestricted, and read. restricted is the default
#If the method is U..., put unrestricted here and the .com files will have in them guess=mix.
#If you put read here, the .com files will contain guess=tcheck, which sometimes makes things faster, sometimes not.
#The use of read requires a specifically defined checkpoint file name using the keyword checkpoint.
method2 restricted
charge 0
multiplicity 1
#oniomchargemult 1 1
processors 3
#*** memory --The following "word" is copied exactly to the gaussian input file after %mem=.
memory 7gb

复制代码

最基本的是修改计算水平，一定要确保和freqinHP一致，否则每一步的能量都与freqinHP中读取的能量相差太大，会被判定为结构不合理而中止运行。

还需要修改的是proganal文件。该文件用于实时分析轨迹。

<font size="2">C2C6=Distance(2,6)
C3C26=Distance(3,26)
C5C7=Distance(5,7)
printf("%s %.3f %s %.3f %s %.3f ","C2C6",C2C6,"C3C26",C3C26,"C5C7",C5C7)
if (runpoint>500) {
print " Too many points. XXXX"
}
if ((C2C6<1.7) && (C5C7<1.7)) {
print "Vinylfuran is dienophile XXXX"
}
if ((C2C6<1.7) && (C3C26<1.7)) {
print "Cyclopentadienone is dienophile XXXX"
}
if (C2C6>2.5) {
print "Returned to SM XXXX"
}
system("date '+%b:%d:%Y %T'")
system("tail -1 Echeck | grep XXXX")
}</font>

复制代码

需要修改的是这一段，即轨迹的终止条件。这里的例子中freqinHP是一个DA反应的过渡态，proganal中测量三个C-C键键长，在键长小于临界值时认为已经运行到产物形成阶段并结束该轨迹。

修改好这些之后就可以提交运算了。首先cd到n400内，新建一个tmp文件夹存放程序的临时文件，并运行

（绝对路径的前项省略）/binall400/progdynstarterHP tmp

如果没有报错的话，至少会在n400里生成以下文件：

g09.com geoPlusVel geoRecord Echeck dynfollowfile runpointnumber isomernumber

其中isomernumber指的是“当前是第几条轨迹”。一条轨迹完成后程序并不会停止运行，而是继续从初始构型开始生成更多的轨迹。geoPlusVel记录了当前的几何构型和各个原子的速度；geoRecord记录了起始构型的几何结构和速度；dynfollowfile则按照proganal指定的输出轨迹信息。画风如下：

根据自己书写的终止条件，可以在dynfollowfile里查找已经结束的轨迹情况。

bad total energy的含义是轨迹的第二步电子能量与起始构型相差太大，被判定为发生了错误而重新开始下一条轨迹。在progdyn.conf中有一个initialdist的选项，指定轨迹的初始化方式。在这里除了过渡态的虚频模式之外，其他振动模式也都给了初始位移（选项4），因此发生bad total energy是很正常的。 initialdist的设置对于轨迹结果至关重要，默认设置为4，是物理意义以及与实验的相符都较好的选择。

由于Gaussian的并行效率惩罚，并不会只运行一个PROGDYN任务。前面展示的progdyn.conf中CPU和内存都设置得很小，正是因为同时运行着许多同类任务。Singleton习惯在一个20核心节点上同时运行6个任务。

得到大量轨迹后最简单的应用就是用于得到产物比例的信息。本例中可能得到两个产物（Vinylfuran和Cyclopentadienone分别作为亲双烯体）。得到的通往两种产物的轨迹比例为1.47，与实验值1.6符合可以接受。

PROGDYN的代码一般在使用过的文章的SI里都有；这里提供下载：链接: https://pan.baidu.com/s/10sUwJgdqZ4v1tUxaP2ogiQ?pwd=munm 提取码: munm

本文首先发表于作者的公众号Rita是美少女上所以图片上经常有水印(x