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

yyyyss11

mutron 发表于 2020-7-18 04:30
试了一下ORCA里面做MD可以固定住一个二面角，所以应该没问题可以去试一下，ORCA的MD有更丰富的选项。至少 ...

我去试试orca的MD功能，谢谢指导。

qinzhong605

社长，我拿您的文件包，xtb最新的6.3.2跑，第一步只有1042帧，不太对呢，不是应该有2000帧吗？对xtb不太熟悉，盼望指点一下，谢谢了。@sobereva

qinzhong605

知道怎么回事了，社长传的文件包里md.in应该为md.inp，初次接触，没注意到这个问题。

mutron

1.老师，第二步用GFN0-xTB优化会导致漏掉一些能量更低的构象。如果跳过第二步直接第三步用GFN2-xTB在GBSA溶剂模型下优化会多出7个构象，能量相对于文中第三步得到的最低构象（-128.462874 a.u.）会再低0.32-1.91 kcal/mol。（我觉得第二步的必要性不是太强，能省出来的计算时间很有限，尤其是结合crest批量优化）



2.crest能够获得电子能量更低的一批构象。我看到老师21楼给出的crest结果“能量最低的至少前50个结构里没有任何一个和真实的能量最低结构有相似性”，是否可能只是crest得出的一批能量更低的构象没有一个和老师的结果有重复？crest得出的一批构象中最低的能量（-128.469685）比文中第三步得到的最低构象（-128.462874）低了4.26 kcal/mol



3.GFN2-xTB的能量精度还是不够，用来排序得保留相对能量比较宽的一批构象（见下），碰上如本例柔性不小的分子下一步DFT优化需要处理的结构太多算不动。所以中间得加上一步比较靠谱又便宜的能量计算，比如老师推荐的B97-3c，或者RI-PWPB95-D3(BJ)/def2-TZVP（-f），获取相对准确的能量排序再筛选

比如crest得到的能量最低构象，其GBSA下的GFN2-xTB能量：-128.469685，B97-3c结合SMD模型的能量：-2320.74009298；用B97-3c重新排序的能量最低构象分别为：-128.465825和-2320.750725。GFN2-xTB能量后者比前者高了2.42 kcal/mol，B97-3c能量则反过来低了6.67 kcal/mol（RI-PWPB95-D3(BJ)/def2-TZVP(-f)则是低了4.92 kcal/mol），误差达到7-9 kcal/mol。由crest通过后续步骤算出来能量最低的构象在GFN2-xTB排序中排第141（B97-3c排序排第2），相对默认的GFN2-xTB能量排序第1位GFN2-xTB能量高了3.78 kcal/mol，也就是说如果按默认的GFN2-xTB排序，后续DFT计算需要算200多个构象才能把最稳定构象包括进去

RI-PWPB95-D3(BJ)/def2-TZVP(-f)和RI-PWPB95-D3(BJ)/def2-QZVPP在本例相对能量相差大概零点几kcal/mol，其和B97-3c对于crest能量最低的17个构象是一样的，只是顺序有变化



4.isostat结构去重的时候，结构去重阈值(埃)选取0.5可能会有点大，导致相当一部分构象被合并了。比如crest经过B97-3c能量排序后最低的17个构象（计算过程用0.2(埃)），如果用0.5(埃)会导致其中4个构象消失，但观看结构会发现烃基扭动方式和取向明显不一样。又比如拿crest得出的所有构象和文中第三步能量最低的四个构象作对比，如果选取0.5(埃)，isostat会给出3个结构有相似，但仔细观看结构会发现差得不少。如果选取选取0.2甚至0.3(埃)则不会有相似结果



5.使用mdopt对MD或者crest结果做结构优化（normal opt或者vtight）、isostat去重后，会得到好几百甚至上千个相对能量在6kcal/mol以内的构象，哪怕用很便宜比较靠谱的方法算一次能量也耗时颇久。此外，第四步用Gaussian优化过后会出现不少相同的结构（比如文中第四步里面的第5个和第2个构象），DFT优化浪费的时间就更多了。如果在xtb优化完加上一步，使用gau_xtb计算力常数再优化一遍能去掉不少相同结构，减轻后面DFT计算的负担（我觉得通过计算力常数收敛到同一个构象，相对于提高阈值判定为相同构象会更稳妥）



6.几种方法最后得到的最低能量结构的自由能（Shermo考虑Grimme的准RRHO模型）：

-2320.82217947   0.00 kcal/mol   N/A        本文结果
-2320.82816084  -3.75                2（a）   本文第三步后使用B97-3c排序
-2320.82719758  -3.15                10         MD_GFN0-xTB_500ps_500K
-2320.83092026  -5.48                10         MD_GFN2-xTB_GBSA_500ps_400K
-2320.82594269  -2.36               N/A        crest_GFN2-xTB_GBSA排序
-2320.82994678  -4.87                20         crest_B97-3c排序
-2320.82986194  -4.82（MD）    20+22    crest能量最低10个结构继续用MD搜索
-2320.83013936  -4.99               20+25    crest能量最低10个结构继续用MTD搜索

（a）：GFN2-xTB+GBSA溶剂模型下得到构象、在RI-PWPB95-D3(BJ)/def2-TZVP(-f)结合SMD模型水平下获得的电子能量对比，比本文最稳定构象能量低的数量。反映最低能量区域附近的采样充分程度

从结果看来，因为老师在26楼谈到的原因，做完crest如果还需要对真正能量最低结构区域采样充分、找到真正最稳定构象，有必要再做一遍MD或者MTD，具体条件（构象附近搜索）需要摸索



7.老师可否考虑在molclus里面加上批量做MD的功能，需要做MD（MTD）的构象一多就开始变得麻烦了



8.molclus 1.9.3版本在计算itask=2、3且没有提供单点任务模板文件时会给出自由能加多一遍Gcorr的错误结果，但最新版1.9.5好像吉布斯自由能G总是=0了，不知道我哪里理解错还是程序问题



暂时找到的最稳定构象的NCI图：

既有图1的内氢键，pi-pi堆积也和图2差不多

sobereva

mutron 发表于 2020-7-27 00:01
1.老师，第二步用GFN0-xTB优化会导致漏掉一些能量更低的构象。如果跳过第二步直接第三步用GFN2-xTB在GBSA溶 ...

前一阵子已经发布了molclus 1.9.6版

本文只是作为演示用，需要更严格的最小点把MD轨迹跑得更长、从而采样更充分即可，简单粗暴但很有效。虽然在细节上做各种改进可能能让当前体系搜索效率更高，但不一定有普适性。

tjuptz

sobereva 发表于 2020-4-7 06:59
先跑一下，产生mdrestart文件，然后再篡改mdrestart文件内容，通过restart=true来读入初速度（PS：ORCA的 ...

请问老师
1 二进制的mdrestart用文本编辑器如notepad++的二进制插件 打开可以嘛？
2 mdrestart里不是已经记录了之前md的速度吗？
3 对于其他体系也是加到最后就可以了吗？还是跟输入坐标一致的？

sobereva

tjuptz 发表于 2020-8-28 15:46
请问老师
1 二进制的mdrestart用文本编辑器如notepad++的二进制插件 打开可以嘛？
2 mdrestart里不是已 ...

xtb的mdrestart是文本文件，直接用文本编辑器打开
是记录了，你可以篡改
里面记录的速度矢量数目得和原子数对应

tjuptz

sobereva 发表于 2020-8-29 00:53
xtb的mdrestart是文本文件，直接用文本编辑器打开
是记录了，你可以篡改
里面记录的速度矢量数目得和原 ...

果然是，之前想错了。谢谢老师

naoki

Sob老师，我想咨询一下xtb做动力学可以用显式溶剂吗，比如水和离子液体，如果一样跑的话，盒子里溶剂一多是不是计算耗时巨大无比

sobereva

naoki 发表于 2020-9-15 16:32
Sob老师，我想咨询一下xtb做动力学可以用显式溶剂吗，比如水和离子液体，如果一样跑的话，盒子里溶剂一多是 ...

虽然可以用，但一方面会显著增大耗时，另一方面xtb自身没法控压，也没法达到理想地模拟溶液状态的效果。非要用显式溶剂应当考虑gromacs跑

naoki

sobereva 发表于 2020-9-16 09:07
虽然可以用，但一方面会显著增大耗时，另一方面xtb自身没法控压，也没法达到理想地模拟溶液状态的效果。 ...

谢谢Sob老师，如果想考察溶液中的聚合反应是不是首选在lammps中用ReaxFF呢

sobereva

naoki 发表于 2020-9-16 09:15
谢谢Sob老师，如果想考察溶液中的聚合反应是不是首选在lammps中用ReaxFF呢

如果模拟尺度很大，是如此
尺度不是很大的话，可以考虑CP2K跑GFN-xTB的动力学，不过我没亲自试过

naoki

sobereva 发表于 2020-9-16 09:24
如果模拟尺度很大，是如此
尺度不是很大的话，可以考虑CP2K跑GFN-xTB的动力学，不过我没亲自试过

非常感谢Sob老师的指点~

wzkchem5

naoki 发表于 2020-9-16 09:15
谢谢Sob老师，如果想考察溶液中的聚合反应是不是首选在lammps中用ReaxFF呢

如果是传统自由基聚合或者某些特别快的阳离子聚合，那可以考虑
如果是可控活性自由基聚合、阴离子聚合、配位聚合，可能1毫秒都未必能聚一个单体，那样就不可能直接模拟了，可以用一条足够长的聚合物链，用metadynamics模拟聚一个单体的情况，然后推而广之到整个聚合过程

naoki

wzkchem5 发表于 2020-9-16 14:12
如果是传统自由基聚合或者某些特别快的阳离子聚合，那可以考虑
如果是可控活性自由基聚合、阴离子聚合、 ...

有道理，谢谢您的指点。我看到之前有文章用ReaxFF做了环氧开环聚合的加速模拟，不知道靠不靠谱