计算化学公社

标题: 求助：在显性溶剂下跑MD，如何导入库里面没有的溶剂，这两种情况应该怎么考虑 [打印本页]

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-16 15:23
标题: 求助：在显性溶剂下跑MD，如何导入库里面没有的溶剂，这两种情况应该怎么考虑
本帖最后由哇卡卡卡卡卡于 2022-11-29 17:08 编辑

各位老师好！我现在遇到一个问题就是，在显性溶剂下跑MD，因为这个溶剂库里面没有，所以需要自己构建。我现在想到有两种情况：第一种就是用MCPB建模得到的tleap.in文件将溶剂水改为二氯甲烷（导入二氯甲烷的mol2文件和frcmod文件），最终生成的pdb文件（TSb2_solv.pdb），溶剂很整齐的排布在溶质的周围（图1）。这里设置的是10埃，不知道为什么总共生成的溶剂分子总共有2000多个，而当溶剂是水时，10埃的溶剂盒子水分子617个，所以我感觉这样直接生硬导入是不是不合理。
然后第二种就是自己先把二氯甲烷的溶剂盒子先平衡，然后再导入，但是我在用packmol构建溶剂盒子的时候，我用15埃的盒子，溶剂分子个数为300，我发现，最终的pdb文件都变成了一团（图2），我应该再减小溶剂分子个数吗，不知道为什么同样大小的盒子，分子数却不一样这两种情况。通过packmol构建的溶剂盒子，我应该怎么生成对应的prmtop文件和inpcrd文件呢
我应该优先采取哪一种，构建盒子的话，盒子大小和溶剂个数应该怎么把握，请各位老师指导！

(, 下载次数 Times of downloads: 53) (, 下载次数 Times of downloads: 50)

作者
Author: Frozen-Penguin 时间: 2022-10-16 18:36
第一种情况，npt模拟做的时间够长且体系没有崩溃的话是可以到达正常状态的，但是如果距离正常状态很远的话体系就可能会崩溃，并且amber用GPU模拟时不允许盒子体积大幅变化，所以应该检查一下密度和理论值相差多少，如果相差太多就不太好。第二种情况，我根据密度算出来边长15埃的盒子里面大约应该有32个二氯甲烷分子，300个太多了，不知道是不是算错了。

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-17 08:06
本帖最后由哇卡卡卡卡卡于 2022-10-17 08:53 编辑

Frozen-Penguin 发表于 2022-10-16 18:36
第一种情况，npt模拟做的时间够长且体系没有崩溃的话是可以到达正常状态的，但是如果距离正常状态很远的话 ...

好的，谢谢老师！老师，我还有个疑问，我刚刚重新试了一下，可能是前面弄错了，15埃的盒子，里面装了300个左右的二氯甲烷分子，您说你计算这个密度和理论值，这个应该怎么计算呀

作者
Author: Frozen-Penguin 时间: 2022-10-17 17:32
本帖最后由 Frozen-Penguin 于 2022-10-17 17:34 编辑

哇卡卡卡卡卡发表于 2022-10-17 08:06
好的，谢谢老师！老师，我还有个疑问，我刚刚重新试了一下，可能是前面弄错了，15埃的盒子，里面装了300 ...

边长15埃的盒子,体积是15^3*10^-30 m^3
二氯甲烷的密度是1.325*10^6 g m^-3，分子量85 g mol^-1
阿伏伽德罗常数6.02*10^23
根据这些就可以算出来应该有多少个分子

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-18 11:06

Frozen-Penguin 发表于 2022-10-17 17:32
边长15埃的盒子,体积是15^3*10^-30 m^3
二氯甲烷的密度是1.325*10^6 g m^-3，分子量85 g mol^-1
阿伏伽 ...

谢谢老师，我自己再去算一下

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-19 14:39

Frozen-Penguin 发表于 2022-10-16 18:36
第一种情况，npt模拟做的时间够长且体系没有崩溃的话是可以到达正常状态的，但是如果距离正常状态很远的话 ...

老师您好，我有问题想请教您一下，您说体系崩溃是指什么状态呢，溶剂跑开算吗？

作者
Author: Frozen-Penguin 时间: 2022-10-19 23:26

哇卡卡卡卡卡发表于 2022-10-19 14:39
老师您好，我有问题想请教您一下，您说体系崩溃是指什么状态呢，溶剂跑开算吗？

溶剂看起来奇怪一般是周期性边界条件问题，如果是体系崩溃一般能量会超过上限，在.out文件里面显示为星号

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-24 09:30

Frozen-Penguin 发表于 2022-10-19 23:26
溶剂看起来奇怪一般是周期性边界条件问题，如果是体系崩溃一般能量会超过上限，在.out文件里面显示为星号

老师，您好！不好意思，又来问您了，对第一种情况下的盒子跑平衡时，出现的一种情况就是一直难以达到平衡状态，计算的RMSD值很大（RMSD值一直在不停的增大，现在RMSD值上升到900多），我网上查阅资料说RMSD值的大小与体系有关，表示为最终平衡的结构与最初结构相差的程度。
我目前的想到的解决办法就是：
1、时间拉长，设置到1000ns，让nc轨迹继续跑，不管RMSD的值，老师这里我想问一下，像第一种情况，平衡的时间要达到多少呢，1fs算长嘛，还是要更长
2、不知道是不是因为溶剂没有预先平衡，导致在对第一种情况下跑平衡的时候，很难达到平衡，根据老师您说的，我看out文件，体系并没有崩溃，所以除了办法1，现在也正在对二氯甲烷溶剂盒子进行平衡
NSTEP =330989000 TIME(PS) =  331288.998  TEMP(K) = 296.58  PRESS = -27.2
Etot =    2050.0020  EKtot =    2277.2624  EPtot    =    -227.2605
BOND =    592.4830  ANGLE =    1205.9269  DIHED    =       7.7898
1-4 NB =       7.4487  1-4 EEL =    -22.2567  VDWAALS =    -2165.0178
EELEC  =       51.6736  EHBOND  =       0.0000  RESTRAINT  =       94.6920
EAMBER (non-restraint)  =    -321.9525
EKCMT  =    508.5422  VIRIAL  =    555.4682  VOLUME    =    80004.0216
                                                Density =       1.0385
Ewald error estimate: 0.2390E-03

老师，关于以上出现的问题，以及自己的思考，希望老师方便的时候可以指导我一下！

作者
Author: Frozen-Penguin 时间: 2022-10-24 10:32

哇卡卡卡卡卡发表于 2022-10-24 09:30
老师，您好！不好意思，又来问您了，对第一种情况下的盒子跑平衡时，出现的一种情况就是一直难以达到平衡 ...

计算RMSD的时候不要考虑溶剂，算溶剂的RMSD意义不大，判断溶剂是否平衡应该看体积。溶质分子的RMSD不太可能达到900埃。另外注意你的表述中使用的时间单位有误。

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-24 11:12

Frozen-Penguin 发表于 2022-10-24 10:32
计算RMSD的时候不要考虑溶剂，算溶剂的RMSD意义不大，判断溶剂是否平衡应该看体积。溶质分子的RMSD不太可 ...

谢谢老师的指导！我刚刚去计算了溶质分子的RMSD，0.04左右波动。然后通过cpptraj中autoimage去查看最后一帧的结构，我感觉溶剂分子就填充在一个立方体盒子里面，应该怎么通过体积去看溶剂是否平衡呢，和加热之后的结构进行对比，除了单个分子位置上移动，整体上都是在一个立方体的盒子里面

作者
Author: Frozen-Penguin 时间: 2022-10-24 12:12

哇卡卡卡卡卡发表于 2022-10-24 11:12
谢谢老师的指导！我刚刚去计算了溶质分子的RMSD，0.04左右波动。然后通过cpptraj中autoimage去查看最后一 ...

把.out文件里面的体积或者密度提取出来，如果持续一段时间变化不大（并且和实际密度接近），就可以认为达到平衡了。由于使用周期性边界条件，autoimage处理溶剂分子肯定是在一个立方体盒子里面的，盒子的形状不变，只有大小会变化。

作者
Author: 哇卡卡卡卡卡 时间: 2022-10-25 10:25

Frozen-Penguin 发表于 2022-10-24 12:12
把.out文件里面的体积或者密度提取出来，如果持续一段时间变化不大（并且和实际密度接近），就可以认为达 ...

感谢老师的指导！我查看了out文件，我发现一开始体积变化很明显，但是密度一直都没有特别大的变化，然后我看体积没什么变化了，我提取出nc轨迹的最后一帧，拿去算QM，溶质分子用QM，溶剂分子用立场uff，找到了过渡态，normal termination。
老师，这里我有一个小问题，要是比如一开始在体系没有完全平衡的情况下去算QM，会出现一种什么情况，会正常结束吗

作者
Author: 哇卡卡卡卡卡 时间: 2022-11-28 08:13

Frozen-Penguin 发表于 2022-10-16 18:36
第一种情况，npt模拟做的时间够长且体系没有崩溃的话是可以到达正常状态的，但是如果距离正常状态很远的话 ...

老师您好，前面根据您说的，构建了溶质溶剂盒子，但是溶剂盒子事先没有平衡好，最终跑MD的时候也能成功平衡，但是这个是系统自定义插入多少个溶剂分子（溶剂分子个数过多）。我现在已经平衡好了溶剂盒子，得到了溶剂盒子的top文件和crd文件，在导入预先平衡好的溶剂盒子时，出现了一些错误，我是生成了预平衡好的溶剂盒子的lib文件，导入的时候但是并不能生成最终溶剂+溶质盒子的top文件，文件上传到了附件，老师，不好意思，希望您方便的时候可以帮我看一下，谢谢老师。

作者
Author: Frozen-Penguin 时间: 2022-11-28 18:40

哇卡卡卡卡卡发表于 2022-11-28 08:13
老师您好，前面根据您说的，构建了溶质溶剂盒子，但是溶剂盒子事先没有平衡好，最终跑MD的时候也能成功平 ...

先准备好单个溶质分子的lib文件和单个溶剂分子的lib文件，溶质分子的pdb和溶剂分子的pdb。然后进入tleap，加载力场参数和lib文件，然后用loadpdb分别读取溶质和溶剂的pdb，例如
mol = loadpdb mol.pdb
sol = loadpdb sol.pdb
然后用solvatebox或者solvateoct加溶剂
solvatebox mol sol 10.0
然后加离子就可以保存top和crd了
如果溶质和溶剂的lib都没问题就可以这样操作，因为单独模拟溶剂可以正常进行所以溶剂的lib应该是没问题的

作者
Author: 哇卡卡卡卡卡 时间: 2022-11-29 18:56

Frozen-Penguin 发表于 2022-11-28 18:40
先准备好单个溶质分子的lib文件和单个溶剂分子的lib文件，溶质分子的pdb和溶剂分子的pdb。然后进入tleap ...

老师，是这样的。我按照您说的这个方法生成的溶剂溶质盒子，溶剂是整齐排列的在溶质周围，如图1所示，而且溶剂分子的个数也是系统自定义的进行填充，solvatebox mol sol 10.0只能改变盒子的大小，盒子大小不变时，不能改变溶剂的个数。现在我遇到的问题就是，溶剂分子填充数有点多，导致盒子过于紧密，希望减少溶剂个数。

作者
Author: Frozen-Penguin 时间: 2022-11-30 08:54

哇卡卡卡卡卡发表于 2022-11-29 18:56
老师，是这样的。我按照您说的这个方法生成的溶剂溶质盒子，溶剂是整齐排列的在溶质周围，如图1所示，而 ...

这个sol应该用类似TIP3PBOX的已经平衡的溶剂体系，而不是单个溶剂分子，用平衡过的溶剂体不应该出现过于致密的情况

作者
Author: 哇卡卡卡卡卡 时间: 2022-11-30 09:20

Frozen-Penguin 发表于 2022-11-30 08:54
这个sol应该用类似TIP3PBOX的已经平衡的溶剂体系，而不是单个溶剂分子，用平衡过的溶剂体不应该出现过于 ...

老师，我也去尝试了用平衡好的溶剂盒子作为sol，导入sol.pdb、sol.lib还有底物的pdb和lib文件。最终是可以形成溶剂溶质盒子的TSb2_solv.pdb文件，但是不知道为什么合成后的pdb文件因为溶剂分子的单元出错了无法生成top文件和crd文件。试了好久也不知道是哪里出现问题了，不知道该怎么改
Writing pdb file: TSb2_solv.pdb
printing CRYST1 record to PDB file with box info
Checking Unit.
FATAL: Atom .R<MOL 3>.A<Cl2 5> does not have a type.
......
/public/software/apps/amber20/bin/teLeap: Fatal Error!
Failed to generate parameters

作者
Author: Frozen-Penguin 时间: 2022-11-30 13:57
lib文件不能用相同的残基名，暂时想不到其他可能的问题，单独能模拟说明lib没问题，应该只是有些地方冲突了

作者
Author: 哇卡卡卡卡卡 时间: 2022-12-1 15:28

Frozen-Penguin 发表于 2022-11-30 13:57
lib文件不能用相同的残基名，暂时想不到其他可能的问题，单独能模拟说明lib没问题，应该只是有些地方冲突了

好的老师，非常感谢老师耐心地指导！

作者
Author: monte 时间: 2022-12-1 21:40
我之前用gromacs做MD，溶剂也是二氯甲烷（DCM），力场是Amber14SB+gaff，然后每次溶剂都会裂开，体积变得巨大，溶剂分子之间间隙变得很大，后来换了oplsaa力场就没事了，以我的经验看是力场的原因，虽然我也不知道是为什么。DCM的电荷我都是用gaussian+Multifwn算的，键长键角都是几乎一样，可以说没有差别，但是Amber力场就是算不了

欢迎光临计算化学公社 (http://bbs.keinsci.com/)