Cu2O几何优化scf不收敛的问题

西麦西麦

各位老师好，计算新手想请教一下Cu2O的晶胞和几何优化问题，我再我想请教一下，我用cp2K做Cu2O的晶胞优化采用OT（CG）的方法，试了所有的preconditions 结构和scf都不收敛，特别是用full_all的时候还出现了自旋污染的情况，不知道是因为什么原因。计算采用的是6.1的版本，Cu加U方法用的是Mulliken（LOWDIN无法做晶胞优化），同时打开了ENFORCE_OCCUPATION（不开DFT+U energy报负），max_scf从100-400都试了，奇怪的是明明不到100步的时候已经小于设定的10-6次方了，为什么scf还是不收敛？另外我的energy gap从10-4试到了10-2次方都是这种情况，而且采用smear的方法也不行，请问各位老师应该如何调整参数使scf和结构收敛？非常感谢

1.out.log中关于scf收敛的部分
Step     Update method      Time    Convergence         Total energy    Change  
110 OT LS       0.30E-03    6.6                           -14340.6341249525
  111 OT CG       0.30E-03   14.3     0.00000989    -14340.6341470311  2.84E-06
  112 OT LS       0.34E-04    6.7                             -14340.6341441943
   113 OT CG       0.34E-04   14.4     0.00000989    -14340.6341467077  3.23E-07
   114 OT LS       0.39E-05    6.7                           -14340.6341463843
   115 OT CG       0.39E-05   14.7     0.00000989    -14340.6341466708  3.69E-08
   116 OT LS       0.44E-06    6.8                           -14340.6341466339
   117 OT CG       0.44E-06   14.3     0.00000989    -14340.6341466666  4.21E-09
   118 OT LS       0.50E-07    6.7                            -14340.6341466624
   119 OT CG       0.50E-07   14.3     0.00000989    -14340.6341466661  4.78E-10
   120 OT LS       0.57E-08    6.7                             -14340.6341466656
   121 OT CG       0.57E-08   14.4     0.00000989    -14340.6341466661  5.46E-11

2. 输入文件
&GLOBAL
   PROJECT_NAME Cu2O
   RUN_TYPE  CELL_OPT
&END GLOBAL
&MOTION
   &CELL_OPT
     OPTIMIZER  BFGS
     MAX_ITER  200
     &BFGS
       TRUST_RADIUS     2.0000000000000001E-01
       USE_MODEL_HESSIAN  F
     &END BFGS
   &END CELL_OPT
   &CONSTRAINT
     &FIXED_ATOMS
     &END FIXED_ATOMS
   &END CONSTRAINT
&END MOTION
&FORCE_EVAL
   METHOD  QS
   STRESS_TENSOR  ANALYTICAL
   &DFT
     BASIS_SET_FILE_NAME CP2K_basis/BASIS_MOLOPT
     POTENTIAL_FILE_NAME CP2K_basis/GTH_POTENTIALS
     UKS  T
     MULTIPLICITY  1
     PLUS_U_METHOD  MULLIKEN
     &SCF
       MAX_SCF  200
       EPS_SCF     1.0000000000000002E-06
       SCF_GUESS  RESTART
       &OT  T
         MINIMIZER  CG
         PRECONDITIONER  FULL_KINETIC
         ENERGY_GAP     2.0000000000000001E-04
       &END OT
       &OUTER_SCF  F
         OPTIMIZER  BISECT
         EPS_SCF     1.0999999999999999E-07
         MAX_SCF  8
       &END OUTER_SCF
       &SMEAR  F
         METHOD  FERMI_DIRAC
         ELECTRONIC_TEMPERATURE     3.0000000000000000E+02
       &END SMEAR
       &MIXING  T
         METHOD  BROYDEN_MIXING
         ALPHA     1.0000000000000001E-01
         BETA     1.2000000000000000E+00
         NBUFFER  60
       &END MIXING
       &PRINT
         &RESTART  ON
         &END RESTART
       &END PRINT
     &END SCF
     &QS
       EPS_DEFAULT     9.9999999999999555E-09
       METHOD  GPW
     &END QS
     &MGRID
       NGRIDS  4
       CUTOFF     4.0000000000000000E+02
       REL_CUTOFF     4.0000000000000000E+01
     &END MGRID
     &XC
       DENSITY_CUTOFF     1.0000000000000000E-10
       GRADIENT_CUTOFF     1.0000000000000000E-10
       TAU_CUTOFF     1.0000000000000000E-10
       &XC_FUNCTIONAL  NO_SHORTCUT
         &PBE  T
         &END PBE
       &END XC_FUNCTIONAL
       &VDW_POTENTIAL
         POTENTIAL_TYPE  PAIR_POTENTIAL
         &PAIR_POTENTIAL
           TYPE  DFTD3
           PARAMETER_FILE_NAME CP2K_basis/dftd3.dat
           REFERENCE_FUNCTIONAL PBE
           &PRINT_DFTD  OFF
           &END PRINT_DFTD
         &END PAIR_POTENTIAL
       &END VDW_POTENTIAL
     &END XC
     &PRINT
       &PLUS_U  SILENT
       &END PLUS_U
     &END PRINT
   &END DFT
   &SUBSYS
     &CELL
       A     1.6877898331949275E+01    0.0000000000000000E+00    0.0000000000000000E+00
       B    -4.8463810916734079E-03    1.6883587422865769E+01    0.0000000000000000E+00
       C     5.4072294791101226E-03    1.6300353401160559E-03    1.6880735154053269E+01
       MULTIPLE_UNIT_CELL  1 1 1
     &END CELL
     &COORD
     略
     &END COORD
     &KIND Cu
       BASIS_SET DZVP-MOLOPT-SR-GTH
       POTENTIAL GTH-PBE-q11
       &DFT_PLUS_U  T
         L  2
         U_MINUS_J     2.4989541273776900E-01
         &ENFORCE_OCCUPATION  T
           NELEC     4.7830000000000004E+00    4.7830000000000004E+00
           ORBITALS  2 -2 1 -1 0
           EPS_SCF     1.0000000000000001E-05
           MAX_SCF  200
           SMEAR  T
         &END ENFORCE_OCCUPATION
       &END DFT_PLUS_U
     &END KIND
     &KIND O
       BASIS_SET DZVP-MOLOPT-SR-GTH
       POTENTIAL GTH-PBE-q6
     &END KIND
     &TOPOLOGY
       NUMBER_OF_ATOMS  384
       MULTIPLE_UNIT_CELL  1 1 1
     &END TOPOLOGY
   &END SUBSYS
&END FORCE_EVAL

西麦西麦

卡开发发 发表于 2022-5-13 06:56
虽然比较忙，不过我还是做一点点提示吧：

Materials Project是个通过高通量计算的数据库，也就是说必 ...

感谢发发老师百忙之中的指导！我按Eudaimonia和发发老师说的diag的试了一下就正常了，虽然OT、smear、diag方法对于不同体系应该怎么选择还是有些疑惑，但是找到尝试的路径了，另附上在cp2k group上一位老师关于scf不收敛的建议，最后感叹一下，计算想完全搞明白真的太难了。。。

卡开发发

西麦西麦 发表于 2022-5-12 22:23
非常感谢老师的指导，豁然开朗，我认真学习一下，学习心得与总结随后附上，再次感谢！

虽然比较忙，不过我还是做一点点提示吧：

Materials Project是个通过高通量计算的数据库，也就是说必然是按照某个规则产生输入或输出，很多磁性材料人都摸不清楚磁性的情况这些规则也就仅供参考了。之所以算出来0.012μB有可能是因为体系被按照自旋极化的情况做初猜，因为收敛标准的问题没有完全被算到0，当然慎重起见可以看看0和1算完的结果到底哪个能量低（展宽k点请慎重测试，对磁矩有影响）。另外要注意，表面还是有可能出现磁矩与块体有差异（比如块体无磁性但表面的表层有局域化的未成对电子）的情况。

Diag鲁棒性要更好，尤其是偏金属性的体系，恰好过渡金属氧化物有可能会出现bulk是半导体但是切完面gap就不见了的情况。CP2K的DFT+U太拉胯了。

西麦西麦

Eudaimonia 发表于 2022-5-11 20:42
1.这个情况我也不太确定到底什么时候该+U

2.这个要确定以往文献中的结构和你是否完全一致，不一致的话 ...

非常感谢老师的指导，豁然开朗，我认真学习一下，学习心得与总结随后附上，再次感谢！

Eudaimonia

西麦西麦 发表于 2022-5-10 13:57
老师你好，感谢指导，因为我原来是做实验的刚接触计算，有些东西理解的不太到位，我把现在的疑惑总结了一下 ...

1.这个情况我也不太确定到底什么时候该+U

2.这个要确定以往文献中的结构和你是否完全一致，不一致的话可能没有比较是否类似的必要，下面是我在MP里面找到的不同结构的Cu2S和Cu，可以看到他们的DOS并不是完全相同的


3.可以，参考[CP2K] cp2k中的磁矩信息在哪里查看呢，带隙应该是计算PDOS吧，计算DOS的时候不能加K点的

4.最好使用Diag方法进行优化，在我的某个帖子中卡开发发老师也认为Diag方法远比OT更加稳健，目前我的认识是除非大到算不动的体系，尽量避免使用OT
磁矩设置在&kind部分加入
       MAGNETIZATION  X  
X为指定磁矩进行设置
是否可以增加alpha最好进行一些测试

5.论坛之前有一个很好的讨论，关于Fe的，可以找一下

6.如果磁矩，k点等等参数设置恰当，重排也不一定能说不对

7.只要你保证自己的所有模型都在一个条件下进行优化，可以更换参数

其他：今天刚好没啥事跑了一下你的Restart文件，报错
用你的Restart文件通过Multiwfn重新创建了inp文件进行测试，都正常收敛了(都是14-15步收敛)，并未遇上SCF收敛困难的问题

并且计算出的能量不加U和磁矩的结果能量最低，其次是+磁矩的，高了大概0.0015Hatree，最高是+U+磁矩的，说明在该计算体系起码+磁矩是多此一举

具体设置细节是用了2021-Jun-3,  Multiwfn 3.8(dev)版本，更新版本的Multiwfn有其他细节改动，在较老版本的cp2k中容易报错，cp2k版本8.2
直接载入，生成输入文件即可，磁矩和+U部分均在&KIND Cu 下进行对应更改，磁矩设置为1，并同时更改自旋多重度为3，+U复制你的restart文件
单点能计算能量收敛限是5 E-6，几何优化会更高一些，设置为3 E-6
更改基组为你设置的基组
从个人经验来看NBROYDEN适当加大总是有利于收敛，改为12

西麦西麦

Eudaimonia 发表于 2022-5-10 09:42
我的意思是+U如果只是用于凑实验带隙(即只有算PDOS等有效)，对几何优化没有影响(或者更确切一点，对你的 ...

老师，我还想请教一下优化的时候可以更换参数吗？比如我算CEL-OPT的时候用DIIS的方法，在GEO-OPT的时候再换成CG的方法，并且调整其他如cutoff alpha之类的参数？

西麦西麦

老师你好，感谢指导，因为我原来是做实验的刚接触计算，有些东西理解的不太到位，我把现在的疑惑总结了一下
1. 我有查阅了一些文献，+U主要应用于高度局域化的情况矫正，Fe/Co/Ni常见（关于d轨道离域特性的一些讨论参考了一下论坛内容（https://exp.newsmth.net/topic/ar ... 9f81b9bf8cb9b732df3）），但是一些文献也讨论了Cu2O的磁性特征，认为Cu2O具有磁性，同时提出了Cu2O中Cu的d轨道和O的p轨道有强相互作用且p轨道出现电子局域化（https://www.sciencedirect.com/sc ... i/S156717391500200X），我有一些理解上的困难，到底Cu2O什么时候加U？考虑O的p轨道的局域化的时候意思是应该O加U，Cu不加U吗？
2.因为后面想讨论一些total和projected DOS的内容，这个时候我的优化是不是要重点考虑优化的结构和以往文献中PDOS的贴合性？
3.与实验特性数据的参考方面，因为cp2k好像不能算磁矩（这个我还不是太确定，还在研究），另外带隙应该是只能用k点算，那么在+U计算的情况下还有没有其他可以对照的数据来判断计算结果和参数设置是否正确？
4.我学习了Eudaimonia老师之前关于加速scf收敛的帖子（http://bbs.keinsci.com/thread-27630-1-1.html，http://bbs.keinsci.com/thread-27115-1-1.html），有一些收获，想请教一下我这个是否适合用Diag方法优化？另外是否需要参考material的磁矩值设置磁矩，可以请老师指导一下怎么设置吗？还有是否可以考虑增加alpha的值来加速收敛？
5.关于磁矩和磁性的问题我一直弄不太明白，请问老师是否可以推荐一些材料学习
6.下图是我用DIIS优化后得到的（111）表面，感觉表面Cu重排的及其扭曲，这个应该不对吧？


再次感谢老师的指导！

Eudaimonia

西麦西麦 发表于 2022-5-9 23:30
老师你好，如上述回复所说，我看很多文献上讨论Cu2O的时候都采用了DFT+U的方式，也有一些文献没有加U，但 ...

我的意思是+U如果只是用于凑实验带隙(即只有算PDOS等有效)，对几何优化没有影响(或者更确切一点，对你的体系的几何优化没有较大影响)的话，可以考虑不+U使用高k点做几何优化，然后放大晶胞，再+U，不用k点算一次单点能得到PDOS结果

表面重构不一定代表不正确，可能确实会发生表面重构

Eudaimonia

西麦西麦 发表于 2022-5-9 23:22
老师你好，我刚刚开始学这个有些地方还不是特别清楚，我看很多文献上讨论Cu2O的时候都采用了DFT+U的方式 ...

形成能只有在你们的软件、泛函、基组等等参数设定一致时才有可比性，这个差距比较大不说明任何问题

这个磁性非常非常弱，而且+U和+磁矩是两码事

+U大部分应用是凑实验带隙，你可以确定一下带隙有多大差别

西麦西麦

输入文件已付

西麦西麦

Eudaimonia 发表于 2022-5-9 13:24
印象里非磁性体系+U的比较少见，你这个是非磁性体系，建议再考虑一下要不要+U
我之前提到我不确定 ‘DFT ...

老师你好，如上述回复所说，我看很多文献上讨论Cu2O的时候都采用了DFT+U的方式，也有一些文献没有加U，但是material project上Cu2O标注了一定的磁矩（0.012 μB），所以我再计算的时候考虑了+U的方法，所以对于这个地方我的确也不是很清楚，但是现在+U的确收敛不了，不知问题在哪里，我感觉DFT+U之后明显scf很难收敛，有一些帖子也说到了这个问题，我之前用DIIS的方法是优化成功了，但是如老师所说，表面Cu重构的厉害，得到了一个不正确的结构（不好意思老师，这个传图我还没有搞清楚，没有上传图片），查阅了一些资料初步考虑是由于DIIS算法不稳定导致的，所以改回了CG方法，另外老师说的这个“如果没有改进尝试几何优化完后换大晶胞+U做一次单点能再进行分析”我没有太理解，能否请老师再指点一二，另输入文件与坐标信息已付，请老师指导，非常感谢！

西麦西麦

jiangning198511 发表于 2022-5-9 13:06
按道理一般后d区的d轨道库伦效应应该不是很强，所以一般LDA+U 主要计算Fe Co Ni等氧化物，你说的对不上是 ...

老师你好，我刚刚开始学这个有些地方还不是特别清楚，我看很多文献上讨论Cu2O的时候都采用了DFT+U的方式，也有一些文献没有加U，但是我看material project上Cu2O也有一定的磁矩（0.012 μB），所以我再计算的时候考虑了+U的方法，我的做了一下（111）面氧空位的形成能算出来是2.2eV，但文献上看的是1.4-1.7eV，感觉差距有些大，另外对比了PDOS，不过之前那一版算的有问题我又重新算了一遍，另外就是我看Cu2O计算中用GGA或者HYBRID的泛函比较多，目前我用的是GGA的泛函，不知道这个有没有什么影响？感谢老师指导！

Eudaimonia

西麦西麦 发表于 2022-5-9 10:57
老师你好，请问用cp2k算Cu2O体系不用加U吗？但是我之前不加U算的和之前他们用vasp算的包括一些实验结果对 ...

印象里非磁性体系+U的比较少见，你这个是非磁性体系，建议再考虑一下要不要+U
我之前提到我不确定 ‘DFT+U对几何优化结构是否有改进作用’ 可以仔细研究一下，如果没有改进尝试几何优化完后换大晶胞+U做一次单点能再进行分析

Cu可能很容易发生表面重构（依个人经验而言）

另外，输入文件和坐标文件建议用附件传，直接贴的这一大坨妨碍阅读

jiangning198511

西麦西麦 发表于 2022-5-9 10:57
老师你好，请问用cp2k算Cu2O体系不用加U吗？但是我之前不加U算的和之前他们用vasp算的包括一些实验结果对 ...

按道理一般后d区的d轨道库伦效应应该不是很强，所以一般LDA+U 主要计算Fe Co Ni等氧化物，你说的对不上是指那些结果对不上 差距多大？

西麦西麦

jiangning198511 发表于 2022-5-9 10:46
看错输入文件了 你这个体系可以考虑小胞加K点，对Cu体系而已 一般不用+U计算

老师你好，请问用cp2k算Cu2O体系不用加U吗？但是我之前不加U算的和之前他们用vasp算的包括一些实验结果对不上啊？另外不用考虑d轨道的强库伦作用和对磁性的影响吗？