LDA+U+SOC converge issue for scf run

1,167 views
Skip to first unread message

Hong Tang

unread,
Oct 1, 2021, 11:23:41 AM10/1/21
to BerkeleyGW Help
Dear BGW masters,

Probably I should post this on qe forum. But I think maybe some experts here can also help.   
I try to run  LAD+U+SOC for a system. But it seems it is very difficult to converge for the scf run.  The structure of the system  was relaxed from vasp with PBE. I also tried in qe with PBE+U+SOC, and it also can not converge.   I checked the structure several times, It should has no problem. 
I already tried change    mixing_mode from 'plain' to  'local-TF',  reduce  mixing_beta to  0.2
 and increase  mixing_ndim to  10,  but it looks do not help. 

Please hlep.

Best,
Hong


Meng Wu

unread,
Oct 1, 2021, 12:39:21 PM10/1/21
to BerkeleyGW Help, Hong Tang
Hi Hong,

I have some suggestions.
1. From my experience, LDA + U should be numerically more stable than PBE + U.
2. Increase energy cutoff.
3. Try 'cg' instead of 'david' for the 'diagonalization'.
4. If this is a metal, use more kpoints.
5. Check your lattice constants and internal coordinates, make sure that the system has the intended symmetry. Symmetrize the crystal structure manually if necessary (e.g., 0.666662222 -> 0.666666667 ).
6. Add constraints for magnetization (see 'constrained_magnetization') if you know the magnetic ground state (from your intuition or experimental data).
7. Check if there are any ghost states in your pseudopotential. Use another pseudopotential.
8. Start with LDA + U + noncollinear + noSOC, and then continue the scf calculation with LDA + U + noncollinear + SOC.

Let me know if you still have problems after trying all the above tricks.

Best!
Meng

Hong Tang

unread,
Oct 1, 2021, 5:00:37 PM10/1/21
to BerkeleyGW Help, wu1m...@berkeley.edu, Hong Tang
Hi Meng,

Thank you so much for helping!  I will try with your suggestions. 
For the pseudopotential,  I used the ones from  www.pseudo-dojo.org, which are PBE and FR. I just set the input_dft='pz'  in qe.  I suppose this is LDA.  Is any other places I can find LDA pseudopotentials? 

Best,
Hong

Meng Wu

unread,
Oct 1, 2021, 5:07:52 PM10/1/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

You are very welcome.

Setting input_dft = 'pz' for PBE pseudopotentials seems strange. I suggest you find true LDA pseudopotentials from this page: https://www.quantum-espresso.org/pseudopotentials. Besides, I am also maintaining a small collection of ONCVPSP (LDA, PBE, PBEsol) pseudopotentials in my Github repo: https://github.com/wu2meng3/ONCVPseudoPack, just for your reference.

Best!
Meng

--
Meng Wu
----------------
Physics Department
UC Berkeley

Hong Tang

unread,
Oct 1, 2021, 6:04:37 PM10/1/21
to BerkeleyGW Help, wu1m...@berkeley.edu, BerkeleyGW Help, Hong Tang

Dear Meng,

Thank you so much for the information. 
Actually, I am doing CrI3 system (ribbon). I know you did this 1L system. Right now, I changed the pps, I hope it will go well. 

Best,
Hong

Hong Tang

unread,
Oct 4, 2021, 5:02:28 PM10/4/21
to BerkeleyGW Help, Hong Tang, wu1m...@berkeley.edu, BerkeleyGW Help
Dear Meng,

I did the 01-scf calculation as you suggested. I did first with LDA+U and converged to conv_thr =1D-10,  then I do LDA+U+SOC, however, I can only converge to  conv_thr =1D-5. 

Then I do 02-wfn, but I get error,  please see below. Please help me.


Hong


___________
error:

.....
      k(   16) = (   0.5000000   0.7142857   4.1687274), wk =   0.0625000

     Dense  grid:  2344729 G-vectors     FFT dimensions: ( 320, 216,  72)

     Estimated max dynamical RAM per process >     909.22 MB

     Estimated total dynamical RAM >     142.07 GB
     Generating pointlists ...
     new r_m :   0.0371 (alat units)  2.1028 (a.u.) for type    1
     new r_m :   0.0371 (alat units)  2.1028 (a.u.) for type    2

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine read_scf (1):
     Reading ldaU ns

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine read_scf (1):
     Reading ldaU ns
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     stopping ...

Meng Wu

unread,
Oct 4, 2021, 6:25:57 PM10/4/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

I have tried it on monolayer CrI3. Just want to double-check with you.

If the scf run, did you set the following?
nspin = 4,
noncolin=.TRUE.,
lspinorb=.FALSE.,

lda_plus_u=.TRUE.,
lda_plus_u_kind=1,
Hubbard_U(1)=1.5D0,
Hubbard_J(1,1)=0.5D0,
Hubbard_J(2,1)=0.0D0,

And in the continued scf run, did you set the following?
nspin = 4,
noncolin=.TRUE.,
lspinorb=.TRUE.,

lda_plus_u=.TRUE.,
lda_plus_u_kind=1,
Hubbard_U(1)=1.5D0,
Hubbard_J(1,1)=0.5D0,
Hubbard_J(2,1)=0.0D0,

startingwfc='file',

Could you also post your input files?

Best!
Meng

Hong Tang

unread,
Oct 4, 2021, 6:37:44 PM10/4/21
to BerkeleyGW Help, wu1m...@berkeley.edu, BerkeleyGW Help, Hong Tang
Dear Meng, 
please see the following,  


scf in 
&control
   prefix = 'CrI3'
   calculation = 'scf'
   restart_mode = 'restart'    ! 'from_scratch'
   wf_collect = .true.      ! .false.
   ! tstress = .true.
   ! tprnfor = .true.
   outdir = './'
   wfcdir = './'
   pseudo_dir = './' ,  max_seconds=12600
/
&system
   ibrav = 0
   nat = 18
   ntyp = 2
   nbnd = 400   ! 77*2  occup.
   ecutwfc = 70.0 ,   ! input_dft='pz', 
   nspin = 4
   noncolin=.TRUE.
   starting_magnetization(1)=1.0 , starting_magnetization(2)=1.0 , constrained_magnetization='total direction', fixed_magnetization(3)=90 ,
   lspinorb=.TRUE.
   angle1(1)=90.0  ,  angle1(2)=90.0 ,
   angle2(1)=90.0  ,  angle2(2)=90.0 ,
   lda_plus_u=.true.
   lda_plus_u_kind = 1
   Hubbard_U(1)=1.5 ,  Hubbard_J(1,1)=0.5 , 
/
&electrons
   electron_maxstep = 300
   conv_thr = 1.0d-5
   mixing_mode = 'local-TF'
   mixing_beta = 0.2
   mixing_ndim = 10
   diagonalization =   'david'    ! 'cg',    !
   diago_david_ndim = 4
   diago_full_acc = .true.
/
CELL_PARAMETERS angstrom
    30.0000000000000000    0.0000000000000000    0.0000000000000000
     0.0000000000000000   21.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    6.9715520833838687
ATOMIC_SPECIES
  Cr   24.000   Cr.upf
  I    53.000    I.upf  
! H    1.000     H.upf
ATOMIC_POSITIONS angstrom
Cr        7.973519498588501619e+00 1.056907996607464639e+01 1.985114235172969810e+00
Cr        1.202810129810380424e+01 1.056955010504596260e+01 1.985202558709478149e+00
Cr        5.959334221702262901e+00 1.056632620656042576e+01 5.470852466435861317e+00
Cr        1.403974277913469848e+01 1.056632744959639325e+01 5.471753972011825340e+00
I        1.000107516870635038e+01 1.2121
......
......



the bands in for 02-wfn
&control
   prefix = 'CrI3'
   calculation = 'bands'
   restart_mode = 'from_scratch'
   wf_collect = .true.
   ! tstress = .false.
   ! tprnfor = .false.
   outdir = './'
   wfcdir = './'
   pseudo_dir = './'
/
&system
   ibrav = 0
   nat = 18
   ntyp = 2
   nbnd = 1500   ! 77*2  occup.
   ecutwfc = 70.0 ,   ! input_dft='pz', 
   nspin = 4
   noncolin=.TRUE.
   starting_magnetization(1)=1.0 , starting_magnetization(2)=1.0 , constrained_magnetization='total direction', fixed_magnetization(3)=90 ,
   lspinorb=.TRUE.
   angle1(1)=90.0  ,  angle1(2)=90.0 ,
   angle2(1)=90.0  ,  angle2(2)=90.0 ,
   lda_plus_u=.true.
   lda_plus_u_kind = 1
   Hubbard_U(1)=1.5 ,  Hubbard_J(1,1)=0.5 , 
/
&electrons
   electron_maxstep = 300
   conv_thr = 1.0d-5
   mixing_mode = 'local-TF'
   mixing_beta = 0.2
   mixing_ndim = 10
   diagonalization = 'david'
   diago_david_ndim = 4
   diago_full_acc = .true.
   startingwfc = 'random'
/
CELL_PARAMETERS angstrom
    30.0000000000000000    0.0000000000000000    0.0000000000000000
     0.0000000000000000   21.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000    6.9715520833838687
ATOMIC_SPECIES
  Cr   24.000   Cr.upf
  I    53.000    I.upf  
! H    1.000     H.upf
ATOMIC_POSITIONS angstrom
Cr        7.973519498588501619e+00 1.056907996607464639e+01 1.985114235172969810e+00
Cr        1.202810129810380424e+01 1.056955010504596260e+01 1.985202558709478149e+00
Cr        5.959334221702262901e+00 1.056632620656042576e+01 5.470852466435861317e+00
Cr        1.403974277913469848e+01 1.056632744959639325e+01 5.471753972011825340e+00
I        1.000107516870635038e+01 1.212148527025875921e+01 2.969539086694043650e+00
I        1.38840762403175705
.....
.....






Hong Tang

unread,
Oct 4, 2021, 6:44:11 PM10/4/21
to BerkeleyGW Help, Hong Tang, wu1m...@berkeley.edu, BerkeleyGW Help
Dear Meng,

I just checked your previous email again. 
In the 1st scf run, I did set the same as you posted.  However, in the continued scf run, I did not set " startingwfc='file',  "   So this maybe the problem, right? 

Best,
Hong

Meng Wu

unread,
Oct 4, 2021, 6:52:44 PM10/4/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

I thought your problem was the convergence of DFT charge and spin density, so I suggested you use two steps of scf run (LDA+U+noncollinear+noSOC, LDA+U+noncollinear+SOC) to converge the charge and spin density. This is because sometimes +U and SOC will lead to bad convergence behavior. But in your input files, it seems you are doing one scf run and a following nscf run, which makes me confused. In your original email, you said 

> But it seems it is very difficult to converge for the scf run.

In this way, you should first converge the scf run (using the two-step approach or tweak other parameters), and after that consider nscf run to generate WFN.

Best!
Meng

Hong Tang

unread,
Oct 4, 2021, 6:54:13 PM10/4/21
to BerkeleyGW Help, Hong Tang, wu1m...@berkeley.edu, BerkeleyGW Help
Dear Meng,

I just check the QE manual, it says, if   restart_mode= 'restart',   The default for
startingwfc and startingpot is set to 'file'. 

so, if in the continued scf run, although I did not set startingwfc='file',   it is still  startingwfc='file'  by default, right? 


Best,
Hong

Hong Tang

unread,
Oct 4, 2021, 7:00:51 PM10/4/21
to BerkeleyGW Help, Hong Tang, wu1m...@berkeley.edu, BerkeleyGW Help
Dear Meng, 

I did do 2 runs for the scf. I did not post you the in file of the 1st scf. What I posted to you is the continued scf and the nscf.   

in the 1st scf, I set lspinorb=.false.     otherwise, I can not converge to conv_thr =1D-10. 
in the continued scf, I set   lspinorb=.true.,  however, I can only   converge to conv_thr =1D-5.    not conv_thr =1D-10.

Best,
Hong

Meng Wu

unread,
Oct 4, 2021, 7:53:49 PM10/4/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

I guess so, but you need to double-check.

Best!
Meng

Meng Wu

unread,
Oct 4, 2021, 7:57:20 PM10/4/21
to Hong Tang, BerkeleyGW Help
I see. Did you use PBE or LDA? LDA might have better convergence.

Meng

Hong Tang

unread,
Oct 5, 2021, 10:15:22 AM10/5/21
to BerkeleyGW Help, wu1m...@berkeley.edu, BerkeleyGW Help, Hong Tang
Hi Meng,
Yes, I used LDA. I used the pseudopotentials from your page. 
 I noticed that I did not set Hubbard_J(2,1)=0.0D0,   and during the running, the  Hubbard_J(2,1) was  recalculated as a non zero value.  Does this matter? 

Best,
Hong

Meng Wu

unread,
Oct 5, 2021, 12:30:41 PM10/5/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

From https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm464, Hubbard_J(2,1) is by default non-zero, but I am not sure if it will affect convergence in your case. 2D CrI3 is pretty easy to converge in my case.

Best!
Meng

Hong Tang

unread,
Oct 5, 2021, 12:37:00 PM10/5/21
to BerkeleyGW Help, wu1m...@berkeley.edu, BerkeleyGW Help, Hong Tang
I just want to report you about this. It is very weird. I am running the scf again. I am doing the 1st scf with lspinorb=.false. , I set    Hubbard_U(1)=1.5D0 ,  Hubbard_J(1,1)=0.5D0 , Hubbard_J(2,1)=0.0D0, 
however, when it runs , 
in the out file, it shows

     atomic species   valence    mass     pseudopotential
        Cr            14.00    24.00000     Cr( 1.00)
        I             17.00    53.00000     I ( 1.00)


     Full LDA+U calculation (l_max = 2) with parameters (eV):
     U(  1) =   1.5000   J(  1) =   0.5000   B(  1) =   0.0574


It becomes B(  1) =   0.0574 , the same value I got in my runs several days before. 
Do you let me know why is that? 

Best,
Hong



Meng Wu

unread,
Oct 6, 2021, 3:20:53 AM10/6/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

If you read the description of Hubbard_J on https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm464 more carefully, you will find the following sentence,

If B or E2 or E3 are not specified or set to 0 they will be
calculated from J using atomic ratios.


I think your set of Hubbard parameters are fine.

Best!
Meng

Hong Tang

unread,
Oct 6, 2021, 10:27:03 AM10/6/21
to BerkeleyGW Help, wu1m...@berkeley.edu, BerkeleyGW Help, Hong Tang
Thank you so much for pointing out this for me.
I am trying  with  diagonalization =   'cg'  , however, it will crash after several iterations,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine  cdiaghg (321):
      problems computing cholesky
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

It looks the only way I can do is  diagonalization =   'david'  .  

Best,
Hong


Meng Wu

unread,
Oct 6, 2021, 2:13:43 PM10/6/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

I cannot guarantee it will work, but you can try diago_full_acc = .false. along with diagonalization = 'cg'. 

If it still has problems of computing Cholesky, you might want to try tweaking the number of bands (e.g., 400 -> 399 or 401), the number of diagonalization groups (e.g., pw.x -nd XXX), or change the number of processes (e.g., pw.x -np XXX) used. Sometimes, this error will go away magically.

This is truly a problem of numerical stability and I don't have a perfect solution to that. If it still won't work, you should post this problem to the QE mail-list.

Meng

-- 
You received this message because you are subscribed to a topic in the Google Groups "BerkeleyGW Help" group.
To unsubscribe from this topic, visit https://groups.google.com/a/berkeleygw.org/d/topic/help/jWCPNrkY1X4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to help+uns...@berkeleygw.org.
To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/0e961713-9838-480e-adce-cdcf7fdddbd3n%40berkeleygw.org.

Hong Tang

unread,
Oct 6, 2021, 2:32:15 PM10/6/21
to BerkeleyGW Help, wu1m...@berkeley.edu, BerkeleyGW Help, Hong Tang
Thank you so much! Meng.
I will try.

Best,
Hong

Hong Tang

unread,
Oct 8, 2021, 5:14:45 PM10/8/21
to BerkeleyGW Help, Hong Tang, wu1m...@berkeley.edu, BerkeleyGW Help
Hi Meng,

I have another question related to spinor BSE.
1. If I do all the folders from 01-scf to 06-wfnq_fi with LAD+U+noncollinear+noSOC, then proceed the spinor BSE with the results.   
2. If I do all the folders from 01-scf to 06-wfnq_fi with LAD+U+noncollinear+SOC, then proceed the spinor BSE with the results.

In BerkelyGE3.0.1, the workflow for the above 1 and 2 are the same, right? They are all the so-called full spinor BSE, right? 

Best,
Hong

Meng Wu

unread,
Oct 8, 2021, 9:35:45 PM10/8/21
to Hong Tang, BerkeleyGW Help
Hi Hong,

Yes, you are right. As long as you use noncol = .true. in QE, you are dealing with spinor wavefunctions. The full-spinor feature in BGW does not have any special treatment for lspinorb = .true. or .false..

Meng

Dibya Prakash Rai

unread,
Apr 16, 2024, 6:18:46 AM4/16/24
to BerkeleyGW Help, Meng Wu, BerkeleyGW Help, Hong Tang
Hi All,

I followed this discussion. As I faced the same problem. I am using QE-7.2
I tried all the possible suggestions mentioned in this discussion but with no success.
Can anyone help me?
*the result is SCF not converged 

My input looks like this.

&CONTROL
  calculation = 'scf'
  etot_conv_thr =   6.0000000000d-05
  forc_conv_thr =   1.0000000000d-04
  outdir = 'temp'
  prefix = 'bulk_soc'
  pseudo_dir = '../../pseudo/'
  tstress = .true.
  tprnfor = .true.
  verbosity = 'high'
/
&SYSTEM
  degauss =   1.4699723600d-02
  ecutrho =   750.0
  ecutwfc =   50.0
  ibrav = 0
    nat= 6,
    ntyp= 3,
    nbnd = 136,
    starting_magnetization(2)=1.0,
    starting_magnetization(3)=-1.0,
    occupations='smearing',
    degauss=0.005,
    smearing='mv',
    nspin = 4,
    noncolin = .true.
    lspinorb = .true.
/
&ELECTRONS
   diagonalization='david'
   conv_thr = 1.0e-8
   mixing_beta = 0.7
/

ATOMIC_SPECIES
O      15.9994   O.rel-pbe-n-kjpaw_psl.1.0.0.UPF
Ru1     101.07   Ru.rel-pbe-spn-kjpaw_psl.1.0.0.UPF
Ru2     101.07   Ru.rel-pbe-spn-kjpaw_psl.1.0.0.UPF

ATOMIC_POSITIONS crystal
Ru1          0.0000000000       0.0000000000       0.0000000000
Ru2          0.5000000000       0.5000000000       0.5000000000
O            0.1943035720       0.8056964280       0.5000000000
O            0.3056964280       0.3056964280       0.0000000000
O            0.6943035720       0.6943035720       0.0000000000
O            0.8056964280       0.1943035720       0.5000000000

CELL_PARAMETERS {angstrom}
      4.5228052431       0.0000000000       0.0000000000
      0.0000000000       4.5228052431       0.0000000000
      0.0000000000       0.0000000000       3.1218966094

K_POINTS AUTOMATIC
6 6 6 1 1 1

HUBBARD {atomic}
U Ru1-4d 2.0
U Ru2-4d 2.0

Robinson Juma Musembi

unread,
Apr 16, 2024, 3:52:24 PM4/16/24
to Dibya Prakash Rai, BerkeleyGW Help, Meng Wu, Hong Tang
Hello
I have tried to run the scf file but it did not converge also, I checked the quantum espresso manual on the usage of nspin = 4 and noncollinear =.true. , the manual says if noncollinear =.true. then nspin =4 can be omitted, I put a comment on nspin =4 and left noncolin=.true., again I ran the scf file and also it did not converge.
I will continue trying. I am using Quantum Espresso 7.2.

Prof. Robinson J.  Musembi
Solid State and Materials Research
Department of Physics, University of Nairobi
P.O. Box 30197 Nairobi
or P.O. Box 24026 GPO Nairobi 00100
Tel: mobile +254 0725 220 788 : Skype at: siafu2001
Tel: office +254 020 444 75 52 (Department of Physics)

Mailtrack Sender notified by
Mailtrack
16/04/24, 22:47:45

You received this message because you are subscribed to the Google Groups "BerkeleyGW Help" group.
To unsubscribe from this group and stop receiving emails from it, send an email to help+uns...@berkeleygw.org.
To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/86283ba8-5a46-4f82-a2b2-f720eda56728n%40berkeleygw.org.

-----------------------------------------
A KEBS 9001:2015 Certified Organization, No. KEBS/QMS/RF:064 Rev. 03

Reply all
Reply to author
Forward
0 new messages