请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 614|回复: 12

[其它程序] 半经典分子动力学程序PROGDYN使用教程

[复制链接]

159

帖子

2

威望

1391

eV
积分
1590

Level 5 (御坂)

发表于 2018-8-14 08:53:50 | 显示全部楼层 |阅读模式
本帖最后由 Accelerator 于 2018-8-14 10:01 编辑

    德州农工的单重态子(D A Singleton)编写的半经典分子动力学程序PROGDYN是研究bifurcation现象以及显式溶剂化等的行业标准。这里介绍一下这个程序(不带ONIOM的版本)的使用。
(典型的bifurcation势能面。这种情况下无法通过TST预测产物比例,只能进行分子动力学模拟)
    PROGDYN目前是基于bash和awk的程序。虽然是公开的,但移植到不同机器上需要修改源代码,并不十分方便。
    在介绍如何使用之前,首先讲解一下PROGDYN的工作流程。该程序需要两个输入文件,freqinHP和progdyn.conf。动力学轨迹的初始结构用Gaussian做频率计算(同时注意写上freq=hpmode保留力常数的更多有效数字),得到的输出文件重命名为freqinHP;progdyn.conf则配置动力学的相关信息。对每一条轨迹,程序从freqinHP中读取初始构型,通过不同的方式初始化后做半经典分子动力学模拟(对每一步求QM的力常数用于下一步),直到轨迹达到终止条件。
2.jpg
  程序的主体位于binall400文件夹内(可以随意重命名);bin里包含很多文件,但progdyn只需要其中的randgen就够了。n400则是与存放输入和输出文件的场所。
1.jpg
    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. echo
  17. echo ORIGDIR at the beginning of run:
  18. echo $origdir
  19. ls $origdir
  20. echo
  21. echo SCRATCHDIR at the beginning of run:
  22. echo $scratchdir
  23. ls $scratchdir
  24. echo
  25. echo PROGGRAMDIR at the beginning of run::
  26. echo $proggramdir
  27. ls $proggramdir

  28. rm -f nogo    # assume that if someone is starting a job, they want it to go.
  29. rm -f diagnostics goingwell tempdone # diagnostics contains extra info from previous runs, other two files are from older versions of progdyn
  30. if (test -s g09.com) then
  31.    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/1fkeddCzzQG4z6frEPAV8uw 密码: ypw9

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


评分

参与人数 6威望 +1 eV +25 收起 理由
agent99 + 5 谢谢分享
ABetaCarw + 5 谢谢分享
kulaomega + 5 谢谢分享
让你变成回忆 + 5 谢谢分享
sobereva + 1
zsu007 + 5 谢谢分享

查看全部评分

1万

帖子

25

威望

1万

eV
积分
35741

管理员

公社社长

发表于 2018-8-14 09:46:34 | 显示全部楼层
有的图片显示不出来
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

159

帖子

2

威望

1391

eV
积分
1590

Level 5 (御坂)

 楼主| 发表于 2018-8-14 10:02:47 | 显示全部楼层
sobereva 发表于 2018-8-14 09:46
有的图片显示不出来

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

1万

帖子

25

威望

1万

eV
积分
35741

管理员

公社社长

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

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

之前的图有的看不到是因为直接复制了微信的图片,微信图片没法外链。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

172

帖子

1

威望

1509

eV
积分
1701

Level 5 (御坂)

那个考了三年没考上的

发表于 2018-8-14 16:26:10 | 显示全部楼层
弱弱问一句
https://thuchemst.github.io
为啥停更了?
恍惚月余,深谙人与人之间的差距。以后还应努力学习,才能与强者比肩。

25

帖子

0

威望

89

eV
积分
114

Level 2 能力者

发表于 2018-8-15 00:57:49 | 显示全部楼层
有点意思

159

帖子

2

威望

1391

eV
积分
1590

Level 5 (御坂)

 楼主| 发表于 2018-8-15 01:37:24 | 显示全部楼层
ABetaCarw 发表于 2018-8-14 16:26
弱弱问一句
https://thuchemst.github.io
为啥停更了?

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

101

帖子

0

威望

879

eV
积分
980

Level 4 (黑子)

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

101

帖子

0

威望

879

eV
积分
980

Level 4 (黑子)

发表于 2018-8-15 07:46:33 | 显示全部楼层
国内好像做quasi classical trajectory做得比较多的是张东辉,他一般是写成准经典轨线。

159

帖子

2

威望

1391

eV
积分
1590

Level 5 (御坂)

 楼主| 发表于 2018-8-15 10:01:35 | 显示全部楼层
Daniel_Arndt 发表于 2018-8-15 07:41
我当初使用的时候,使用的Progdyn还不依赖openbabel。当初使用的时候,有个大问题,就是输出的xyz文件中的 ...

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

159

帖子

2

威望

1391

eV
积分
1590

Level 5 (御坂)

 楼主| 发表于 2018-8-15 10:04:57 | 显示全部楼层
本帖最后由 Accelerator 于 2018-8-15 10:12 编辑

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

101

帖子

0

威望

879

eV
积分
980

Level 4 (黑子)

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

159

帖子

2

威望

1391

eV
积分
1590

Level 5 (御坂)

 楼主| 发表于 2018-8-26 04:34:35 | 显示全部楼层
Daniel_Arndt 发表于 2018-8-26 00:24
这会想起来当时还有一个问题。我当时在progdyn.conf里设置的processors是8,然后我考虑到一个节点是24核, ...

还有这种事情.jpg
不过我也是在同一个节点运行很多任务,还没有遇到过这种情况。
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2018-11-17 07:21 , Processed in 0.122310 second(s), 28 queries .

快速回复 返回顶部 返回列表