1! 2! Copyright (C) 2004-2011 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! 8MODULE uspp_param 9 ! 10 ! ... Ultrasoft and Norm-Conserving pseudopotential parameters 11 ! 12 USE upf_kinds, ONLY : DP 13 USE upf_params, ONLY : npsx 14 USE pseudo_types, ONLY : pseudo_upf 15 ! 16 SAVE 17 PRIVATE :: randy 18 ! 19 TYPE (pseudo_upf), ALLOCATABLE, TARGET :: upf(:) 20 21 INTEGER :: & 22 nh(npsx), &! number of beta functions per atomic type 23 nhm, &! max number of different beta functions per atom 24 nbetam, &! max number of beta functions 25 iver(3,npsx) ! version of the atomic code 26 INTEGER :: & 27 lmaxkb, &! max angular momentum 28 lmaxq ! max angular momentum + 1 for Q functions 29 INTEGER :: & 30 nvb, &! number of species with Vanderbilt PPs (CPV) 31 ish(npsx) ! for each specie the index of the first beta 32 ! function: ish(1)=1, ish(i)=1+SUM(nh(1:i-1)) 33CONTAINS 34 ! 35 !---------------------------------------------------------------------------- 36 FUNCTION n_atom_wfc( nat, ityp, noncolin ) 37 !---------------------------------------------------------------------------- 38 ! 39 ! ... Find number of starting atomic orbitals 40 ! 41 IMPLICIT NONE 42 ! 43 INTEGER, INTENT(IN) :: nat, ityp(nat) 44 LOGICAL, INTENT(IN), OPTIONAL :: noncolin 45 INTEGER :: n_atom_wfc 46 ! 47 INTEGER :: na, nt, n 48 LOGICAL :: non_col 49 ! 50 ! 51 non_col = .FALSE. 52 IF ( PRESENT (noncolin) ) non_col=noncolin 53 n_atom_wfc = 0 54 ! 55 DO na = 1, nat 56 ! 57 nt = ityp(na) 58 ! 59 DO n = 1, upf(nt)%nwfc 60 ! 61 IF ( upf(nt)%oc(n) >= 0.D0 ) THEN 62 ! 63 IF ( non_col ) THEN 64 ! 65 IF ( upf(nt)%has_so ) THEN 66 ! 67 n_atom_wfc = n_atom_wfc + 2 * upf(nt)%lchi(n) 68 ! 69 IF ( ABS( upf(nt)%jchi(n)-upf(nt)%lchi(n) - 0.5D0 ) < 1.D-6 ) & 70 n_atom_wfc = n_atom_wfc + 2 71 ! 72 ELSE 73 ! 74 n_atom_wfc = n_atom_wfc + 2 * ( 2 * upf(nt)%lchi(n) + 1 ) 75 ! 76 END IF 77 ! 78 ELSE 79 ! 80 n_atom_wfc = n_atom_wfc + 2 * upf(nt)%lchi(n) + 1 81 ! 82 END IF 83 END IF 84 END DO 85 END DO 86 ! 87 RETURN 88 ! 89 END FUNCTION n_atom_wfc 90END MODULE uspp_param 91 92! <<<<<<<<<<<<<<<~~~~<<<<<<<<<<<<<<<<----------------- 93 94MODULE uspp 95 ! 96 ! Ultrasoft PPs: 97 ! - Clebsch-Gordan coefficients "ap", auxiliary variables "lpx", "lpl" 98 ! - beta and q functions of the solid 99 ! 100 USE upf_kinds, ONLY: DP 101 USE upf_params, ONLY: lmaxx, lqmax 102 IMPLICIT NONE 103 PRIVATE 104 SAVE 105 PUBLIC :: nlx, lpx, lpl, ap, aainit, indv, nhtol, nhtolm, indv_ijkb0, & 106 nkb, nkbus, vkb, dvan, deeq, qq_at, qq_nt, nhtoj, ijtoh, beta, & 107 becsum, ebecsum, deallocate_uspp 108 109 PUBLIC :: okvan, nlcc_any 110 PUBLIC :: qq_so, dvan_so, deeq_nc 111 PUBLIC :: dbeta 112 INTEGER, PARAMETER :: & 113 nlx = (lmaxx+1)**2, &! maximum number of combined angular momentum 114 mx = 2*lqmax-1 ! maximum magnetic angular momentum of Q 115 INTEGER :: &! for each pair of combined momenta lm(1),lm(2): 116 lpx(nlx,nlx), &! maximum combined angular momentum LM 117 lpl(nlx,nlx,mx) ! list of combined angular momenta LM 118 REAL(DP) :: & 119 ap(lqmax*lqmax,nlx,nlx) 120 ! Clebsch-Gordan coefficients for spherical harmonics 121 ! 122 INTEGER :: nkb, &! total number of beta functions, with struct.fact. 123 nkbus ! as above, for US-PP only 124 ! 125 INTEGER, ALLOCATABLE ::& 126 indv(:,:), &! indes linking atomic beta's to beta's in the solid 127 nhtol(:,:), &! correspondence n <-> angular momentum l 128 nhtolm(:,:), &! correspondence n <-> combined lm index for (l,m) 129 ijtoh(:,:,:), &! correspondence beta indexes ih,jh -> composite index ijh 130 indv_ijkb0(:) ! first beta (index in the solid) for each atom 131 ! 132 LOGICAL :: & 133 okvan = .FALSE.,& ! if .TRUE. at least one pseudo is Vanderbilt 134 nlcc_any=.FALSE. ! if .TRUE. at least one pseudo has core corrections 135 ! 136 COMPLEX(DP), ALLOCATABLE, TARGET :: & 137 vkb(:,:) ! all beta functions in reciprocal space 138 REAL(DP), ALLOCATABLE :: & 139 becsum(:,:,:) ! \sum_i f(i) <psi(i)|beta_l><beta_m|psi(i)> 140 REAL(DP), ALLOCATABLE :: & 141 ebecsum(:,:,:) ! \sum_i f(i) et(i) <psi(i)|beta_l><beta_m|psi(i)> 142 REAL(DP), ALLOCATABLE :: & 143 dvan(:,:,:), &! the D functions of the solid 144 deeq(:,:,:,:), &! the integral of V_eff and Q_{nm} 145 qq_nt(:,:,:), &! the integral of q functions in the solid (ONE PER NTYP) used to be the qq array 146 qq_at(:,:,:), &! the integral of q functions in the solid (ONE PER ATOM !!!!) 147 nhtoj(:,:) ! correspondence n <-> total angular momentum 148 ! 149 COMPLEX(DP), ALLOCATABLE :: & ! variables for spin-orbit/noncolinear case: 150 qq_so(:,:,:,:), &! Q_{nm} 151 dvan_so(:,:,:,:), &! D_{nm} 152 deeq_nc(:,:,:,:) ! \int V_{eff}(r) Q_{nm}(r) dr 153 ! 154 ! spin-orbit coupling: qq and dvan are complex, qq has additional spin index 155 ! noncolinear magnetism: deeq is complex (even in absence of spin-orbit) 156 ! 157 REAL(DP), ALLOCATABLE :: & 158 beta(:,:,:) ! beta functions for CP (without struct.factor) 159 REAL(DP), ALLOCATABLE :: & 160 dbeta(:,:,:,:,:) ! derivative of beta functions w.r.t. cell for CP (without struct.factor) 161 ! 162CONTAINS 163 ! 164 !----------------------------------------------------------------------- 165 subroutine aainit(lli) 166 !----------------------------------------------------------------------- 167 ! 168 ! this routine computes the coefficients of the expansion of the product 169 ! of two real spherical harmonics into real spherical harmonics. 170 ! 171 ! Y_limi(r) * Y_ljmj(r) = \sum_LM ap(LM,limi,ljmj) Y_LM(r) 172 ! 173 ! On output: 174 ! ap the expansion coefficients 175 ! lpx for each input limi,ljmj is the number of LM in the sum 176 ! lpl for each input limi,ljmj points to the allowed LM 177 ! 178 ! The indices limi,ljmj and LM assume the order for real spherical 179 ! harmonics given in routine ylmr2 180 ! 181 USE upf_invmat 182 implicit none 183 ! 184 ! input: the maximum li considered 185 ! 186 integer :: lli 187 ! 188 ! local variables 189 ! 190 integer :: llx, l, li, lj 191 real(DP) , allocatable :: r(:,:), rr(:), ylm(:,:), mly(:,:) 192 ! an array of random vectors: r(3,llx) 193 ! the norm of r: rr(llx) 194 ! the real spherical harmonics for array r: ylm(llx,llx) 195 ! the inverse of ylm considered as a matrix: mly(llx,llx) 196 ! 197 if (lli < 0) call upf_error('aainit','lli not allowed',lli) 198 199 if (lli*lli > nlx) call upf_error('aainit','nlx is too small ',lli*lli) 200 201 llx = (2*lli-1)**2 202 if (2*lli-1 > lqmax) & 203 call upf_error('aainit','ap leading dimension is too small',llx) 204 205 allocate (r( 3, llx )) 206 allocate (rr( llx )) 207 allocate (ylm( llx, llx )) 208 allocate (mly( llx, llx )) 209 210 r(:,:) = 0.0_DP 211 ylm(:,:) = 0.0_DP 212 mly(:,:) = 0.0_DP 213 ap(:,:,:)= 0.0_DP 214 215 ! - generate an array of random vectors (uniform deviate on unitary sphere) 216 217 call gen_rndm_r(llx,r,rr) 218 219 ! - generate the real spherical harmonics for the array: ylm(ir,lm) 220 221 call ylmr2(llx,llx,r,rr,ylm) 222 223 !- store the inverse of ylm(ir,lm) in mly(lm,ir) 224 225 call invmat(llx, ylm, mly) 226 227 !- for each li,lj compute ap(l,li,lj) and the indices, lpx and lpl 228 do li = 1, lli*lli 229 do lj = 1, lli*lli 230 lpx(li,lj)=0 231 do l = 1, llx 232 ap(l,li,lj) = compute_ap(l,li,lj,llx,ylm,mly) 233 if (abs(ap(l,li,lj)) > 1.d-3) then 234 lpx(li,lj) = lpx(li,lj) + 1 235 if (lpx(li,lj) > mx) & 236 call upf_error('aainit','mx dimension too small', lpx(li,lj)) 237 lpl(li,lj,lpx(li,lj)) = l 238 end if 239 end do 240 end do 241 end do 242 243 deallocate(mly) 244 deallocate(ylm) 245 deallocate(rr) 246 deallocate(r) 247 248 return 249 end subroutine aainit 250 ! 251 !----------------------------------------------------------------------- 252 subroutine gen_rndm_r(llx,r,rr) 253 !----------------------------------------------------------------------- 254 ! - generate an array of random vectors (uniform deviate on unitary sphere) 255 ! 256 USE upf_const, ONLY: tpi 257 implicit none 258 ! 259 ! first the I/O variables 260 ! 261 integer :: llx ! input: the dimension of r and rr 262 real(DP) :: & 263 r(3,llx), &! output: an array of random vectors 264 rr(llx) ! output: the norm of r 265 ! 266 ! here the local variables 267 ! 268 integer :: ir 269 real(DP) :: costheta, sintheta, phi 270 271 do ir = 1, llx 272 costheta = 2.0_DP * randy() - 1.0_DP 273 sintheta = SQRT ( 1.0_DP - costheta*costheta) 274 phi = tpi * randy() 275 r (1,ir) = sintheta * cos(phi) 276 r (2,ir) = sintheta * sin(phi) 277 r (3,ir) = costheta 278 rr(ir) = 1.0_DP 279 end do 280 281 return 282 end subroutine gen_rndm_r 283! 284!------------------------------------------------------------------------ 285 FUNCTION randy ( irand ) 286 !------------------------------------------------------------------------ 287 ! 288 ! x=randy(n): reseed with initial seed idum=n ( 0 <= n <= ic, see below) 289 ! if randy is not explicitly initialized, it will be 290 ! initialized with seed idum=0 the first time it is called 291 ! x=randy() : generate uniform real(DP) numbers x in [0,1] 292 ! 293 use upf_kinds, only : DP 294 implicit none 295 ! 296 REAL(DP) :: randy 297 INTEGER, optional :: irand 298 ! 299 INTEGER , PARAMETER :: m = 714025, & 300 ia = 1366, & 301 ic = 150889, & 302 ntab = 97 303 REAL(DP), PARAMETER :: rm = 1.0_DP / m 304 INTEGER :: j 305 INTEGER, SAVE :: ir(ntab), iy, idum=0 306 LOGICAL, SAVE :: first=.true. 307 ! 308 IF ( present(irand) ) THEN 309 idum = MIN( ABS(irand), ic) 310 first=.true. 311 END IF 312 313 IF ( first ) THEN 314 ! 315 first = .false. 316 idum = MOD( ic - idum, m ) 317 ! 318 DO j=1,ntab 319 idum=mod(ia*idum+ic,m) 320 ir(j)=idum 321 END DO 322 idum=mod(ia*idum+ic,m) 323 iy=idum 324 END IF 325 j=1+(ntab*iy)/m 326 IF( j > ntab .OR. j < 1 ) call upf_error('randy','j out of range',ABS(j)+1) 327 iy=ir(j) 328 randy=iy*rm 329 idum=mod(ia*idum+ic,m) 330 ir(j)=idum 331 ! 332 RETURN 333 ! 334 END FUNCTION randy 335 ! 336 !----------------------------------------------------------------------- 337 function compute_ap(l,li,lj,llx,ylm,mly) 338 !----------------------------------------------------------------------- 339 !- given an l and a li,lj pair compute ap(l,li,lj) 340 implicit none 341 ! 342 ! first the I/O variables 343 ! 344 integer :: & 345 llx, &! the dimension of ylm and mly 346 l,li,lj ! the arguments of the array ap 347 348 real(DP) :: & 349 compute_ap, &! this function 350 ylm(llx,llx),&! the real spherical harmonics for array r 351 mly(llx,llx) ! the inverse of ylm considered as a matrix 352 ! 353 ! here the local variables 354 ! 355 integer :: ir 356 357 compute_ap = 0.0_DP 358 do ir = 1,llx 359 compute_ap = compute_ap + mly(l,ir)*ylm(ir,li)*ylm(ir,lj) 360 end do 361 362 return 363 end function compute_ap 364 ! 365 !----------------------------------------------------------------------- 366 SUBROUTINE deallocate_uspp() 367 !----------------------------------------------------------------------- 368 ! 369 IF( ALLOCATED( nhtol ) ) DEALLOCATE( nhtol ) 370 IF( ALLOCATED( indv ) ) DEALLOCATE( indv ) 371 IF( ALLOCATED( nhtolm ) ) DEALLOCATE( nhtolm ) 372 IF( ALLOCATED( nhtoj ) ) DEALLOCATE( nhtoj ) 373 IF( ALLOCATED( indv_ijkb0 ) ) DEALLOCATE( indv_ijkb0 ) 374 IF( ALLOCATED( ijtoh ) ) DEALLOCATE( ijtoh ) 375 IF( ALLOCATED( vkb ) ) DEALLOCATE( vkb ) 376 IF( ALLOCATED( becsum ) ) DEALLOCATE( becsum ) 377 IF( ALLOCATED( ebecsum ) ) DEALLOCATE( ebecsum ) 378 IF( ALLOCATED( qq_at ) ) DEALLOCATE( qq_at ) 379 IF( ALLOCATED( qq_nt ) ) DEALLOCATE( qq_nt ) 380 IF( ALLOCATED( dvan ) ) DEALLOCATE( dvan ) 381 IF( ALLOCATED( deeq ) ) DEALLOCATE( deeq ) 382 IF( ALLOCATED( qq_so ) ) DEALLOCATE( qq_so ) 383 IF( ALLOCATED( dvan_so ) ) DEALLOCATE( dvan_so ) 384 IF( ALLOCATED( deeq_nc ) ) DEALLOCATE( deeq_nc ) 385 IF( ALLOCATED( beta ) ) DEALLOCATE( beta ) 386 IF( ALLOCATED( dbeta ) ) DEALLOCATE( dbeta ) 387 ! 388 END SUBROUTINE deallocate_uspp 389 ! 390END MODULE uspp 391 392