Encount with problem ------ SCF not NOT converged

[Y Han]

Dear all,

Good morning, I'm a beginer with CP2K, and I'm writting to ask a question and seeking your kind assistance.

I tried to run a GEO_OPT for Fe2O3 bulk with DZVP-MOLOPT-SR-GTH basis set and PBE potential. The DFT+U was applied by the method mulliken.

However, The SCF calculation has been non-convergent (as figure below), and the energy has been oscillating. Attached are my input file, CIF file, and the output structure file.

Could you please point out any issues in these files or settings, thank you very much for your time and kind assistance.

Best regards,

Yun Han

[Y Han]

Here is the copied input setting:

&GLOBAL
  PROJECT Fe2O3Bulk
  PRINT_LEVEL LOW
  RUN_TYPE GEO_OPT
  PREFERRED_DIAG_LIBRARY SL  #
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep
  &SUBSYS
    &CELL
      CELL_FILE_FORMAT CIF
      CELL_FILE_NAME Fe2O3.cif
      PERIODIC XYZ #Direction of applied PBC (geometry aspect)
    &END CELL
    &TOPOLOGY
      COORD_FILE_FORMAT CIF
      COORD_FILE_NAME Fe2O3.cif
    &END TOPOLOGY
    &KIND Fe
      ELEMENT Fe
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q16
      &DFT_PLUS_U
        L 2
        U_MINUS_J 0.1223
      &END DFT_PLUS_U
&END KIND
    &KIND O    
      ELEMENT O
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-PBE-q6
    &END KIND
  &END SUBSYS

  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  POTENTIAL
#   WFN_RESTART_FILE_NAME CuO-111-Vo-RESTART.wfn
    CHARGE    0 #Net charge
#    MULTIPLICITY    1 #Spin multiplicity
    RELAX_MULTIPLICITY 0.001
    UKS TRUE
    PLUS_U_METHOD MULLIKEN
#    &KPOINTS
#      SCHEME GAMMA
#    &END KPOINTS
    &QS
      EPS_DEFAULT 1.0E-10 #Set all EPS_xxx to values such that the energy will be correct up to this value
    &END QS
    &POISSON
      PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
      PSOLVER PERIODIC #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(BJ)
          REFERENCE_FUNCTIONAL PBE
          #CALCULATE_C9_TERM T #Calculate C9-related three-body term, more accurate for large system
        &END PAIR_POTENTIAL
      &END VDW_POTENTIAL
    &END XC
    &MGRID
  COMMENSURATE TRUE
      CUTOFF 400
      REL_CUTOFF 55
    &END MGRID
    &SCF
      MAX_SCF 100
      EPS_SCF 1.0E-05 #Convergence threshold of density matrix of inner 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 4 #Default is 4. Number of previous steps stored for the actual mixing scheme
      &END MIXING
      &SMEAR
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE 300 #Electronic temperature of Fermi-Dirac smearing in K
      &END SMEAR
      ADDED_MOS  50 50 #Number of virtual MOs to solve
      &PRINT
        &RESTART #Note: Use "&RESTART OFF" can prevent generating wfn file
          BACKUP_COPIES 0 #Maximum number of backup copies of wfn file. 0 means never
        &END RESTART
      &END PRINT
    &END SCF
&PRINT  
  &PLUS_U MEDIUM
    FILENAME PLUS_U
  &END PLUS_U
&END PRINT
  &END DFT
&END FORCE_EVAL

&MOTION
  &GEO_OPT
    TYPE MINIMIZATION #Search for minimum
    KEEP_SPACE_GROUP F #If T, then space group will be detected and preserved
    OPTIMIZER CG #Can also be CG (more robust for difficult cases) or LBFGS
     &CG
      MAX_STEEP_STEPS  0
      RESTART_LIMIT 9.0E-01
     &END CG
#    &BFGS
#     TRUST_RADIUS 0.2 #Trust radius (maximum stepsize) in Angstrom
#     RESTART_HESSIAN T #If read initial Hessian, uncomment this line and specify the file in the next line
#     RESTART_FILE_NAME to_be_specified
#    &END BFGS
    MAX_ITER 1000  #Maximum number of geometry optimization
    #The following thresholds of geometry convergence are the default ones
    MAX_DR 1E-3 #Maximum geometry change
    RMS_DR 1E-3 #RMS geometry change
    MAX_FORCE 2E-3 #Maximum force
    RMS_FORCE 2E-3 #RMS force
  &END GEO_OPT
#  &CONSTRAINT
#    &FIXED_ATOMS
#      COMPONENTS_TO_FIX XYZ
#      LIST 1 2 4 5 7 8 9 10 11 13 14 16 17 19 20 22 23 25  26  28  29  31  32  33  34  35  36  37  38  39  40  41  43  44  46  47  49  50  52  53  55  56  58  59  61  62  64  65  67  68  70  71  73  74  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  95  96  98  99 101 102 104 105 106 107 108 110 111 112 113 114 115 116 117 118 119 120 122 123 125 126 128 129 131 132 133 134 135 136 137 138 139 140 141 143 144 146 147 149 150 152 153 155 156 158 159 161 162 164 165 166
#    &END FIXED_ATOMS
#  &END CONSTRAINT
  &PRINT
    &TRAJECTORY
      FORMAT xyz
    &END TRAJECTORY
    &RESTART
      BACKUP_COPIES 0 #Maximum number of backing up restart file, 0 means never
    &END RESTART
  &END PRINT
&END MOTION

[Dragoon]

I think the alpha for mixing (and default value of beta) are generally too large. I don't have any experience with RELAX_MULTIPLICITY keyword but I think you don't need it, as the the spin are not fixed anyway with the smearing. (I think RELAX_MULTIPLICITY is meant for OT and you can search this keyword in this group to find out more information). Also I think it's generally better to use a more restricted convergence criteria, e.g., EPS_DEFAULT=1e-12, EPS_SCF=1e-7, EPS_PGF_ORB=1e-7, even though I used 1e-6 for my own case. Also I would recommend to use OT + FULL_KINETIC + ENERGY_GAP=0.002 with fixed non-ferromagnetic coupling multiplicity, which again seems to work for my own case (transition metal cluster within protein).