SCF Convergence Issues with OT and Diagonalization for CZTS System (Vacuum, Sulfur Layer, and Hydrogen Passivation

[L Heidarizadeh]

Hello CP2K community,

I am running molecular dynamics (MD) simulations on a Cu2ZnSnS4 (CZTS) system using DFT in CP2K. Below is a detailed description of my system and the modifications I applied, followed by the SCF convergence issue I am facing.

Cu2ZnSnS4 (CZTS) system modeled in a periodic box.
Unit cell dimensions: 10.8 × 10.8 × 10.8 Å.
The goal is to study surface interactions and electronic properties with a vacuum layer.

A 20 Å vacuum layer was added in the Z direction to simulate surface effects: 10.8 × 10.8 × 30.8 Å

A layer of sulfur (S) atoms was added to the surface to stabilize the system and account for surface states.

I attempted hydrogen passivation by capping the dangling bonds with H atoms to further stabilize the surface.
I tried running the SCF loop with and without hydrogen passivation, but both cases failed to converge.

SCF Settings and Methods Tried:
Orbital Transformation (OT):
MINIMIZER: DIIS
PRECONDITIONER: FULL_SINGLE_INVERSE
ENERGY_GAP: 0.001
N_HISTORY_VEC: 7
Diagonalization:
I disabled the OT section and enabled diagonalization as a fallback method, but the SCF still did not converge. ( I tried different parameters setting)
SCF Parameters:
SCF_GUESS: ATOMIC
EPS_SCF: 1.0E-6
MAX_SCF: 100

The SCF loop exits after a few minutes, failing to converge under both OT and diagonalization methods.
Are there specific SCF settings or preconditioners that can improve convergence for systems with large vacuum gaps?
Are there alternative strategies for handling surfaces and vacuum layers that could make the system more stable for electronic structure calculations?
Has anyone successfully applied hydrogen passivation to stabilize surfaces and improve SCF convergence in CP2K?

Any suggestions or advice would be greatly appreciated!

Thank you for your help and support.

[Marcella Iannuzzi]

Hi ..

Maybe it simply needs to run for more iterations to converge.

With the information you provide it is hard to guess.

Is the electronic structure calculation of the bulk working fine?

Can you reproduce with your settings (BS, PP, XC etc) the known bulk properties?

Regards

Marcella

[L Heidarizadeh]

Hi,
  Thank you for the suggestion! I’ve already tested the bulk CZTS system with the same basis sets, pseudopotentials, and exchange-correlation functional, and it converged successfully. However, the convergence problem arises only after introducing the 20 Å vacuum layer and surface modifications. I've Adjusted MAX_SCF values (up to 500 iterations) and tried various mixing parameters, including Broyden mixing, but no improvement. I have attached my input file for more clarification (before switching to OT).

Thank you again for your help and suggestions!

Best regards,

Layla

[Marcella Iannuzzi]

Hi ...

It is good that the bulk system converges. Do you also obtain the correct electronic structure? The energy cutoff seems very low. 

Are you using k-points for the bulk?  Is the PBE functional good enough for this type of systems?  12 grids are too many, just use 4 or 5

I suppose that depending on how you cleave the bulk, there might be dangling bonds that should be saturated.

Maybe the surface needs to go through a reconstruction (just guessing), in this case it might help to adjust the coordinates to avoid too many dangling bonds.

Diagonalization is recommended if the energy gap is very small, which can be the case if you have unrelaxed dangling bonds at the surface.

In this case, a smaller mixing ALPHA parameter might help, like 0.005.

[L Heidarizadeh]

Hi Marcella,

Thank you for your feedback and suggestions. I wanted to update you on what I’ve tried so far and ask for some additional guidance.
Since I didn’t include a &KPOINTS block, CP2K defaulted to Γ-point-only sampling. For my 2x2x1 periodic system, I realize this might not be sufficient, especially for accurate TDDFT calculations. should I add a 2 × 2 × 1 k-point grid to sample the Brillouin zone more effectively?

For MD simulations, I’ve been using Γ-point-only sampling since the focus is on atomic forces. Would you recommend defining a small k-point grid (e.g., 2 × 2 × 1) for MD runs as well, or is Γ-point-only sufficient in this case?

I increased the energy cutoff to 500 Ry but still faced convergence issues with the surface model. I also adjusted the mixing parameter ALPHA down to 0.005, but the issue persisted.
I also applied DFT+U  However, I’m still encountering issues with the surface convergence.

I tried adding a sulfur layer and hydrogen passivation to address dangling bonds, but the SCF still failed to converge. You mentioned surface reconstruction—would you suggest running a geometry optimization on the surface before attempting SCF calculations? Also, are there any specific techniques or guidelines for identifying and capping dangling bonds effectively to stabilize the surface?

I really appreciate your insights so far. If you have any further recommendations, I would be grateful.

Best regards,

Layla

[Marcella Iannuzzi]

Hi Layla

It should not be too difficult to verify whether your bulk calculations are accurate enough. 

You can compare some properties with the literature, for instance band gap and density of states. 

The k-point sampling can also be replaced by adding more replicas of the unit cell.

If the 2x2x1 is not sufficient in bulk calculations, it is also not sufficient for slab models,  and this for sure would affect the forces.

But if I correctly understood, you did not manage to calculate forces yet. 

A geometry optimisation before running MD is for sure meaningful in this case, 

but this would not solve the convergence problem, I fear, since the SCF is exactly the same. 

So the problem is still the single point calculation for the set of coordinates of your model. 

Maybe you can share the output of the SCF, such that we get an idea of how bad the problem is. 

Regards

Marcella

[L Heidarizadeh]

Hi Marcella,

Thanks so much again for the feedback . (I believe hit the reply wrongly when I respond to this message ,sorry for that).  I obtained a band gap of about 1 eV, whereas the theoretical and experimental values are around 1.5 eV. My current system is a 2 × 2 × 1 supercell with 64 atoms.  Do you think increasing the size of the supercell would help improve the accuracy, and I should try a larger system? And yes, i didn't calculate the forces.  Thank you again.

Best regards,

Layla