计算化学公社

标题: 求助：跑AIMD时小分子与表面相距越来越远如何解决？ [打印本页]

作者
Author: joeson 时间: 2022-4-24 22:44
标题: 求助：跑AIMD时小分子与表面相距越来越远如何解决？
本帖最后由 joeson 于 2022-4-24 22:59 编辑

各位专家、老师，我为探究一个小分子(SiCl4)在Si表面上的解离吸附的多种路径，用CP2K做了AIMD，但结果是SiCl4离表面越来越远，请问这种情况该如何解决？望各位老师解惑。

作者
Author: sobereva 时间: 2022-4-25 02:13
如果之前没做过优化，先做优化
并且先在低温下跑看看什么情况。并且热浴建议改用CSVR

作者
Author: 丁越 时间: 2022-4-25 08:52

sobereva 发表于 2022-4-25 02:13
如果之前没做过优化，先做优化
并且先在低温下跑看看什么情况。并且热浴建议改用CSVR

我之前也碰到了这个问题，我是在CP2K_INPUT/FORCE_EVEL/EXTERNAL_POTENTIAL中添加球形限制势想来解决这个问题。我查了一下xTB的手册，里面有两种形式的限制势，logfermi和polynomial。前者的好处是势函数的梯度在一定半径内就衰减为零，所以不会造成人为的压缩分子。后者虽然函数形式简单，但是存在弊端就是势能梯度在球半径内不会衰减为零，会造成压缩分子。所以我设置了如下的函数形式，跑一段时间后分子确实被束缚在一定的球半径内。但我不知道如下两种形式是否正确，所以请教下sob老师：

1. &EXTERNAL_POTENTIAL
      ATOMS_LIST 129..138
      FUNCTION A*(sqrt((X-a)^2+(Y-b)^2+(Z-c)^2)/R)^8
      PARAMETERS A R a b c
      VALUES    1.0 5 9.2 10.8 5.6
      UNITS  eV angstrom angstrom angstrom angstrom
  &END EXTERNAL_POTENTIAL

2. &EXTERNAL_POTENTIAL
      ATOMS_LIST 129..138
      FUNCTION A*log10(1+exp(8*sqrt((X-a)^2+(Y-b)^2+(Z-c)^2)-R))
      PARAMETERS A R a b c
      VALUES    1.0E-3 6 9.2 10.8 2
      UNITS  eV angstrom angstrom angstrom angstrom
  &END EXTERNAL_POTENTIAL

(, 下载次数 Times of downloads: 28)

(, 下载次数 Times of downloads: 26)

(, 下载次数 Times of downloads: 22)

作者
Author: joeson 时间: 2022-4-25 09:46

sobereva 发表于 2022-4-25 02:13
如果之前没做过优化，先做优化
并且先在低温下跑看看什么情况。并且热浴建议改用CSVR

感谢Sob老师百忙之中的答疑解惑。

作者
Author: sobereva 时间: 2022-4-25 09:51

丁越发表于 2022-4-25 08:52
我之前也碰到了这个问题，我是在CP2K_INPUT/FORCE_EVEL/EXTERNAL_POTENTIAL中添加球形限制势想来解决这个 ...

我只用过xtb程序的这种限制势，CP2K里我没试过，因此不做过多评论
不过从你的图来看，我感觉限制得有点极端。若半径设得很大，外势搞成大圆盖似的扣在石墨烯上，尽量减小人为影响，我认为更好

作者
Author: 丁越 时间: 2022-4-25 10:08

sobereva 发表于 2022-4-25 09:51
我只用过xtb程序的这种限制势，CP2K里我没试过，因此不做过多评论
不过从你的图来看，我感觉限制得有点 ...

谢谢老师答疑，我知道了。

作者
Author: www1 时间: 2022-8-31 12:29

丁越发表于 2022-4-25 08:52
我之前也碰到了这个问题，我是在CP2K_INPUT/FORCE_EVEL/EXTERNAL_POTENTIAL中添加球形限制势想来解决这个 ...

请问您这个限制势的function里的 A R abc分别代表什么含义呢。这几个参数如何根据体系去设定呢，只是限制水的话大约加大的限制势就可以了呢。
感谢~

作者
Author: 丁越 时间: 2022-8-31 14:50

www1 发表于 2022-8-31 12:29
请问您这个限制势的function里的 A R abc分别代表什么含义呢。这几个参数如何根据体系去设定呢，只是限制 ...

A,R就是下面球形限制势公式中的常量，a，b，c是把球形限制势放置的中心位置

作者
Author: www1 时间: 2022-8-31 16:36

丁越发表于 2022-8-31 14:50
A,R就是下面球形限制势公式中的常量，a，b，c是把球形限制势放置的中心位置

非常感谢您的回复~
那也就是说abc是球形限制势起始的位置，那么这个限制势的范围是哪个参数设定的呢，比如说这个球形限制势的半径。
还有一个问题就是，设置的大小如何去考虑呢，还是需要自己去test一下~

作者
Author: 丁越 时间: 2022-8-31 18:51

www1 发表于 2022-8-31 16:36
非常感谢您的回复~
那也就是说abc是球形限制势起始的位置，那么这个限制势的范围是哪个参数设定的呢， ...

范围即公式中的R.

这个需要自己设置不同大小参数测试一下。

作者
Author: www1 时间: 2022-8-31 20:08

丁越发表于 2022-8-31 18:51
范围即公式中的R.

这个需要自己设置不同大小参数测试一下。

感谢感谢！

作者
Author: ene 时间: 2022-8-31 20:51
直接跑动力学吸附不上，加了外部限制跑动力学得到一个能吸附结果，这有什么意义么……

作者
Author: chands 时间: 2022-8-31 21:19
本帖最后由 chands 于 2022-8-31 21:21 编辑

ene 发表于 2022-8-31 20:51
直接跑动力学吸附不上，加了外部限制跑动力学得到一个能吸附结果，这有什么意义么……

@丁越这个问题我想过。
（1）如果这个外部势是metadynamics的一个collective variable，就好办，但是我不知怎么设（如果表面是刚性的，比如金属的表面，可以将collective variable设成点到面的距离，但是气液界面的情况我不知怎么解决。）
（2）如果不设外部势，设一个wall，好像没什么问题，但是如果算吸附自由能，不管是用PMF还是metadynamics，算平均的时候可能有一部分数据把wall算进去，我想不出什么办法。

作者
Author: 丁越 时间: 2022-8-31 21:53
（1）对于MTD，用PLUMED设置WALL会更方便一点，将你的这个距离CV设置为UPPER_WALL的ARG，比如这篇文章的SI(J. Phys. Chem. C 2021, 125, 15283−15291)可以参考一下。
（2）我不太清楚，你要不咨询一下其他的人吧，我也是正在学习中

作者
Author: chands 时间: 2022-8-31 22:09

丁越发表于 2022-8-31 21:53
（1）对于MTD，用PLUMED设置WALL会更方便一点，将你的这个距离CV设置为UPPER_WALL的ARG，比如这篇文章的SI( ...

好滴我去瞄一下。

作者
Author: scprosper 时间: 2023-3-25 18:41

丁越发表于 2022-4-25 08:52
我之前也碰到了这个问题，我是在CP2K_INPUT/FORCE_EVEL/EXTERNAL_POTENTIAL中添加球形限制势想来解决这个 ...

老师您好，我也想采用这种添加外部势的方式给盒子的上壁添加一个LJ势的排斥项，从而防止盒子里面的氢原子跑出去，输入文件该部分如下
&EXTERNAL_POTENTIAL
ATOMS_LIST 118..225 #H
FUNCTION 4*e*((s/(Z-25.4688))^12) #LJ势的12阶排斥项
PARAMETERS e s
VALUES    0.001908 2.5711
UNITS  eV angstrom
&END EXTERNAL_POTENTIAL
程序能正常开始运行，但是跑到第13步左右就死机不动，看轨迹文件第2个离子步就有一些氢原子飞到天上去了，out文件里面有warning如下：
*** WARNING in external_potential_methods.F:124 :: ASSERTION (cond)  ***
*** failed at line 125                   Error (0.294764E-06) in ***
*** computing numerical derivatives larger then(0.100000E-11) .    ***
请问您知道这个是什么原因导致的吗？如何解决呢？

作者
Author: 丁越 时间: 2023-3-26 21:16

scprosper 发表于 2023-3-25 18:41
老师您好，我也想采用这种添加外部势的方式给盒子的上壁添加一个LJ势的排斥项，从而防止盒子里面的氢原子 ...

如果没有优化结构就跑动力学，那么先优化了再跑。
这种限制势能没啥特殊目的一般不要加，你可以试试降低模拟温度防止H原子跑飞（但我不清楚你模拟的啥）

作者
Author: ustbhh@126.com 时间: 2023-5-11 14:34

丁越发表于 2022-4-25 08:52
我之前也碰到了这个问题，我是在CP2K_INPUT/FORCE_EVEL/EXTERNAL_POTENTIAL中添加球形限制势想来解决这个 ...

请问xTB是什么？我想加一个正方体势限制原子运动，可以实现吗

作者
Author: 丁越 时间: 2023-5-11 18:59

ustbhh@126.com 发表于 2023-5-11 14:34
请问xTB是什么？我想加一个正方体势限制原子运动，可以实现吗

半经验的理论方法。
我不清楚你这样限制是要模拟啥

作者
Author: 朱陈 时间: 2025-12-30 17:00

丁越发表于 2022-4-25 08:52
我之前也碰到了这个问题，我是在CP2K_INPUT/FORCE_EVEL/EXTERNAL_POTENTIAL中添加球形限制势想来解决这个 ...

丁老师，我也遇到了类似的问题，请问您哪一个方法更好呢？

作者
Author: 丁越 时间: 2025-12-30 19:24

朱陈发表于 2025-12-30 17:00
丁老师，我也遇到了类似的问题，请问您哪一个方法更好呢？

把模拟的温度降低，然后升温速率一定要慢一点。

作者
Author: 朱陈 时间: 2025-12-31 09:33

丁越发表于 2025-12-30 19:24
把模拟的温度降低，然后升温速率一定要慢一点。

谢谢您的回复，我是在400 K下跑MTD的，要是不加限制势的话采样空间太大了。我也试过加wall，但是有个CV是配位数，不好加。还试过设置第三个CV限制分子运动，填高斯峰的时候宽度设置0，但是跑了一段时间发现设置第三个CV的两个原子会一起运动。所以现在用您说的这个方法了，我还有个疑问是如果给分子中的一个原子限制，是不是也可以限制分子的运动呀。

作者
Author: 丁越 时间: 2025-12-31 09:45

朱陈发表于 2025-12-31 09:33
谢谢您的回复，我是在400 K下跑MTD的，要是不加限制势的话采样空间太大了。我也试过加wall，但是有个CV是 ...

要是限制采样的空间大小的话就不是添加球形限制势的事儿了，当然在实际模拟中也不要用这种限制势能，会导致模拟失真。推荐用PLUMED进行MTD模拟，可以通过UPPER_WALLS或者LOWER_WALLS关键词限制CVs的范围。

作者
Author: 朱陈 时间: 2025-12-31 10:15

丁越发表于 2025-12-31 09:45
要是限制采样的空间大小的话就不是添加球形限制势的事儿了，当然在实际模拟中也不要用这种限制势能，会导 ...

老师，请问这个和cp2k里wall_plus和wall_minus是一回事吗？还是plumed的更好用

作者
Author: 丁越 时间: 2025-12-31 16:09
本帖最后由丁越于 2025-12-31 16:15 编辑

朱陈发表于 2025-12-31 10:15
老师，请问这个和cp2k里wall_plus和wall_minus是一回事吗？还是plumed的更好用

是。

CP2K定义一个MTD输入文件太繁琐了，而且CVs的支持的类型、后处理分析的模块都远没有PLUMED丰富和好用，所以我一直用的都是PLUMED。

作者
Author: 朱陈 时间: 2026-1-4 20:30

丁越发表于 2025-12-31 16:09
是。

CP2K定义一个MTD输入文件太繁琐了，而且CVs的支持的类型、后处理分析的模块都远没有PLUMED丰富和 ...

好嘞，谢谢您，我也去研究一下。您觉得我之前设置第三个CV限制分子运动，填高斯峰的时候宽度设置0，是不是不合理，我是从cp2k论坛里的学的。

作者
Author: 丁越 时间: 2026-1-4 20:46

朱陈发表于 2026-1-4 20:30
好嘞，谢谢您，我也去研究一下。您觉得我之前设置第三个CV限制分子运动，填高斯峰的时候宽度设置0，是不 ...

没有意义，限制CVs范围就按之前说的直接设置wall就行了。

作者
Author: 朱陈 时间: 2026-1-4 22:27

丁越发表于 2026-1-4 20:46
没有意义，限制CVs范围就按之前说的直接设置wall就行了。

好嘞，谢谢您

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