1!
2! Copyright (C) 2004 PWSCF 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!
9!---------------------------------------------------------------
10SUBROUTINE run_pseudo
11  !---------------------------------------------------------------
12  !
13  !     this routine is a driver to a pseudopotential calculation
14  !     with the parameters given in input
15  !
16  !
17  USE kinds, ONLY : dp
18  USE radial_grids, ONLY : ndmx
19  USE ld1_parameters, ONLY : nwfsx
20  USE ld1inc, ONLY : enl, lpaw, nlcc, lsd, latt, pawsetup, &
21                     nstoaets, grid, nspin, iter, rhos, rhoc, &
22                     nwfts, enlts, llts, jjts, iswts, octs, phits, &
23                     vxt, enne, vh, vpsloc, file_potscf, beta, tr2,  &
24                     eps0, file_recon, deld, vpstot, nbeta, ddd, etots, &
25                     paw_energy, iswitch, lgipaw_reconstruction, use_paw_as_gipaw
26  USE atomic_paw, ONLY : new_paw_hamiltonian
27  IMPLICIT NONE
28
29  INTEGER :: &
30       ns, &    ! counter on pseudowavefunctions
31       n,  &    ! counter on mesh
32       is       ! counter on spin
33
34  real(DP) :: &
35       vnew(ndmx,2)   ! the potential
36
37  INTEGER :: &
38       iunps, ios
39
40  LOGICAL :: &
41       conv       ! if true convergence reached
42
43  real(DP) :: &
44       dddnew(nwfsx,nwfsx,2),   & ! the new D coefficients
45       vd(2*(ndmx+nwfsx+nwfsx)), & ! Vloc and D in one array for mixing
46       vdnew(2*(ndmx+nwfsx+nwfsx)) ! the new vd array
47  INTEGER :: &
48       nerr                       ! error message
49
50  real(DP), PARAMETER :: thresh=1.e-10_dp
51  INTEGER, PARAMETER :: itmax=200
52  CHARACTER(len=256) :: nomefile
53
54  !
55  !     initial estimate of the eigenvalues
56  !
57  DO ns=1,nwfts
58     enlts(ns)=enl(nstoaets(ns))
59  ENDDO
60  !
61  !    compute an initial estimate of the potential
62  !
63  CALL guess_initial_wfc()
64  IF (.not.lpaw) THEN
65     CALL start_potps ( )
66  ELSE
67     CALL new_paw_hamiltonian (vpstot, ddd, etots,pawsetup, nwfts, &
68                               llts, jjts, nspin, iswts, octs, phits, enlts)
69     DO is=1,nspin
70        vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)-pawsetup%psloc(1:grid%mesh)
71     ENDDO
72     CALL vdpack (grid%mesh, ndmx, nbeta, nwfsx, nspin, vpstot, ddd, vd, "PACK")
73     DO is=1,nspin
74        vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)+pawsetup%psloc(1:grid%mesh)
75     ENDDO
76  ENDIF
77  !
78  !     iterate to self-consistency
79  !
80  DO iter=1,itmax
81     CALL ascheqps_drv(vpstot, nspin, thresh, .false., nerr)
82
83     IF (.not.lpaw) THEN
84        !
85        CALL chargeps(rhos,phits,nwfts,llts,jjts,octs,iswts)
86        CALL new_potential(ndmx,grid%mesh,grid,0.0_dp,vxt,lsd,&
87             nlcc,latt,enne,rhoc,rhos,vh,vnew,1)
88
89        DO is=1,nspin
90           vpstot(:,is)=vpstot(:,is)-vpsloc(:)
91        ENDDO
92
93        IF (file_potscf/=' ') THEN
94           IF (iter<10) THEN
95              WRITE(nomefile,'(a,"_",i1)') trim(file_potscf), iter
96           ELSEIF(iter<100) THEN
97              WRITE(nomefile,'(a,"_",i2)') trim(file_potscf), iter
98           ELSEIF(iter<1000) THEN
99              WRITE(nomefile,'(a,"_",i3)') trim(file_potscf), iter
100           ELSE
101              CALL errore('run_pseudo','problem with iteration',1)
102           ENDIF
103           OPEN(unit=18,file=trim(nomefile), status='unknown', &
104                                             err=100, iostat=ios)
105100        CALL errore('run_pseudo','opening file' // nomefile,abs(ios))
106           IF (lsd==1) THEN
107              DO n=1,grid%mesh
108                 WRITE(18,'(5e20.12)') grid%r(n),vnew(n,1)-vpstot(n,1), &
109                         vnew(n,1), vnew(n,2)-vpstot(n,2), vnew(n,2)
110              ENDDO
111           ELSE
112              DO n=1,grid%mesh
113                 WRITE(18,'(3e26.15)') grid%r(n),vnew(n,1)-vpstot(n,1), &
114                          vnew(n,1)
115              ENDDO
116           ENDIF
117           CLOSE(18)
118        ENDIF
119
120        CALL vpack(grid%mesh,ndmx,nspin,vnew,vpstot,1)
121        CALL dmixp(grid%mesh*nspin,vnew,vpstot,beta,tr2,iter,3,eps0,conv,itmax)
122        CALL vpack(grid%mesh,ndmx,nspin,vnew,vpstot,-1)
123
124        DO is=1,nspin
125           DO n=1,grid%mesh
126              vpstot(n,is)=vpstot(n,is)+vpsloc(n)
127           ENDDO
128        ENDDO
129        CALL newd_at
130        !
131     ELSE
132        !
133        CALL new_paw_hamiltonian (vnew, dddnew, etots, &
134             pawsetup, nwfts, llts, jjts, nspin, iswts, octs, phits, enlts)
135        DO is=1,nspin
136           vnew(1:grid%mesh,is)=vnew(1:grid%mesh,is)-pawsetup%psloc(1:grid%mesh)
137        ENDDO
138        CALL vdpack (grid%mesh, ndmx, nbeta, nwfsx, nspin, vnew, dddnew, vdnew, "PACK")
139        CALL dmixp((grid%mesh+nbeta*nbeta)*nspin,vdnew,vd,beta,tr2,iter,3,eps0,conv,itmax)
140        CALL vdpack (grid%mesh, ndmx, nbeta, nwfsx, nspin, vpstot, ddd, vd, "UNDO")
141        DO is=1,nspin
142           vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)+pawsetup%psloc(1:grid%mesh)
143        ENDDO
144        !
145     ENDIF
146
147     IF (conv) THEN
148        IF (nerr /= 0) THEN
149           IF (iswitch==2) THEN
150              CALL infomsg ('run_pseudo','BEWARE! Errors in PS-KS equations')
151           ELSE
152              CALL errore ('run_pseudo','Errors in PS-KS equation', 1)
153           ENDIF
154        ENDIF
155        GOTO 900
156     ENDIF
157  ENDDO
158  CALL infomsg('run_pseudo','Warning: convergence not achieved')
159  !
160  !    final calculation with all states
161  !
162900 CONTINUE
163
164  CALL ascheqps_drv(vpstot, nspin, thresh, .true., nerr)
165
166  IF (.not.lpaw) THEN
167     CALL elsdps ( )
168  ELSE
169     CALL new_paw_hamiltonian (vnew, dddnew, etots, pawsetup, nwfts, &
170                llts, jjts, nspin, iswts, octs, phits, enlts, paw_energy)
171     CALL elsdps_paw()
172  ENDIF
173
174  IF ( lgipaw_reconstruction.and.(.not.use_paw_as_gipaw) ) &
175       CALL calculate_gipaw_orbitals()
176
177  IF (file_recon/=' ') CALL write_paw_recon ( )
178  !
179  !    compute logarithmic derivatives
180  !
181  IF ( deld > 0.0_dp) CALL lderivps ( )
182  !
183  !    compute expansion in partial waves
184  !
185  IF ( deld > 0.0_dp) CALL partial_wave_expansion ( )
186
187  RETURN
188END SUBROUTINE run_pseudo
189