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