Cp2k跑NVT动力学，常温下收敛，升温就不收敛的问题

吞木木

各位老师好：
      最近有在用Cp2k练习模拟AIMD，根据文献中的说法是先在跑一个298.15k的NVT下的平衡，然后在平衡相下选择一个相对合适的结构去跑3842K的NVT的动力学，以上是整体的模拟思路。刚开始接触Cp2k和分子动力学，假如有很基础的常识性错误，还请老师们原谅，同时很期待各位老师们能够解决我的问题。
      在298.15k下跑的时候（步数15000，步长1fs），没有任何的警告和报错，跑出来的结果也是很不错的，下面有298.15k的.inp文件、.ener文件以及.pdb格式的坐标文件（.out文件太大上传空间不够了）。
      第一次尝试跑3842k的NVT的时候，我选择了第13494步的结构作为初始结构（格式为pdb），这个时候首先进行进行WARNING的部分是：
GENERATE|  Preliminary Number of Bonds generated:                            72
GENERATE| WARNING in connectivity generation!
GENERATE| Two molecules/residues named (MOL1) have different  number of atoms.
GENERATE| Molecule starting at position (19) has Nr. <6> of atoms.
GENERATE| while the other same molecules have Nr. <17> of atoms!
GENERATE| Increasing bondparm_factor by 1.05.. An error was found in the generated
GENERATE| connectivity. Retry...
GENERATE| Present value of BONDPARM_FACTOR (   1.547810 ).

      
       虽然有警告，但是.out文件中还是终归还是出现了下面的信息：


GENERATE|  Preliminary Number of Bonds generated:                           145
GENERATE|  Achieved consistency in connectivity generation.
GENERATE|  Number of Bonds generated:                                       145
GENERATE|  Preliminary Number of Bends generated:                           740
GENERATE|  Number of Bends generated:                                       740
GENERATE|  Number of UB generated:                                          740
GENERATE|  Preliminary Number of Torsions generated:                       3508
GENERATE|  Number of Torsions generated:                                   3508
GENERATE|  Number of Impropers generated:                                     1
GENERATE|  Number of 1-4 interactions generated:                            472

      根据3842k下跑出来的.ener文件，可以看到整个体系的温度和能量在1385步之前都看着挺正常的，但是在1386步体系的温度骤增，控制不住了，升高的很快，这个时候就出现了收敛的报错：

  Leaving inner SCF loop after reaching   200 steps.
  Electronic density on regular grids:       -264.8236188371      -20.8236188371
  Core density on regular grids:              243.9999999897       -0.0000000103
  Total charge density on r-space grids:      -20.8236188474
  Total charge density g-space grids:         -20.8236188474
  Overlap energy of the core charge distribution:               0.00000587383592
  Self energy of the core charge distribution:              -1229.61887969338568
  Core Hamiltonian energy:                                    307.23924582209207
  Hartree energy:                                             520.53050115421161
  Exchange-correlation energy:                               -135.50625408044229
  Dispersion energy:                                           -0.15667569146005
  Total energy:                                              -537.51205661514837
  outer SCF iter =   21 RMS gradient =   0.11E-01 energy =       -537.5120566151
  outer SCF loop FAILED to converge after   21 iterations or 4200 steps
*** WARNING in qs_scf.F:576 :: SCF run NOT converged ***

      接下来，我说一下我解决以上问题的思路，有错误的地方还请老师们指出，首先根据（Warning in connectivity generation (google.com)）谷歌论坛里的提示，先解决第一个报错信息，发现可能是pdb格式文件作为初始结构的问题，然后我就把pdb格式的文件转换成.xyz格式的文件，虽然第一条WARNING不再显示，但是之后依旧会出现不收敛，温度突然之间升高的问题。然后我就觉得这么改只是一种换汤不换药的做法，是不可取的。然后我接下来的思路就是，既然在298.15k的初始构型能够完美的运行NVT，那我就直接用最原始的构型去跑3842k的NVT不就好了？当我直接在298.15k的输入文件中直接改掉温度继续跑的时候（依旧是pdb格式的坐标文件），在跑到2000多步的时候依旧是出现不收敛的情况。
     接下来，我有尝试把3842k降低一点，换成2000k，这个时候整个体系还是在4000步左右就不收敛了，跑的步数比3842k多了1000多步。
  
      还请各位老师能够指出我的问题，万分感激。


2022/03/08更新以下内容：

    首先按照卢天老师的建议，通过关键词（annealing）我对我之前的体系进行了慢慢的升温（其中所有的步长均设为0.5fs，缩小了一半），从10k开始升温到3800K左右，然后通过vmd选取相对合理的结构进行接下来的NVT平衡，依旧会遇到与上边相同的不收敛的问题。


    然后为了排除是体系和初始结构构建的问题，我换了一个体系，并且从CCDC晶体库中直接下载了别人做实验弄好的cif文件（如附件HMX.cif），作为初始构型（但是直接下载的cif文件里边直接通过关键词&topology&coord_file_name会有未识别的文字进行报错，所以我通过Multiwfn转换了一下文件格式为cif，如附件hmx.cif），由于是直接下载的cif文件，所以我并没有进行几何优化或者晶胞优化。直接去跑了298.15k下的nvt（在输入文件中通过关键词multiple对盒子进行了扩胞处理2*1*1），输入文件如附件nvt_hmx_298.inp，结果The number of warnings for this run is : 0，而且通过vmd查看轨迹文件的时候基本没有发现分子产生大幅度的移动现象，然后虽然体系的温度并不能直接判断整个体系到达了平衡状态，但是体系达到平衡状态后，温度一定不会有太大的漂移，因此我把最容易判断体系平衡的温度作为评判标准，298.15k下的温度如附件298.15k.png，可以看到在5ps左右的时候温度基本处于相对稳定状态。
    紧接着，我进行了以下几种情况的处理
    （1）缩短步长为0.5fs，初始构型为原来的cif文件（hmx.cif），直接更改刚才的输入文件nvt_hmx_298.inp中的温度为3500k，继续跑10000步，其跑出来的结果同情况（3），温度升不上去，一直处于1000+k
    （2）缩短步长为0.5fs，初始构型为原来的cif文件（hmx.cif），通过关键词annealing对整个体系进行升温处理，由于对升温的步数没有掌控好，直接最后升温到了几万k，当到达6000k左右的时候已经开始不收敛了，还在收敛的时候大概（1500步到1800步），因此我选取了大概3500k下的温度作为我下一步的初始构型（体系温度如附件升温.png），跑3500k下的nvt（体系温度变化如附件升温_nvt_3500k.png），温度一直上不去（跑出来的.ener文件如附件nvt_hmx-1.ener），这一点和之前的“症状有所不同”。   
      报错信息如下（附件out文件的报错信息_部分.out）
  Electronic density on regular grids:       -447.9999999511        0.0000000489
  Core density on regular grids:              447.9999999821       -0.0000000179
  Total charge density on r-space grids:        0.0000000310
  Total charge density g-space grids:           0.0000000310

  Overlap energy of the core charge distribution:               0.00003518178974
  Self energy of the core charge distribution:              -2382.55372291188041
  Core Hamiltonian energy:                                    693.21446767731504
  Hartree energy:                                             982.58206925632578
  Exchange-correlation energy:                               -230.83928148553358
  Dispersion energy:                                           -0.27085260828278

  Total energy:                                              -937.86728489026609
*** WARNING in qs_scf.F:576 :: SCF run NOT converged ***


    （3）缩短步长为0.5fs，初始构型为跑完298.15k下的平衡之后，选取第9731步的构型为初始构型（如附件9731.pdb），没有进行升温处理直接跑3500k下的nvt，输出的能量文件如下（nvt_hmx-2.ener），做出体系的温度的折线图（如图3500k_298平衡.png），报错情况同（2）在2000步左右就开始偶尔报错不收敛了，之所以并没有选取常温下跑完平衡后的构型（9731步）进行升温处理，然后再去跑3500k下的nvt，是因为我觉得跑出来的结果会和情况（2）类似。

     注：所有的输入文件都是根据如附件nvt_hmx_298.inp进行更改条件得来的。

     还请各位老师，帮帮忙，看一下问题到底出在什么地方，谢谢啦！

      
      

sobereva

以优化过的结构开始，从0K逐渐升温到目标温度，这是最稳妥的方法
一下子把温度拉得很高容易由于某些原子速度过大导致瞬间结构变形厉害、出现不合理的接触等问题导致模拟崩溃或者SCF难收敛
另外，高温时原子运动很快，必要时候尝试减小步长

吞木木

sobereva 发表于 2022-2-14 04:58
以优化过的结构开始，从0K逐渐升温到目标温度，这是最稳妥的方法
一下子把温度拉得很高容易由于某些原子速 ...

谢谢卢天老师的指点，我去试试老师的方法

吞木木

sobereva 发表于 2022-2-14 04:58
以优化过的结构开始，从0K逐渐升温到目标温度，这是最稳妥的方法
一下子把温度拉得很高容易由于某些原子速 ...

老师老师，我更新了一下我的问题还麻烦您有时间给解答一下，然后对于我新选的体系，别人也发表过相关的论文，不过人家是用MS跑的aimd，我用cp2k想要模拟一下，结果就是感觉进行不下去，就在报错，给老师添麻烦了，也给大家添麻烦了