1!=========================================================================
2!
3! Utilities:
4!
5! (1) eps0sym     Originally by MJ      Last Modified 1/12/2010 (MJ)
6!
7!     This utility symmetrizes eps0mat file.
8!
9!=========================================================================
10
11#include "f_defs.h"
12
13program eps0sym
14
15  use global_m
16  use epsread_hdf5_m
17  use epswrite_hdf5_m
18#ifdef HDF5
19  use hdf5
20#endif
21  use write_matrix_m
22
23  implicit none
24
25  character :: ajname*6,adate*11,outfile*80,infile*80
26  real(DP)  :: ecuts
27  real(DP), allocatable :: dFreqGrid(:), epsdiag(:,:,:)
28  complex(DPC), allocatable :: dFreqBrd(:)
29  integer   :: i,ig,j, &
30    ng,ngq,nmtx, &
31    ii,qgrid(3),freq_dep,nFreq,nfreq_imag,&
32    jj,kk,ll,nargs,imap,iout, &
33    nq, gx, gy, gz, gmax, &
34    ig0, igx, igy, igz, jgx, jgy, jgz, jout
35  real(DP), allocatable :: ekin(:)
36  real(DP) :: q(3,1), qk(3,1), errorMax
37  SCALAR, allocatable :: eps(:,:), tempeps(:,:)
38  SCALAR :: errorMaxTemp
39  integer, allocatable :: isort(:),isorti(:),kx(:),ky(:),kz(:)
40  integer, allocatable :: minusgidx(:), map(:), old(:,:)
41  integer :: nmtx0_of_q(1), isize, error
42  logical :: use_hdf5
43#ifdef HDF5
44  integer(HID_T) :: file_id
45#endif
46
47  write(6,*) MYFLAVOR // ' version is used to symmetrize file'
48
49!------------------
50! Get file names from command-line arguments
51
52  nargs = iargc()
53
54  if (nargs .ne. 2) then
55    call die('Usage: eps0sym eps0mat_in eps0mat_out')
56  endif
57
58  call getarg(1,infile)
59  call getarg(2,outfile)
60
61  use_hdf5 = .false.
62#ifdef HDF5
63
64  call h5open_f(error)
65
66  call h5fopen_f(trim(infile), H5F_ACC_RDONLY_F, file_id, error)
67  if(error == 0) use_hdf5 = .true.
68
69  if(use_hdf5) then
70    call h5fclose_f(file_id, error)
71    call system("cp "//trim(infile)//" "//trim(outfile))
72    call h5fopen_f(trim(outfile), H5F_ACC_RDWR_F, file_id, error)
73    call read_eps_grid_sizes_hdf5(ngq, nq, ecuts, nfreq, nfreq_imag, nmtx, qgrid, freq_dep, infile)
74    if (freq_dep.ne.0) then
75      call die('eps0sym: Full frequency not supported')
76    endif
77    ng = ngq
78
79    call read_eps_qgrid_hdf5(nq,q,nmtx0_of_q,infile)
80    qk = q
81  else
82
83#endif
84
85    call open_file(unit=11,file=TRUNC(infile),form='unformatted',status='old')
86
87!--------------------------
88!  Read in the eps0mat file
89!
90
91    read(11) ajname,adate
92    read(11) freq_dep,nFreq
93    if (freq_dep.ne.0) then
94      call die('eps0sym: Full frequency not supported')
95    endif
96    read(11) (qgrid(ii),ii=1,3)
97    if (freq_dep .eq. 2) then
98      SAFE_ALLOCATE(dFreqGrid,(nFreq))
99      SAFE_ALLOCATE(dFreqBrd,(nFreq))
100      read(11) (dFreqGrid(i),i=1,nFreq),(dFreqBrd(i),i=1,nFreq)
101    else
102      read(11)
103    endif
104    read(11)
105    read(11)
106    read(11) ecuts
107    read(11) nq,(q(j,1),j=1,3)
108    !if (nq .ne. 1) then
109    !  call die('This only works for the q->0 point.')
110    !endif
111    read(11) ng
112    read(11) ngq,nmtx
113
114    call close_file(11)
115
116#ifdef HDF5
117  endif
118#endif
119
120  SAFE_ALLOCATE(kx, (ng))
121  SAFE_ALLOCATE(ky, (ng))
122  SAFE_ALLOCATE(kz, (ng))
123  SAFE_ALLOCATE(isort, (ngq))
124  SAFE_ALLOCATE(isorti, (ngq))
125  SAFE_ALLOCATE(ekin, (ngq))
126  SAFE_ALLOCATE(eps, (nmtx,nmtx))
127
128#ifdef HDF5
129  if(use_hdf5) then
130
131    SAFE_ALLOCATE(old,(3,ngq))
132    call read_eps_old_gvecs_hdf5(ngq,old,infile)
133    kx(:)=old(1,:)
134    ky(:)=old(2,:)
135    kz(:)=old(3,:)
136
137    call read_eps_gvecsofq_hdf5(ngq,isort,isorti,ekin,1,infile)
138
139    call read_eps_matrix_ser_hdf5(eps, nmtx, 1, 1, infile)
140
141  else
142#endif
143
144    call open_file(unit=11,file=TRUNC(infile),form='unformatted',status='old')
145    read(11)
146    read(11)
147    read(11)
148    read(11)
149    read(11)
150    read(11)
151    read(11)
152    read(11)
153    read(11) ng,(kx(i),ky(i),kz(i),i=1,ng)
154    read(11) ngq,nmtx,(isort(ig),isorti(ig),ig=1,ngq)
155    read(11) (ekin(ig),ig=1,ngq)
156    read(11) (qk(j,1),j=1,3)
157    do jj = 1, nmtx
158      read(11) (eps(ii,jj),ii=1,nmtx)
159    enddo
160    call close_file(11)
161
162#ifdef HDF5
163  endif
164#endif
165
166  ! Since we want the q=0 dielectric function but we have the
167  ! q0<>0 but small dielectric, we can average over q0 and -q0
168  ! to get a better dielectric (linear terms in q0 will be canceled)
169
170! Calculate the maximum
171  gmax = 0
172  do ii = 1,ng
173    if (abs(kx(ii)) > gmax) gmax = abs(kx(ii))
174    if (abs(ky(ii)) > gmax) gmax = abs(ky(ii))
175    if (abs(kz(ii)) > gmax) gmax = abs(kz(ii))
176  enddo
177! Create a map
178!  write(6,*) gmax
179  SAFE_ALLOCATE(map, ((2*gmax+1)*(2*gmax+1)*(2*gmax+1)))
180  map = 0
181  do ii = 1, ngq
182    iout = isort(ii)
183    if (iout.eq.0) cycle
184    gx = kx(iout)
185    gy = ky(iout)
186    gz = kz(iout)
187    imap = ((gx+gmax)*(2*gmax+1)+gy+gmax)*(2*gmax+1)+gz+gmax+1
188    map(imap) = ii
189  enddo
190
191! For each g, find -g in the gvector list
192! and also mark which ii corresponds to G=0
193
194  SAFE_ALLOCATE(minusgidx, (nmtx))
195  minusgidx = 0
196  do ii=1,nmtx
197    iout = isort(ii)
198    gx = kx(iout)
199    gy = ky(iout)
200    gz = kz(iout)
201    imap = ((-gx+gmax)*(2*gmax+1)-gy+gmax)*(2*gmax+1)-gz+gmax+1
202    minusgidx(ii) = map(imap)
203    if (gx .eq. 0 .and. gy .eq. 0 .and. gz .eq. 0) ig0 = ii
204  enddo
205
206!  do ii = 1, min(100,nmtx)
207!    iout = isort(ii)
208!    gx = kx(iout)
209!    gy = ky(iout)
210!    gz = kz(iout)
211!    mgx = kx(isort(minusgidx(ii)))
212!    mgy = ky(isort(minusgidx(ii)))
213!    mgz = kz(isort(minusgidx(ii)))
214!    write(6,'(7i4)') ii, gx, gy , gz, mgx, mgy, mgz
215!  enddo
216
217! Set the wings to zero
218! This is as per Baldereschi and Tosatti, PRB 17, 4710 (1978)
219! This is perhaps not the correct thing to do. One should
220! still symmetrize epsilon - but the wings should come out
221! whatever they need to be automatically. What is given in that
222! paper by Baldereschi and Tosatti is at q=0 and the dielectric
223! function is weird at that point. What is needed in the GW
224! code is the average in the minibz...
225  !do ii=1,nmtx
226  !  do jj=1,nmtx
227  !    if (ii .eq. ig0 .and. jj .eq. ig0) cycle
228  !    if (ii .eq. ig0) eps(ii,jj) = 0.0d0
229  !    if (jj .eq. ig0) eps(ii,jj) = 0.0d0
230  !  enddo
231  !enddo
232
233! Copy eps(q0) into a temporary
234
235  SAFE_ALLOCATE(tempeps, (nmtx,nmtx))
236  tempeps = eps
237
238  errorMax=0D0
239  errorMaxTemp=0D0
240
241  write(6,'(5a4,3x,3a4,1x)',advance='no') "ig", "ig'", "Gx", "Gy", "Gz", "G'x", "G'y", "G'z"
242  ! handle the fact that the two components of the complex numbers will be written in the complex case
243#ifdef CPLX
244  write(6,'(6a16)') "Re eps(G,G')", "Im eps(G,G')", "Re eps(-G,-G')", "-Im eps(-G,-G')", "Re diff", "Im diff"
245#else
246  write(6,'(3a16)') "eps(G,G')", "eps(-G,-G')*", "difference"
247#endif
248
249!  do ii = 1, min(100,nmtx)
250!    do jj = 1, min(100,nmtx)
251  do ii = 1, nmtx
252    do jj = 1, nmtx
253      iout = isort(ii)
254      igx = kx(iout)
255      igy = ky(iout)
256      igz = kz(iout)
257      jout = isort(jj)
258      jgx = kx(jout)
259      jgy = ky(jout)
260      jgz = kz(jout)
261      kk = minusgidx(ii)
262      ll = minusgidx(jj)
263      if ((kk .le. nmtx) .and. (ll .le. nmtx)) then
264        errorMaxTemp=eps(ii,jj)-MYCONJG(eps(kk,ll))
265        if (abs(errorMaxTemp) .gt. errorMax) errorMax = abs(errorMaxTemp)
266        if (ii .lt. 20 .and. jj .lt. 20) then
267          write(6,'(5i4,3x,3i4,1x,6E16.5)') ii,jj,igx,igy,igz,jgx,jgy,jgz,eps(ii,jj),MYCONJG(eps(kk,ll)),errorMaxTemp
268        endif
269      endif
270    enddo
271  enddo
272
273  write(6,*) "Error = eps(G,G')-eps*(-G,-G')"
274  write(6,'("The max error in your matrix is",E12.5)') errorMax
275  write(6,*) "Symmetrizing the matrix"
276
277! Now add in contribution from -q0 which means conjg(eps(-g,-gp))
278
279  do ii=1,nmtx
280    do jj=1,nmtx
281
282      kk = minusgidx(ii)
283      ll = minusgidx(jj)
284      if ((kk .le. nmtx) .and. (ll .le. nmtx)) then
285        tempeps(ii,jj) = tempeps(ii,jj) + MYCONJG(eps(kk,ll))
286      endif
287    enddo
288  enddo
289! Average over q0 and -q0 and put back into eps
290
291  eps = 0.5d0*tempeps
292  SAFE_DEALLOCATE(tempeps)
293  SAFE_DEALLOCATE(minusgidx)
294
295#ifdef HDF5
296  if(use_hdf5) then
297    SAFE_DEALLOCATE(old)
298    call write_gvec_indices_hdf(ng,isort,isorti,ekin,1,outfile)
299
300    isize=SCALARSIZE
301
302    call write_matrix_ser_hdf(eps, nmtx, 1, 1, outfile)
303
304    SAFE_ALLOCATE(epsdiag,(isize,nmtx,1))
305
306    do ii = 1, nmtx
307      epsdiag(1,ii,1) = dble(eps(ii,ii))
308#ifdef CPLX
309      epsdiag(2,ii,1) = IMAG(eps(ii,ii))
310#endif
311    enddo
312
313    call write_matrix_diagonal_hdf(epsdiag, nmtx, 1, isize, outfile)
314
315    SAFE_DEALLOCATE(epsdiag)
316
317    call h5close_f(error)
318  else
319#endif
320
321    call open_file(unit=20,file=TRUNC(outfile),form='unformatted',status='replace')
322
323    ajname='chiGG0'
324    write(20) ajname,adate
325    write(20) freq_dep,nFreq
326    write(20) (qgrid(ii),ii=1,3)
327    if (freq_dep .eq. 2) then
328      write(20) (dFreqGrid(i),i=1,nFreq),(dFreqBrd(i),i=1,nFreq)
329    else
330      write(20)
331    endif
332    write(20)
333    write(20)
334    write(20) ecuts
335    write(20) nq,(q(j,1),j=1,3)
336    write(20) ng,(kx(ig),ky(ig),kz(ig),ig=1,ng)
337    write(20) ngq,nmtx,(isort(ig),isorti(ig),ig=1,ngq)
338    write(20) (ekin(ig),ig=1,ngq)
339    write(20) (qk(j,1),j=1,3)
340    do jj = 1, nmtx
341      write(20) (eps(ii,jj),ii=1,nmtx)
342    enddo
343
344    call close_file(20)
345
346#ifdef HDF5
347  endif
348#endif
349
350end program eps0sym
351