CP2K跑AIMD求助，目的：计算溶液体系Z轴方向的氢键数目分布

1285050231

各位同学、老师：
我想使用CP2K做固-液界面的AIMD模拟，整体思路如下：
1.使用MS做晶胞切面，例如（111）面，使用redefine将其改成正交构型，添加真空层（厚度大概25A），使用supercell进行扩胞，使其XY平面满足25A*25A的最小尺寸要求。
2.利用MS自带的AC模块构建溶液体系（溶液组成为：KOH、正丙醇、水分子）。PS：最开始溶液相只画了K+、正丙醇、水分子。今天请教了sob老师，要添加平衡的OH-根离子，后续会加入，使其溶液相电荷总量为0，即charge=0.
3.使用MS的build layer模块组合固-液两相，保证XY平面尺寸满足25A*25A的最低要求。
3.利用MS直接导出cif文件，放入Multiwfn进行CP2K输入文件的编写。目前是直接使用NPT系综进行计算，发现.out文件跑一半任务就中断了。使用.restart文件续算，也是一样的报错结果。
NPT系综参数如下：
&MD
    ENSEMBLE NPT_F
    STEPS 10000 #Number of steps to run
    TIMESTEP 1.0 #Step size in fs. Decrease it properly for high temperature simulation
    TEMPERATURE 298.15 #Initial and maintained temperature (K)
#   COMVEL_TOL 0 #Uncomment this can remove translation motion of center-of-mass every step
    &THERMOSTAT
      TYPE CSVR
      &CSVR
        TIMECON 100 #Time constant in fs. Smaller/larger results in stronger/weaker temperature coupling
      &END CSVR
    &END THERMOSTAT
    &BAROSTAT
      PRESSURE 1.01325 #Initial and maintained pressure (bar)
      TIMECON 250 #Barostat time constant (fs)
      VIRIAL Z #Relax the cell along which cartesian axes
    &END BAROSTAT
    &PRINT
      &PROGRAM_RUN_INFO
        &EACH
          MD     1 #Output frequency of MD information, 0 means never
        &END EACH
      &END PROGRAM_RUN_INFO
    &END PRINT
  &END MD


疑问如下：
1.整个思路是否存在不合理的地方。例如，买了CP2K的教程，在236页有案例《TiO2与水界面的模拟》，其中提到“由于初始结构偏离平衡状态较为显著，因此先做2 ps的NVT动力学。”针对我所求解的这个体系，是否也应该效仿一下，先跑2ps NVT再进行NPT。
2.对于建构的模型：①XY方向是否一定要满足25A*25A的最小要求，因为画的越大，所需要的分子越多，算的越慢②最后建模完的模型是否需要考虑尽量保证其无限接近立方体（在建模时候可以调整真空层厚度进行微调）。
3.我在进行NVT/NPT前是否需要进行结构优化呢？PS：老师找代算把丙醇换成另一种同分异构体，代算给的最后说明是：先做结构优化，再做NVT。针对代算这一思路我表示怀疑的，因为代算给的代码和他的说明对不上。
4.针对关于NVT/NPT的选择问题，针对计算Z轴方向的氢键数目或Z轴方向的水分子分布。是否有一个比较固定的思路。例如，①NVT→NPT②直接NPT③NPT→NVT。我在类似的固-液界面例子上三种计算思路都见过，不确定哪一种是最合理或者说最稳妥？通用性最好，就是尽可能保证其合理 不被质疑。

sobereva

1285050231 发表于 2024-9-11 09:25
1.“发现.out文件跑一半任务就中断了”，.out文件写到一半就不写了。这是.out文件最后那一部分，突然就中 ...

1 严格按下文编译
CP2K第一性原理程序在CentOS中的简易安装方法
http://sobereva.com/586（http://bbs.keinsci.com/thread-21608-1-1.html）

2 这种说法不实，绝对不要随便相信网上的小视频
DFT计算对范德华作用根本没有截断半径一说。即便是DFT-D3色散校正的截断半径，也是允许超过盒子最短边长一半的。从物理上来说，范德华作用到10埃就已经衰减得颇为充分了。无论从什么角度，都没有理由盒子下限必须到25埃

3 MD没有几何优化那种层面的“收敛”。用词必须准确

4 NPT_F。理应如此

1285050231

sobereva 发表于 2024-9-11 01:36
“发现.out文件跑一半任务就中断了” 指代不明，说清楚怎么断的、报错是什么

不存在“25A*25A的最小要求 ...

1.“发现.out文件跑一半任务就中断了”，.out文件写到一半就不写了。这是.out文件最后那一部分，突然就中断了。
  ----------------------------------- OT ---------------------------------------
  Minimizer      : DIIS                : direct inversion
                                         in the iterative subspace
                                         using   7 DIIS vectors
                                         safer DIIS on
  Preconditioner : FULL_KINETIC        : inversion of T + eS
  Precond_solver : DEFAULT
  stepsize       :    0.15000000                  energy_gap     :    0.20000000
  eps_taylor     :   0.10000E-15                  max_taylor     :             4
  ----------------------------------- OT ---------------------------------------

  Step     Update method      Time    Convergence         Total energy    Change
  ------------------------------------------------------------------------------
     1 OT DIIS     0.15E+00   43.9     0.00002402    -74203.5090124383  4.80E-01
     2 OT SD       0.15E+00   15.0     0.00084341    -74243.4686953126 -4.00E+01
     3 OT SD       0.15E+00   17.0     0.00032434    -74250.4201344706 -6.95E+00
     4 OT SD       0.15E+00   16.8     0.00028082    -74251.7556910641 -1.34E+00
     5 OT SD       0.15E+00   16.8     0.00041374    -74252.4957362250 -7.40E-01
     6 OT DIIS     0.15E+00   15.1     0.00086388    -74251.8863354406  6.09E-01
     7 OT DIIS     0.15E+00   14.8     0.00085655    -74243.1139618725  8.77E+00
     8 OT DIIS     0.15E+00   16.4     0.00084633    -74243.3455312088 -2.32E-01
     9 OT DIIS     0.15E+00   16.8     0.00023922    -74253.0525640467 -9.71E+00
    10 OT SD       0.15E+00   16.7     0.00049087    -74254.1101209751 -1.06E+00
    11 OT DIIS     0.15E+00   16.7     0.00053035    -74255.4415508160 -1.33E+00
    12 OT SD       0.15E+00   16.8     0.00029771    -74252.9173701697  2.52E+00
    13 OT SD       0.15E+00   16.8     0.00031419    -74253.4693215293 -5.52E-01
    14 OT DIIS     0.15E+00   16.8     0.00056504    -74253.6472361066 -1.78E-01
    15 OT SD       0.15E+00   16.8     0.00088073    -74248.5703772239  5.08E+00
    16 OT DIIS     0.15E+00   16.8     0.00049381    -74251.9277156837 -3.36E+00
    17 OT DIIS     0.15E+00   16.8     0.00035888    -74252.7865259087 -8.59E-01
    18 OT SD       0.15E+00   16.7     0.00026900    -74256.6042538035 -3.82E+00
    19 OT SD       0.15E+00   16.7     0.00039414    -74254.0591334464  2.55E+00
    20 OT SD       0.15E+00   16.8     0.00118470    -74240.3496510453  1.37E+01
    21 OT DIIS     0.15E+00   16.7     0.00130744    -74246.9227949049 -6.57E+00

这是服务器的报错信息：
m3co1706:302857:0:302857] Caught signal 11 (Segmentation fault: address not mapped to object at address (nil))
==== backtrace ====
    0  /usr/lib64/libucs.so.0(+0x17970) [0x2b2e481ee970]
    1  /usr/lib64/libucs.so.0(+0x17b22) [0x2b2e481eeb22]
    2  cp2k.popt() [0x36ba332]
    3  cp2k.popt() [0x36afe9a]
    4  /public1/soft/oneAPI/2022.1/compiler/latest/linux/compiler/lib/intel64/libiomp5.so(__kmp_invoke_microtask+0x93) [0x2b2a197b9bb3]
    5  /public1/soft/oneAPI/2022.1/compiler/latest/linux/compiler/lib/intel64/libiomp5.so(__kmp_fork_call+0x3dc) [0x2b2a19735fac]
    6  /public1/soft/oneAPI/2022.1/compiler/latest/linux/compiler/lib/intel64/libiomp5.so(__kmpc_fork_call+0x1a5) [0x2b2a196f7cb5]
    7  cp2k.popt() [0x36b6584]
    8  cp2k.popt() [0x35cc3d2]
    9  cp2k.popt() [0x36de8c5]
   10  cp2k.popt() [0x1bd6957]
   11  cp2k.popt() [0x1bd5ff9]
   12  cp2k.popt() [0x2256abf]
   13  cp2k.popt() [0x1edfc0c]
   14  cp2k.popt() [0x1edec0e]
   15  cp2k.popt() [0x10d6d25]
   16  cp2k.popt() [0x10cc68b]
   17  cp2k.popt() [0xe12239]
   18  cp2k.popt() [0x17ebddb]
   19  cp2k.popt() [0xf67abe]
   20  cp2k.popt() [0x912204]
   21  cp2k.popt() [0x84f0d3]
   22  cp2k.popt() [0x7628d3]
   23  cp2k.popt() [0x7614a9]
   24  cp2k.popt() [0x615c55]
   25  cp2k.popt() [0x6146e4]
   26  cp2k.popt() [0x612a6d]
   27  cp2k.popt() [0x610cb2]
   28  /lib64/libc.so.6(__libc_start_main+0xf5) [0x2b2a19abd555]
   29  cp2k.popt() [0x610bc0]
===================

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   RANK 0 PID 302852 RUNNING AT m3co1706
=   KILLED BY SIGNAL: 9 (Killed)
===================================================================================

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   RANK 1 PID 302853 RUNNING AT m3co1706
=   KILLED BY SIGNAL: 9 (Killed)
===================================================================================

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   RANK 2 PID 302854 RUNNING AT m3co1706
=   KILLED BY SIGNAL: 9 (Killed)
===================================================================================

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   RANK 3 PID 302855 RUNNING AT m3co1706
=   KILLED BY SIGNAL: 9 (Killed)
===================================================================================
…………
2.关于25A*25A的最小要求，这个是从B站学习MS建模构建液相盒子时学到的知识，就是盒子尺寸要满足“大于常用的范德华力截断半径（12.5A）的两倍，即25A”。不确定这句话是否我的体系也需要考虑？
体系是固相-Ni(111)单晶与液相-KOH、丙醇、水构建的溶液盒子模型。不知道sob老师说的体系特征是不是这个？
3.“MD只要跑起来不崩”这句话的含义是：只要收敛不报错就可以不进行结构优化，我理解的对吗？
4.做NPT条件盒子让密度自发变得合理：关于这个具体的参数设计，请教一下sob老师，我应该采用NPT_F还是NPT_I？参考《TiO2与水界面的模拟》案例倾向于使用NPT_F模型，即各项异性。且同时仿照案例参数将VIRIAL设为Z使得模拟过程中只有盒子Z方向运行改变（液相在Z轴上方/固相在Z轴下方），这样是否可行？

sobereva

“发现.out文件跑一半任务就中断了” 指代不明，说清楚怎么断的、报错是什么

不存在“25A*25A的最小要求”这种说法。有的时候边长十几埃都可以，没体系特征没法说

MD只要跑起来不崩，就说明优化不是必须的。通常只要初始结构没有明显不合理的地方（如过近的接触），结构优化可以省

但凡液体部分的初始密度不合理，就必须得做NPT调节盒子让密度自发变得合理。之后统计氢键，是NPT还是NVT无所谓