1!
2! Copyright (C) 2004i-2010 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!---------------------------------------------------------------
9SUBROUTINE scf(ic)
10  !---------------------------------------------------------------
11  !
12  !   this routine performs the atomic self-consistent procedure
13  !   self-interaction-correction allowed
14  !
15  USE kinds, ONLY : dp
16  USE funct, ONLY : dft_is_meta
17  USE radial_grids, ONLY : ndmx
18  USE constants, ONLY: e2
19  USE ld1inc, ONLY : grid, zed, psi, isic, vpot, vh, vxt, rho, iter, &
20                     lsd, rel, latt, enne, beta, nspin, tr2, eps0, &
21                     nwf, nn, ll, jj, enl, oc, isw, core_state, frozen_core, &
22                     tau, vtau, vsic, vsicnew, vhn1, egc, relpert, noscf
23  IMPLICIT NONE
24
25  INTEGER, INTENT(in) :: ic
26
27  LOGICAL:: meta, conv
28  INTEGER:: nerr, nstop, n, i, is, id, nin
29  real(DP) ::  vnew(ndmx,2), vtaunew(ndmx), rhoc1(ndmx), ze2
30  INTEGER, PARAMETER :: maxter=200
31  real(DP), PARAMETER :: thresh=1.0e-10_dp
32  !
33  !
34  meta = dft_is_meta()
35  ze2 = - zed * e2
36  rhoc1=0.0_dp
37  IF (.not.frozen_core.or.ic==1) psi=0.0_dp
38  DO iter = 1,maxter
39     nerr=0
40     vnew=vpot
41     vtaunew=vtau
42     DO n=1,nwf
43        IF (oc(n) >= 0.0_dp) THEN
44           IF (ic==1.or..not.frozen_core.or..not.core_state(n)) THEN
45              is=isw(n)
46              IF (isic /= 0 .and. iter > 1) vnew(:,is)=vpot(:,is)-vsic(:,n)
47              IF (rel == 0) THEN
48                 ! nonrelativistic calculation
49                 IF ( meta ) THEN
50                   ! Meta-GGA version of lschps
51                    CALL lschps_meta (2, zed, thresh, grid, nin, nn(n), ll(n),&
52                         enl(n), vnew(1,is), vtaunew, psi(1,1,n), nstop)
53                 ELSE
54                !     print  *, "Solving nonrelativistic nonmeta equation"
55                    CALL ascheq (nn(n),ll(n),enl(n),grid%mesh,grid,&
56                      vnew(1,is), & ! potential
57                      ze2,thresh,psi(1,1,n),nstop)
58                 END IF
59              ELSEIF (rel == 1) THEN
60                 ! relativistic scalar calculation
61                 IF ( meta ) THEN
62                    CALL lschps_meta (1, zed, thresh, grid, nin, nn(n), ll(n),&
63                         enl(n), vnew(1,is), vtaunew, psi(1,1,n), nstop)
64                 ELSE
65                    CALL lschps (1, zed, thresh, grid, nin, nn(n), ll(n),&
66                         enl(n), vnew(1,is), psi(1,1,n), nstop)
67                 END IF
68                 IF (nstop>0.and.oc(n)<1.e-10_DP) nstop=0
69              ELSEIF (rel == 2) THEN
70                 CALL dirsol (ndmx,grid%mesh,nn(n),ll(n),jj(n),iter,enl(n), &
71                      thresh,grid,psi(1,1,n),vnew(1,is),nstop)
72              ELSE
73                 CALL errore('scf','relativistic not programmed',1)
74              ENDIF
75              !      write(6,*) nn(n),ll(n),enl(n)
76              ! if (nstop /= 0) write(6,'(4i6)') iter,nn(n),ll(n),nstop
77              nerr=nerr+nstop
78           ENDIF
79        ELSE
80           enl(n)=0.0_dp
81           psi(:,:,n)=0.0_dp
82        ENDIF
83     ENDDO
84
85     !
86     ! calculate charge density (spherical approximation)
87     !
88     rho=0.0_dp
89
90     IF (noscf) GOTO 500
91     DO n=1,nwf
92        rho(1:grid%mesh,isw(n))=rho(1:grid%mesh,isw(n)) + &
93        oc(n)*(psi(1:grid%mesh,1,n)**2+psi(1:grid%mesh,2,n)**2)
94     ENDDO
95
96
97     !
98     ! calculate kinetc energy density (spherical approximation)
99     !
100     IF ( meta ) CALL kin_e_density (ndmx, grid%mesh, nwf, &
101         ll, oc, psi, grid%r, grid%r2, grid%dx, tau)
102     !
103     ! calculate new potential
104     !
105
106     CALL new_potential ( ndmx, grid%mesh, grid, zed, vxt, &
107          lsd, .false., latt, enne, rhoc1, rho, vh, vnew, 1 )
108     !
109     ! calculate SIC correction potential (if present)
110     !
111     IF (isic /= 0) THEN
112        DO n=1,nwf
113           IF (oc(n) >= 0.0_dp) THEN
114              is=isw(n)
115              CALL sic_correction(n,vhn1,vsicnew,egc)
116              !
117              ! use simple mixing for SIC correction
118              !
119              vsic(:,n) = (1.0_dp-beta)*vsic(:,n)+beta*vsicnew(:)
120           ENDIF
121        ENDDO
122     ENDIF
123     !
124     ! mix old and new potential
125     !
126     id = 3
127     IF (isic /= 0 .and. relpert)  id=1
128     !
129     CALL vpack(grid%mesh,ndmx,nspin,vnew,vpot,1)
130     CALL dmixp(grid%mesh*nspin,vnew,vpot,beta,tr2,iter,id,eps0,conv,maxter)
131     CALL vpack(grid%mesh,ndmx,nspin,vnew,vpot,-1)
132!        write(6,*) iter, eps0
133     !
134     ! mix old and new metaGGA potential - use simple mixing
135     !
136     IF ( meta ) vtau(:) = (1.0_dp-beta)*vtaunew(:)+beta*vtau(:)
137     !
138500  IF (noscf) THEN
139        conv=.true.
140        eps0 = 0.0_DP
141     ENDIF
142     IF (conv) THEN
143        IF (nerr /= 0) CALL infomsg ('scf','warning: at least one error in KS equations')
144        EXIT ! exit cycle
145     ENDIF
146  ENDDO
147  IF ( .not. conv ) CALL infomsg('scf','warning: convergence not achieved')
148
149END SUBROUTINE scf
150