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