1 2! Copyright (C) 2001-2007 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! 9!----------------------------------------------------------------------------------- 10SUBROUTINE init_us_0 11 !--------------------------------------------------------------------------------- 12 !! This routine performs the following task: for each uspp or paw pseudopotential 13 !! the l-dependent aumentation charge \(\text{ q_nb_mb_l}\)(r), stored in 14 !! \(\text{qfuncl}\)(ir,nmb,l), is: 15 ! 16 !! * transformed in reciprocal space by bessel transform up to qmax = sqrt(ecutrho); 17 !! * smoothed by multiplying with a filter function \(\textrm{filter}\)(q/qmax,a,nn); 18 !! * brought back in real space, 19 ! 20 !! where it overwrites the original array. The filter function is: 21 !! \[ \text{filter}(x,a,\text{nn}) = e^{-\text{axx}} \sum_{k=0,\text{nn}} 22 !! \frac{\text{axx}^k}{k!}\ . \] 23 ! 24 USE kinds, ONLY: DP 25 USE gvect, ONLY: ecutrho 26 USE io_global, ONLY: stdout 27 USE constants, ONLY: fpi, sqrt2, eps8, eps6 28 USE atom, ONLY: rgrid 29 USE ions_base, ONLY: ntyp => nsp 30 USE cell_base, ONLY: omega, tpiba 31 USE us, ONLY: dq 32 USE uspp_param, ONLY: upf, lmaxq, nbetam 33 USE mp_bands, ONLY: intra_bgrp_comm 34 USE mp, ONLY: mp_sum 35 ! 36 IMPLICIT NONE 37 ! 38 ! ... local variables 39 ! 40 ! sdg 41 ! LOGICAL, PARAMETER :: tprint=.true. ! whether the q_l(r) and its relatives are printed or not 42 LOGICAL, PARAMETER :: tprint=.FALSE. ! whether the q_l(r) and its relatives are printed or not 43 INTEGER, PARAMETER :: nn=16 ! smoothing parameter, order of the polynomial inverse gaussian 44 ! approximant 45 REAL(DP), PARAMETER:: a=22.0_DP ! smoothing parameter, exponent of the gaussian decaying factor 46 ! a=0.0 ; nn=0.0 would be no smoothing. 47 ! 48 INTEGER :: nt, ih, jh, nb, mb, ijv, l, m, ir, ir0, iq, is, startq, lastq, ilast, ndm, ia 49 ! various counters 50 INTEGER :: nqxq 51 REAL(DP), ALLOCATABLE :: qrad_q(:,:,:), qrad_r(:,:,:), qrad_rs(:,:,:), ffrr(:) 52 REAL(DP), ALLOCATABLE :: power_0(:,:), power_r(:,:), power_q(:,:), power_rs(:,:), power_qs(:,:) 53 REAL(DP), ALLOCATABLE :: aux(:), aux1(:) 54 ! various work space 55 REAL(DP) :: q, qi 56 ! the modulus of g for each shell 57 ! q-point grid for interpolation 58 REAL(DP), ALLOCATABLE :: ylmk0(:) 59 ! the spherical harmonics 60 INTEGER, EXTERNAL :: sph_ind 61 INTEGER :: lnb, lmb 62 REAL(DP) :: qmax, rcut, drcut 63 REAL(DP) :: target_ratio, ratio, ratio_s, fac 64 ! 65 CHARACTER(LEN=6) :: filename 66 ! 67 CALL start_clock( 'init_us_0' ) 68 ! 69 ! ... Initialization of the variables 70 ! 71 drcut = ABS(LOG(eps8))/SQRT(ecutrho) 72 ! 73 DO nt = 1, ntyp 74 ! 75 rcut = rgrid(nt)%r(upf(nt)%kkbeta) 76 WRITE (stdout,*) 'RCUT:', rcut, drcut, rcut+drcut 77 DO ir = upf(nt)%kkbeta, upf(nt)%mesh 78 IF ( rgrid(nt)%r(ir) < rcut+drcut ) upf(nt)%kkbeta=ir 79 ENDDO 80 ! 81 ENDDO 82 ! 83 ndm = MAXVAL( upf(:)%kkbeta ) 84 nqxq = INT( SQRT(ecutrho) / dq + 4 ) 85 ! 86 IF (tprint) THEN 87 WRITE (stdout,*) " PSEUDOPOTENTIAL REPORT " 88 WRITE (stdout,*) ' NDM :', ndm, ' ', upf(1:ntyp)%kkbeta 89 WRITE (stdout,*) ' LMAXQ :', lmaxq, ' NBETAM :', nbetam 90 ENDIF 91 ! 92 ALLOCATE( aux(ndm), aux1(ndm) ) 93 ALLOCATE( qrad_q(nqxq, nbetam*(nbetam+1)/2,lmaxq) ) 94 ALLOCATE( qrad_r(ndm , nbetam*(nbetam+1)/2,lmaxq) ) 95 ALLOCATE( qrad_rs(ndm, nbetam*(nbetam+1)/2,lmaxq) ) 96 ALLOCATE( ffrr(ndm) ) 97 ALLOCATE( power_0(nbetam*(nbetam+1)/2 ,0:lmaxq) ) 98 ALLOCATE( power_r(nbetam*(nbetam+1)/2 ,0:lmaxq), power_q(nbetam*(nbetam+1)/2 ,0:lmaxq) ) 99 ALLOCATE( power_rs(nbetam*(nbetam+1)/2,0:lmaxq), power_qs(nbetam*(nbetam+1)/2,0:lmaxq) ) 100 ALLOCATE( ylmk0(lmaxq*lmaxq) ) 101 ! 102 ! ... the following prevents an out-of-bound error: upf(nt)%nqlc=2*lmax+1 103 ! but in some versions of the PP files lmax is not set to the maximum 104 ! l of the beta functions but includes the l of the local potential 105 ! 106 DO nt = 1, ntyp 107 ! 108 upf(nt)%nqlc = MIN( upf(nt)%nqlc, lmaxq ) 109 IF (upf(nt)%nqlc < 0) upf(nt)%nqlc = 0 110 ! 111 ENDDO 112 ! 113 ! ... Here for the US types we compute the Fourier transform of the 114 ! Q functions. 115 ! 116 CALL divide( intra_bgrp_comm, nqxq, startq, lastq ) 117 ! 118 qmax = SQRT(ecutrho) 119 WRITE (stdout, *) ' qmax : sqrt(ecutrho) =', SQRT(ecutrho), dq*nqxq 120 WRITE (stdout,'(a,f6.2,a,i4,4(a,f11.8))') 'FILTER : a=',a,', nn=',nn, & 121 ', filter(1.1d0)=', filter(1.1d0,a,nn), & 122 ', filter(1.0d0)=', filter(1.0d0,a,nn), & 123 ', filter(0.9d0)=', filter(0.9d0,a,nn), & 124 ', filter(0.8d0)=', filter(0.8d0,a,nn) 125 ! 126 DO nt = 1, ntyp 127 WRITE (stdout,*) ' NT = ', nt 128 ! 129 IF ( upf(nt)%tvanp ) THEN 130 !- 131 IF (tprint) THEN 132 filename = 'qr_' ! the radial q_l(r) as defined in the pseudopotential 133 ! 134 DO nb = 1, upf(nt)%nbeta 135 DO mb = nb, upf(nt)%nbeta 136 ! 137 ijv = mb * (mb - 1) / 2 + nb 138 lnb = upf(nt)%lll(nb) ; lmb = upf(nt)%lll(mb) 139 WRITE (filename(4:4),'(i1)') nb 140 WRITE (filename(5:5),'(i1)') mb 141 WRITE (filename(6:6),'(i1)') nt 142 OPEN (4,file=filename,form='formatted', status='unknown') 143 WRITE (4, *) '# the radial q_l(r) as defined in the pseudopotential' 144 WRITE (4, *) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, & 145 ' kkbeta :',upf(nt)%kkbeta 146 ! 147 DO ir =1,upf(nt)%kkbeta 148 write (4,'(12f16.10)') rgrid(nt)%r(ir), (upf(nt)%qfuncl(ir,ijv,l), l=0, lnb+lmb) 149 ENDDO 150 ! 151 CLOSE (4) 152 ! 153 ENDDO 154 ENDDO 155 ! 156 ENDIF 157 !- 158 qrad_q(:,:,:) = 0.0_DP ; qrad_r(:,:,:) = 0.0_DP 159 power_0(:,:) = 0.0_DP ; power_r(:,:) = 0.0_DP ; power_q(:,:) = 0.0_DP 160 ! 161 DO l = 0, upf(nt)%nqlc-1 162 ! 163 ! 1) first of all compute the integrated power spectrum of the Qs in real space. 164 ! 165 DO nb = 1, upf(nt)%nbeta 166 DO mb = nb, upf(nt)%nbeta 167 ! 168 ijv = mb * (mb - 1) / 2 + nb 169 aux1(1) = 0.0_DP 170 ! 171 ir0 = 1 172 IF (rgrid(nt)%r(ir0) < eps8) ir0 = 2 173 ! 174 DO ir = ir0, upf(nt)%kkbeta 175 aux1(ir) = (upf(nt)%qfuncl(ir,ijv,l)/rgrid(nt)%r(ir))**2 176 ENDDO 177 ! 178 CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, power_0(ijv,l+1) ) 179 ! 180 ENDDO 181 ENDDO 182 ! 183 ! 2) compute the fourier transform of the Qs and their integrated power spectum in 184 ! reciprocal space - FIXME: use routine compute_qrad in init_us_1 185 ! 186 DO iq = startq, lastq 187 ! 188 q = (iq - 1) * dq 189 ! 190 ! ... here we compute the spherical bessel function for the given |q| 191 ! 192 CALL sph_bes( upf(nt)%kkbeta, rgrid(nt)%r, q, l, aux ) 193 ! 194 ! .. and here we integrate with all the Q functions 195 ! 196 DO nb = 1, upf(nt)%nbeta 197 DO mb = nb, upf(nt)%nbeta 198 ! 199 ijv = mb * (mb - 1) / 2 + nb 200 lnb = upf(nt)%lll(nb) 201 lmb = upf(nt)%lll(mb) 202 ! 203 IF ( (l >= ABS(lnb-lmb)) .AND. (l <= lnb+lmb) .AND. (MOD(l+lnb+lmb,2)==0) ) THEN 204 ! 205 DO ir = 1, upf(nt)%kkbeta 206 aux1 (ir) = aux (ir) * upf(nt)%qfuncl(ir,ijv,l) 207 ENDDO 208 ! 209 ! ... computes the Fourier transform of q_nb_mb_l(r) up to qmax 210 ! 211 CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, qrad_q(iq,ijv,l+1) ) 212 ! 213 ! ... and update the integrated power spectrum in reciprocal space 214 ! 215 power_q(ijv,l+1) = power_q(ijv,l+1) + q*q * dq * qrad_q(iq,ijv,l+1)**2 216 ! 217 ENDIF 218 ! 219 ENDDO 220 ENDDO 221 ! 222 ! 3) back-fourier transform of the Qs to real space. 223 ! 224 DO ir =1, upf(nt)%kkbeta 225 ! 226 ! ... q_nb_mb_l(r) from the back fourier transform up to qmax of q_nb_mb_l(q) 227 ! 228 qrad_r(ir,1:nbetam*(nbetam+1)/2,l+1) = qrad_r(ir,1:nbetam*(nbetam+1)/2,l+1) + & 229 q*q * dq * aux(ir) * rgrid(nt)%r(ir)**2 * & 230 qrad_q(iq,1:nbetam*(nbetam+1)/2,l+1) 231 ENDDO 232 ! 233 ENDDO 234 ! 235 ENDDO 236 ! 237 CALL mp_sum( qrad_q , intra_bgrp_comm ) 238 CALL mp_sum( power_q, intra_bgrp_comm ) ; power_q (:,:) = power_q (:,:) * 8.0_DP/fpi 239 CALL mp_sum( qrad_r , intra_bgrp_comm ) ; qrad_r(:,:,:) = qrad_r(:,:,:) * 8.0_DP/fpi 240 ! 241 ! 4) compute integrated power spectrum of the Qs in real space (completeness check). 242 ! 243 DO l = 0, upf(nt)%nqlc-1 244 ! 245 DO nb = 1, upf(nt)%nbeta 246 DO mb = nb, upf(nt)%nbeta 247 ! 248 ijv = mb * (mb - 1)/2 + nb 249 aux1(1) = 0.0_DP 250 ! 251 ir0 = 1 252 IF (rgrid(nt)%r(ir0) < eps8) ir0 = 2 253 ! 254 DO ir = ir0, upf(nt)%kkbeta 255 aux1(ir) = (qrad_r(ir,ijv,l+1)/rgrid(nt)%r(ir))**2 256 ENDDO 257 ! 258 CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, power_r(ijv,l+1) ) 259 ! 260 ENDDO 261 ENDDO 262 ! 263 ENDDO 264 ! 265 ! ... computes the ratio of power_0 and power_q to see how complete is the reciprocal 266 ! space expansion. 267 ! 268 power_0(:,0) = 0.0_DP 269 power_q(:,0) = 0.0_DP 270 power_r(:,0) = 0.0_DP 271 target_ratio = 1.0_DP 272 ! 273 DO nb = 1, upf(nt)%nbeta 274 DO mb = nb, upf(nt)%nbeta 275 ! 276 ijv = mb * (mb-1) / 2 + nb 277 lnb = upf(nt)%lll(nb) 278 lmb = upf(nt)%lll(mb) 279 ! 280 DO l = 0, lnb+lmb 281 power_0(ijv,0) = power_0(ijv,0) + power_0(ijv,l+1) 282 power_q(ijv,0) = power_q(ijv,0) + power_q(ijv,l+1) 283 power_r(ijv,0) = power_r(ijv,0) + power_r(ijv,l+1) 284 ENDDO 285 ! 286 !write (stdout, *) ' nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb 287 !write (stdout,'(a,12f16.10)') 'power_0 ',(power_0 (ijv,l+1), l=0,lnb+lmb) 288 !write (stdout,'(a,12f16.10)') 'power_r ',(power_r (ijv,l+1), l=0,lnb+lmb) 289 !write (stdout,'(a,12f16.10)') 'power_q ',(power_q (ijv,l+1), l=0,lnb+lmb) 290 !write (stdout,*) 'ratio ',1.d0-(power_r (ijv,0)/power_q (ijv,0)), & 291 ! 1.d0-(power_r (ijv,0)/power_0 (ijv,0)) 292 IF (power_0(ijv,0)>eps8) target_ratio = MIN(target_ratio, power_r(ijv,0)/power_0(ijv,0)) 293 ! 294 ENDDO 295 ENDDO 296 ! 297 WRITE (stdout,*) ' TARGET Qs SPILLOVER : 1.d0-target_ratio, target_ratio ', & 298 1.d0-target_ratio, target_ratio 299 ! 300 fac = 1.2_DP 301 ! 302 99 CONTINUE 303 ! 304 ! ... smooth the Fourier transform of the Qs and compute its integrated power spectrum 305 ! 306 power_qs(:,:) = 0.0_DP 307 ! 308 DO nb = 1, upf(nt)%nbeta 309 DO mb = nb, upf(nt)%nbeta 310 ! 311 ijv = mb * (mb - 1) / 2 + nb 312 lnb = upf(nt)%lll(nb) 313 lmb = upf(nt)%lll(mb) 314 ! 315 DO l = 0, upf(nt)%nqlc-1 316 ! 317 IF ( (l >= ABS(lnb-lmb)) .AND. (l <= lnb+lmb) .AND. (MOD(l+lnb+lmb,2)==0) ) THEN 318 DO iq = startq, lastq 319 q = (iq - 1) * dq 320 power_qs(ijv,l+1) = power_qs(ijv,l+1) + q*q * dq * & 321 (qrad_q(iq,ijv,l+1)*filter(fac*q/qmax,a,nn))**2 322 ENDDO 323 ENDIF 324 ! 325 power_qs(ijv,0) = power_qs(ijv,0) + power_qs(ijv,l+1) 326 ! 327 ENDDO 328 ! 329 ENDDO 330 ENDDO 331 ! 332 power_qs(:,:) = power_qs(:,:) * 8.0_DP/fpi 333 ! 334 CALL mp_sum( power_qs, intra_bgrp_comm ) 335 ! 336 ratio = 1.0_DP ; ratio_s = 1.0_DP 337 ! 338 DO nb = 1, upf(nt)%nbeta 339 DO mb = nb, upf(nt)%nbeta 340 ijv = mb * (mb - 1) / 2 + nb 341 IF (power_q(ijv,0)>eps8) ratio = MIN(ratio, power_qs(ijv,0)/power_q(ijv,0)) 342 ENDDO 343 ENDDO 344 ! 345 ! WRITE (stdout,*) ' filter factor and power_qs/power_q ratio:', fac, ratio 346 ! 347 ! ... with a given fac the smoothed Qs qre built 348 ! 349 qrad_rs(:,:,:) = 0.0_DP ; ffrr(:) = 0.0_DP 350 ! 351 DO l = 0, upf(nt)%nqlc-1 352 ! 353 DO iq = startq, lastq 354 ! 355 q = (iq - 1) * dq 356 ! 357 ! ... here we compute the spherical bessel function for the given |q| ... 358 ! 359 CALL sph_bes( upf(nt)%kkbeta, rgrid(nt)%r, q, l, aux ) 360 ! 361 ! ... and here we integrate with all the Q functions 362 ! 363 ! 5) back-fourier transform of the smoothed Qs to real space. 364 ! 365 DO ir = 1, upf(nt)%kkbeta 366 ! 367 ! ... q_nb_mb_l(r) from the back fourier transform up to qmax of q_nb_mb_l(q) 368 ! 369 qrad_rs(ir,1:nbetam*(nbetam+1)/2,l+1) = qrad_rs(ir,1:nbetam*(nbetam+1)/2,l+1) & 370 + aux(ir) * q*q * dq * rgrid(nt)%r(ir)**2 & 371 * qrad_q(iq,1:nbetam*(nbetam+1)/2,l+1) * filter(fac*q/qmax,a,nn) 372 ! 373 ! ... build the filter function in real space from the back fourier transform up 374 ! to qmax 375 ! 376 IF (l==0) ffrr(ir) = ffrr(ir) + q*q * dq * aux(ir) * rgrid(nt)%r(ir)**2 & 377 * filter(fac*q/qmax,a,nn) 378 ! 379 ENDDO 380 ! 381 ENDDO 382 ! 383 ENDDO 384 ! 385 CALL mp_sum( qrad_rs, intra_bgrp_comm ) ; qrad_rs(:,:,:) = qrad_rs(:,:,:) * 8.0_DP/fpi 386 CALL mp_sum( ffrr, intra_bgrp_comm ) ; ffrr(:) = ffrr(:) * 8.0_DP/fpi 387 ! 388 ! 6) compute intergrated power spectrum of the Qs in real space (completeness check). 389 ! 390 DO l = 0, upf(nt)%nqlc-1 391 ! 392 DO nb = 1, upf(nt)%nbeta 393 DO mb = nb, upf(nt)%nbeta 394 ! 395 ijv = mb * (mb-1)/2 + nb 396 aux1(1) = 0.0_DP 397 ! 398 ir0 = 1 399 IF (rgrid(nt)%r(ir0) < eps8) ir0 = 2 400 ! 401 DO ir = ir0, upf(nt)%kkbeta 402 aux1(ir) = (qrad_rs(ir,ijv,l+1)/rgrid(nt)%r(ir))**2 403 ENDDO 404 ! 405 CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, power_rs(ijv,l+1) ) 406 ! 407 ENDDO 408 ENDDO 409 ! 410 ENDDO 411 ! 412 power_rs(:,0) = 0.0_DP 413 ! 414 DO nb = 1, upf(nt)%nbeta 415 DO mb = nb, upf(nt)%nbeta 416 ! 417 ijv = mb * (mb-1)/2 + nb 418 lnb = upf(nt)%lll(nb) 419 lmb = upf(nt)%lll(mb) 420 ! 421 DO l = 0, lnb+lmb 422 power_rs(ijv,0) = power_rs(ijv,0) + power_rs(ijv,l+1) 423 ENDDO 424 ! 425 !write (stdout, *) ' nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb 426 !write (stdout,'(a,12f16.10)') 'power_rs',(power_rs(ijv,l+1), l=0,lnb+lmb) 427 !write (stdout,'(a,12f16.10)') 'power_qs',(power_qs(ijv,l+1), l=0,lnb+lmb) 428 !write (stdout,*) 'ratio ',1.d0-(power_rs(ijv,0)/power_qs(ijv,0)), & 429 ! 1.d0-(power_qs(ijv,0)/power_q (ijv,0)) 430 IF (power_qs(ijv,0) > eps8) ratio_s = MIN(ratio_s, power_rs(ijv,0)/power_qs(ijv,0)) 431 ! 432 ENDDO 433 ENDDO 434 ! 435 WRITE (stdout,*) ' filter factor, 1.d0-power_rs/power_qs and power_qs/power_q ratio :', & 436 fac, 1.d0-ratio_s, ratio 437 fac = fac - 0.05_DP 438 !sdg 439 IF (ratio < target_ratio .AND. (1.d0-ratio_s) < eps8) GOTO 99 440 ! 441 ! IF (ratio < target_ratio .AND. (1.d0-ratio_s) < eps6) GOTO 99 442 ! 443 fac = fac + 0.05_DP ! reset the last successful value of fac 444 ! 445 !- save the smoothed real space qs in qfuncl 446 ! 447 DO nb = 1, upf(nt)%nbeta 448 DO mb = nb, upf(nt)%nbeta 449 ! 450 ijv = mb * (mb-1)/2 + nb 451 lnb = upf(nt)%lll(nb) 452 lmb = upf(nt)%lll(mb) 453 ! 454 DO l = 0, upf(nt)%nqlc-1 455 ! 456 IF ( (l>=ABS(lnb-lmb)) .AND. (l<=lnb+lmb) .AND. (MOD(l+lnb+lmb,2)==0) ) THEN 457 ! 458 upf(nt)%qfuncl(1:upf(nt)%kkbeta,ijv,l) = qrad_rs(1:upf(nt)%kkbeta,ijv,l+1) 459 ! 460 ENDIF 461 ! 462 ENDDO 463 ! 464 ENDDO 465 ENDDO 466 ! 467 ENDIF 468 ! 469 IF (tprint) THEN 470 !- 471 filename = 'ffff' 472 OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown') 473 WRITE (4, *) '# filter function : a=',a,', nn=',nn,', fac=', fac 474 WRITE (4, *) '# nqxq :', nqxq,' dq :',dq, ' qmax :',qmax 475 DO iq = 1, nqxq 476 q = (iq-1)*dq 477 WRITE (4,'(2f16.10)') q, filter( fac*q/qmax, a, nn ) 478 ENDDO 479 CLOSE (4) 480 !- 481 filename = 'ffrr' 482 OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown') 483 WRITE (4, *) '# filter function : a=',a,', nn=',nn,', fac=', fac 484 WRITE (4, *) '# kkbeta :', upf(nt)%kkbeta 485 DO ir = 1, upf(nt)%kkbeta 486 write (4,'(2f16.10)') rgrid(nt)%r(ir), ffrr(ir) 487 ENDDO 488 CLOSE (4) 489 !- 490 DO nb = 1, upf(nt)%nbeta 491 DO mb = nb, upf(nt)%nbeta 492 ! 493 ijv = mb * (mb - 1) / 2 + nb 494 lnb = upf(nt)%lll(nb) ; lmb = upf(nt)%lll(mb) 495 WRITE (filename(4:4),'(i1)') nb 496 WRITE (filename(5:5),'(i1)') mb 497 WRITE (filename(6:6),'(i1)') nt 498 !- 499 filename(1:3) = 'qq_' ! the radial fourier transform of q_l in reciprocal space 500 OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown') 501 WRITE (4,*) '# the radial fourier transform of q_l in reciprocal space' 502 WRITE (4,*) '# nb :', nb, lnb,' mb :', mb, lmb,' lmax :', lnb+lmb, ' nqxq :', nqxq 503 DO iq=1,nqxq 504 q = (iq-1)*dq 505 WRITE (4,'(12f16.10)') q, (qrad_q(iq,ijv,l+1), l=0,lnb+lmb ) 506 ENDDO 507 CLOSE (4) 508 !- 509 filename(1:3) = 'qqs' ! the smoothed radial fourier transform of q_l in reciprocal space 510 OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown') 511 WRITE (4,*) '# the smoothed radial fourier transform of q_l in reciprocal space' 512 WRITE (4,*) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, ' nqxq :',nqxq 513 DO iq = 1, nqxq 514 q = (iq-1)*dq 515 WRITE (4,'(12f16.10)') q,(qrad_q(iq,ijv,l+1)*filter(fac*q/qmax,a,nn), l=0,lnb+lmb ) 516 ENDDO 517 CLOSE (4) 518 !- 519 filename(1:3) = 'qrq' ! the radial q_l(r) as obtained back-transforming the q_l(q) 520 OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown') 521 WRITE (4,*) '# the radial q_l(r) as obtained back-transforming the q_l(q)' 522 WRITE (4,*) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, ' kkbeta :',upf(nt)%kkbeta 523 DO ir = 1, upf(nt)%kkbeta 524 WRITE (4,'(12f16.10)') rgrid(nt)%r(ir),(qrad_r(ir,ijv,l+1), l=0,lnb+lmb ) 525 ENDDO 526 CLOSE (4) 527 !- 528 filename(1:3) = 'qrs' ! the radial q_l(r) as obtained back-transforming the smoothed q_l(q) 529 OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown') 530 WRITE (4,*) '# the radial q_l(r) as obtained back-transforming the smoothed q_l(q)' 531 WRITE (4,*) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, ' kkbeta :',upf(nt)%kkbeta 532 DO ir = 1, upf(nt)%kkbeta 533 WRITE (4,'(12f16.10)') rgrid(nt)%r(ir),(qrad_rs(ir,ijv,l+1), l=0,lnb+lmb ) 534 ENDDO 535 CLOSE (4) 536 !- 537 ENDDO 538 ENDDO 539 ! 540 ENDIF 541 ! ntyp 542 ENDDO 543 ! 544 DEALLOCATE( ylmk0 ) 545 DEALLOCATE( power_0, power_r, power_q, power_rs, power_qs ) 546 DEALLOCATE( qrad_q, qrad_r, qrad_rs, ffrr ) 547 DEALLOCATE( aux1, aux ) 548 ! 549 CALL stop_clock( 'init_us_0' ) 550 ! 551 RETURN 552 ! 553 ! 554 ! 555 CONTAINS 556 ! 557 REAL(DP) FUNCTION filter( x, a, n ) 558 ! 559 IMPLICIT NONE 560 ! 561 REAL(DP), INTENT(IN) :: x, a 562 INTEGER, INTENT(IN) :: n 563 REAL(DP) :: axx, ff 564 INTEGER :: k 565 ! 566 axx = a * x * x 567 ff = 1.0_DP 568 ! 569 DO k = n, 1, -1 570 ff = (1._DP + axx/DBLE(k)*ff) 571 ENDDO 572 ! 573 filter = ff * EXP(-axx) 574 ! 575 RETURN 576 ! 577 END FUNCTION filter 578 ! 579 ! 580END SUBROUTINE init_us_0 581