1! 2! Copyright (C) 2001-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! 8SUBROUTINE compute_qdipol( dpqq ) 9 !! This routine computes the term dpqq, i.e. the dipole moment of the 10 !! augmentation charge. The output is given on cartesian coordinates. 11 ! 12 USE kinds, ONLY: DP 13 USE constants, ONLY: fpi 14 USE atom, ONLY: rgrid 15 USE ions_base, ONLY: ntyp => nsp 16 USE uspp, ONLY: nhtol, nhtolm, indv, nlx, ap 17 USE uspp_param, ONLY: upf, nbetam, nh, nhm 18 ! 19 IMPLICIT NONE 20 ! 21 REAL(DP) :: dpqq(nhm, nhm, 3, ntyp) 22 !! Dipole moment of augmentation charge 23 REAL(DP), ALLOCATABLE :: qrad2(:,:,:), qtot(:,:,:), aux(:) 24 REAL(DP) :: fact 25 INTEGER :: nt, l, ir, nb, mb, ijv, ilast, ipol, ih, ivl, jh, jvl, lp, ndm 26 ! 27 CALL start_clock( 'cmpt_qdipol' ) 28 ! 29 ndm = MAXVAL( upf(1:ntyp)%kkbeta ) 30 ALLOCATE( qrad2( nbetam , nbetam, ntyp) ) 31 ALLOCATE( aux( ndm) ) 32 ALLOCATE( qtot( ndm, nbetam, nbetam) ) 33 34 qrad2(:,:,:)=0.d0 35 dpqq=0.d0 36 37 DO nt = 1, ntyp 38 IF ( upf(nt)%tvanp ) THEN 39 l=1 40 ! 41 ! Only l=1 terms enter in the dipole of Q 42 ! 43 DO nb = 1, upf(nt)%nbeta 44 DO mb = nb, upf(nt)%nbeta 45 ijv = mb * (mb-1) /2 + nb 46 IF ( ( l >= ABS(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .AND. & 47 ( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .AND. & 48 (MOD(l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) THEN 49 qtot(1:upf(nt)%kkbeta,nb,mb) = upf(nt)%qfuncl(1:upf(nt)%kkbeta,ijv,l) 50 ENDIF 51 ENDDO 52 ENDDO 53 do nb=1, upf(nt)%nbeta 54 ! 55 ! the Q are symmetric with respect to indices 56 ! 57 do mb=nb, upf(nt)%nbeta 58 if ( ( l >= ABS(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .AND. & 59 ( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .AND. & 60 (MOD(l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) THEN 61 do ir = 1, upf(nt)%kkbeta 62 aux(ir)=rgrid(nt)%r(ir)*qtot(ir, nb, mb) 63 ENDDO 64 call simpson ( upf(nt)%kkbeta, aux, rgrid(nt)%rab, & 65 qrad2(nb,mb,nt) ) 66 ENDIF 67 ENDDO 68 ENDDO 69 ENDIF 70 ! ntyp 71 ENDDO 72 73 DO ipol = 1, 3 74 fact = -SQRT(fpi/3.d0) 75 IF (ipol == 1) lp=3 76 IF (ipol == 2) lp=4 77 IF (ipol == 3) THEN 78 lp=2 79 fact=-fact 80 ENDIF 81 DO nt = 1, ntyp 82 IF ( upf(nt)%tvanp ) THEN 83 DO ih = 1, nh(nt) 84 ivl = nhtolm(ih, nt) 85 mb = indv(ih, nt) 86 DO jh = ih, nh(nt) 87 jvl = nhtolm(jh, nt) 88 nb = indv(jh,nt) 89 IF (ivl > nlx) CALL errore( 'compute_qdipol',' ivl > nlx', ivl ) 90 IF (jvl > nlx) CALL errore( 'compute_qdipol',' jvl > nlx', jvl ) 91 IF (nb > nbetam) & 92 CALL errore( 'compute_qdipol',' nb out of bounds', nb ) 93 IF (mb > nbetam) & 94 CALL errore( 'compute_qdipol',' mb out of bounds', mb ) 95 IF (mb > nb) CALL errore( 'compute_qdipol',' mb > nb', 1 ) 96 dpqq(ih,jh,ipol,nt) = fact * ap(lp,ivl,jvl) * qrad2(mb,nb,nt) 97 dpqq(jh,ih,ipol,nt) = dpqq(ih,jh,ipol,nt) 98 ! WRITE( stdout,'(3i5,2f15.9)') ih,jh,ipol,dpqq(ih,jh,ipol,nt) 99 ENDDO 100 ENDDO 101 ENDIF 102 ENDDO 103 ENDDO 104 ! 105 DEALLOCATE( qtot ) 106 DEALLOCATE( aux ) 107 DEALLOCATE( qrad2 ) 108 ! 109 CALL stop_clock( 'cmpt_qdipol' ) 110 ! 111 RETURN 112 ! 113END SUBROUTINE compute_qdipol 114