1! 2! Copyright (C) 2001-2013 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 calculate_wing(n_set, orthonorm) 11 12!this subroutine calculate the terms 13!\Sum_G <G|\tilde{w^P_i}>\epsilon(G'=0, G, iw) 14! it requires the file .e_head 15 USE io_global, ONLY : stdout, ionode, ionode_id 16 USE io_files, ONLY : prefix, tmp_dir, diropn 17 USE kinds, ONLY : DP 18 USE wannier_gw 19 USE mp, ONLY : mp_bcast, mp_sum 20 USE gvect, ONLY : mill, ngm, gstart,g,ngm_g, ig_l2g 21 USE cell_base, ONLY : tpiba 22 USE mp_wave, ONLY : mergewf,splitwf 23 USE mp_pools, ONLY : intra_pool_comm 24 USE mp_world, ONLY : mpime, nproc, world_comm 25 USE wvfct, ONLY : npwx, npw 26 USE cell_base, ONLY : at,bg 27 28 implicit none 29 30 INTEGER, EXTERNAL :: find_free_unit 31 32 INTEGER, INTENT(in) :: n_set !defines the number of states to be read from disk at the same time 33 INTEGER, INTENT(in) :: orthonorm!if ==1 opens orthonormalized products of wannier file, if ==2 reduced one 34 35 36 INTEGER iun, iungprod 37 INTEGER :: n_g, ngm_k 38 INTEGER :: ig, igg, iw, iiw, i, ii, ipol 39 REAL(kind=DP) :: omega_g 40 REAL(kind=DP), ALLOCATABLE :: freqs(:) 41 COMPLEX(kind=DP), ALLOCATABLE :: e_head(:,:,:), e_head_g0(:) 42 LOGICAL :: exst 43 COMPLEX(kind=DP), ALLOCATABLE :: tmpspacei(:,:) 44 REAL(kind=DP), ALLOCATABLE :: wing(:,:,:), wing_c(:,:,:) 45 REAL(kind=DP) :: sca 46 REAL(kind=DP), ALLOCATABLE :: fact(:) 47 REAL(kind=DP) :: qq 48 INTEGER :: npwx_g 49 INTEGER, ALLOCATABLE :: k2g_ig_l2g(:) 50 51 npwx_g=npwx 52 call mp_sum(npwx_g,world_comm) 53 54!read file .e_head 55 56 write(stdout,*) 'Routine calculate_wing' 57 FLUSH(stdout) 58 allocate(fact(ngm)) 59 allocate(k2g_ig_l2g(ngm)) 60 61 fact(:)=0.d0 62 if(gstart==2) fact(1)=0.d0 63 do ig=gstart,npw 64! qq = g(1,ig)**2.d0 + g(2,ig)**2.d0 + g(3,ig)**2.d0 65! fact(ig)=1.d0/tpiba/dsqrt(qq) 66 fact(ig)=dsqrt(vg_q(ig)) 67 end do 68 69 70 call ktogamma_ig_l2g ( k2g_ig_l2g, at, bg ) 71 72 write(stdout,*) 'ATT0.1' 73 FLUSH(stdout) 74 75 76 77 if(ionode) then 78 iun = find_free_unit() 79 open( unit= iun, file=trim(tmp_dir)//'/_ph0/'//trim(prefix)//'.e_head', status='old',form='unformatted') 80 read(iun) n_g 81 read(iun) omega_g 82 endif 83 84 85 call mp_bcast(n_g, ionode_id,world_comm) 86 call mp_bcast(omega_g, ionode_id,world_comm) 87 allocate(freqs(n_g+1)) 88 89 if(ionode) then 90 read(iun) freqs(1:n_g+1) 91 read(iun) ngm_k 92 endif 93 94 write(stdout,*) 'ATT0.2' 95 FLUSH(stdout) 96 97 98 call mp_bcast(freqs(:), ionode_id,world_comm) 99 call mp_bcast(ngm_k, ionode_id,world_comm) 100 101 allocate(e_head_g0(ngm_k)) 102 allocate(e_head(npw, n_g+1,3)) 103 e_head(:,:,:) = (0.d0,0.d0) 104 do ii=1,n_g+1 105 do ipol=1,3 106 e_head_g0(:)=(0.d0,0.d0) 107 if(ionode) read(iun) e_head_g0(1:ngm_k) 108 call splitwf(e_head(:, ii,ipol),e_head_g0,npw,k2g_ig_l2g,mpime,nproc,ionode_id,intra_pool_comm) 109 enddo 110 enddo 111 if(ionode) close(iun) 112 113 write(stdout,*) 'ATT1' 114 FLUSH(stdout) 115 116 deallocate(e_head_g0) 117!loop on n_set groups 118 write(stdout,*) 'ATT2' 119 FLUSH(stdout) 120 121 122 allocate(tmpspacei(max_ngm,n_set)) 123 iungprod = find_free_unit() 124 125 if(orthonorm==0) then 126 CALL diropn( iungprod, 'wiwjwfc', max_ngm*2, exst ) 127 else if(orthonorm==1) then 128 CALL diropn( iungprod, 'wiwjwfc_on', max_ngm*2, exst ) 129 else 130 CALL diropn( iungprod, 'wiwjwfc_red', max_ngm*2, exst ) 131 endif 132 allocate(wing(numw_prod, n_g+1,3)) 133 wing(:,:,:)=0.d0 134 !allocate(wing_c(numw_prod, n_g+1,3)) 135 !wing_c(:,:,:)=0.d0 136 137 138 139 do iiw=1,ceiling(real(numw_prod)/real(n_set)) 140!read states 141 do iw=(iiw-1)*n_set+1,min(iiw*n_set,numw_prod) 142 CALL davcio(tmpspacei(:,iw-(iiw-1)*n_set),max_ngm*2,iungprod,iw,-1) 143 enddo 144 145 146 write(stdout,*) 'ATT3' 147 FLUSH(stdout) 148 149 150 !loop on states 151 do iw=(iiw-1)*n_set+1,min(iiw*n_set,numw_prod) 152 do ipol=1,3 153 do i=1, n_g+1 154 sca=0.d0 155 do ig=1,max_ngm 156! sca=sca+2.d0*real(tmpspacei(ig,iw-(iiw-1)*n_set)*conjg(e_head(ig,i,ipol)))*fact(ig) 157 sca=sca+2.d0*dble((tmpspacei(ig,iw-(iiw-1)*n_set))*conjg(e_head(ig,i,ipol)))!*fact(ig)!ATTENZIONE 158 enddo 159 call mp_sum(sca,world_comm) 160 wing(iw,i,ipol)=sca 161 162 163 enddo 164 enddo 165 166 enddo 167 enddo 168 169 write(stdout,*) 'ATT4' 170 FLUSH(stdout) 171 172 173!write terms on file 174 175 if(ionode) then 176 iun = find_free_unit() 177 open( unit= iun, file=trim(tmp_dir)//trim(prefix)//'.wing', status='unknown',form='unformatted') 178 write(iun) n_g 179 write(iun) omega_g 180 write(iun) numw_prod 181 do ipol=1,3 182 do i=1,n_g+1 183 write(iun) wing(1:numw_prod,i,ipol) 184 enddo 185! do i=1,n_g+1 186! write(iun) wing_c(1:numw_prod,i,ipol) 187! enddo 188 enddo 189 close(iun) 190 endif 191 192 deallocate(tmpspacei) 193 close(iungprod) 194 deallocate(fact) 195 deallocate(freqs) 196 ! if(ionode) deallocate (e_head_g) 197 deallocate(e_head) 198 deallocate(wing) 199 deallocate(k2g_ig_l2g) 200 return 201end subroutine calculate_wing 202 203 204 !----------------------------------------------------------------------- 205 SUBROUTINE ktogamma_ig_l2g ( k2g_ig_l2g, at, bg ) 206 !---------------------------------------------------------------------- 207 ! 208 ! This routine generates all the reciprocal lattice vectors 209 ! contained in the sphere of radius gcutm. Furthermore it 210 ! computes the indices nl which give the correspondence 211 ! between the fft mesh points and the array of g vectors. 212 ! 213 USE gvect, ONLY : ig_l2g, g, gg, ngm, ngm_g, gcutm, & 214 mill, gstart 215 USE fft_base, ONLY : dfftp, dffts 216! 217 USE kinds, ONLY : DP 218 USE constants, ONLY : eps8 219 220 221 222 223 224 IMPLICIT NONE 225 ! 226 INTEGER, INTENT(out) :: k2g_ig_l2g(ngm) 227 REAL(DP), INTENT(IN) :: at(3,3), bg(3,3) 228 ! here a few local variables 229 ! 230 REAL(DP) :: t (3), tt 231 INTEGER :: ngm_, n1, n2, n3, n1s, n2s, n3s 232 ! 233 REAL(DP), ALLOCATABLE :: g2sort_g(:) 234 ! array containing all g vectors, on all processors: replicated data 235 INTEGER, ALLOCATABLE :: mill_g(:,:), mill_unsorted(:,:) 236 ! array containing all g vectors generators, on all processors: 237 ! replicated data 238 INTEGER, ALLOCATABLE :: igsrt(:) 239 ! 240 INTEGER :: ni, nj, nk, i, j, k, ipol, ng, igl, indsw,ig 241 ! 242 ! counters 243 ! 244 ! set the total number of fft mesh points and and initial value of gg 245 ! The choice of gcutm is due to the fact that we have to order the 246 ! vectors after computing them. 247 ! 248 ! 249 ! set d vector for unique ordering 250 ! 251 ! and computes all the g vectors inside a sphere 252 ! 253 ALLOCATE( mill_g( 3, ngm_g*3 ),mill_unsorted( 3, ngm_g*3 ) ) 254 ALLOCATE( igsrt( ngm_g*3 ) ) 255 ALLOCATE( g2sort_g( ngm_g*3 ) ) 256 g2sort_g(:) = 1.0d20 257 ! 258 ! save present value of ngm 259 ! 260 261 ! 262 ngm_ = 0 263 264 ! 265 ! max miller indices (same convention as in module stick_set) 266 ! 267 ni = (dfftp%nr1-1)/2 268 nj = (dfftp%nr2-1)/2 269 nk = (dfftp%nr3-1)/2 270 ! 271 iloop: DO i = -ni, ni 272 ! 273 274 jloop: DO j = -nj, nj 275 ! 276 277 kloop: DO k = -nk, nk 278 ! 279 280 t(:) = i * bg (:,1) + j * bg (:,2) + k * bg (:,3) 281 tt = sum(t(:)**2) 282 IF (tt <= gcutm) THEN 283 ngm_ = ngm_ + 1 284 mill_unsorted( :, ngm_ ) = (/ i,j,k /) 285 IF ( tt > eps8 ) THEN 286 g2sort_g(ngm_) = tt 287 ELSE 288 g2sort_g(ngm_) = 0.d0 289 ENDIF 290 ENDIF 291 ENDDO kloop 292 ENDDO jloop 293 ENDDO iloop 294 295 296 igsrt(1) = 0 297 CALL hpsort_eps( ngm_, g2sort_g, igsrt, eps8 ) 298 mill_g(1,1:ngm_) = mill_unsorted(1,igsrt(1:ngm_)) 299 mill_g(2,1:ngm_) = mill_unsorted(2,igsrt(1:ngm_)) 300 mill_g(3,1:ngm_) = mill_unsorted(3,igsrt(1:ngm_)) 301 DEALLOCATE( g2sort_g, igsrt, mill_unsorted ) 302 303 304 305 do ig=1,ngm 306 do ng=1,ngm_ 307 if(mill_g(1,ng)==mill(1,ig).and.mill_g(2,ng)==mill(2,ig).and.mill_g(3,ng)==mill(3,ig)) then 308 k2g_ig_l2g(ig)=ng 309 exit 310 endif 311 enddo 312 enddo 313 314 ! 315 DEALLOCATE( mill_g ) 316 317 318 319 END SUBROUTINE ktogamma_ig_l2g 320