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