半经典分子动力学程序PROGDYN使用教程

Accelerator

    德州农工的单重态子（D A Singleton）编写的半经典分子动力学程序PROGDYN是研究bifurcation现象以及显式溶剂化等的行业标准。这里介绍一下这个程序（不带ONIOM的版本）的使用。

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

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

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

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

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

1. #export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:'/home/ymma/mopac':'~/bin/lib'
   
2. export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:'/home/accelerator/progdyn/bin/lib'
   
3. #export BABEL_LIBDIR='/scratch/user/ymma/babel/build/lib'
   
4. #export BABEL_DATADIR='/scratch/user/ymma/babel/openbabel-2.4.1/data'
   
5. export LC_ALL=C
   
6. echo $1
   
7. scratchdir=$1
   
8. export g09root=/home/accelerator
   
9. . $g09root/g09/bsd/g09.profile
   
10. origdir=`pwd`
    
11. cd $origdir
    
12. logfile=ymmalog
    
13. randdir=/home/accelerator/progdyn/bin
    
14. proggramdir=/home/accelerator/progdyn/binall400
    
15. freqfile=/home/accelerator/progdyn/binall400/freqinHP
    
16. 
    
17. echo
    
18. echo ORIGDIR at the beginning of run:
    
19. echo $origdir
    
20. ls $origdir
    
21. echo
    
22. echo SCRATCHDIR at the beginning of run:
    
23. echo $scratchdir
    
24. ls $scratchdir
    
25. echo
    
26. echo PROGGRAMDIR at the beginning of run::
    
27. echo $proggramdir
    
28. ls $proggramdir
    
29. 
    
30. rm -f nogo    # assume that if someone is starting a job, they want it to go.
    
31. rm -f diagnostics goingwell tempdone # diagnostics contains extra info from previous runs, other two files are from older versions of progdyn
    
32. if (test -s g09.com) then
    
33.    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等许多重要信息。其中关于每一条设置的含义都有详细的注释。画风是这样的：

1. #***method, charge, multiplicity, memory, processors, title
   
2. #*** method --The following word is copied exactly to the gaussian input file.
   
3. method MP2/cc-pvdz
   
4. #To do a nonstandard route, make nonstandard 1.  For normal calcs, use nonstandard 0 or else leave it out.
   
5. #Then make a file called "nonstandard" containing the nonstandard route with no extra lines.
   
6. nonstandard 0
   
7. # 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
   
8. #NMRtype 1
   
9. #NMRmethod2 B97D/6-31G*
   
10. #NMRmethod LC-wPBE/6-31G*
    
11. #NMRmethod3 B3LYP/cc-pvtz
    
12. #NMRevery 4
    
13. #NMRrand 1
    
14. #NMRcc 1
    
15. #loadlimit 10.0
    
16. #geometry linear
    
17. rotationmode 1
    
18. #*** method2 --The options here are restricted, unrestricted, and read. restricted is the default
    
19. #If the method is U..., put unrestricted here and the .com files will have in them guess=mix.
    
20. #If you put read here, the .com files will contain guess=tcheck, which sometimes makes things faster, sometimes not.
    
21. #The use of read requires a specifically defined checkpoint file name using the keyword checkpoint.
    
22. method2 restricted
    
23. charge 0
    
24. multiplicity 1
    
25. #oniomchargemult 1 1
    
26. processors 3
    
27. #*** memory --The following "word" is copied exactly to the gaussian input file after %mem=.
    
28. memory 7gb

复制代码

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

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

1. <font size="2">C2C6=Distance(2,6)
   
2. C3C26=Distance(3,26)
   
3. C5C7=Distance(5,7)
   
4. printf("%s %.3f   %s %.3f   %s %.3f  ","C2C6",C2C6,"C3C26",C3C26,"C5C7",C5C7)
   
5.    if (runpoint>500) {
   
6.       print "   Too many points.  XXXX"
   
7.       }
   
8.    if ((C2C6<1.7) && (C5C7<1.7)) {
   
9.       print "Vinylfuran is dienophile XXXX"
   
10.       }
    
11.    if ((C2C6<1.7) && (C3C26<1.7)) {
    
12.       print "Cyclopentadienone is dienophile XXXX"
    
13.       }
    
14.    if (C2C6>2.5) {
    
15.       print "Returned to SM XXXX"
    
16.       }
    
17.    system("date '+%b:%d:%Y %T'")
    
18.    system("tail -1 Echeck | grep XXXX")
    
19.    }</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

sobereva

有的图片显示不出来

Accelerator

sobereva 发表于 2018-8-14 09:46
有的图片显示不出来

我这里看正常。刚刚重新上传了一下 一共是5个图
发帖和编辑都要审核好不方便...

sobereva

Accelerator 发表于 2018-8-14 10:02
我这里看正常。刚刚重新上传了一下 一共是5个图
发帖和编辑都要审核好不方便...

说明内容有敏感词
在论坛里发山寨广告的太多，之前还有机器人大规模轰炸过本论坛，不得不设了很多敏感词，碰到敏感词的时候就需要管理员审核

之前的图有的看不到是因为直接复制了微信的图片，微信图片没法外链。

ABetaCarw

弱弱问一句
https://thuchemst.github.io
为啥停更了？

pyscf

有点意思

Accelerator

ABetaCarw 发表于 2018-8-14 16:26
弱弱问一句
https://thuchemst.github.io
为啥停更了？

因为去年的我懒 然后今年的我不再负责这个了(
ps 如果你关注了清化科协的话就会发现清化科协也停更了(不是

Daniel_Arndt

我当初使用的时候，使用的Progdyn还不依赖openbabel。当初使用的时候，有个大问题，就是输出的xyz文件中的电子能有些错位。具体来说，就是第一帧的电子能出现在第二帧中，第二帧的电子能出现在第三帧中，依此类推。还有个小问题，就是原子数目超过50时，高斯的log文件的格式有点变化，然后需要自己手动修改Progdyn的代码，否则出错。不知道这两个问题现在有没有解决。

Daniel_Arndt

国内好像做quasi classical trajectory做得比较多的是张东辉，他一般是写成准经典轨线。

Accelerator

Daniel_Arndt 发表于 2018-8-15 07:41
我当初使用的时候，使用的Progdyn还不依赖openbabel。当初使用的时候，有个大问题，就是输出的xyz文件中的 ...

感谢提醒。检查了一下，输出的轨迹里单点能仍然是错位的。超过50个原子的体系还没有尝试过。
查看了一下代码，应该只是在书写traj的时候把旧的单点和progdynb生成的新的构型写到了一起，并不会影响轨迹的正确性。前面在计算动能时还有一句注释“This is for the point prior to the current point”，看来Singleton是知道这件事情的。

Accelerator

借楼赞美一下Singleton其人；LZ现在在这里做暑研，快要70岁了还每天亲自写代码，手把手教学生，到现在已经拿到5个MOOC证书了。组会就像放课后茶会一样，感觉这里的学生工作真的是出于爱去做的，并不是因为怕老板什么的。我超喜欢这里的。

Daniel_Arndt

这会想起来当时还有一个问题。我当时在progdyn.conf里设置的processors是8，然后我考虑到一个节点是24核，就写了个很简单的bash脚本，在一个节点上同时跑三个轨线。结果，三个文件夹里出现的第一个xyz文件完全相同。我后来发现是randgen在同一秒中内生成的随机数完全一样，awk在同一秒内的随机数也完全一样。其实也不是什么大问题，就是处理数据时删掉两个文件即可。我后来自己给randgen加了个锁。

Accelerator

Daniel_Arndt 发表于 2018-8-26 00:24
这会想起来当时还有一个问题。我当时在progdyn.conf里设置的processors是8，然后我考虑到一个节点是24核， ...

还有这种事情.jpg
不过我也是在同一个节点运行很多任务，还没有遇到过这种情况。

chemqzx

网盘分享没有了，方便的话可以上传一下吗？谢谢，十分感谢

pyscf

请问这个程序的公开链接在哪？好想找不到了
网盘链接也失效了。。。。