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