About scf error: rho0(ir) < 0

[hitliaomq]

Hello, developers,

I'm trying to use the code to calculate the equation of state for a system. In some volume, it works fine, but in some volume, the scf error (rho0(ir)<0) occurred. How to fix this error?

The following is the last rows:

Charge Density :: Spherical Val Charge in IS   =    10.10566437837118
Charge Density :: Spherical Val Charge in VP   =     7.61140077508768
Charge Density :: Sph. DeepCore Charge in MT-Sphere =     9.99999994575232
Charge Density :: Sph. DeepCore Charge in WS-Sphere =     9.99999994575232
Charge Density :: Sph. DeepCore Charge in EX-Sphere =     9.99999994575232
Charge Density :: Sph. DeepCore Charge in Insc. Sph.=     9.99999994575232
Charge Density :: Sph. DeepCore Charge in VP        =     9.99999994575232
Charge Density :: Sph. SemiCore Charge in MT-Sphere =     7.99922668305437
Charge Density :: Sph. SemiCore Charge in WS-Sphere =     7.99999999540989
Charge Density :: Sph. SemiCore Charge in EX-Sphere =     7.99999999540989
Charge Density :: Sph. SemiCore Charge in Insc. Sph.=     7.99922668198491
Charge Density :: Sph. SemiCore Charge in VP        =     7.99969479201574
For id =     1, ir =   993, r(ir) =   2.38411532
ERROR: rho0(ir) < 0
VALUE:     -0.9781096D-04
STOP AT constructChargeDensity

The i_new file is generated by the ase_must.py, as follows:

================================================================================
                           Input Parameter Data File
================================================================================
                  * * *              * * *              * * *
================================================================================
                           System Related Parameters
================================================================================
No. Atoms in System (> 0)  ::  1
================================================================================
                  * * *              * * *              * * *
================================================================================
Default Potential Input File Name  ::  Cu_ss_pot Zn_ss_pot Sn_ss_pot Al_ss_pot
Default Potential Output File Name  ::  CuZnSnAl_mt_v
Default Potential Input File Form  ::  0
No. Iterations (> 0)  ::  100
Method of SCF Calculation  ::  3
Potential Type (>= 0)  ::  0
Exch-Corr. LDA Type (>= 0)  ::  0
Spin Index Param (>= 1)  ::  1
Kx, Ky, Kz Division (> 0)  ::  12, 12, 12
Mixing algorithm  ::  2
Mixing quantity type  ::  1
Symmetrize BZ Integration  ::  1

The position.dat is as follows:

1.0

0.         3.65246291 3.65246291
3.65246291 0.         3.65246291
3.65246291 3.65246291 0.

CPA  0. 0. 0.  Cu 0.7222537109795818  Zn 0.26683691213263494  Sn 0.010885418756503172  Al 2.3958131279963568e-05

Thank you very much.

[wangy2014]

This error may happen if the CPA calculation is not properly converged. You may check your output file to make sure if the the CPA convergence has indeed achieved at each energy point. Usually increasing the number of k-points and/or decreasing the mixing parameter for CPA iteration will improve the CPA convergence.
Another cause of the problem could be related to the potential being divergent. You may use simple mixing scheme to run 20 SCF iterations to stabilize the potential change from iteration to iteration, and then switch to more aggressive Broyden mixing scheme to speed up the SCF convergence.
The starting potential you use may also be ill conditioned for the unit cell volume you are calculating. You may use the converged potential from the calculation for the other unit cell volume as the starting potential for the volume you are trying to calculate. So long as the lattice constant difference between the two cases is less than 0.2 atomic units, this trick usually works quite well.

[hitliaomq]

@wangy2014 Prof.Wang, Thank you very much for your suggestions. I will try those solutions.