1! 2! Copyright (C) 2002-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 10 11 SUBROUTINE compute_dvan_x() 12 ! 13 ! calculate array dvan(iv,jv,is) 14 ! 15 ! rw**2 * vrps = [ ( Vpsnl(r) - Vpsloc(r) )* Rps(r) * r^2 ] 16 ! = [ DVpsnl(r) * Rps(r) * r^2 ] 17 ! dion = (2l+1) / < Rps(r) | DVpsnl(r) | Rps(r) > 18 19 USE kinds, ONLY: DP 20 use uspp, only: dvan, nhtolm, indv 21 use uspp_param, only: upf, nhm, nh 22 use ions_base, only: nsp 23 ! 24 implicit none 25 ! 26 integer :: is, iv, jv 27 real(DP) :: fac 28 ! 29 if( allocated( dvan ) ) deallocate( dvan ) 30 allocate( dvan( nhm, nhm, nsp ) ) 31 dvan(:,:,:) =0.d0 32 ! 33 do is = 1, nsp 34 ! fac converts ry to hartree 35 fac = 0.5d0 36 do iv=1,nh(is) 37 do jv=1,nh(is) 38 if ( nhtolm(iv,is) == nhtolm(jv,is) ) then 39 dvan( iv, jv, is ) = fac * upf(is)%dion( indv(iv,is), indv(jv,is) ) 40 endif 41 end do 42 end do 43 end do 44 RETURN 45 END SUBROUTINE compute_dvan_x 46 47 48 49!------------------------------------------------------------------------------! 50 51 52 53 SUBROUTINE pseudopotential_indexes_x( ) 54 55 use upf_params, only: lmaxx ! 56 use ions_base, only: nsp, & ! number of specie 57 na, & ! number of atoms for each specie 58 nat, & ! total number of atom 59 ityp ! the atomi specie for each atom 60 use uspp, only: nkb, & ! 61 nkbus ! 62 use uspp_param, only: ish, &! 63 upf, &! 64 lmaxkb, &! 65 nhm, &! 66 nbetam, &! 67 nh, &! 68 lmaxq ! 69 use uspp, only: nhtol, &! 70 nhtolm, &! 71 indv, &! 72 ijtoh, &! 73 indv_ijkb0 ! 74 75 IMPLICIT NONE 76 77 ! 78 INTEGER :: is, iv, ind, il, lm, ih, jh, ijv, ijkb0, ia 79 ! ------------------------------------------------------------------ 80 ! find number of beta functions per species, max dimensions, 81 ! total number of beta functions (all and Vanderbilt only) 82 ! ------------------------------------------------------------------ 83 lmaxkb = -1 84 nkb = 0 85 nkbus = 0 86 ! 87 do is = 1, nsp 88 ind = 0 89 do iv = 1, upf(is)%nbeta 90 lmaxkb = max( lmaxkb, upf(is)%lll( iv ) ) 91 ind = ind + 2 * upf(is)%lll( iv ) + 1 92 end do 93 nh(is) = ind 94 ish(is)=nkb 95 nkb = nkb + na(is) * nh(is) 96 if( upf(is)%tvanp ) nkbus = nkbus + na(is) * nh(is) 97 end do 98 nhm = MAXVAL( nh(1:nsp) ) 99 nbetam = MAXVAL( upf(1:nsp)%nbeta ) 100 if (lmaxkb > lmaxx) call errore(' pseudopotential_indexes ',' l > lmax ',lmaxkb) 101 lmaxq = 2*lmaxkb + 1 102 ! 103 ! the following prevents an out-of-bound error: nqlc(is)=2*lmax+1 104 ! but in some versions of the PP files lmax is not set to the maximum 105 ! l of the beta functions but includes the l of the local potential 106 ! 107 do is=1,nsp 108 upf(is)%nqlc = MIN ( upf(is)%nqlc, lmaxq ) 109 IF ( upf(is)%nqlc < 0 ) upf(is)%nqlc = 0 110 end do 111 if (nkb <= 0) call errore(' pseudopotential_indexes ',' not implemented ?',nkb) 112 113 if( allocated( nhtol ) ) deallocate( nhtol ) 114 if( allocated( indv ) ) deallocate( indv ) 115 if( allocated( nhtolm ) ) deallocate( nhtolm ) 116 if( allocated( ijtoh ) ) deallocate( ijtoh ) 117 if( allocated( indv_ijkb0 ) ) deallocate( indv_ijkb0 ) 118 ! 119 allocate(nhtol(nhm,nsp)) 120 allocate(indv (nhm,nsp)) 121 allocate(nhtolm(nhm,nsp)) 122 allocate(ijtoh(nhm,nhm,nsp)) 123 allocate(indv_ijkb0(nat)) 124 125 ! ------------------------------------------------------------------ 126 ! definition of indices nhtol, indv, nhtolm 127 ! ------------------------------------------------------------------ 128 ! 129 ijkb0 = 0 130 do is = 1, nsp 131 ind = 0 132 do iv = 1, upf(is)%nbeta 133 lm = upf(is)%lll(iv)**2 134 do il = 1, 2* upf(is)%lll( iv ) + 1 135 lm = lm + 1 136 ind = ind + 1 137 nhtolm( ind, is ) = lm 138 nhtol( ind, is ) = upf(is)%lll( iv ) 139 indv( ind, is ) = iv 140 end do 141 end do 142 ! 143 ! ijtoh map augmentation channel indexes ih and jh to composite 144 ! "triangular" index ijh 145 ijtoh(:,:,is) = -1 146 ijv = 0 147 do ih = 1,nh(is) 148 do jh = ih,nh(is) 149 ijv = ijv+1 150 ijtoh(ih,jh,is) = ijv 151 ijtoh(jh,ih,is) = ijv 152 end do 153 end do 154 ! 155 ! ijkb0 is just before the first beta "in the solid" for atom ia 156 ! i.e. ijkb0+1,.. ijkb0+nh(ityp(ia)) are the nh beta functions of 157 ! atom ia in the global list of beta functions 158 do ia = 1,nat 159 IF ( ityp(ia) == is ) THEN 160 indv_ijkb0(ia) = ijkb0 161 ijkb0 = ijkb0 + nh(is) 162 END IF 163 end do 164 165 end do 166 167 RETURN 168 END SUBROUTINE pseudopotential_indexes_x 169 170 171!------------------------------------------------------------------------------! 172 173 174 SUBROUTINE compute_xgtab_x( xgmin, xgmax ) 175 ! 176 USE kinds, ONLY : DP 177 USE cell_base, ONLY : tpiba, tpiba2 178 USE mp, ONLY : mp_max 179 USE mp_global, ONLY : intra_bgrp_comm 180 USE gvect, ONLY : gg 181 USE pseudopotential, ONLY : xgtab 182 USE betax, ONLY : mmx, refg 183 ! 184 IMPLICIT NONE 185 ! 186 REAL(DP), INTENT(OUT) :: xgmax, xgmin 187 ! 188 INTEGER :: ig 189 REAL(DP) :: xg, dxg, res 190 ! 191 IF( ALLOCATED( xgtab ) ) & 192 DEALLOCATE( xgtab ) 193 194 ALLOCATE( xgtab( mmx ) ) 195 ! 196 xgmin = 0.0d0 197 xgmax = SQRT( refg * mmx ) 198 dxg = (xgmax - xgmin) / DBLE( mmx - 1 ) 199 ! 200 DO ig = 1, SIZE( xgtab ) 201 xgtab(ig) = xgmin + DBLE(ig-1) * dxg 202 END DO 203 ! 204 xgtab = xgtab**2 / tpiba**2 205 ! 206 RETURN 207 END SUBROUTINE compute_xgtab_x 208 209 210!------------------------------------------------------------------------------! 211 212 213 SUBROUTINE build_pstab_x( ) 214 215 USE kinds, ONLY : DP 216 USE atom, ONLY : rgrid 217 USE ions_base, ONLY : nsp, rcmax, zv 218 USE cell_base, ONLY : tpiba, tpiba2 219 USE splines, ONLY : init_spline, allocate_spline, kill_spline, nullify_spline 220 USE pseudo_base, ONLY : formfn, formfa 221 USE uspp_param, only : upf 222 USE control_flags, only : tpre 223 use gvect, ONLY : gg, gstart 224 USE cp_interfaces, ONLY : compute_xgtab 225 USE pseudopotential, ONLY : vps_sp, dvps_sp, xgtab 226 USE local_pseudo, ONLY : vps0 227 USE betax, ONLY : mmx 228 229 IMPLICIT NONE 230 231 INTEGER :: is, ig 232 REAL(DP) :: xgmax, xgmin 233 LOGICAL :: compute_tab 234 ! 235 IF( .NOT. ALLOCATED( rgrid ) ) & 236 CALL errore( ' build_pstab_x ', ' rgrid not allocated ', 1 ) 237 IF( .NOT. ALLOCATED( upf ) ) & 238 CALL errore( ' build_pstab_x ', ' upf not allocated ', 1 ) 239 ! 240 IF( ALLOCATED( vps_sp ) .AND. ALLOCATED( dvps_sp ) ) THEN 241 ! 242 DO is = 1, nsp 243 CALL kill_spline( vps_sp(is), 'a' ) 244 CALL kill_spline(dvps_sp(is),'a') 245 END DO 246 DEALLOCATE( vps_sp ) 247 DEALLOCATE(dvps_sp) 248 ! 249 END IF 250 ! 251 IF( ALLOCATED( vps_sp ) .OR. ALLOCATED( dvps_sp ) ) THEN 252 CALL errore( ' build_pstab_x ', ' inconsistent allocation ', 1 ) 253 END IF 254 ! 255 CALL compute_xgtab( xgmin, xgmax ) 256 ! 257 ALLOCATE( vps_sp(nsp)) 258 ALLOCATE( dvps_sp(nsp)) 259 ! 260 DO is = 1, nsp 261 262 CALL nullify_spline( vps_sp( is ) ) 263 CALL nullify_spline( dvps_sp( is ) ) 264 265 CALL allocate_spline( vps_sp(is), mmx, xgmin, xgmax ) 266 CALL allocate_spline( dvps_sp(is), mmx, xgmin, xgmax ) 267 268 call formfn( rgrid(is)%r, rgrid(is)%rab, & 269 upf(is)%vloc(1:rgrid(is)%mesh), zv(is), rcmax(is), & 270 xgtab, 1.0d0, tpiba2, rgrid(is)%mesh, mmx, & 271 tpre, vps_sp(is)%y, vps0(is), dvps_sp(is)%y ) 272 ! obsolete BHS format 273 !call formfa( vps_sp(is)%y, dvps_sp(is)%y, rc1(is), rc2(is), & 274 ! wrc1(is), wrc2(is), rcl(:,is,lloc(is)), & 275 ! al(:,is,lloc(is)), bl(:,is,lloc(is)), zv(is), & 276 ! rcmax(is), xgtab, 1.0d0, tpiba2, mmx, 2 , tpre ) 277 278 ! WRITE( 13, "(3D16.8)" ) ( xgtab(ig), vps_sp(is)%y(ig), dvps_sp(is)%y(ig), ig = 1, mmx ) 279 280 CALL init_spline( vps_sp(is) ) 281 CALL init_spline( dvps_sp(is) ) 282 283 END DO 284 285 RETURN 286 END SUBROUTINE build_pstab_x 287 288 289!------------------------------------------------------------------------------! 290 291 292 SUBROUTINE build_cctab_x( ) 293 294 USE kinds, ONLY : DP 295 USE atom, ONLY : rgrid 296 USE uspp_param, ONLY : upf 297 USE ions_base, ONLY : nsp, rcmax 298 USE cell_base, ONLY : tpiba, tpiba2 299 USE splines, ONLY : init_spline, allocate_spline, kill_spline, nullify_spline 300 USE pseudo_base, ONLY : compute_rhocg 301 USE gvect, ONLY : gg, gstart 302 USE cp_interfaces, ONLY : compute_xgtab 303 USE pseudopotential, ONLY : rhoc1_sp, rhocp_sp, xgtab 304 USE betax, ONLY : mmx 305 306 IMPLICIT NONE 307 308 INTEGER :: is 309 REAL(DP) :: xgmax, xgmin 310 LOGICAL :: compute_tab 311 ! 312 IF( .NOT. ALLOCATED( rgrid ) ) & 313 CALL errore( ' build_cctab_x ', ' rgrid not allocated ', 1 ) 314 IF( .NOT. ALLOCATED( upf ) ) & 315 CALL errore( ' build_cctab_x ', ' upf not allocated ', 1 ) 316 ! 317 IF( ALLOCATED( rhoc1_sp ) .AND. ALLOCATED( rhocp_sp ) ) THEN 318 ! 319 DO is = 1, nsp 320 CALL kill_spline(rhoc1_sp(is),'a') 321 CALL kill_spline(rhocp_sp(is),'a') 322 END DO 323 DEALLOCATE(rhoc1_sp) 324 DEALLOCATE(rhocp_sp) 325 ! 326 END IF 327 ! 328 IF( ALLOCATED( rhoc1_sp ) .OR. ALLOCATED( rhocp_sp ) ) THEN 329 CALL errore( ' build_cctab_x ', ' inconsistent allocation ', 1 ) 330 END IF 331 ! 332 CALL compute_xgtab( xgmin, xgmax ) 333 ! 334 ALLOCATE( rhoc1_sp(nsp)) 335 ALLOCATE( rhocp_sp(nsp)) 336 ! 337 DO is = 1, nsp 338 339 CALL nullify_spline( rhoc1_sp( is ) ) 340 CALL nullify_spline( rhocp_sp( is ) ) 341 342 IF( upf(is)%nlcc ) THEN 343 ! 344 CALL allocate_spline( rhoc1_sp(is), mmx, xgmin, xgmax ) 345 CALL allocate_spline( rhocp_sp(is), mmx, xgmin, xgmax ) 346 ! 347 CALL compute_rhocg( rhoc1_sp(is)%y, rhocp_sp(is)%y, rgrid(is)%r, & 348 rgrid(is)%rab, upf(is)%rho_atc(:), xgtab, 1.0d0, tpiba2, & 349 rgrid(is)%mesh, mmx, 1 ) 350 ! 351 CALL init_spline( rhoc1_sp(is) ) 352 CALL init_spline( rhocp_sp(is) ) 353 ! 354 END IF 355 356 END DO 357 358 RETURN 359 END SUBROUTINE build_cctab_x 360 361 362!------------------------------------------------------------------------------! 363 364 365 SUBROUTINE compute_betagx_x( tpre ) 366 ! 367 ! calculation of array betagx(ig,iv,is) 368 ! 369 USE kinds, ONLY : DP 370 USE ions_base, ONLY : nsp 371 USE uspp_param, ONLY : upf, nh, nhm 372 USE atom, ONLY : rgrid 373 USE uspp, ONLY : nhtol, indv 374 USE betax, only : refg, betagx, mmx, dbetagx 375 ! 376 IMPLICIT NONE 377 ! 378 LOGICAL, INTENT(IN) :: tpre 379 ! 380 INTEGER :: is, iv, l, il, ir, nr 381 REAL(DP), ALLOCATABLE :: dfint(:), djl(:), fint(:), jl(:) 382 REAL(DP) :: xg 383 ! 384 CALL start_clock('betagx') 385 ! 386 IF( .NOT. ALLOCATED( rgrid ) ) & 387 CALL errore( ' compute_betagx_x ', ' rgrid not allocated ', 1 ) 388 IF( .NOT. ALLOCATED( upf ) ) & 389 CALL errore( ' compute_betagx_x ', ' upf not allocated ', 1 ) 390 ! 391 IF( ALLOCATED( betagx ) ) DEALLOCATE( betagx ) 392 IF( ALLOCATED( dbetagx ) ) DEALLOCATE( dbetagx ) 393 ! 394 ALLOCATE( betagx ( mmx, nhm, nsp ) ) 395 IF ( tpre ) ALLOCATE( dbetagx( mmx, nhm, nsp ) ) 396 ! 397 do is = 1, nsp 398 ! 399 nr = upf(is)%kkbeta 400 ! 401 do iv = 1, nh(is) 402 ! 403 l = nhtol(iv,is) 404 ! 405!$omp parallel default(none), private( dfint, djl, fint, jl, il, xg, ir ), & 406!$omp shared( tpre, nr, mmx, refg, l, is, rgrid, upf, indv, iv, betagx, dbetagx ) 407 if ( tpre ) then 408 allocate( dfint( nr ) ) 409 allocate( djl ( nr ) ) 410 end if 411 ! 412 allocate( fint ( nr ) ) 413 allocate( jl ( nr ) ) 414 ! 415!$omp do 416 interp_tab : do il = 1, mmx 417 ! 418 xg = sqrt( refg * (il-1) ) 419 call sph_bes ( nr, rgrid(is)%r, xg, l, jl ) 420! 421 if( tpre )then 422 ! 423 call sph_dbes1 ( nr, rgrid(is)%r, xg, l, jl, djl) 424 ! 425 endif 426 ! 427 ! beta(ir)=r*beta(r) 428 ! 429 do ir = 1, nr 430 fint(ir) = rgrid(is)%r(ir) * jl(ir) * & 431 upf(is)%beta( ir, indv(iv,is) ) 432 end do 433 call simpson_cp90(nr,fint,rgrid(is)%rab,betagx(il,iv,is)) 434 ! 435 if(tpre) then 436 do ir = 1, nr 437 dfint(ir) = rgrid(is)%r(ir) * djl(ir) * & 438 upf(is)%beta( ir, indv(iv,is) ) 439 end do 440 call simpson_cp90(nr,dfint,rgrid(is)%rab,dbetagx(il,iv,is)) 441 endif 442 ! 443 end do interp_tab 444!$omp end do 445 ! 446 deallocate(jl) 447 deallocate(fint) 448 ! 449 if (tpre) then 450 deallocate(djl) 451 deallocate(dfint) 452 end if 453 ! 454!$omp end parallel 455 ! 456 end do 457 ! 458 end do 459 CALL stop_clock('betagx') 460 RETURN 461 END SUBROUTINE compute_betagx_x 462 463 464!------------------------------------------------------------------------------! 465 466 467 SUBROUTINE compute_qradx_x( tpre ) 468 ! 469 ! calculation of array qradx(igb,iv,jv,is) for interpolation table 470 ! (symmetric wrt exchange of iv and jv: a single index ijv is used) 471 ! 472 ! qradx(ig,l,k,is) = 4pi/omega int_0^r dr r^2 j_l(qr) q(r,l,k,is) 473 ! 474 ! 475 ! 476 USE kinds, ONLY : DP 477 use io_global, only : stdout 478 USE ions_base, ONLY : nsp 479 USE uspp_param, ONLY : upf, nh, nhm, nbetam, lmaxq, ish 480 USE atom, ONLY : rgrid 481 USE uspp, ONLY : indv 482 USE betax, only : refg, qradx, mmx, dqradx 483 use smallbox_gvec, only : ngb 484 USE cp_interfaces, ONLY : fill_qrl 485 ! 486 IMPLICIT NONE 487 ! 488 LOGICAL, INTENT(IN) :: tpre 489 ! 490 INTEGER :: is, iv, l, il, ir, jv, ijv, ierr 491 INTEGER :: nr 492 REAL(DP), ALLOCATABLE :: dfint(:), djl(:), fint(:), jl(:), qrl(:,:,:) 493 REAL(DP) :: xg 494 495 496 CALL start_clock('qradx') 497 498 IF( .NOT. ALLOCATED( rgrid ) ) & 499 CALL errore( ' compute_qradx_x ', ' rgrid not allocated ', 1 ) 500 IF( .NOT. ALLOCATED( upf ) ) & 501 CALL errore( ' compute_qradx_x ', ' upf not allocated ', 1 ) 502 503 IF( ALLOCATED( qradx ) ) DEALLOCATE( qradx ) 504 IF( ALLOCATED( dqradx ) ) DEALLOCATE( dqradx ) 505 ! 506 ALLOCATE( qradx( mmx, nbetam*(nbetam+1)/2, lmaxq, nsp ) ) 507 ! 508 IF ( tpre ) ALLOCATE( dqradx( mmx, nbetam*(nbetam+1)/2, lmaxq, nsp ) ) 509 510 DO is = 1, nsp 511 ! 512 IF( .NOT. upf(is)%tvanp ) CYCLE 513 ! 514 ! qqq and beta are now indexed and taken in the same order 515 ! as vanderbilts ppot-code prints them out 516 ! 517 WRITE( stdout,'(A,5I5)') 'is, nh(is), ngb, kkbeta, lmaxq = ', & 518 & is, nh(is), ngb, upf(is)%kkbeta, upf(is)%nqlc 519 ! 520 nr = upf(is)%kkbeta 521 ! 522 ALLOCATE( qrl( nr, upf(is)%nbeta*(upf(is)%nbeta+1)/2, upf(is)%nqlc) ) 523 ! 524 call fill_qrl ( is, qrl ) 525 ! 526 do l = 1, upf(is)%nqlc 527 ! 528!$omp parallel default(none), private( djl, dfint, fint, jl, il, iv, jv, ijv, xg, ir ), & 529!$omp shared( tpre, nr, mmx, refg, rgrid, l, upf, qrl, qradx, dqradx, is ) 530 IF ( tpre ) THEN 531 ALLOCATE( djl ( nr ) ) 532 ALLOCATE( dfint( nr ) ) 533 END IF 534 ! 535 ALLOCATE( fint( nr ) ) 536 ALLOCATE( jl ( nr ) ) 537!$omp do 538 interp_tab : do il = 1, mmx 539 ! 540 xg = sqrt( refg * DBLE(il-1) ) 541 ! 542 call sph_bes ( nr, rgrid(is)%r, xg, l-1, jl(1) ) 543 ! 544 if( tpre ) then 545 ! 546 call sph_dbes1 ( nr, rgrid(is)%r, xg, l-1, jl, djl) 547 ! 548 endif 549 ! 550 ! 551 do iv = 1, upf(is)%nbeta 552 do jv = iv, upf(is)%nbeta 553 ijv = jv * ( jv - 1 ) / 2 + iv 554 ! 555 ! note qrl(r)=r^2*q(r) 556 ! 557 do ir = 1, nr 558 fint( ir ) = qrl( ir, ijv, l ) * jl( ir ) 559 end do 560 call simpson_cp90 & 561 (nr,fint(1),rgrid(is)%rab,qradx(il,ijv,l,is)) 562 ! 563 if( tpre ) then 564 do ir = 1, nr 565 dfint(ir) = qrl(ir,ijv,l) * djl(ir) 566 end do 567 call simpson_cp90 & 568 (nr,dfint(1),rgrid(is)%rab,dqradx(il,ijv,l,is)) 569 end if 570 ! 571 end do 572 end do 573 ! 574 ! 575 end do interp_tab 576!$omp end do 577 ! 578 DEALLOCATE ( jl ) 579 DEALLOCATE ( fint ) 580 ! 581 if ( tpre ) then 582 DEALLOCATE(djl) 583 DEALLOCATE ( dfint ) 584 end if 585 ! 586!$omp end parallel 587 ! 588 end do 589 ! 590 DEALLOCATE ( qrl ) 591 WRITE( stdout,*) 592 WRITE( stdout,'(20x,a)') ' qqq ' 593 ! 594 do iv=1,upf(is)%nbeta 595 WRITE( stdout,'(8f9.4)') (upf(is)%qqq(iv,jv),jv=1,upf(is)%nbeta) 596 end do 597 WRITE( stdout,*) 598 ! 599 end do 600 601 CALL stop_clock('qradx') 602 603 RETURN 604 END SUBROUTINE compute_qradx_x 605 606!------------------------------------------------------------------------------! 607 608 SUBROUTINE exact_qradb_x( tpre ) 609 ! 610 USE kinds, ONLY : DP 611 use io_global, only: stdout 612 USE ions_base, ONLY: nsp 613 USE uspp_param, ONLY: upf, nh, nhm, nbetam, lmaxq 614 use uspp_param, only: lmaxkb, ish 615 USE atom, ONLY: rgrid 616 USE uspp, ONLY: indv 617 use uspp, only: qq_nt, beta 618 USE betax, only: refg, qradx, mmx, dqradx 619 use smallbox_gvec, only: ngb 620 use control_flags, only: iprint, iverbosity 621 use cell_base, only: ainv 622 use constants, only: pi, fpi 623 use qgb_mod, only: qgb, dqgb 624 use smallbox_gvec, only: gb, gxb 625 use small_box, only: omegab, tpibab 626 USE cp_interfaces, ONLY: fill_qrl 627 ! 628 IMPLICIT NONE 629 ! 630 LOGICAL, INTENT(IN) :: tpre 631 ! 632 INTEGER :: is, iv, l, il, ir, jv, ijv, ierr 633 INTEGER :: ig, i,j, jj, nr 634 REAL(DP), ALLOCATABLE :: dfint(:), djl(:), fint(:), jl(:), qrl(:,:,:) 635 REAL(DP) :: xg, c, betagl, dbetagl 636 REAL(DP), ALLOCATABLE :: dqradb(:,:,:,:) 637 REAL(DP), ALLOCATABLE :: dqrad( :, :, :, :, :, : ) 638 REAL(DP), ALLOCATABLE :: qradb(:,:,:,:) 639 REAL(DP), ALLOCATABLE :: ylmb(:,:), dylmb(:,:,:,:) 640 COMPLEX(DP), ALLOCATABLE :: dqgbs(:,:,:) 641 LOGICAL :: tvanp 642 643 tvanp = .FALSE. 644 DO is = 1, nsp 645 tvanp = tvanp .OR. upf(is)%tvanp 646 END DO 647 IF( .NOT. tvanp ) & 648 return 649 650 IF( .NOT. ALLOCATED( rgrid ) ) & 651 CALL errore( ' exact_qradb_x ', ' rgrid not allocated ', 1 ) 652 IF( .NOT. ALLOCATED( upf ) ) & 653 CALL errore( ' exact_qradb_x ', ' upf not allocated ', 1 ) 654 655 IF( ALLOCATED( qradx ) ) DEALLOCATE( qradx ) 656 IF( ALLOCATED( dqradx ) ) DEALLOCATE( dqradx ) 657 ! 658 ALLOCATE( qradx( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp ) ) 659 ! 660 IF ( tpre ) ALLOCATE( dqradx( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp ) ) 661 662 ALLOCATE( qradb( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp ) ) 663 qradb(:,:,:,:) = 0.d0 664 665 666 DO is = 1, nsp 667 668 IF( .NOT. upf(is)%tvanp ) CYCLE 669 ! 670 ! qqq and beta are now indexed and taken in the same order 671 ! as vanderbilts ppot-code prints them out 672 ! 673 WRITE( stdout,*) ' nlinit nh(is), ngb, is, kkbeta, lmaxq = ', & 674 & nh(is), ngb, is, upf(is)%kkbeta, upf(is)%nqlc 675 ! 676 nr = upf(is)%kkbeta 677 ! 678 IF ( tpre ) THEN 679 ALLOCATE( djl ( nr ) ) 680 ALLOCATE( dfint( nr ) ) 681 END IF 682 ! 683 ALLOCATE( fint( nr ) ) 684 ALLOCATE( jl ( nr ) ) 685 ALLOCATE( qrl( nr, upf(is)%nbeta*(upf(is)%nbeta+1)/2, upf(is)%nqlc) ) 686 ! 687 call fill_qrl ( is, qrl ) 688 ! qrl = 0.0d0 689 ! 690 do l = 1, upf(is)%nqlc 691 ! 692 do il = 1, ngb 693 ! 694 xg = sqrt( gb( il ) * tpibab * tpibab ) 695 ! 696 call sph_bes ( nr, rgrid(is)%r, xg, l-1, jl(1) ) 697 ! 698 if( tpre ) then 699 ! 700 call sph_dbes1 ( nr, rgrid(is)%r, xg, l-1, jl, djl) 701 ! 702 endif 703 ! 704 ! 705 do iv = 1, upf(is)%nbeta 706 do jv = iv, upf(is)%nbeta 707 ijv = jv * ( jv - 1 ) / 2 + iv 708 ! 709 ! note qrl(r)=r^2*q(r) 710 ! 711 do ir = 1, nr 712 fint( ir ) = qrl( ir, ijv, l ) * jl( ir ) 713 end do 714 call simpson_cp90 & 715 (nr,fint(1),rgrid(is)%rab,qradx(il,ijv,l,is)) 716 ! 717 if( tpre ) then 718 do ir = 1, nr 719 dfint(ir) = qrl(ir,ijv,l) * djl(ir) 720 end do 721 call simpson_cp90 & 722 (nr,dfint(1),rgrid(is)%rab,dqradx(il,ijv,l,is)) 723 end if 724 ! 725 end do 726 end do 727 ! 728 ! 729 end do 730 end do 731 ! 732 DEALLOCATE ( jl ) 733 DEALLOCATE ( qrl ) 734 DEALLOCATE ( fint ) 735 ! 736 if ( tpre ) then 737 DEALLOCATE(djl) 738 DEALLOCATE ( dfint ) 739 end if 740 ! 741 WRITE( stdout,*) 742 WRITE( stdout,'(20x,a)') ' qqq ' 743 ! 744 do iv=1, upf(is)%nbeta 745 WRITE( stdout,'(8f9.4)') (upf(is)%qqq(iv,jv),jv=1, upf(is)%nbeta) 746 end do 747 WRITE( stdout,*) 748 ! 749 end do 750 751 allocate( ylmb( ngb, lmaxq*lmaxq ), STAT=ierr ) 752 IF( ierr /= 0 ) & 753 CALL errore(' exact_qradb ', ' cannot allocate ylmb ', 1 ) 754! 755 call ylmr2 (lmaxq*lmaxq, ngb, gxb, gb, ylmb) 756 757 do is = 1, nsp 758 759 IF( .NOT. upf(is)%tvanp ) CYCLE 760 ! 761 ! calculation of array qradb(igb,iv,jv,is) 762 ! 763 ! qradb(ig,l,k,is) = 4pi/omega int_0^r dr r^2 j_l(qr) q(r,l,k,is) 764 ! 765 if( iverbosity > 2 ) WRITE( stdout,*) ' qradb ' 766 ! 767 c = fpi / omegab 768 ! 769 do iv= 1, upf(is)%nbeta 770 do jv = iv, upf(is)%nbeta 771 ijv = jv*(jv-1)/2 + iv 772 do ig=1,ngb 773 do l=1,upf(is)%nqlc 774 qradb(ig,ijv,l,is)= c*qradx(ig,ijv,l,is) 775 enddo 776 enddo 777 enddo 778 enddo 779 ! 780 ! --------------------------------------------------------------- 781 ! stocking of qgb(igb,ijv,is) and of qq(iv,jv,is) 782 ! --------------------------------------------------------------- 783 ! 784 do iv= 1,nh(is) 785 do jv=iv,nh(is) 786 ! 787 ! compact indices because qgb is symmetric 788 ! 789 ijv = jv*(jv-1)/2 + iv 790 call qvan2b(ngb,iv,jv,is,ylmb,qgb(1,ijv,is),qradb ) 791! 792 qq_nt(iv,jv,is)=omegab*DBLE(qgb(1,ijv,is)) 793 qq_nt(jv,iv,is)=qq_nt(iv,jv,is) 794! 795 end do 796 end do 797 798 end do 799! 800 if (tpre) then 801! --------------------------------------------------------------- 802! arrays required for stress calculation, variable-cell dynamics 803! --------------------------------------------------------------- 804 allocate(dqradb(ngb,nbetam*(nbetam+1)/2,lmaxq,nsp)) 805 allocate(dylmb(ngb,lmaxq*lmaxq,3,3)) 806 allocate(dqgbs(ngb,3,3)) 807 allocate( dqrad( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp, 3, 3 ) ) 808 dqrad(:,:,:,:,:,:) = 0.d0 809 ! 810 call dylmr2_(lmaxq*lmaxq, ngb, gxb, gb, ainv, dylmb) 811 ! 812 do is=1,nsp 813 ! 814 IF( .NOT. upf(is)%tvanp ) CYCLE 815 ! 816 do iv= 1, upf(is)%nbeta 817 do jv=iv, upf(is)%nbeta 818 ijv = jv*(jv-1)/2 + iv 819 do l=1,upf(is)%nqlc 820 do ig=1,ngb 821 dqradb(ig,ijv,l,is) = dqradx(ig,ijv,l,is) 822 enddo 823 do i=1,3 824 do j=1,3 825 dqrad(1,ijv,l,is,i,j) = & 826 -qradb(1,ijv,l,is) * ainv(j,i) 827 do ig=2,ngb 828 dqrad(ig,ijv,l,is,i,j) = & 829 & -qradb(ig,ijv,l,is)*ainv(j,i) & 830 & -c*dqradb(ig,ijv,l,is)* & 831 & gxb(i,ig)/gb(ig)* & 832 & (gxb(1,ig)*ainv(j,1)+ & 833 & gxb(2,ig)*ainv(j,2)+ & 834 & gxb(3,ig)*ainv(j,3)) 835 enddo 836 enddo 837 enddo 838 end do 839 enddo 840 enddo 841 ! 842 do iv= 1,nh(is) 843 do jv=iv,nh(is) 844 ! 845 ! compact indices because qgb is symmetric 846 ! 847 ijv = jv*(jv-1)/2 + iv 848 call dqvan2b(ngb,iv,jv,is,ylmb,dylmb,dqgbs,dqrad,qradb ) 849 do i=1,3 850 do j=1,3 851 do ig=1,ngb 852 dqgb(ig,ijv,is,i,j)=dqgbs(ig,i,j) 853 enddo 854 enddo 855 enddo 856 end do 857 end do 858 end do 859 deallocate(dqrad) 860 deallocate(dqgbs) 861 deallocate(dylmb) 862 deallocate(dqradb) 863 end if 864 865 deallocate( ylmb ) 866 deallocate( qradb ) 867 868 IF( ALLOCATED( qradx ) ) DEALLOCATE( qradx ) 869 IF( ALLOCATED( dqradx ) ) DEALLOCATE( dqradx ) 870 871 RETURN 872 END SUBROUTINE exact_qradb_x 873 874 875!------------------------------------------------------------------------------! 876 877 878 LOGICAL FUNCTION check_tables_x( gmax ) 879 ! 880 ! check table size against cell variations 881 ! 882 ! 883 USE kinds, ONLY : DP 884 USE betax, ONLY : refg, mmx 885 USE mp, ONLY : mp_max 886 USE mp_global, ONLY : intra_bgrp_comm 887 USE gvecw, ONLY : ngw 888 USE cell_base, ONLY : tpiba2 889 USE small_box, ONLY : tpibab 890 USE smallbox_gvec, ONLY : gb, ngb 891 USE gvect, ONLY : gg 892 USE fft_base, ONLY : dfftp 893 ! 894 IMPLICIT NONE 895 ! 896 REAL(DP), INTENT(OUT) :: gmax 897 REAL(DP) :: g2, g2b 898 ! 899 g2 = MAXVAL( gg( 1:dfftp%ngm ) ) 900 ! 901 g2 = g2 * tpiba2 / refg 902 ! 903 IF( ALLOCATED( gb ) ) THEN 904 ! 905 g2b = MAXVAL( gb( 1:ngb ) ) 906 g2b = g2b * tpibab * tpibab / refg 907 gmax = MAX( g2, g2b ) 908 ! 909 ELSE 910 ! 911 gmax = g2 912 ! 913 END IF 914 ! 915 CALL mp_max( gmax, intra_bgrp_comm ) 916 ! 917 check_tables_x = .FALSE. 918 IF( INT( gmax ) + 2 >= mmx ) check_tables_x = .TRUE. 919 ! 920 RETURN 921 END FUNCTION check_tables_x 922 923 924!------------------------------------------------------------------------------! 925 926 927 SUBROUTINE interpolate_beta_x( tpre ) 928 ! 929 ! interpolate array beta(ig,iv,is) 930 ! 931 ! 932 USE kinds, ONLY : DP 933 USE control_flags, only: iverbosity 934 USE constants, only: pi, fpi 935 USE io_global, only: stdout 936 USE gvecw, only: ngw 937 USE ions_base, only: nsp 938 USE gvect, only: gg, g, gstart 939 USE uspp_param, only: upf, lmaxq, lmaxkb, nh 940 USE uspp, only: qq_nt, nhtolm, beta, dbeta 941 USE cell_base, only: ainv, omega, tpiba2, tpiba 942 USE betax, ONLY : refg, betagx, dbetagx 943 944 IMPLICIT NONE 945 946 LOGICAL, INTENT(IN) :: tpre 947 948 REAL(DP), ALLOCATABLE :: ylm(:,:), dylm(:,:,:,:) 949 REAL(DP) :: c, g2, betagl, dbetagl 950 INTEGER :: is, iv, lp, ig, jj, i, j 951 952 ALLOCATE( ylm( ngw, (lmaxkb+1)**2 ) ) 953 CALL ylmr2 ( (lmaxkb+1)**2, ngw, g, gg, ylm) 954 ! 955 ! 956 do is = 1, nsp 957 ! 958 ! calculation of array beta(ig,iv,is) 959 ! 960 if( iverbosity > 2 ) WRITE( stdout,*) ' beta ' 961 c = fpi / sqrt(omega) 962 do iv = 1, nh(is) 963 lp = nhtolm( iv, is ) 964 do ig = gstart, ngw 965 g2 = gg( ig ) * tpiba * tpiba / refg 966 jj = int( g2 ) + 1 967 betagl = betagx( jj+1, iv, is ) * ( g2 - DBLE(jj-1) ) + betagx( jj, iv, is ) * ( DBLE(jj) - g2 ) 968 beta( ig, iv, is ) = c * ylm( ig, lp ) * betagl 969 end do 970 if( gstart == 2 ) then 971 beta( 1, iv, is ) = c * ylm( 1, lp ) * betagx( 1, iv, is ) 972 end if 973 end do 974 end do 975 976 if (tpre) then 977 ! 978 ! calculation of array dbeta required for stress, variable-cell 979 ! 980 allocate( dylm( ngw, (lmaxkb+1)**2, 3, 3 ) ) 981 ! 982 call dylmr2_( (lmaxkb+1)**2, ngw, g, gg, ainv, dylm ) 983 ! 984 do is = 1, nsp 985 if( iverbosity > 2 ) WRITE( stdout,*) ' dbeta ' 986 c = fpi / sqrt(omega) 987 do iv = 1, nh(is) 988 lp = nhtolm(iv,is) 989 if( ngw > 0 ) then 990 betagl = betagx(1,iv,is) 991 do i=1,3 992 do j=1,3 993 dbeta( 1, iv, is, i, j ) = -0.5d0 * beta( 1, iv, is ) * ainv( j, i ) & 994 & - c * dylm( 1, lp, i, j ) * betagl ! SEGNO 995 enddo 996 enddo 997 end if 998 do ig = gstart, ngw 999 g2 = gg(ig) * tpiba * tpiba / refg 1000 jj=int(g2)+1 1001 betagl = betagx( jj+1, iv, is ) * ( g2 - DBLE(jj-1) ) + & 1002 & betagx( jj , iv, is ) * ( DBLE(jj) - g2 ) 1003 dbetagl= dbetagx( jj+1, iv, is ) * ( g2 - DBLE(jj-1) ) + & 1004 & dbetagx( jj , iv, is ) * ( DBLE(jj) - g2 ) 1005 do i=1,3 1006 do j=1,3 1007 dbeta( ig, iv, is, i, j ) = & 1008 & - 0.5d0 * beta( ig, iv, is ) * ainv( j, i ) & 1009 & - c * dylm( ig, lp, i, j ) * betagl & ! SEGNO 1010 & - c * ylm ( ig, lp ) *dbetagl * g(i,ig)/gg(ig)& 1011 & * ( g( 1, ig ) * ainv( j, 1 ) + g( 2, ig ) * ainv( j, 2 ) + g( 3, ig ) * ainv( j, 3 ) ) 1012 end do 1013 end do 1014 end do 1015 end do 1016 end do 1017 ! 1018 deallocate(dylm) 1019 ! 1020 end if 1021 ! 1022 deallocate(ylm) 1023 1024 RETURN 1025 END SUBROUTINE interpolate_beta_x 1026 1027 1028!------------------------------------------------------------------------------! 1029 1030 1031 SUBROUTINE interpolate_qradb_x( tpre ) 1032 ! 1033 ! interpolate array qradb(ig,iv,is) 1034 ! 1035 ! 1036 USE kinds, ONLY : DP 1037 use control_flags, only: iprint, iverbosity 1038 use io_global, only: stdout 1039 use gvecw, only: ngw 1040 use cell_base, only: ainv 1041 use uspp, only: qq_nt, nhtolm, beta 1042 use constants, only: pi, fpi 1043 use ions_base, only: nsp 1044 use uspp_param, only: upf, lmaxq, lmaxkb, nbetam, nh 1045 use qgb_mod, only: qgb, dqgb 1046 use smallbox_gvec, only: gb, gxb, ngb 1047 use small_box, only: omegab, tpibab 1048 USE betax, ONLY: qradx, dqradx, refg, mmx 1049! 1050 implicit none 1051 1052 LOGICAL, INTENT(IN) :: tpre 1053 1054 integer is, l, ig, ir, iv, jv, ijv, i,j, jj, ierr 1055 real(dp), allocatable:: fint(:), jl(:), dqradb(:,:,:,:) 1056 real(dp), allocatable:: ylmb(:,:), dylmb(:,:,:,:) 1057 REAL(DP), ALLOCATABLE :: dqrad( :, :, :, :, :, : ) 1058 REAL(DP), ALLOCATABLE :: qradb( :, :, :, : ) 1059 complex(dp), allocatable:: dqgbs(:,:,:) 1060 real(dp) xg, c, betagl, dbetagl, g2 1061 LOGICAL :: tvanp 1062 1063 tvanp = .FALSE. 1064 DO is = 1, nsp 1065 tvanp = tvanp .OR. upf(is)%tvanp 1066 END DO 1067 IF( .NOT. tvanp ) & 1068 return 1069 1070 allocate( qradb( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp ), STAT=ierr ) 1071 IF( ierr /= 0 ) & 1072 CALL errore(' interpolate_qradb ', ' cannot allocate qradb ', 1 ) 1073 ! 1074 qradb(:,:,:,:) = 0.d0 1075! 1076 allocate( ylmb( ngb, lmaxq*lmaxq ), STAT=ierr ) 1077 IF( ierr /= 0 ) & 1078 CALL errore(' interpolate_qradb ', ' cannot allocate ylmb ', 1 ) 1079! 1080 call ylmr2 (lmaxq*lmaxq, ngb, gxb, gb, ylmb) 1081 1082 do is = 1, nsp 1083 ! 1084 IF( .NOT. upf(is)%tvanp ) CYCLE 1085 ! 1086 ! calculation of array qradb(igb,iv,jv,is) 1087 ! 1088 ! qradb(ig,l,k,is) = 4pi/omega int_0^r dr r^2 j_l(qr) q(r,l,k,is) 1089 ! 1090 if( iverbosity > 2 ) WRITE( stdout,*) ' qradb ' 1091 ! 1092 c = fpi / omegab 1093 ! 1094 do iv= 1, upf(is)%nbeta 1095 do jv = iv, upf(is)%nbeta 1096 ijv = jv*(jv-1)/2 + iv 1097 do l=1, upf(is)%nqlc 1098 qradb(1,ijv,l,is) = c * qradx(1,ijv,l,is) 1099 end do 1100 do ig=2,ngb 1101 g2=gb(ig)*tpibab*tpibab/refg 1102 jj=int(g2)+1 1103 do l=1,upf(is)%nqlc 1104 if(jj.ge.mmx) then 1105 qradb(ig,ijv,l,is)=0.d0 1106 else 1107 qradb(ig,ijv,l,is)= & 1108 & c*qradx(jj+1,ijv,l,is)*(g2-DBLE(jj-1))+ & 1109 & c*qradx(jj,ijv,l,is)*(DBLE(jj)-g2) 1110 endif 1111 enddo 1112 enddo 1113 enddo 1114 enddo 1115! 1116! --------------------------------------------------------------- 1117! stocking of qgb(igb,ijv,is) and of qq(iv,jv,is) 1118! --------------------------------------------------------------- 1119 do iv= 1,nh(is) 1120 do jv=iv,nh(is) 1121! 1122! compact indices because qgb is symmetric 1123! 1124 ijv = jv*(jv-1)/2 + iv 1125 call qvan2b(ngb,iv,jv,is,ylmb,qgb(1,ijv,is),qradb ) 1126! 1127 qq_nt(iv,jv,is)=omegab*DBLE(qgb(1,ijv,is)) 1128 qq_nt(jv,iv,is)=qq_nt(iv,jv,is) 1129! 1130 end do 1131 end do 1132 1133 end do 1134! 1135 if (tpre) then 1136! --------------------------------------------------------------- 1137! arrays required for stress calculation, variable-cell dynamics 1138! --------------------------------------------------------------- 1139 allocate(dqradb(ngb,nbetam*(nbetam+1)/2,lmaxq,nsp)) 1140 allocate(dylmb(ngb,lmaxq*lmaxq,3,3)) 1141 allocate(dqgbs(ngb,3,3)) 1142 allocate( dqrad( ngb, nbetam*(nbetam+1)/2, lmaxq, nsp, 3, 3 ) ) 1143 dqrad(:,:,:,:,:,:) = 0.d0 1144 ! 1145 call dylmr2_( lmaxq*lmaxq, ngb, gxb, gb, ainv, dylmb ) 1146 ! 1147 do is=1,nsp 1148 ! 1149 IF( .NOT. upf(is)%tvanp ) CYCLE 1150 ! 1151 do iv= 1, upf(is)%nbeta 1152 do jv=iv, upf(is)%nbeta 1153 ijv = jv*(jv-1)/2 + iv 1154 do l=1,upf(is)%nqlc 1155 dqradb(1,ijv,l,is) = dqradx(1,ijv,l,is) 1156 do ig=2,ngb 1157 g2=gb(ig)*tpibab*tpibab/refg 1158 jj=int(g2)+1 1159 if(jj.ge.mmx) then 1160 dqradb(ig,ijv,l,is) = 0.d0 1161 else 1162 dqradb(ig,ijv,l,is) = & 1163 dqradx(jj+1,ijv,l,is)*(g2-DBLE(jj-1)) + & 1164 dqradx(jj,ijv,l,is)*(DBLE(jj)-g2) 1165 endif 1166 enddo 1167 do i=1,3 1168 do j=1,3 1169 dqrad(1,ijv,l,is,i,j) = - qradb(1,ijv,l,is) * ainv(j,i) 1170 do ig=2,ngb 1171 dqrad(ig,ijv,l,is,i,j) = & 1172 & - qradb(ig,ijv,l,is)*ainv(j,i) & 1173 & - c * dqradb(ig,ijv,l,is)* & 1174 & gxb(i,ig)/gb(ig)* & 1175 & (gxb(1,ig)*ainv(j,1)+ & 1176 & gxb(2,ig)*ainv(j,2)+ & 1177 & gxb(3,ig)*ainv(j,3)) 1178 enddo 1179 enddo 1180 enddo 1181 end do 1182 enddo 1183 enddo 1184 ! 1185 do iv= 1,nh(is) 1186 do jv=iv,nh(is) 1187 ! 1188 ! compact indices because qgb is symmetric 1189 ! 1190 ijv = jv*(jv-1)/2 + iv 1191 call dqvan2b(ngb,iv,jv,is,ylmb,dylmb,dqgbs,dqrad,qradb ) 1192 do i=1,3 1193 do j=1,3 1194 do ig=1,ngb 1195 dqgb(ig,ijv,is,i,j)=dqgbs(ig,i,j) 1196 enddo 1197 enddo 1198 enddo 1199 end do 1200 end do 1201 end do 1202 deallocate(dqrad) 1203 deallocate(dqgbs) 1204 deallocate(dylmb) 1205 deallocate(dqradb) 1206 end if 1207 deallocate(ylmb) 1208 deallocate(qradb) 1209 1210 RETURN 1211 END SUBROUTINE interpolate_qradb_x 1212 1213 1214!------------------------------------------------------------------------------! 1215 1216 1217 SUBROUTINE exact_beta_x( tpre ) 1218 ! 1219 ! compute array beta without interpolation 1220 ! 1221 ! 1222 USE control_flags, only : iverbosity 1223 USE kinds, ONLY : DP 1224 USE constants, only : pi, fpi 1225 USE io_global, only : stdout 1226 USE gvecw, only : ngw 1227 USE ions_base, only : nsp 1228 USE uspp_param, only : upf, lmaxq, lmaxkb, nh, nhm 1229 USE uspp, only : qq_nt, nhtolm, beta, nhtol, indv, dbeta 1230 USE cell_base, only : ainv, omega, tpiba2, tpiba 1231 USE atom, ONLY : rgrid 1232 USE gvect, only : gg, g, gstart 1233 1234 IMPLICIT NONE 1235 1236 LOGICAL, INTENT(IN) :: tpre 1237 1238 REAL(DP), ALLOCATABLE :: ylm(:,:), dylm(:,:,:,:) 1239 REAL(DP) :: c, g2, betagl, dbetagl 1240 INTEGER :: is, iv, lp, ig, jj, i, j, nr 1241 INTEGER :: l, il, ir 1242 REAL(DP), ALLOCATABLE :: dfint(:), djl(:), fint(:), jl(:) 1243 REAL(DP), ALLOCATABLE :: betagx ( :, :, : ), dbetagx( :, :, : ) 1244 REAL(DP) :: xg 1245 1246 IF( .NOT. ALLOCATED( rgrid ) ) & 1247 CALL errore( ' exact_beta_x ', ' rgrid not allocated ', 1 ) 1248 IF( .NOT. ALLOCATED( upf ) ) & 1249 CALL errore( ' exact_beta_x ', ' upf not allocated ', 1 ) 1250 1251 ALLOCATE( ylm( ngw, (lmaxkb+1)**2 ) ) 1252 ALLOCATE( betagx ( ngw, nhm, nsp ) ) 1253 IF (tpre) ALLOCATE( dbetagx( ngw, nhm, nsp ) ) 1254 1255 CALL ylmr2 ( (lmaxkb+1)**2, ngw, g, gg, ylm) 1256 1257 ! 1258 do is = 1, nsp 1259 ! 1260 nr = upf(is)%kkbeta 1261 ! 1262 if ( tpre ) then 1263 allocate( dfint( nr ) ) 1264 allocate( djl ( nr ) ) 1265 end if 1266 ! 1267 allocate( fint ( nr ) ) 1268 allocate( jl ( nr ) ) 1269 ! 1270 do iv = 1, nh(is) 1271 ! 1272 l = nhtol(iv,is) 1273 ! 1274 do il = 1, ngw 1275 ! 1276 xg = sqrt( gg( il ) * tpiba * tpiba ) 1277 call sph_bes (nr, rgrid(is)%r, xg, l, jl ) 1278 ! 1279 if( tpre )then 1280 ! 1281 call sph_dbes1 ( nr, rgrid(is)%r, xg, l, jl, djl) 1282 ! 1283 endif 1284 ! 1285 ! beta(ir)=r*beta(r) 1286 ! 1287 do ir = 1, nr 1288 fint(ir) = rgrid(is)%r(ir) * jl(ir) * & 1289 upf(is)%beta( ir, indv(iv,is) ) 1290 end do 1291 call simpson_cp90(nr,fint,rgrid(is)%rab,betagx(il,iv,is)) 1292 ! 1293 if(tpre) then 1294 do ir = 1, nr 1295 dfint(ir) = rgrid(is)%r(ir) * djl(ir) * & 1296 upf(is)%beta( ir, indv(iv,is) ) 1297 end do 1298 call simpson_cp90(nr,dfint,rgrid(is)%rab,dbetagx(il,iv,is)) 1299 endif 1300 ! 1301 end do 1302 end do 1303! 1304 deallocate(jl) 1305 deallocate(fint) 1306 ! 1307 if (tpre) then 1308 deallocate(djl) 1309 deallocate(dfint) 1310 end if 1311 ! 1312 end do 1313 ! 1314 do is = 1, nsp 1315 ! 1316 ! calculation of array beta(ig,iv,is) 1317 ! 1318 if( iverbosity > 2 ) WRITE( stdout,*) ' beta ' 1319 c = fpi / sqrt(omega) 1320 do iv = 1, nh(is) 1321 lp = nhtolm( iv, is ) 1322 do ig = 1, ngw 1323 betagl = betagx( ig, iv, is ) 1324 beta( ig, iv, is ) = c * ylm( ig, lp ) * betagl 1325 end do 1326 end do 1327 end do 1328 1329 if (tpre) then 1330 ! 1331 ! calculation of array dbeta required for stress, variable-cell 1332 ! 1333 allocate( dylm( ngw, (lmaxkb+1)**2, 3, 3 ) ) 1334 ! 1335 call dylmr2_( (lmaxkb+1)**2, ngw, g, gg, ainv, dylm ) 1336 ! 1337 do is = 1, nsp 1338 if( iverbosity > 2 ) WRITE( stdout,*) ' dbeta ' 1339 c = fpi / sqrt(omega) 1340 do iv = 1, nh(is) 1341 lp = nhtolm(iv,is) 1342 betagl = betagx(1,iv,is) 1343 do i=1,3 1344 do j=1,3 1345 dbeta(1,iv,is,i,j)=-0.5d0*beta(1,iv,is)*ainv(j,i) & 1346 & -c*dylm(1,lp,i,j)*betagl ! SEGNO 1347 enddo 1348 enddo 1349 do ig=gstart,ngw 1350 betagl = betagx(ig,iv,is) 1351 dbetagl= dbetagx(ig,iv,is) 1352 do i=1,3 1353 do j=1,3 1354 dbeta(ig,iv,is,i,j)= & 1355 & -0.5d0*beta(ig,iv,is)*ainv(j,i) & 1356 & -c*dylm(ig,lp,i,j)*betagl & ! SEGNO 1357 & -c*ylm (ig,lp)*dbetagl*g(i,ig)/gg(ig) & 1358 & *(g(1,ig)*ainv(j,1)+ & 1359 & g(2,ig)*ainv(j,2)+ & 1360 & g(3,ig)*ainv(j,3)) 1361 end do 1362 end do 1363 end do 1364 end do 1365 end do 1366 ! 1367 deallocate(dylm) 1368 ! 1369 end if 1370 ! 1371 deallocate(ylm) 1372 IF( ALLOCATED( betagx ) ) DEALLOCATE( betagx ) 1373 IF( ALLOCATED( dbetagx ) ) DEALLOCATE( dbetagx ) 1374 1375 RETURN 1376 END SUBROUTINE exact_beta_x 1377! 1378! 1379!------------------------------------------------------------------------------! 1380! 1381! 1382 SUBROUTINE fill_qrl_x( is, qrl ) 1383 ! 1384 ! fill l-components of Q(r) as in Vanderbilt's approach 1385 ! 1386 USE uspp_param, ONLY: upf 1387 USE atom, ONLY: rgrid 1388 USE kinds, ONLY: DP 1389 USE io_global, ONLY: stdout 1390 ! 1391 IMPLICIT NONE 1392 ! 1393 INTEGER, INTENT(IN) :: is 1394 REAL(DP), INTENT(OUT) :: qrl( :, :, : ) 1395 ! 1396 INTEGER :: iv, jv, ijv, lmin, lmax, l, ir, i 1397 INTEGER :: dim1, dim2, dim3 1398 ! 1399 IF( .NOT. ALLOCATED( rgrid ) ) & 1400 CALL errore( ' fill_qrl_x ', ' rgrid not allocated ', 1 ) 1401 IF( .NOT. ALLOCATED( upf ) ) & 1402 CALL errore( ' fill_qrl_x ', ' upf not allocated ', 1 ) 1403 1404 dim1 = SIZE( qrl, 1 ) 1405 dim2 = SIZE( qrl, 2 ) 1406 dim3 = SIZE( qrl, 3 ) 1407 ! 1408 IF ( upf(is)%kkbeta > dim1 ) & 1409 CALL errore ('fill_qrl', 'bad 1st dimension for array qrl', 1) 1410 ! 1411 qrl = 0.0d0 1412 ! 1413 do iv = 1, upf(is)%nbeta 1414 ! 1415 do jv = iv, upf(is)%nbeta 1416 ! 1417 ijv = (jv-1)*jv/2 + iv 1418 ! 1419 IF ( ijv > dim2) & 1420 CALL errore ('fill_qrl', 'bad 2nd dimension for array qrl', 2) 1421 1422 ! notice that L runs from 1 to Lmax+1 1423 1424 lmin = ABS (upf(is)%lll(jv) - upf(is)%lll(iv)) + 1 1425 lmax = upf(is)%lll(jv) + upf(is)%lll(iv) + 1 1426 1427 ! WRITE( stdout, * ) 'QRL is, jv, iv = ', is, jv, iv 1428 ! WRITE( stdout, * ) 'QRL lll jv, iv = ', upf(is)%lll(jv), upf(is)%lll(iv) 1429 ! WRITE( stdout, * ) 'QRL lmin, lmax = ', lmin, lmax 1430 ! WRITE( stdout, * ) '---------------- ' 1431 1432 IF ( lmin < 1 .OR. lmax > dim3) THEN 1433 WRITE( stdout, * ) ' lmin, lmax = ', lmin, lmax 1434 CALL errore ('fill_qrl', 'bad 3rd dimension for array qrl', 3) 1435 END IF 1436 1437 1438 do l = lmin, lmax 1439 do ir = 1, upf(is)%kkbeta 1440 IF( upf(is)%tvanp ) THEN 1441 ! BEWARE: index l in upf%qfuncl(l) runs from 0 to lmax, not from 1 to lmax+1 1442 qrl(ir,ijv,l)=upf(is)%qfuncl(ir,ijv,l-1) 1443 ENDIF 1444 end do 1445 end do 1446 end do 1447 end do 1448 RETURN 1449 END SUBROUTINE fill_qrl_x 1450