半经典分子动力学程序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

黑胡桃

楼主好，我最近用到这个程序，但是出现了报错。我看了一下，应该是初始的构象读取出错了，所有原子都在一个位置。但是我看了freqinHP文件是没有问题的，所有的路径应该也是对的不知道哪里出问题了

youxx189

Accelerator 发表于 2019-3-21 13:29
抱歉才看到，网盘链接已经在ls更新

您好  请问现在有安装包嘛？求！

baobaoxiaoqizi

Accelerator 发表于 2023-3-30 22:09
只需一个。初始速度和初始结构是自动生成的

十分感谢您的回复，理解了。
还有一个问题，我尝试您的算例跑动力学，但是发现总是报错，看了log文件的输出信息，发现内存不够:
Out-of-memory error in routine RdGeom-AllZAt (IEnd=      21250090 MxCore=      20000000)
Use %mem=23MW to provide the minimum amount of memory required to complete this step.
Error termination via Lnk1e in /data/home/baoly/soft/g16/l101.exe at Mon Apr 10 20:09:43 2023.

而且命令行:
%nproc=1
Will use up to    1 processors via shared memory.
%mem=20000000，
这与progdyn.conf文件所设置的内存数和核数不一致，不知道是哪里出现的问题，想向您请教。
[网盘链接：https://pan.baidu.com/s/17_3u5n1NzAG_qpXc1UAeqA  提取码：i43t]

Accelerator

baobaoxiaoqizi 发表于 2023-3-30 20:19
想请教一下，初始的freqinHP频率文件是只包含了一个结构吗？还是要先给过渡态不同的初始速度，生成多个起始 ...

只需一个。初始速度和初始结构是自动生成的

baobaoxiaoqizi

想请教一下，初始的freqinHP频率文件是只包含了一个结构吗？还是要先给过渡态不同的初始速度，生成多个起始结构呢？这部分是怎么操作的呢？

Kikyou

Accelerator 发表于 2022-5-18 17:11
链接: https://pan.baidu.com/s/10sUwJgdqZ4v1tUxaP2ogiQ?pwd=munm 提取码: munm

收到，谢谢老师！

Accelerator

Kikyou 发表于 2022-5-16 18:46
求一下软件的下载链接

链接: https://pan.baidu.com/s/10sUwJgdqZ4v1tUxaP2ogiQ?pwd=munm 提取码: munm

Kikyou

求一下软件的下载链接

sobereva

星斗如盘 发表于 2021-8-10 19:43
使用IOp(1/44=N)，N设置为多少比较合适呢？

随意
这只是个随机数种子，不懂什么叫随机数种子的话google一下便知

星斗如盘

sobereva 发表于 2021-8-10 19:34
ADMP和BOMD的情况完全不同，我前面说的是BOMD
Gaussian里的ADMP没法等效地实现

使用IOp(1/44=N)，N设置为多少比较合适呢？

sobereva

星斗如盘 发表于 2021-8-10 19:20
这个具体怎么操作呢？我使用ADMP，通过shell循环了100次，结果100次的结果都一样。

ADMP和BOMD的情况完全不同，我前面说的是BOMD
Gaussian里的ADMP没法等效地实现

星斗如盘

sobereva 发表于 2020-4-16 13:49
其实用Gaussian的BOMD也可以做同样的准经典动力学模拟，感觉用起来更清爽和省事。只需通过shell循环就可以 ...

这个具体怎么操作呢？我使用ADMP，通过shell循环了100次，结果100次的结果都一样。
以前尝试过用xtb跑动力学，循环100次是可以得到不同的结果，不过使用xtb不能更改计算水平。现在想在b3lyp-D3/6-31g(d)水平下跑一下过渡态

星斗如盘

有PROGDYN编译后的安装包吗？

zwy2020

Daniel_Arndt 发表于 2021-8-7 12:29
程序作者Daniel A. Singleton的使用了progdyn的文章里面会在Supporting Information部分附上代码。

嗯嗯，找到了，还有一些问题，希望请教您