计算化学公社

标题: 硅酸盐熔体体系，CP2K收敛慢于VASP，请教如何改进？ [打印本页]

作者
Author: wth1219 时间: 2023-8-18 13:38
标题: 硅酸盐熔体体系，CP2K收敛慢于VASP，请教如何改进？
卢老师、各位同行，大家好。我想用CP2K做硅酸盐熔体体系的FPMD模拟，其间参考卢老师的文章《解决CP2K中SCF不收敛的方法》对SCF部分的设置做了多次调整，但CP2K的计算均慢于VASP，其原因为CP2K达到SCF收敛所用的圈数显著大于VASP。谨此向卢老师和各位同行请教如何对CP2K的设置做进一步的优化？计算的详情如下，还请大家不吝赐教！

体系：组成为Fe4 Mg30 Na Ca2 Al3 Si24 O89的硅酸盐熔体，其中Fe均为+2价、视作闭壳层，具体坐标见CP2K输入文件的&COORD部分；计算均在Gamma点下开展；
泛函：PBE；
模拟温度：6000 K（NVT系综）；
CPU：Intel(R) Xeon(R) Gold 6248（80核），所有计算均在同一个计算节点上开展；

1. CP2K计算
版本：2022.1；
输入文件如附件CP2K-1.inp、CP2K-2.inp、CP2K-3.inp、CP2K-4.inp所示；
基本设置：DZVP基组，CUTOFF = 800 Ry，EPS_SCF = 1.0E-6 a.u.（EPS_DEFAULT = 1.0E-12 a.u.），SCF部分的设置见下；

1.1. CP2K-1.inp
      100步用时：5810 s；
      SCF部分设置：
         &SCF
            SCF_GUESS RESTART
            EPS_SCF 1.0E-6
            MAX_SCF 50
            &OUTER_SCF
            EPS_SCF 1.0E-6
            MAX_SCF 5
            &END OUTER_SCF
            &OT
            MINIMIZER DIIS
            &END OT
         &END SCF

1.2. CP2K-2.inp
      100步用时：7448 s；
      SCF部分设置：
         &SCF
            SCF_GUESS RESTART
            EPS_SCF 1.0E-6
            MAX_SCF 50
            &OUTER_SCF
            EPS_SCF 1.0E-6
            MAX_SCF 5
            &END OUTER_SCF
            &OT
            MINIMIZER DIIS
            ALGORITHM IRAC
            &END OT
         &END SCF

1.3. CP2K-3.inp
      因前10步SCF均未收敛，故而终止计算；
      SCF部分设置：
         &SCF
            SCF_GUESS RESTART
            EPS_SCF 1.0E-6
            MAX_SCF 50
            &OUTER_SCF
            EPS_SCF 1.0E-6
            MAX_SCF 5
            &END OUTER_SCF
            &OT
            MINIMIZER DIIS
            PRECONDITIONER FULL_ALL
            ENERGY_GAP 0.002
            &END OT
         &END SCF

1.4. CP2K-4.inp
      100步用时：3854 s；
      SCF部分设置：
         &SCF
            SCF_GUESS RESTART
            EPS_SCF 1.0E-6
            MAX_SCF 300
            &DIAGONALIZATION
            ALGORITHM STANDARD
            &END DIAGONALIZATION
            &MIXING
            METHOD BROYDEN_MIXING
            ALPHA 0.1
            NBROYDEN 12
            &END MIXING
            &SMEAR
            METHOD FERMI_DIRAC
            ELECTRONIC_TEMPERATURE 300
            &END SMEAR
            ADDED_MOS 20
         &END SCF

2. VASP计算
版本：6.4.1；
输入文件如附件INCAR所示；
基本设置：ENCUT = 600 eV，EDIFF = 1E-5 eV，ALGO = FAST，ISMEAR = 0，SIGMA = 0.05；
赝势中Na、Mg价电子数与CP2K不一致（CP2K中为9和12、VASP中为7和8），其余一致；
100步用时：2169 s；

使用CP2K和VASP时每步（离子步）FPMD模拟的用时如下图所示，可以看到VASP每步的用时显著低于CP2K，其原因为CP2K达到SCF收敛所用的圈数显著大于VASP（而且从CP2K的收敛情况来看，似乎这个熔体已经与金属接近了？）。不知道这些输入文件还能进行哪些优化，使CP2K的计算速度接近或大于VASP？不胜感激！

（CP2K和VASP的SCF收敛限在一个数量级上（VASP还更严格一些）且均为精度较高的计算（CP2K的CUTOFF = 800 Ry，VASP的ENCUT = 600 eV），赝势中Na、Mg价电子数的差异应该不会造成收敛速度的差异。）

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

(, 下载次数 Times of downloads: 14) (, 下载次数 Times of downloads: 2)
(, 下载次数 Times of downloads: 2) (, 下载次数 Times of downloads: 6)
(, 下载次数 Times of downloads: 4)

作者
Author: sobereva 时间: 2023-8-18 15:12
没具体结构文件，别人不好进行测试
没必要开smearing
EPS_SCF用1E-5一般就够了
另外，为节约时间，可以考虑用GAPW，对Na、Mg用pob-DZVP-rev2全电子基组，其它的还用之前的，这样cutoff用400就够了，可能速度快不少

作者
Author: wth1219 时间: 2023-8-18 15:25

sobereva 发表于 2023-8-18 15:12
没具体结构文件，别人不好进行测试
没必要开smearing
EPS_SCF用1E-5一般就够了

卢老师好，非常感谢您的回复！熔体结构见于附件CP2K-n.inp文件中的&COORD部分。

SCF的收敛限10-6是我们课题组通用的，可能不太好改。您提出的这些方案学生晚上回去都试一下。

作者
Author: sobereva 时间: 2023-8-18 15:40
可以先用1E-5跑一阵子等结构弛豫比较充分、SCF比最开始的结构更容易收敛了，再切换到1E-6

作者
Author: wth1219 时间: 2023-8-18 23:12

sobereva 发表于 2023-8-18 15:40
可以先用1E-5跑一阵子等结构弛豫比较充分、SCF比最开始的结构更容易收敛了，再切换到1E-6

卢老师好，学生测试了您提出的方案，结果不甚理想，谨作说明如下。所有计算均以原贴中的CP2K-4.inp输入文件为基础修改。

1. 删除输入文件中的&SMEAR部分和ADDED_MOS关键词后，30步（离子步）FPMD模拟中只有3~13步共计11步能够收敛，其余步SCF均不收敛，而且SCF不收敛时每步SCF的Convergence几乎没有小于0.01的时候。这个体系可能需要使用&SMEAR保证收敛；

2. 将Na和Mg的基组按照您的建议修改、赝势改为ALL、使用GAPW方法时，无论输入文件中是否包含&SMEAR（及ADDED_MOS），SCF均无法收敛且Convergence高达几十上百，导致FPMD模拟在两步之内崩坏（温度升至数万K）；

您建议的“先用1E-5跑一阵子”学生会继续测试，但如下面粘贴的上述CP2K-4.inp文件最后一步的SCF收敛情况所示，改为10-5之后大约会少5圈左右的SCF、大约5 s，而且因为这个模拟温度太高了，每一步硅酸盐熔体的结构均伴随着化学键的断裂、可能会发生结构的剧烈变化，导致SCF难以收敛。（学生用的这个初始结构是他人文章中在4000 K下用VASP跑成的。）

不知除了这些之外是否还有其它的优化方案？

SCF WAVEFUNCTION OPTIMIZATION

  Step    Update method    Time Convergence       Total energy Change
  ------------------------------------------------------------------------------
   1 Broy./Diag. 0.10E+00 0.6    0.71796802    -4045.6454464221 -4.05E+03
   2 Broy./Diag. 0.10E+00 0.8    1.62566911    -4042.3488118599  3.30E+00
   3 Broy./Diag. 0.10E+00 0.8    1.57627141    -4054.3316716558 -1.20E+01
   4 Broy./Diag. 0.10E+00 0.9    1.05438199    -4050.6302016778  3.70E+00
   5 Broy./Diag. 0.10E+00 0.8    0.81842120    -4042.0889755874  8.54E+00
   6 Broy./Diag. 0.10E+00 0.9    1.05621236    -4045.5201921228 -3.43E+00
   7 Broy./Diag. 0.10E+00 0.9    0.43237197    -4046.2294194020 -7.09E-01
   8 Broy./Diag. 0.10E+00 0.9    0.41543703    -4044.9093263496  1.32E+00
   9 Broy./Diag. 0.10E+00 0.9    0.06233185    -4046.4838416825 -1.57E+00
10 Broy./Diag. 0.10E+00 0.9    0.07575778    -4046.0200066948  4.64E-01
11 Broy./Diag. 0.10E+00 0.9    0.08663157    -4045.4641620875  5.56E-01
12 Broy./Diag. 0.10E+00 1.0    0.05793022    -4045.7191835564 -2.55E-01
13 Broy./Diag. 0.10E+00 0.9    0.03410296    -4045.8390869861 -1.20E-01
14 Broy./Diag. 0.10E+00 0.9    0.03551187    -4045.6349273443  2.04E-01
15 Broy./Diag. 0.10E+00 0.9    0.03015356    -4045.8082945056 -1.73E-01
16 Broy./Diag. 0.10E+00 0.9    0.00644314    -4045.7149369847  9.34E-02
17 Broy./Diag. 0.10E+00 1.1    0.00318764    -4045.6646175427  5.03E-02
18 Broy./Diag. 0.10E+00 0.9    0.00295663    -4045.6702183130 -5.60E-03
19 Broy./Diag. 0.10E+00 0.9    0.00213015    -4045.6439173030  2.63E-02
20 Broy./Diag. 0.10E+00 0.9    0.00145262    -4045.6486850038 -4.77E-03
21 Broy./Diag. 0.10E+00 0.9    0.00119654    -4045.6475310848  1.15E-03
22 Broy./Diag. 0.10E+00 1.0    0.00075674    -4045.6465525036  9.79E-04
23 Broy./Diag. 0.10E+00 0.9    0.00063684    -4045.6495983606 -3.05E-03
24 Broy./Diag. 0.10E+00 0.9    0.00049605    -4045.6466308865  2.97E-03
25 Broy./Diag. 0.10E+00 0.9    0.00026266    -4045.6473517677 -7.21E-04
26 Broy./Diag. 0.10E+00 1.0    0.00029518    -4045.6462291764  1.12E-03
27 Broy./Diag. 0.10E+00 0.9    0.00022718    -4045.6468837599 -6.55E-04
28 Broy./Diag. 0.10E+00 1.0    0.00016916    -4045.6466996428  1.84E-04
29 Broy./Diag. 0.10E+00 0.9    0.00020867    -4045.6473470367 -6.47E-04
30 Broy./Diag. 0.10E+00 0.9    0.00012179    -4045.6462627433  1.08E-03
31 Broy./Diag. 0.10E+00 0.9    0.00001162    -4045.6469816214 -7.19E-04
32 Broy./Diag. 0.10E+00 0.9    0.00002282    -4045.6469636009  1.80E-05
33 Broy./Diag. 0.10E+00 1.0    0.00001041    -4045.6470090546 -4.55E-05
34 Broy./Diag. 0.10E+00 0.9    0.00001564    -4045.6470972373 -8.82E-05
35 Broy./Diag. 0.10E+00 0.9    0.00001830    -4045.6469373483  1.60E-04
36 Broy./Diag. 0.10E+00 0.9    0.00001317    -4045.6469508707 -1.35E-05
37 Broy./Diag. 0.10E+00 0.9    0.00001251    -4045.6470047684 -5.39E-05
38 Broy./Diag. 0.10E+00 0.9    0.00000546    -4045.6469395685  6.52E-05
39 Broy./Diag. 0.10E+00 1.0    0.00000448    -4045.6469560892 -1.65E-05
40 Broy./Diag. 0.10E+00 0.9    0.00000254    -4045.6469525956  3.49E-06
41 Broy./Diag. 0.10E+00 0.9    0.00000207    -4045.6469618810 -9.29E-06
42 Broy./Diag. 0.10E+00 0.9    0.00000252    -4045.6469546552  7.23E-06
43 Broy./Diag. 0.10E+00 0.9    0.00000147    -4045.6469625592 -7.90E-06
44 Broy./Diag. 0.10E+00 1.0    0.00000018    -4045.6469603861  2.17E-06

  *** SCF run converged in 44 steps ***

作者
Author: ms860309 时间: 2023-8-18 23:39
高溫直覺應該要降time step
當然這樣整體時間又會更長

作者
Author: 乐平 时间: 2023-8-19 00:25
CP2K 的 CUTOFF = 800 Ry 是不是太大了？
如果我没弄错的话，1 Ry = 13.6 eV

而你在 VASP 里设置的 ENCUT = 600 eV，远远小于 CP2K 里的 800 Ry。会不会是这原因呢？

另外，VASP 6.4.1 还可以在跑 MD 的时候训练机器学习势场，这样还能跑得更快。后续还能直接调用训练好的势场，进一步加速 MD。

作者
Author: sobereva 时间: 2023-8-19 04:47

乐平发表于 2023-8-19 00:25
CP2K 的 CUTOFF = 800 Ry 是不是太大了？
如果我没弄错的话，1 Ry = 13.6 eV

由于两个程序用的赝势不同所以情况不同
Na用小核GTH赝势的时候即便600 Ry也可能有一点虚假的动力学行为。这是为什么前面我提到可以考虑对Na靠GAPW结合全电子基组来降低对CUTOFF的要求

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

作者
Author: sobereva 时间: 2023-8-19 04:51

wth1219 发表于 2023-8-18 23:12
卢老师好，学生测试了您提出的方案，结果不甚理想，谨作说明如下。所有计算均以原贴中的CP2K-4.inp输入文 ...

没你的GAPW时候的输入文件没法判断正确性

可以步长适当改小，跑个500步看看此时平均每步多少轮SCF能收敛

alpha设那么小反倒可能增加SCF达到收敛的平均圈数

作者
Author: 卡开发发 时间: 2023-8-19 12:33

乐平发表于 2023-8-19 00:25
CP2K 的 CUTOFF = 800 Ry 是不是太大了？
如果我没弄错的话，1 Ry = 13.6 eV

vasp的ENCUT并不对应cp2k的CUTOFF，两者不是一回事，应当对照ENAUG和CUTOFF，虽然最终很可能ENAUG还是会比较小，这个和赝势构造以及计算框架两方面都有关系。另外，GAPW并非只能对全电子基组，即便如此也可以以上述类似的方式降低网格使用需求。

作者
Author: wth1219 时间: 2023-8-31 18:43

sobereva 发表于 2023-8-19 04:51
没你的GAPW时候的输入文件没法判断正确性

可以步长适当改小，跑个500步看看此时平均每步多少轮SCF能收 ...

卢老师好，开&SMEAR（及ADDED_MOS）时GAPW计算的输入文件已附件发送作为示例（学生同时测试了删除&SMEAR及ADDED_MOS，使用XC_DERIV SPLINE2这两种情况。），烦请您审阅、指正；

以原文中收敛情况最好的CP2K-4.inp（100步用时3854 s）为基础，将时间步长改为0.1 fs、跑500步的用时为16300 s（= 5*3260 s），SCF收敛略有改进；

上述CP2K-4.inp中ALPHA设置为0.1是经过测试的，当ALPHA设置为0.2、0.3、0.4时，100步的用时分别为4272 s、4899 s、6401 s。

再次感谢您的指教！

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

作者
Author: wth1219 时间: 2023-9-24 01:36
报告两个测试进展：

一是在原贴VASP输入文件的基础上，将ISMEAR设置为-1（Fermi smearing）、SIGMA设置为k(B)T（0.52 eV），计算用时进一步缩短为了1673 s。

二是如果参考上一点中VASP的设置，将CP2K-4.inp中的ELECTRONIC_TEMPERATURE直接设置为6000 K（虽然不知道这是否合理。），计算用时也会缩短（2786 s）但仍高于VASP的计算用时。

作者
Author: 759884279 时间: 2024-5-13 09:07

sobereva 发表于 2023-8-19 04:47
由于两个程序用的赝势不同所以情况不同
Na用小核GTH赝势的时候即便600 Ry也可能有一点虚假的动力学行为 ...

Sob老师，请问如果体系含有Be-F-Na，是否考虑这三种元素都使用GAPW结合全电子基组来计算呢？

作者
Author: sobereva 时间: 2024-5-13 22:17

759884279 发表于 2024-5-13 09:07
Sob老师，请问如果体系含有Be-F-Na，是否考虑这三种元素都使用GAPW结合全电子基组来计算呢？

可以如此

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