1! 2! Copyright (C) 2003 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! 9!----------------------------------------------------------------------- 10SUBROUTINE dielec(do_zstar) 11 !----------------------------------------------------------------------- 12 ! 13 ! calculates the dielectric tensor and effective charges 14 ! 15 USE constants, ONLY : fpi 16 USE cell_base, ONLY : omega 17 USE ions_base, ONLY : nat, zv, ityp 18 USE mp_pools, ONLY : intra_pool_comm 19 USE mp, ONLY : mp_sum 20 USE io_files, ONLY : seqopn 21 USE klist, ONLY : wk, ngk 22 USE wvfct, ONLY: nbnd, npwx 23 USE cgcom 24 25 IMPLICIT NONE 26 LOGICAL :: do_zstar 27 ! 28 INTEGER :: npw,ibnd,ipol,jpol,na,nu,ik 29 CHARACTER(len=7) :: filbar, fildwf 30 real(DP) :: w, weight 31 real(DP), ALLOCATABLE :: work(:,:) 32 COMPLEX(DP), ALLOCATABLE :: dpsi2(:,:), dpsi3(:,:) 33 LOGICAL :: done 34 ! 35 CALL start_clock('dielec') 36 ! 37 ALLOCATE (dpsi2( npwx, nbnd)) 38 ALLOCATE (dpsi3( npwx, nbnd)) 39 ALLOCATE (work( nbnd, 3)) 40 ! 41 epsilon0(:,:) = 0.d0 42 IF (do_zstar) zstar (:,:,:) = 0.d0 43 ! do ik=1,nks 44 ik = 1 45 npw= ngk(ik) 46 weight = wk(ik) 47 w = fpi/omega * weight 48 ! 49 !** calculate Effective Charges (<DeltaV*psi(ion)|DeltaPsi(E)>) 50 ! 51 ! read DeltaPsi(E) 52 ! pol. 1 53 ipol=1 54 iudwf=10+ipol 55 WRITE(fildwf,'("fildwx",i1)') ipol 56 CALL seqopn (iudwf,fildwf,'unformatted',done) 57 READ (iudwf) dpsi 58 CLOSE(unit=iudwf) 59 ! pol. 2 60 ipol=2 61 iudwf=10+ipol 62 WRITE(fildwf,'("fildwx",i1)') ipol 63 CALL seqopn (iudwf,fildwf,'unformatted',done) 64 READ (iudwf) dpsi2 65 CLOSE(unit=iudwf) 66 ! pol. 3 67 ipol=3 68 iudwf=10+ipol 69 WRITE(fildwf,'("fildwx",i1)') ipol 70 CALL seqopn (iudwf,fildwf,'unformatted',done) 71 READ (iudwf) dpsi3 72 CLOSE(unit=iudwf) 73 ! 74 IF (.not.do_zstar) GOTO 10 75 ! 76 DO nu = 1,nmodes 77 na = (nu-1)/3+1 78 IF (has_equivalent(na)==0) THEN 79 ! DeltaV*psi(ion) for mode nu is recalculated 80 CALL dvpsi_kb(ik,nu) 81 ! 82 jpol= mod(nu-1,3)+1 83 ! work is the real part of <DeltaV*Psi(ion)|DeltaPsi(E)> 84 CALL pw_dot('N',npw,nbnd,dvpsi,npwx,dpsi ,npwx,work(1,1)) 85 CALL pw_dot('N',npw,nbnd,dvpsi,npwx,dpsi2,npwx,work(1,2)) 86 CALL pw_dot('N',npw,nbnd,dvpsi,npwx,dpsi3,npwx,work(1,3)) 87 DO ipol = 1,3 88 DO ibnd = 1,nbnd 89 zstar(ipol,jpol,na) = zstar(ipol,jpol,na) + 2.0d0*weight*work(ibnd,ipol) 90 ENDDO 91 ENDDO 92 ENDIF 93 ENDDO 9410 CONTINUE 95 !** calculate Dielectric Tensor (<DeltaV*psi(E)\DeltaPsi(E)>) 96 ! 97 DO jpol=1,3 98 ! read DeltaV*Psi(elec) for polarization jpol 99 iubar=jpol 100 WRITE(filbar,'("filbar",i1)') iubar 101 CALL seqopn (iubar,filbar,'unformatted',done) 102 READ (iubar) dvpsi 103 CLOSE(iubar) 104 ! now work is the real part of <DeltaV*psi(E)|DeltaPsi(E)> 105 CALL pw_dot('N',npw,nbnd,dvpsi,npwx,dpsi ,npwx,work(1,1)) 106 CALL pw_dot('N',npw,nbnd,dvpsi,npwx,dpsi2,npwx,work(1,2)) 107 CALL pw_dot('N',npw,nbnd,dvpsi,npwx,dpsi3,npwx,work(1,3)) 108 DO ipol = 1,3 109 DO ibnd = 1,nbnd 110 epsilon0(ipol,jpol) = epsilon0(ipol,jpol) + 4.0d0*w*work(ibnd,ipol) 111 ENDDO 112 ENDDO 113 ENDDO 114 ! end do 115#if defined(__MPI) 116 IF (do_zstar) CALL mp_sum( zstar, intra_pool_comm ) 117 CALL mp_sum( epsilon0, intra_pool_comm ) 118#endif 119 DEALLOCATE(work) 120 DEALLOCATE(dpsi3) 121 DEALLOCATE(dpsi2) 122 ! 123 ! add the diagonal part 124 ! 125 DO ipol=1,3 126 epsilon0(ipol,ipol) = epsilon0(ipol,ipol) + 1.0d0 127 IF (do_zstar) THEN 128 DO na=1,nat 129 zstar(ipol,ipol,na) = zstar(ipol,ipol,na) + zv(ityp(na)) 130 ENDDO 131 ENDIF 132 ENDDO 133 ! 134 CALL stop_clock('dielec') 135 ! 136 RETURN 137END SUBROUTINE dielec 138