1!======================================================================= 2! 3! setup_bse_subsample Originally By DYQ Last Modified 02/09/2015 (DYQ) 4! 5! 6! Given a WFN_co and WFN files and a desired fine k-grid nk_fi_x nk_fi_y nk_fi_z, automatically find 7! k-point and q-points required for clustered subsampling. 8! Input files: 9! WFN_co: file used to generate BSE kernel on coarse grid 10! WFN: WFN file usedfor epsilon, usually the same as WFN_co 11! 12! Output files: 13! kpoints_sub_*.dat: Cluster of subsampled points around each coarse points. One file 14! is written for each coarse point on the full grid. Use to generate 15! WFN_sub wavefunctions. Run kernel for each WFN_sub. 16! epsilon_q0s.inp: Part of epsilon.inp file used to generate epsmat compatible with subsampled 17! k-points. 18! kpoints_wfnq.dat: Contains list of k-points needed to generate WFNq for all q points 19! kpoints_wfnq_*.dat: Separated lists of k points needed to generate WFNq for each q point in epsilon_q0s.inp 20! kpoints_wfnq.dat should contain combined list of all kpoints in all kpoints_wfnq_*.dat files 21! subsample.inp: Part of subsample.inp 22! 23! 24! TODO: 25! - Support HDF5 26! - Support subsampling along arbitrary directions (currently y and xy only) 27! 28!=============================================================================== 29 30#include "f_defs.h" 31 32program setup_bse_subsample 33 34 use global_m 35 use wfn_rho_vxc_io_m 36 use random_m 37 use irrbz_m 38 use fullbz_m 39 use checkbz_m 40 implicit none 41 42 integer :: narg, nk_fi(3),nk_sub, nk_old, nk_new,nk_fold,nk_tot,use_syms,nsub_fact, dir 43 character(len=256) :: fname_wfn_co,fname_wfn, fname_kpts,tmpstr 44 character(len=5) file_fmt 45 46 real(DP) :: dq(3), delta_co(3),delta_fi(3),qshift(3), delta_min 47 real(DP) :: sub_len 48 real(DP), allocatable :: kp_sub(:,:),kpts_new(:,:),kpts_all(:,:),qshifts(:,:) 49 integer :: ii,ik,iq,iq2,nk_sub_xy,nsub_max,ik_sub 50 type(mf_header_t) :: mf 51 logical :: skip_checkbz 52 integer, allocatable :: neq(:), indrk(:),neq_all(:) 53 type(grid) :: gr,gr_co 54 55 PUSH_SUB(setup_bse_subsample) 56 57 ! Input args 58 narg = iargc() 59 if (narg.lt.8 .or. narg.gt.12) then 60 write(0,*) 'Usage: setup_bse_subsample.x ASCII|BIN|HDF5 WFN_co WFN nk_fi_x nk_fi_y nk_fi_z nsub_factor direction' 61 write(0,*) ' [qshift_x qshift_y qshift_z use_syms]' 62 write(0,*) 63 write(0,*) 'Required arguments:' 64 write(0,*) ' ASCII|BIN|HDF5: format of the WFN file' 65 write(0,*) ' WFN_co: WFN_co file to used in the Kernel calculation' 66 write(0,*) ' WFN: WFN file used in the Epsilon/Sigma calculation (typically the same as WFN_co)' 67 write(0,*) ' nl_fi_x nk_fi_y nk_fi_z: the fine k-grid desired for absorption' 68 write(0,*) ' nsub_factor: number of subsampled points = nk_fi/nk_co*nsub_factor' 69 write(0,*) ' direction= 1 (b1-axis), 2 (b2-axis), 3 (diagonal)' 70 write(0,*) 71 write(0,*) 'Optional arguments:' 72 write(0,*) ' qshift_x qshift_y qshift_z: shift applied to subsampled kpoints' 73 write(0,*) ' use_syms = 1 (default), to unfold WFN_co using symmetries' 74 write(0,*) ' = 0 do not unfold WFN_co using symmetries' 75 write(0,*) 76 stop 77 endif 78 call getarg(1, file_fmt) 79 call getarg(2, fname_wfn_co) 80 call getarg(3, fname_wfn) 81 call getarg(4, tmpstr) 82 read(tmpstr,*) nk_fi(1) 83 call getarg(5, tmpstr) 84 read(tmpstr,*) nk_fi(2) 85 call getarg(6, tmpstr) 86 read(tmpstr,*) nk_fi(3) 87 call getarg(7, tmpstr) 88 read(tmpstr,*) nsub_fact 89 call getarg(8, tmpstr) 90 read(tmpstr,*) dir 91 qshift(:) = 0.d0 92 use_syms=1 93 if (narg.eq.9) then 94 call getarg(9, tmpstr) 95 read(tmpstr,*) use_syms 96 elseif (narg.gt.9) then 97 call getarg(9, tmpstr) 98 read(tmpstr,*) qshift(1) 99 call getarg(10, tmpstr) 100 read(tmpstr,*) qshift(2) 101 call getarg(11, tmpstr) 102 read(tmpstr,*) qshift(3) 103 if (narg.eq.12) then 104 call getarg(12, tmpstr) 105 read(tmpstr,*) use_syms 106 endif 107 write(6,*) "qshift = ",qshift 108 endif 109 110 ! Read header of WFN_co 111 if (file_fmt=='ASCII') then 112 call open_file(unit=11, file=TRIM(fname_wfn_co), form='formatted', status='old') 113 call read_mf_header(11, mf) 114 call close_file(11) 115 elseif (file_fmt=='BIN ') then 116 call open_file(unit=11, file=TRIM(fname_wfn_co), form='unformatted', status='old') 117 call read_mf_header(11, mf) 118 call close_file(11) 119 elseif (file_fmt=='HDF5 ') then 120 call die('HDF5 not implemented yet') 121 else 122 call die('Unknown format "'//TRIM(file_fmt)//'". Must be either ASCII, BIN or HDF5.') 123 endif 124 write(6,'(a,3(1x,i0))') 'Read WFN_co k-grid:', mf%kp%kgrid 125 126 ! Determine cluster size and sampling fineness from k-grid 127 do ii=1,3 128 delta_co(ii) = 1.0/mf%kp%kgrid(ii) 129 delta_fi(ii) = 1.0/nk_fi(ii) / nsub_fact 130 enddo 131 delta_min=minval(delta_co)-minval(delta_fi) 132 nsub_max = minval(nk_fi)/delta_min * nsub_fact * 2. 133 134 ! Unfold WFN_co 135 gr_co%nr = mf%kp%nrk 136 SAFE_ALLOCATE(gr_co%r, (3, gr_co%nr)) 137 gr_co%r = mf%kp%rk 138 if (use_syms.eq.1) then 139 call fullbz(mf%crys, mf%syms, gr_co, mf%syms%ntran, skip_checkbz, wigner_seitz=.false., paranoid=.true.) 140 else 141 call fullbz(mf%crys, mf%syms, gr_co, 1, skip_checkbz, wigner_seitz=.false., paranoid=.true.) 142 endif 143 if (.not.skip_checkbz) call checkbz(gr_co%nf, gr_co%f, mf%kp%kgrid, mf%kp%shift, & 144 mf%crys%bdot, TRUNC(fname_wfn), 'k', .false., .false., .false.) 145 146 ! Print out list of cluster points 147 SAFE_ALLOCATE(kp_sub,(3,nsub_max)) 148 SAFE_ALLOCATE(qshifts,(3,nsub_max)) 149 do ik=1,gr_co%nf 150 kp_sub = 0.0 151 qshifts = 0.0 152 sub_len=0.0 153 nk_sub=0 154 do while (sub_len<delta_min) 155 nk_sub=nk_sub+1 156 kp_sub(:,nk_sub) = gr_co%f(:,ik) 157 if (dir.eq.3) then 158 kp_sub(2,nk_sub) = kp_sub(2,nk_sub) + delta_fi(2)*nk_sub 159 kp_sub(1,nk_sub) = kp_sub(1,nk_sub) + delta_fi(1)*nk_sub 160 qshifts(2,nk_sub) = delta_fi(2)*nk_sub 161 qshifts(1,nk_sub) = delta_fi(1)*nk_sub 162 else 163 kp_sub(dir,nk_sub) = kp_sub(dir,nk_sub) + delta_fi(dir)*nk_sub 164 qshifts(dir,nk_sub) = delta_fi(dir)*nk_sub 165 endif 166 sub_len=sub_len+minval(delta_fi) 167 enddo 168 write(fname_kpts,'(a,i4.4,a)') 'kpoints_sub_', ik, '.dat' 169 call open_file(unit=12,file=fname_kpts,form='formatted',status='replace') 170 write(12,'(a)') 'K_POINTS crystal' 171 write(12,'(i8)') nk_sub+1 172 write(12, '(3(f13.9),f6.1)') gr_co%f(:,ik) + qshift(:), 1.0 173 do ik_sub=1,nk_sub 174 write(12, '(3(f13.9),f6.1)') kp_sub(:,ik_sub) + qshift(:), 1.0 175 enddo 176 call close_file(12) 177 enddo ! loop over ik_co 178 write(6,*) "Number of subsampled points per coarse point= ", nk_sub 179 180 call open_file(unit=15,file='subsample.inp',form='formatted',status='replace') 181 write(15,'(i8)') nk_sub+1 182 write(15,'(i8)') gr_co%nf 183 do ik=1,gr_co%nf 184 write(15, '(3(f13.9))') gr_co%f(:,ik) 185 enddo 186 write(15,'(a)') "! List subsampled bsemat.h5 files here" 187 write(15,'(a)') "! List subsampled WFN files here" 188 write(15,'(a)') "! If using finite Q, list subsampled WFNq files here. Otherwise delete this line." 189 call close_file(15) 190 191 ! Read header for WFN file 192 if (file_fmt=='ASCII') then 193 call open_file(unit=11, file=TRIM(fname_wfn), form='formatted', status='old') 194 call read_mf_header(11, mf) 195 call close_file(11) 196 elseif (file_fmt=='BIN ') then 197 call open_file(unit=11, file=TRIM(fname_wfn), form='unformatted', status='old') 198 call read_mf_header(11, mf) 199 call close_file(11) 200 elseif (file_fmt=='HDF5 ') then 201 call die('HDF5 not implemented yet') 202 else 203 call die('Unknown format "'//TRIM(file_fmt)//'". Must be either ASCII, BIN or HDF5.') 204 endif 205 write(6,'(a,3(1x,i0))') 'Read WFN k-grid:', mf%kp%kgrid 206 207 208 ! Find q-points 209 call open_file(unit=13,file='epsilon_q0s.inp',form='formatted',status='replace') 210 call open_file(unit=14,file='kpoints_wfnq.dat',form='formatted',status='replace') 211 write(13,*) 212 write(13,'(a)') 'subsample' 213 write(13,'(a,i0)') 'number_qpoints ', nk_sub*2 214 write(13,'(a)') 'begin qpoints' 215 216 ! Find k-points for WFNq 217 gr%nr = mf%kp%nrk 218 SAFE_ALLOCATE(gr%r, (3, gr%nr)) 219 gr%r = mf%kp%rk 220 call fullbz(mf%crys, mf%syms, gr, mf%syms%ntran, skip_checkbz, wigner_seitz=.false., paranoid=.true.) 221 if (.not.skip_checkbz) call checkbz(gr%nf, gr%f, mf%kp%kgrid, mf%kp%shift, & 222 mf%crys%bdot, TRUNC(fname_wfn), 'k', .false., .false., .false.) 223 224 nk_new = gr%nf*(nk_sub)*(nk_sub+1) 225 SAFE_ALLOCATE(kpts_new, (3,gr%nf)) 226 SAFE_ALLOCATE(indrk, (gr%nf)) 227 SAFE_ALLOCATE(neq, (gr%nf)) 228 SAFE_ALLOCATE(kpts_all, (3,nk_new*2)) 229 SAFE_ALLOCATE(neq_all, (nk_new*2)) 230 nk_new = 0 231 nk_old = 0 232 nk_tot = 0 233 dq = 0 234 write(6,*) "number of subsampled points: ", nk_sub 235 do iq=1,nk_sub 236 if (dir.eq.3) then 237 dq(1) = delta_fi(1)*iq 238 dq(2) = delta_fi(2)*iq 239 else 240 dq(dir) = delta_fi(dir)*iq 241 endif 242 write(13,'(3(f13.9),a)') dq(:), ' 1 1' 243 244 ! Find subgroup that leaves dq invariant 245 call subgrp(dq,mf%syms) 246 forall(ik=1:gr%nf) kpts_new(1:3,ik) = gr%f(1:3, ik) + dq 247 nk_old = nk_new+1 248 nk_new = nk_new + gr%nf 249 ! Fold shifted points to irreducible wedge using subgroup symmetry 250 call irrbz(mf%syms, gr%nf, kpts_new, nk_fold, neq, indrk) 251 kpts_all(1:3,nk_tot+1:nk_tot+nk_fold) = kpts_new(1:3,indrk(1:nk_fold)) 252 neq_all(nk_tot+1:nk_tot+nk_fold) = neq(1:nk_fold) 253 nk_tot = nk_tot + nk_fold 254 255 ! write a list of kpoints corresponding to each q-point 256 write(fname_kpts,'(a,i4.4,a)') 'kpoints_wfnq_', iq*2-1, '.dat' 257 call open_file(unit=17,file=fname_kpts,form='formatted',status='replace') 258 write(17,'(a)') 'K_POINTS crystal' 259 write(17,'(i0)') nk_fold 260 do ik=1,nk_fold 261 write(17, '(3(f13.9),f6.1)') kpts_new(:,indrk(ik)), dble(neq(ik)) 262 enddo 263 call close_file(17) 264 265 ! Time-reversed q point 266 if (dir.eq.3) THEN 267 dq(1) = 1.0d0-delta_fi(1)*iq 268 dq(2) = 1.0d0-delta_fi(2)*iq 269 else 270 dq(dir) = 1.0d0-delta_fi(dir)*iq 271 endif 272 write(13,'(3(f13.9),a)') dq(:), ' 1 1' 273 ! Find subgroup that leaves dq invariant 274 call subgrp(dq,mf%syms) 275 forall(ik=1:gr%nf) kpts_new(1:3,ik) = gr%f(1:3, ik) + dq 276 nk_old = nk_new+1 277 nk_new = nk_new + gr%nf 278 ! Fold shifted points to irreducible wedge using subgroup symmetry 279 call irrbz(mf%syms, gr%nf, kpts_new, nk_fold, neq, indrk) 280 kpts_all(1:3,nk_tot+1:nk_tot+nk_fold) = kpts_new(1:3,indrk(1:nk_fold)) 281 neq_all(nk_tot+1:nk_tot+nk_fold) = neq(1:nk_fold) 282 nk_tot = nk_tot + nk_fold 283 284 write(fname_kpts,'(a,i4.4,a)') 'kpoints_wfnq_', iq*2, '.dat' 285 call open_file(unit=17,file=fname_kpts,form='formatted',status='replace') 286 write(17,'(a)') 'K_POINTS crystal' 287 write(17,'(i0)') nk_fold 288 do ik=1,nk_fold 289 write(17, '(3(f13.9),f6.1)') kpts_new(:,indrk(ik)), dble(neq(ik)) 290 enddo 291 call close_file(17) 292 293 enddo 294 write(13,'(a)') 'end' 295 write(13,*) 296 call close_file(13) 297 298 write(6,*) "Total kpoints in WFNq: ",nk_tot 299 write(14,'(a)') 'K_POINTS crystal' 300 write(14,'(i0)') nk_tot 301 do ik=1,nk_tot 302 write(14, '(3(f13.9),f6.1)') kpts_all(:,ik), dble(neq_all(ik)) 303 enddo 304 call close_file(14) 305 306 SAFE_DEALLOCATE(kpts_new) 307 SAFE_DEALLOCATE(kpts_all) 308 SAFE_DEALLOCATE(neq) 309 SAFE_DEALLOCATE(neq_all) 310 SAFE_DEALLOCATE(indrk) 311 SAFE_DEALLOCATE(kp_sub) 312 SAFE_DEALLOCATE(qshifts) 313 314 POP_SUB(setup_bse_subsample) 315 316end program setup_bse_subsample 317