Ce原子在H2O体系下的AIMD模拟，SCF逐渐不收敛

GEEK

各位老师好，想模拟稀土元素的氧化过程，在AIMD计算说SCF收敛特别慢，麻烦各位老师指导一下,波函数读取的结构优化后的文件。有几点疑问：
1. 在模拟的过程设置LSD是否合理；
2. 是否要设置自旋多重度，MULTIPLICITY ，尝试设置1时，同样会在MD过程中计算速度降慢；
麻烦各位老师看看哪些参数不合适


计算过程中变得不收敛：

1. #     Step Nr.          Time[fs]        Kin.[a.u.]          Temp[K]            Pot.[a.u.]        Cons Qty[a.u.]        UsedTime[s]
   
2.          0            0.000000         0.040376896       500.000000000      -141.860553840      -141.819385241         0.000000000
   
3.          1            1.000000         0.035911946       444.709102620      -141.855512124      -141.818797830        14.187150002
   
4.          2            2.000000         0.028191444       349.103654805      -141.847159051      -141.818156322         4.645514011
   
5.          3            3.000000         0.025454664       315.213233417      -141.844539323      -141.818265909         9.491217852
   
6.          4            4.000000         0.027034235       334.773571567      -141.846523827      -141.818663556         5.566816092
   
7.          5            5.000000         0.027475308       340.235518656      -141.847135191      -141.818826299         5.553884983
   
8.          6            6.000000         0.024820317       307.357917669      -141.844205149      -141.818544023        16.323983908
   
9.          7            7.000000         0.024190519       299.558928097      -141.860595659      -141.835557578       822.669039965
   
10.          8            8.000000         0.029083721       360.153013342      -141.914218770      -141.884280166       105.358803988
    
11.          9            9.000000         0.032945994       407.980768582      -141.855310599      -141.821501215      1500.823399067
    

复制代码

输入文件：

1. #Generated by Multiwfn
   
2. &GLOBAL
   
3.   PROJECT Ce-1
   
4.   PRINT_LEVEL LOW
   
5.   RUN_TYPE MD
   
6. &END GLOBAL
   
7. 
   
8. &FORCE_EVAL
   
9.   METHOD Quickstep
   
10.   &SUBSYS
    
11.     &TOPOLOGY
    
12.       COORDINATE XYZ
    
13.       COORD_FILE_NAME Ce-1.xyz
    
14.       &CENTER_COORDINATES
    
15.       &END CENTER_COORDINATES
    
16.     &END TOPOLOGY
    
17. 
    
18.     &CELL
    
19.       ABC    18.000    18.000    18.000
    
20.       PERIODIC NONE #Direction of applied PBC (geometry aspect)
    
21.     &END CELL
    
22. #   &VELOCITY #You can set initial atomic velocities in this section
    
23. #   &END VELOCITY
    
24.     &KIND Ce
    
25.       ELEMENT Ce
    
26.       BASIS_SET DZVP-MOLOPT-SR-GTH-q12
    
27.       POTENTIAL GTH-PBE
    
28.       &DFT_PLUS_U T
    
29.         L 3 #Quantum number of angular momentum the atomic orbitals to +U. 2=d, 3=f
    
30.         U_MINUS_J [eV] 7.0 #Effective on-site Coulomb interaction parameter U(eff) = U - J
    
31.       &END DFT_PLUS_U
    
32.     &END KIND
    
33.     &KIND O
    
34.       ELEMENT O
    
35.       BASIS_SET DZVP-MOLOPT-SR-GTH-q6
    
36.       POTENTIAL GTH-PBE
    
37.     &END KIND
    
38.     &KIND H
    
39.       ELEMENT H
    
40.       BASIS_SET DZVP-MOLOPT-SR-GTH-q1
    
41.       POTENTIAL GTH-PBE
    
42.     &END KIND
    
43.   &END SUBSYS
    
44. 
    
45.   &DFT
    
46.     LSD
    
47.     BASIS_SET_FILE_NAME  BASIS_MOLOPT
    
48.     BASIS_SET_FILE_NAME  BASIS_MOLOPT_UCL
    
49.     POTENTIAL_FILE_NAME  POTENTIAL
    
50. #   WFN_RESTART_FILE_NAME Ce-1-pos-1-RESTART.wfn
    
51.     CHARGE    0 #Net charge
    
52.     #MULTIPLICITY    1 #Spin multiplicity
    
53.     PLUS_U_METHOD MULLIKEN #The method used in DFT+U. Can also be Lowdin
    
54.     &QS
    
55.       EPS_DEFAULT 1.0E-10 #Set all EPS_xxx to values such that the energy will be correct up to this value
    
56.       EXTRAPOLATION ASPC #Extrapolation for wavefunction during e.g. MD. ASPC is default, PS also be used
    
57.       EXTRAPOLATION_ORDER 3 #Order for PS or ASPC extrapolation. 3 is default
    
58.     &END QS
    
59.     &POISSON
    
60.       PERIODIC NONE #Direction(s) of PBC for calculating electrostatics
    
61.       POISSON_SOLVER WAVELET #The way to solve Poisson equation
    
62.     &END POISSON
    
63.     &XC
    
64.       &XC_FUNCTIONAL PBE
    
65.       &END XC_FUNCTIONAL
    
66.       &VDW_POTENTIAL
    
67.         POTENTIAL_TYPE PAIR_POTENTIAL
    
68.         &PAIR_POTENTIAL
    
69.           PARAMETER_FILE_NAME dftd3.dat
    
70.           TYPE DFTD3
    
71.           REFERENCE_FUNCTIONAL PBE
    
72.           #CALCULATE_C9_TERM T #Calculate C9-related three-body term, more accurate for large system
    
73.         &END PAIR_POTENTIAL
    
74.       &END VDW_POTENTIAL
    
75.     &END XC
    
76.     &MGRID
    
77.       CUTOFF 400
    
78.       REL_CUTOFF 40
    
79.     &END MGRID
    
80.     &SCF
    
81.       MAX_SCF 30
    
82.       EPS_SCF 5.0E-04 #Convergence threshold of density matrix of inner SCF
    
83.       SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
    
84.       &DIAGONALIZATION
    
85.         ALGORITHM STANDARD #Algorithm for diagonalization. DAVIDSON is faster for large systems
    
86.       &END DIAGONALIZATION
    
87.       &MIXING #How to mix old and new density matrices
    
88.         METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
    
89.         ALPHA 0.35 #Default. Mixing 40% of new density matrix with the old one
    
90.         NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
    
91.         BETA 0.5
    
92.       &END MIXING
    
93.       &OUTER_SCF
    
94.         MAX_SCF 50 #Maximum number of steps of outer SCF
    
95.         EPS_SCF 5.0E-04 #Convergence threshold of outer SCF
    
96.       &END OUTER_SCF
    
97. 
    
98.       &PRINT
    
99.         &RESTART OFF #Do not generate wfn file to suppress meaningless I/O cost
    
100.         &END RESTART
     
101.       &END PRINT
     
102.     &END SCF
     
103.   &END DFT
     
104. &END FORCE_EVAL
     
105. 
     
106. &MOTION
     
107.   &MD
     
108.     ENSEMBLE NVT
     
109.     STEPS 200 #Number of steps to run
     
110.     TIMESTEP 1.0 #Step size in fs. Decrease it properly for high temperature simulation
     
111.     TEMPERATURE 500  #298.15 Initial and maintained temperature (K)
     
112. #   COMVEL_TOL 0 #Uncomment this can remove translation motion of center-of-mass every step
     
113. #   ANGVEL_TOL 0 #Uncomment this can remove overall rotation every step
     
114.     ANGVEL_ZERO T #Eliminate overall rotation component from initial velocity
     
115.     &THERMOSTAT
     
116.       TYPE NOSE
     
117.     &END THERMOSTAT
     
118.     &PRINT
     
119.       &PROGRAM_RUN_INFO
     
120.         &EACH
     
121.           MD     1 #Output frequency of MD information, 0 means never
     
122.         &END EACH
     
123.       &END PROGRAM_RUN_INFO
     
124.     &END PRINT
     
125.   &END MD
     
126.   &PRINT
     
127.     &TRAJECTORY
     
128.       &EACH
     
129.         MD   1 #Output frequency of coordinates, 0 means never
     
130.       &END EACH
     
131.       FORMAT xyz
     
132.     &END TRAJECTORY
     
133.     &VELOCITIES
     
134.       &EACH
     
135.         MD     0 #Output frequency of velocities, 0 means never
     
136.       &END EACH
     
137.     &END VELOCITIES
     
138.     &FORCES
     
139.       &EACH
     
140.         MD     0 #Output frequency of forces, 0 means never
     
141.       &END EACH
     
142.     &END FORCES
     
143.     &RESTART
     
144.       BACKUP_COPIES 0 #Maximum number of backing up restart file, 0 means never
     
145.       &EACH
     
146.         MD  1 #Frequency of updating last restart file, 0 means never
     
147.       &END EACH
     
148.     &END RESTART
     
149.     &RESTART_HISTORY OFF
     
150.     &END RESTART_HISTORY
     
151.   &END PRINT
     
152. &END MOTION
     

复制代码

结构：

1. 19
   
2. i =        1, E =      -141.8136644967
   
3. Ce         9.2727741648        9.2649868834        8.9737126548
   
4.   O         7.0202624684        7.5082225438        6.7754137217
   
5.   H         7.6467180677        7.3574166154        6.0363402815
   
6.   H         7.6846840817        7.6912282453        7.4992920324
   
7.   O        11.4584594845       11.3528849374        6.3901993215
   
8.   H        11.1475990598       12.2326209529        6.6717893547
   
9.   H        10.6736097071       10.7747053184        6.5243508957
   
10.   O        10.4561912835        8.7389562777       11.8804532920
    
11.   H        11.2217431039        9.3519882259       11.8799051798
    
12.   H        10.8563655193        7.8554270797       12.0106969656
    
13.   O         8.3116752548       11.8426386457       10.8468607084
    
14.   H         8.0306650705       10.9863176099       11.2357496406
    
15.   H         9.2868175770       11.7139853903       10.8064870900
    
16.   O         7.2375696829       10.9059560683        6.6971539253
    
17.   H         6.5606159176       11.1083637890        6.0220268884
    
18.   H         7.3780656170        9.9412794665        6.5540795124
    
19.   O         9.6309334221        6.6157108219        9.4693031951
    
20.   H         9.2874696748        5.7819670166        9.0871507074
    
21.   H        10.0758401949        6.3472620382       10.3013536206
    

复制代码

sobereva

对于这类分子体系，强烈建议改用ORCA做AIMD，比CP2K适合得多
参考
使用ORCA做从头算动力学(AIMD)的简单例子
http://sobereva.com/576（http://bbs.keinsci.com/thread-20800-1-1.html）
自旋多重度必须根据实际情况设定

wudiazhu

sobereva 发表于 2022-11-23 00:36
对于这类分子体系，强烈建议改用ORCA做AIMD，比CP2K适合得多
参考
使用ORCA做从头算动力学(AIMD)的简单例 ...

老师，请问CP2K和ORCA在AIMD上有什么区别？

GEEK

sobereva 发表于 2022-11-23 00:36
对于这类分子体系，强烈建议改用ORCA做AIMD，比CP2K适合得多
参考
使用ORCA做从头算动力学(AIMD)的简单例 ...

谢谢社长，看了您的帖子，单个稀土原子用ORCA确实更合适，还想请教一下，如果把Ce原子换成Ce团簇，CP2K和ORCA哪个更合适呀？

ChemG

那要看你要模拟哪种尺度的，如果团簇非常大可能更接近于体相，适用cp2k表现周期性特征
如果团簇很小，那肯定是更接近于孤立体系，适用orca。

wzkchem5

你确定你要算单个Ce原子和水的反应？不是算体相Ce或Ce团簇/纳米粒子？
我比较好奇单个Ce原子放在500K的水里，实验上是怎么做的？要气相沉积的话水挥发没了，不依赖高真空的实验手段又不好搞出单个Ce原子。

sobereva

GEEK 发表于 2022-11-23 13:52
谢谢社长，看了您的帖子，单个稀土原子用ORCA确实更合适，还想请教一下，如果把Ce原子换成Ce团簇，CP2K和 ...

我说的是Ce配合物，不是单个原子。你之前贴的输入文件里又不是单个Ce原子

CP2K计算需要定义盒子，其中真空区越大用CP2K越吃亏

GEEK

wzkchem5 发表于 2022-11-23 15:11
你确定你要算单个Ce原子和水的反应？不是算体相Ce或Ce团簇/纳米粒子？
我比较好奇单个Ce原子放在500K的水 ...

谢谢老师，是为了计算团簇或者晶体在水蒸汽下的反应，以水的气相密度建模，H2O分子太稀疏，就把水分子加密了，先拿单个Ce原子测试CP2K的参数

GEEK

sobereva 发表于 2022-11-24 01:01
我说的是Ce配合物，不是单个原子。你之前贴的输入文件里又不是单个Ce原子

CP2K计算需要定义盒子，其中 ...

感谢社长答疑解惑，非常感谢。还有一个问题，在Ce晶体表面，500K下建模还是把H2O分子按照液体的密度建模可以吗？