计算化学公社

标题: Ce原子在H2O体系下的AIMD模拟，SCF逐渐不收敛 [打印本页]
作者
Author: GEEK    时间: 2022-11-22 16:51
标题: Ce原子在H2O体系下的AIMD模拟，SCF逐渐不收敛

各位老师好，想模拟稀土元素的氧化过程，在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

  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

  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

  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

  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

  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
复制代码



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

老师，请问CP2K和ORCA在AIMD上有什么区别？
作者
Author: GEEK    时间: 2022-11-23 13:52
sobereva 发表于 2022-11-23 00:36
对于这类分子体系，强烈建议改用ORCA做AIMD，比CP2K适合得多
参考
使用ORCA做从头算动力学(AIMD)的简单例 ...

谢谢社长，看了您的帖子，单个稀土原子用ORCA确实更合适，还想请教一下，如果把Ce原子换成Ce团簇，CP2K和ORCA哪个更合适呀？
作者
Author: ChemG    时间: 2022-11-23 14:40
那要看你要模拟哪种尺度的，如果团簇非常大可能更接近于体相，适用cp2k表现周期性特征
如果团簇很小，那肯定是更接近于孤立体系，适用orca。
作者
Author: wzkchem5    时间: 2022-11-23 15:11
你确定你要算单个Ce原子和水的反应？不是算体相Ce或Ce团簇/纳米粒子？
我比较好奇单个Ce原子放在500K的水里，实验上是怎么做的？要气相沉积的话水挥发没了，不依赖高真空的实验手段又不好搞出单个Ce原子。
作者
Author: sobereva    时间: 2022-11-24 01:01
GEEK 发表于 2022-11-23 13:52
谢谢社长，看了您的帖子，单个稀土原子用ORCA确实更合适，还想请教一下，如果把Ce原子换成Ce团簇，CP2K和 ...

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

CP2K计算需要定义盒子，其中真空区越大用CP2K越吃亏
作者
Author: GEEK    时间: 2022-11-24 08:38
wzkchem5 发表于 2022-11-23 15:11
你确定你要算单个Ce原子和水的反应？不是算体相Ce或Ce团簇/纳米粒子？
我比较好奇单个Ce原子放在500K的水 ...

谢谢老师，是为了计算团簇或者晶体在水蒸汽下的反应，以水的气相密度建模，H2O分子太稀疏，就把水分子加密了，先拿单个Ce原子测试CP2K的参数
作者
Author: GEEK    时间: 2022-11-24 08:41
sobereva 发表于 2022-11-24 01:01
我说的是Ce配合物，不是单个原子。你之前贴的输入文件里又不是单个Ce原子

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

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3