硅酸盐熔体体系，CP2K收敛慢于VASP，请教如何改进？

wth1219

卢老师、各位同行，大家好。我想用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价电子数的差异应该不会造成收敛速度的差异。）




CP2K-1.inp (10.26 KB, 下载次数 Times of downloads: 14)

2023-8-18 13:38 上传 Uploaded

点击下载
Click to download

CP2K-2.inp (10.29 KB, 下载次数 Times of downloads: 2)

2023-8-18 13:38 上传 Uploaded

点击下载
Click to download

CP2K-3.inp (10.32 KB, 下载次数 Times of downloads: 2)

2023-8-18 13:38 上传 Uploaded

点击下载
Click to download

CP2K-4.inp (10.43 KB, 下载次数 Times of downloads: 6)

2023-8-18 13:38 上传 Uploaded

点击下载
Click to download

INCAR (407 Bytes, 下载次数 Times of downloads: 4)

2023-8-18 13:38 上传 Uploaded

点击下载
Click to download

sobereva

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

wth1219

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

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

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

sobereva

可以先用1E-5跑一阵子等结构弛豫比较充分、SCF比最开始的结构更容易收敛了，再切换到1E-6

wth1219

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 ***

ms860309

高溫直覺應該要降time step
當然這樣整體時間又會更長

乐平

CP2K 的 CUTOFF = 800 Ry 是不是太大了？
如果我没弄错的话 ，1 Ry = 13.6 eV

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

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

sobereva

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

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

sobereva

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

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

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

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

卡开发发

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

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

wth1219

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。

再次感谢您的指教！

6000_Huang-2022_PBE DZVP-UZH-mix_400Ry_test_GAPW.inp (10.45 KB, 下载次数 Times of downloads: 7)

2023-8-31 18:42 上传 Uploaded

点击下载
Click to download

wth1219

报告两个测试进展：

一是在原贴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的计算用时。

759884279

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

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

sobereva

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

可以如此