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!--------------------------------------------------------------
9SUBROUTINE calculate_gipaw_orbitals
10  !------------------------------------------------------------
11
12  USE ld1_parameters, ONLY : nwfsx
13  USE ld1inc,         ONLY : nld, file_logder, zed, grid, &
14       wfc_ae_recon, wfc_ps_recon, nspin, nwftsc, lltsc, enltsc, &
15       nwfs, eltsc, el, nwf, enl, rel, vpot, nbeta, pseudotype, &
16       vpstot, nstoaets, enlts, vnl, lls, betas, ddd, qq, jjs, els, &
17       rcutus, ikk, nwfts, verbosity
18  USE io_global,     ONLY : stdout
19  USE kinds,         ONLY : dp
20  USE radial_grids,  ONLY : ndmx, series
21  IMPLICIT NONE
22
23  !<apsi> from lderiv.f90
24  INTEGER ::        &
25       lam,    &   ! the angular momentum
26       nc,     &   ! counter on logarithmic derivatives
27       is,     &   ! counter on spin
28       ierr,   &   ! used for allocation control
29       ios,    &   ! used for I/O control
30       n,ie        ! generic counter
31
32  real(DP) ::            &
33       aux(ndmx),         & ! the square of the wavefunction
34       aux_dir(ndmx,2),   & ! the square of the wavefunction
35       ze2,              & ! the nuclear charge in Ry units
36       e,                & ! the eigenvalue
37       j                   ! total angular momentum for log_der
38
39  INTEGER  ::       &
40       nbf           ! number of b functions
41
42  real(DP) ::  &
43       jam,     &    ! the total angular momentum
44       lamsq,            & ! combined angular momentum
45       b(0:3),c(4),      & ! used for starting guess of the solution
46       b0e, rr1,rr2,     & ! auxiliary
47       xl1, x4l6, ddx12, &
48       x6l12, x8l20
49
50  real(DP) :: &
51       vaux(ndmx),     &  ! auxiliary: the potential
52       al(ndmx)           ! the known part of the differential equation
53
54  real(DP), EXTERNAL :: int_0_inf_dr
55
56  INTEGER :: &
57       ib,jb,iib,jjb, &  ! counters on beta functions
58       nst,nstop,     &  ! auxiliary for integrals
59       ind               ! counters on index
60
61  REAL ( dp ) :: r1, r2, thrdum = 0.0_dp
62
63  INTEGER :: ns, outer_r_radial_index, n_idx, ik, i, r_write_max_index
64  INTEGER :: ik_stop
65  REAL ( dp ) :: factor, rcut_match, max_val
66  REAL ( dp ) :: wfc_ae_recon2(ndmx,nwfsx)
67  REAL ( dp ) :: wfc_ps_recon2(ndmx,nwfsx)
68  REAL ( dp ) :: de
69  REAL ( dp ), ALLOCATABLE :: f_ae(:), f_ps(:)
70
71  IF ( verbosity == 'high' ) THEN
72     WRITE ( stdout, '( 3A )' ) &
73          "     ------------------------", &
74          " (GI)PAW reconstruction ", &
75          "-------------------------"
76  ENDIF
77  IF ( nld > nwfsx ) &
78       CALL errore ( 'calculate_gipaw_orbitals', 'nld is too large', 1 )
79
80  vaux(:) = 0.0_dp
81  ze2=-zed*2.0_dp
82
83  ! Choose a radius for the normalisation of all-electron wave functions
84  ! In principle the value for radius is arbitrary [apsi]
85  ik_stop = 10
86  DO n = 1, grid%mesh
87     IF ( grid%r(n) < 1.5 ) ik_stop = n
88  ENDDO
89
90  DO is=1,nspin
91
92     DO ns = 1, nwftsc(1)
93        lam = lltsc(ns,1)
94        j = 0.0_dp
95
96        DO n = 1, grid%mesh
97           IF ( grid%r(n) < 15.0 ) THEN
98              outer_r_radial_index = n
99           ENDIF
100        ENDDO
101
102        n_idx = -1
103        IF ( abs ( enltsc(ns,1) ) > 1e-8 ) THEN
104           e = enltsc(ns,1)
105           DO n = 1, nwf
106              IF ( eltsc(ns,1) == el(n) ) THEN
107                 n_idx = n
108              ENDIF
109           ENDDO
110        ELSE
111           DO n = 1, nwf
112              IF ( eltsc(ns,1) == el(n) ) THEN
113                 e = enl(n)
114                 n_idx = n
115              ENDIF
116           ENDDO
117        ENDIF
118
119        IF ( verbosity == 'high' ) THEN
120           WRITE ( stdout, '( /, 5X, 3A )') "========================= ", &
121                trim(el(n_idx)), &
122                " ========================="
123           WRITE ( stdout, * ) "     AE: e(ref), l, n_idx ", e, lam, n_idx
124        ENDIF
125
126           !
127           !    integrate outward up to ik_stop
128           !
129           IF (rel == 1) THEN
130              CALL lschps(3, zed, thrdum, grid, outer_r_radial_index, &
131                   1, lam, e, vpot(1,is), aux, nstop)
132           ELSEIF (rel == 2) THEN
133              CALL dir_outward(ndmx,outer_r_radial_index,lam,j,e,grid%dx,&
134                   aux_dir,grid%r,grid%rab,vpot(1,is))
135              aux(:)=aux_dir(:,1)
136           ELSE
137              CALL intref(lam,e,outer_r_radial_index,grid,vpot(1,is),ze2,aux)
138           ENDIF
139
140           wfc_ae_recon(:grid%mesh,n_idx) = aux(:grid%mesh)
141
142           ! Set the maximum to be +-1
143           wfc_ae_recon(:grid%mesh,n_idx) = wfc_ae_recon(:grid%mesh,n_idx) &
144                / maxval ( abs ( wfc_ae_recon(:ik_stop,n_idx) ) )
145
146     ENDDO
147
148  ENDDO
149  !</apsi> from lderiv.f90
150
151  !<apsi> from lderivps.f90
152
153  DO is = 1, nspin
154
155     IF ( .not. rel < 2 ) THEN
156        CALL errore ( 'calculate_gipaw_orbitals', &
157             'not implemented for rel >= 2', rel )
158     ENDIF
159
160     DO ns = 1, nwftsc(1)
161        lam = lltsc(ns,1)
162        jam=0.0_dp
163
164        xl1=lam+1
165        x4l6=4*lam+6
166        x6l12=6*lam+12
167        x8l20=8*lam+20
168        ddx12=grid%dx*grid%dx/12.0_dp
169        nst=(lam+1)*2
170        nbf=nbeta
171        IF (pseudotype == 1) THEN
172           IF (rel < 2 .or. lam == 0 .or. abs(jam-lam+0.5_dp) < 0.001_dp) THEN
173              ind=1
174           ELSEIF (rel==2 .and. lam>0 .and. abs(jam-lam-0.5_dp)<0.001_dp) THEN
175              ind=2
176           ENDIF
177           DO n=1,grid%mesh
178              vaux(n) = vpstot(n,is) + vnl(n,lam,ind)
179           ENDDO
180           nbf=0.0
181        ELSE
182           DO n=1,grid%mesh
183              vaux(n) = vpstot(n,is)
184           ENDDO
185        ENDIF
186
187        DO n=1,4
188           al(n)=vaux(n)-ze2/grid%r(n)
189        ENDDO
190        CALL series(al,grid%r,grid%r2,b)
191
192        IF ( abs ( enltsc(ns,1) ) > 1e-8 ) THEN
193           e = enltsc(ns,1)
194        ELSE
195           e = enlts(ns)
196        ENDIF
197
198        DO n = 1, nwftsc(1)
199           IF ( eltsc(n,1) == el(nstoaets(ns)) ) THEN
200              n_idx = n
201           ENDIF
202        ENDDO
203        IF ( verbosity == 'high' ) THEN
204           WRITE ( stdout, '( /, 5X, 3A )' ) "========================= ", &
205                trim(eltsc(ns,1)), &
206                " ========================="
207           WRITE ( stdout, * ) "     PS: e(ref), e(eig), n_idx ", &
208                e, enlts(ns), n_idx
209        ENDIF
210
211        rcut_match = -1.0_dp
212        DO n = 1, nwfs
213           ! If this one has already a core radius...
214           IF ( els(n) == el(nstoaets(ns)) ) THEN
215              rcut_match = rcutus(n)
216           ENDIF
217        ENDDO
218        IF ( rcut_match < 0.0_dp ) THEN
219           DO n = 1, nwfs
220              ! If there is one with the same l...
221              IF ( lltsc(ns,1) == lls(n) ) THEN
222                 rcut_match = rcutus(n)
223              ENDIF
224           ENDDO
225        ENDIF
226        IF ( rcut_match < 0.0_dp ) THEN
227           max_val = -1.0_dp
228           DO n = 1, grid%mesh
229              IF ( grid%r(n) < 5.0_dp &
230                   .and. abs ( wfc_ae_recon(n,n_idx) ) > max_val ) THEN
231                 rcut_match = grid%r(n)
232                 max_val = abs ( wfc_ae_recon(n,n_idx) )
233              ENDIF
234           ENDDO
235           ! In case over divergence
236           rcut_match = min ( rcut_match, 5.0_dp )
237        ENDIF
238        ik = 10
239        DO n = 1, grid%mesh
240           IF ( grid%r(n) < rcut_match ) ik = n
241        ENDDO
242        IF ( verbosity == 'high' ) THEN
243           WRITE ( stdout, * ) "     r(cut): ", rcut_match, ik, outer_r_radial_index
244        ENDIF
245
246        DO n = 1, grid%mesh
247           IF ( grid%r(n) < 15.0 ) THEN
248              outer_r_radial_index = n
249           ENDIF
250        ENDDO
251
252           lamsq=(lam+0.5_dp)**2
253           !
254           !     b) find the value of solution s in the first two points
255           !
256           b0e=b(0)-e
257           c(1)=0.5_dp*ze2/xl1
258           c(2)=(c(1)*ze2+b0e)/x4l6
259           c(3)=(c(2)*ze2+c(1)*b0e+b(1))/x6l12
260           c(4)=(c(3)*ze2+c(2)*b0e+c(1)*b(1)+b(2))/x8l20
261           r1 = grid%r(1)
262           r2 = grid%r(2)
263           rr1=(1.0_dp+r1*(c(1)+r1*(c(2)+r1*(c(3)+r1*c(4)))))*r1**(lam+1)
264           rr2=(1.0_dp+r2*(c(1)+r2*(c(2)+r2*(c(3)+r2*c(4)))))*r2**(lam+1)
265           aux(1)=rr1/grid%sqr(1)
266           aux(2)=rr2/grid%sqr(2)
267
268           DO n=1,grid%mesh
269              al(n)=( (vaux(n)-e)*grid%r2(n) + lamsq )*ddx12
270              al(n)=1.0_dp-al(n)
271           ENDDO
272
273           CALL integrate_outward (lam,jam,e,grid%mesh,ndmx,grid,al,b,aux, &
274                betas,ddd,qq,nbf,nwfsx,lls,jjs,ikk,outer_r_radial_index)
275           wfc_ps_recon(:grid%mesh,n_idx) = aux(:grid%mesh) &
276                * sqrt ( grid%r(:grid%mesh) )
277
278           DO n = 1, grid%mesh
279              IF ( grid%r(n) > 5.0 ) exit
280           ENDDO
281
282           IF ( abs ( wfc_ps_recon(ik,n_idx) ) < 1e-5 ) THEN
283              WRITE ( stdout, * ) "     Warning: ", wfc_ps_recon(ik,n_idx), ns
284              CALL errore ( "calculate_gipaw_orbitals", &
285                   "safer to stop here...", ik )
286           ENDIF
287
288           factor = wfc_ae_recon(ik,nstoaets(n_idx)) / wfc_ps_recon(ik,n_idx)
289           wfc_ps_recon(:grid%mesh,n_idx) = wfc_ps_recon(:grid%mesh,n_idx) &
290                * factor
291
292           IF ( verbosity == 'high' ) THEN
293              WRITE ( stdout, * ) "     SCALE: ", &
294                   wfc_ae_recon(ik,nstoaets(n_idx)) &
295                   / wfc_ps_recon(ik,n_idx), grid%r(ik), factor
296              WRITE ( stdout, * ) "     SCALE: ", &
297                   wfc_ae_recon(ik+5,nstoaets(n_idx)) &
298                   / wfc_ps_recon(ik+5,n_idx),grid%r(ik+5)
299
300              ALLOCATE ( f_ae(grid%mesh), f_ps(grid%mesh) )
301
302              f_ae = wfc_ae_recon(:grid%mesh,nstoaets(n_idx)) ** 2
303              f_ps = wfc_ps_recon(:grid%mesh,n_idx) ** 2
304
305              ! Test the norm
306              nst = ( lam + 1 ) * 2
307              IF ( verbosity == 'high' ) THEN
308                 WRITE ( stdout, '(A,3F12.8)' ) "     NORM: ", &
309                      int_0_inf_dr ( f_ae, grid, ik, nst ), &
310                      int_0_inf_dr ( f_ps, grid, ik, nst )
311              ENDIF
312
313              DEALLOCATE ( f_ae, f_ps )
314           ENDIF
315
316        ENDDO
317
318  ENDDO
319  !</apsi> from lderivps.f90
320
321  IF ( verbosity == 'high' ) THEN
322     WRITE ( stdout, '( 3A )' ) &
323          "     ---------------------", &
324          " End of (GI)PAW reconstruction ", &
325          "---------------------"
326  ENDIF
327
328END SUBROUTINE calculate_gipaw_orbitals
329