使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索

Caroline2022

按照Sob老师方法，今天终于找到了能量最低构象，还学会了做旋转动图，很开心。
Sob老师的帖子真的需要反复阅读，踏实地一步步做下去，终能得到想要的结果。在群里或论坛求助，老师有时间就会回复，再不像以前那么苦恼地一个人想到脑壳疼，连个可探讨的人都没有。每天都进步一点点，真的很有成就感。

Caroline2022

sobereva 发表于 2023-3-11 23:12
当时写这个帖子的时候molclus还不支持调用Shermo
如今强烈建议调用Shermo，对这样的柔性体系用Shermo的Q ...

可是，Sob老师我还有个疑问，在您写这个帖子的时候Molclus还不支持调用Shermo，那么这个高精度自由能（即Gaussian计算得到的自由能热校正量+ORCA计算得到的高精度单点能）就是Molclus直接从Gaussian和ORCA的输出文件中读取的么？那么当我模仿您这个例子做计算时（使用软件版本为Molclus 1.9.9.9, Gaussian 16 C.01,  ORCA 5.0.3， Shermo 2.3.5)，虽然显示Molclus调用了Shermo，但是它所做的工作也就仅仅是原来Molclus自己所作的（读取数据并相加）工作，不知道我这样理解对么？

sobereva

Caroline2022 发表于 2023-3-24 08:23
可是，Sob老师我还有个疑问，在您写这个帖子的时候Molclus还不支持调用Shermo，那么这个高精度自由能（即 ...

“自己所作的（读取数据并相加）工作”语义不明
此时molclus是将电子能量与Shermo算出来的热力学校正量相加，Gaussian或ORCA输出文件里的原有的热力学量完全不被使用

Caroline2022

sobereva 发表于 2023-3-24 08:45
“自己所作的（读取数据并相加）工作”语义不明
此时molclus是将电子能量与Shermo算出来的热力学校正量 ...

对不起老师，应该将“自己”两个字删去，即为“Molclus所作的工作”
那么看来我对第二个问题的理解不对。我再认真学习一下Shermo。

Caroline2022

sobereva 发表于 2023-3-11 23:12
当时写这个帖子的时候molclus还不支持调用Shermo
如今强烈建议调用Shermo，对这样的柔性体系用Shermo的Q ...

Sob老师，我按照您的方法，在setting.ini中加入了如下调用Shermo的语句：
--- Below for Shermo ---
  Shermo_path= "/sob/Shermo_2.3.5/Shermo"  // Command for invoking Shermo
  T= 298.15
  P= 1.0
  sclZPE= 0.9806
  ilowfreq= 2
正常运行结束后看到您帖子中还写道下面的文字：
“如果你要考察其它温度下的Boltzmann分布，不仅isostat里要输入相应的温度，在template.gjf里也应当加上“temperature=[温度]”关键词来让Gaussian输出的自由能热校正量对应相应的温度；或者settings.ini里用ibkchk=1，使得每个opt freq任务的chk都被保留下来，便于之后自行用freqchk工具得到其它温度下的自由能热校正量。单点能从Molclus的ibkout选项备份出来的ORCA输出文件里直接就可以读到。”
由于我计算时忘记在template.gjf里加上“temperature=[温度]”关键词，所以想尝试用freqchk工具得到其它温度下的自由能热校正量。但是当我读完您的这两篇博文http://sobereva.com/221和http://sobereva.com/552后，我觉得Shermo在这里应该是已经达到了使用freqchk工具的效果。不知道我这样理解是否正确？

sobereva

Caroline2022 发表于 2023-3-27 08:57
Sob老师，我按照您的方法，在setting.ini中加入了如下调用Shermo的语句：
--- Below for Shermo ---
   ...

用Shermo算热力学量的时候不需要在gjf文件里写temperature，写了也白写

Caroline2022

sobereva 发表于 2023-3-28 04:28
用Shermo算热力学量的时候不需要在gjf文件里写temperature，写了也白写

谢谢Sob老师，我明白了

Caroline2022

Sob老师，我仿照您的方法计算了一个含Si,B,C,H,O的高分子体系的构象（47个原子），按照保留相对能量<=3 kcal/mol的标准也是有6个构象，对这6个构象计算了高精度自由能后用isostat进行了排序（标准0.5 kcal/mol, 0.5 Angstrom)，构象3的能量由于接近被归纳进了构象1，因此只剩下了5个构象。屏幕输出如下：

[root@bogon molclus_1.9.9.9_Linux]# ./isostat
Isostat: A statistical analysis tool for a batch of isomers (64 bit)
Released as a component of Molclus, last update of isostat: 2022-May-28
Developed by Tian Lu, contact: sobereva@sina.com
Beijing Kein Research Center for Natural Sciences (http://www.keinsci.com)

Number of CPU cores to be used:    6

Input the path of input file, e.g. C:\sob\lovelive\nico.xyz
If press ENTER button directly, isomers.xyz in current folder will be loaded
   To change number of cores used in calculation, input value now, e.g. 18


Input the energy threshold for distinguishing different clusters
e.g. 0.5  (in kcal/mol)
If press ENTER button directly, 0.25 kcal/mol will be used
0.5

Input the geometry threshold for distinguishing different clusters
e.g. 0.25  (in Angstrom)
If press ENTER button directly, 0.1 Angstrom will be used
0.5

Analyzing...

There are    47 atoms
There are     6 isomers

Energy of isomer     1 :  -2030.65956700 a.u.,  new cluster
Energy of isomer     2 :  -2030.65779000 a.u.,  new cluster
Energy of isomer     3 :  -2030.65918400 a.u.,  attributed to cluster     1
Energy of isomer     4 :  -2030.65920000 a.u.,  new cluster
Energy of isomer     5 :  -2030.66040900 a.u.,  new cluster
Energy of isomer     6 :  -2030.65848800 a.u.,  new cluster

The lowest energy is  -2030.66040900 a.u.

Nimag information was found, no isomer has imaginary frequency

Sorting clusters according to energy...
Calculation took up wall clock time         0 s

In below output, DE is energy difference with respect to the lowest one. DGmin quantifies minimal geometry difference with respect to all other clusters

#    1 Count:    1  E=  -2030.660409 a.u.  DGmin=   0.27  DE=    0.00 kcal/mol
#    2 Count:    2  E=  -2030.659567 a.u.  DGmin=   0.20  DE=    0.53 kcal/mol
#    3 Count:    1  E=  -2030.659200 a.u.  DGmin=   0.79  DE=    0.76 kcal/mol
#    4 Count:    1  E=  -2030.658488 a.u.  DGmin=   0.27  DE=    1.21 kcal/mol
#    5 Count:    1  E=  -2030.657790 a.u.  DGmin=   0.20  DE=    1.64 kcal/mol

Representative geometry of each cluster along with their energies have been outputted to cluster.xyz in current folder

If you want to calculate Boltzmann distribution, input temperature in K
If you want to exit, directly press ENTER button
298.15
#    1   Count:   1   Ratio:  53.16%
#    2   Count:   2   Ratio:  21.79%
#    3   Count:   1   Ratio:  14.77%
#    4   Count:   1   Ratio:   6.95%
#    5   Count:   1   Ratio:   3.32%
Press ENTER button to exit
  
下面我想用IRI方法看一下能量最低的3个构象的化学键和分子内相互作用，在转换完chk文件后，我想到一个问题：排序后能量最低的构象（E=  -2030.660409 a.u.），对应的是原来的5号构象（已分别在上面标红），那么我分析的是不是也应该是gau00005.fchk文件呢？

sobereva

Caroline2022 发表于 2023-4-1 10:33
Sob老师，我仿照您的方法计算了一个含Si,B,C,H,O的高分子体系的构象（47个原子），按照保留相对能量

是

zcang1995

Sob老师您好，如图比如研究对象是一系列的连续反应，很难确定每一步的生成物，如果每一步先做高温MD，再用本贴的方法您觉得可以更加省时间合理的做生成物结构的猜想吗？

sobereva

zcang1995 发表于 2023-4-12 13:10
Sob老师您好，如图比如研究对象是一系列的连续反应，很难确定每一步的生成物，如果每一步先做高温MD，再用 ...

这跟本帖没什么关系，这是单纯的从头算动力学的事

跑AIMD并不会比自己尝试计算各种生成物的情况节约时间。
如果你只是用GFN-xTB这种很便宜的级别跑动力学，倒也可以尝试写个脚本，以不同方向不同速度让两个分子碰撞，看都能得到什么

zcang1995

sobereva 发表于 2023-4-12 14:20
这跟本帖没什么关系，这是单纯的从头算动力学的事

跑AIMD并不会比自己尝试计算各种生成物的情况节约时 ...

谢谢老师回复，AIMD想跑可是有的样本Ea比较大，（ > 50 kcal/mol)，看了其他帖子好像AIMD算不了，就想看看用XTB试下...

sobereva

zcang1995 发表于 2023-4-12 14:34
谢谢老师回复，AIMD想跑可是有的样本Ea比较大，（ > 50 kcal/mol)，看了其他帖子好像AIMD算不了，就想看 ...

你没搞清楚内在逻辑

xtb做动力学都是一样的BOMD算法，只不过具体用的理论方法不一样而已。如果xtb里的GFN-xTB的合理的话，别的程序做AIMD跑不出来的xtb做动力学照样跑不出来，顶多是GFN-xTB由于便宜可以跑更长的时间尺度而已，因而有可能在同样模拟耗时下增加得到反应出现的概率，但势垒那么高的情况如果不用加速采样技术照样跑不出来

zcang1995

sobereva 发表于 2023-4-12 15:01
你没搞清楚内在逻辑

xtb做动力学都是一样的BOMD算法，只不过具体用的理论方法不一样而已。如果xtb里的 ...

好的，感谢sob老师!

Caroline2022

温馨提示：
oRCA从5.0后不能出现grid4 gridx4关键词，需要改为defgrid3等关键词