求助团簇AIMD计算无法得到正确的几何结构和温度

王一伊

请教各位老师  在利用CP2K计算带有配体的金属团簇时，采用的是Gaussian优化好的结构进行气相中的AIMD模拟。不知道为什么总是会出现升温严重且配体结构变化严重的现象。
左图是初始构型（优化好的）      右图是动力学一段时间后的结构
         
以下是ener文件中部分内容。。。温度就好离谱
#     Step Nr.          Time[fs]        Kin.[a.u.]          Temp[K]            Pot.[a.u.]        Cons Qty[a.u.]        UsedTime

0            0.000000         0.142981713       350.000000000      -849.113780810      -848.970244905         0.000000000
         1            0.500000         0.267722637       655.349003554      -849.239484787      -848.971195166       114.021744962
         2            1.000000         0.703372666      1721.761674600      -849.680986093      -848.977015866        43.312583999
         3            1.500000         1.359945295      3328.963151583      -850.348293183      -848.987683117        43.300933295
         4            2.000000         2.088812106      5113.131065513      -851.088381559      -848.998784870        43.286104960
         5            2.500000         2.763576127      6764.862623562      -851.764177692      -848.999631990        55.444407645
         6            3.000000         3.378855053      8270.982670749      -852.376400763      -848.996313956        55.388999823
         7            3.500000         3.978498917      9738.830190563      -852.976794987      -848.996708690        64.579657502
         8            4.000000         4.570904217     11188.958679626      -853.568060658      -848.995098492        64.690229601
         9            4.500000         5.156715417     12622.945699668      -854.156695073      -848.997310987        61.625999871
        10            5.000000         5.718730596     13998.683257736      -854.722693571      -849.000516282        67.752503497
        11            5.500000         6.227109198     15243.125694133      -855.235616912      -849.004090217        67.584877040
        12            6.000000         6.646192047     16268.984136615      -855.660899998      -849.009107812        67.356619552
        13            6.500000         6.934939290     16975.798546819      -855.956890683      -849.014949907        70.295927319
        14            7.000000         7.012146971     17164.792564424      -856.045718029      -849.024966098        91.843523374
        15            7.500000         6.801568205     16649.324070384      -855.831739671      -849.019814166        95.022219928
        16            8.000000         6.407656617     15685.081488327      -855.418686133      -848.998843818       100.825603608
        17            8.500000         6.072743524     14865.259288113      -855.076945067      -848.990151308       109.675689472
        18            9.000000         5.951678124     14568.907473348      -854.954927782      -848.987277241       131.281410974
        19            9.500000         6.022009625     14741.069526866      -855.048278584      -849.008259897       143.629240992
        20           10.000000         6.057098992     14826.963576389      -855.097884639      -849.020598073       116.071622038
        21           10.500000         5.942687156     14546.898793284      -854.963660351      -848.998498172       103.871716673

部分输入文件如下
  &DFT
    CHARGE   0  
    MULTIPLICITY     1
    &QS
      EPS_DEFAULT 1E-10 #This is default. Set all EPS_xxx to values such that the energy will be correct up to this value
      EXTRAPOLATION ASPC #Extrapolation for wavefunction during e.g. MD. ASPC is default, PS also be used
      EXTRAPOLATION_ORDER 3 #Order for PS or ASPC extrapolation. 3 is default
    &END QS
    &POISSON
      PERIODIC NONE #Direction(s) of PBC for calculating electrostatics
      PSOLVER ANALYTIC #The way to solve Poisson equation
    &END POISSON
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
      &VDW_POTENTIAL
        POTENTIAL_TYPE PAIR_POTENTIAL
        &PAIR_POTENTIAL
          PARAMETER_FILE_NAME dftd3.dat
          TYPE DFTD3
          REFERENCE_FUNCTIONAL PBE
        &END PAIR_POTENTIAL
      &END VDW_POTENTIAL
    &END XC
    &MGRID
      CUTOFF 450
     NGRIDS 5
      REL_CUTOFF 60
    &END MGRID
    &SCF
      MAX_SCF 128
      EPS_SCF 1.0E-05 #Convergence threshold of density matrix during SCF
#     SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
      &DIAGONALIZATION
        ALGORITHM STANDARD #Algorithm for diagonalization. DAVIDSON is faster for large systems
      &END DIAGONALIZATION
      &MIXING #How to mix old and new density matrices
        METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
        ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
        NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
      &END MIXING
      &SMEAR
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE 350 #Electronic temperature of Fermi-Dirac smearing in K
      &END SMEAR
      ADDED_MOS 30 #Number of virtual MOs to be solved
      &PRINT
        &RESTART #Use "&RESTART OFF" can prevent generating wfn file
          BACKUP_COPIES 0 #Maximum number of backup copies of wfn file
        &END RESTART
      &END PRINT
    &END SCF
&END DFT
&END FORCE_EVAL

&MOTION
  &MD
    ENSEMBLE NVT
    STEPS 10000
    TIMESTEP 0.5 #fs. Decrease it properly for high temperature simulation
    TEMPERATURE 350 #Initial and maintained temperature (K)
    &THERMOSTAT
      TYPE NOSE
    &END THERMOSTAT
  &END MD
  &PRINT
    &TRAJECTORY
      &EACH
        MD     1 #Output frequency of geometry
      &END EACH
      FORMAT xyz
    &END TRAJECTORY
    &VELOCITIES
      &EACH
        MD     1 #Output frequency of velocity
      &END EACH
    &END VELOCITIES
    &RESTART
      BACKUP_COPIES 0 #Maximum number of backing up restart file
      &EACH
        MD 10 #Frequency of updating last restart file
      &END EACH
    &END RESTART
  &END PRINT
&END MOTION

由于未加周期性所以认为不存在相邻团簇之间的影响问题，目前还是不明白哪里出了问题，还请各位老师指点~

Quantum198907

我感觉结构变化严重是温度升高太快且升的太高引起的，是不是先考虑优化一下初始结构，再做动力学计算，或者换个热浴方法，比如GLE......下图是我最近做的一个测试（49个水分子体系），个人看法，仅供参考。

王一伊

Quantum198907 发表于 2021-11-12 11:03
我感觉结构变化严重是温度升高太快且升的太高引起的，是不是先考虑优化一下初始结构，再做动力学计算，或者 ...

谢谢老师回复
目前我测试了NOSE和CSVR两种热浴，我再试试您推荐的，感谢~

Quantum198907

王一伊 发表于 2021-11-12 11:23
谢谢老师回复
目前我测试了NOSE和CSVR两种热浴，我再试试您推荐的，感谢~

我也是在学习中，您有什么新的结论或者想法，欢迎分享，共同进步

喵星大佬

你先修正一下轨迹这都什么玩意

王一伊

喵星大佬 发表于 2021-11-12 12:20
你先修正一下轨迹这都什么玩意

冒昧问一句，我搜索了“修正轨迹”的帖子发现大多是原子或分子跑到设定的盒子之外之后的操作，此处的设置没有周期性，盒子里面的团簇是一个整体，请问这儿修正轨迹的意义是什么呢？谢谢老师~

Quantum198907

王一伊 发表于 2021-11-12 15:02
冒昧问一句，我搜索了“修正轨迹”的帖子发现大多是原子或分子跑到设定的盒子之外之后的操作，此处的设置 ...

我的理解是：VMD里的“Drawing Method”改为“DynamicBonds”，设置合适的截断距离，看上去就好多了

喵星大佬

王一伊 发表于 2021-11-12 15:02
冒昧问一句，我搜索了“修正轨迹”的帖子发现大多是原子或分子跑到设定的盒子之外之后的操作，此处的设置 ...

就是让你改成dynamic bond，不要让键画得乱七八糟的


再就是如果不要周期性你为啥要用cp2k去做AIMD，ORCA或者瓜丝不好么，至于温度，就这么几个原子温度肯定控制不好的，有H的体系先用SHAKE约束一下H，跑一阵子之后再放开。

王一伊

Quantum198907 发表于 2021-11-12 15:33
我的理解是：VMD里的“Drawing Method”改为“DynamicBonds”，设置合适的截断距离，看上去就好多了 ...

嗷嗷这样~

王一伊

喵星大佬 发表于 2021-11-12 15:34
就是让你改成dynamic bond，不要让键画得乱七八糟的

好的，谢谢老师的建议~