TABLE OF CONTENTS
ABINIT/mrgscr [ Programs ]
NAME
mrgscr
FUNCTION
This code reads partial (SCR|SUSC) files for different q points creating a single file that can be used to perform a sigma calculation.
COPYRIGHT
Copyright (C) 2005-2024 ABINIT group (RS, MG, MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
NOTES
If the number of SCR files to be merged is equal to 1, the program checks the integrity of the file reporting the list of missing q-points. Note that the list of required q-points depends on the k-mesh used during the calculation of the KSS file. We assume indeed that the same k-mesh is used during the calculation of the matrix elements of sigma.
INPUTS
(Main program)
OUTPUT
Only checking and writing
SOURCE
31 #if defined HAVE_CONFIG_H 32 #include "config.h" 33 #endif 34 35 #include "abi_common.h" 36 37 program mrgscr 38 39 use defs_basis 40 use m_xmpi 41 use m_abicore 42 use m_errors 43 use m_nctk 44 use netcdf 45 use m_hdr 46 use m_crystal 47 use m_pawrhoij 48 use m_dtset 49 50 use defs_abitypes, only : MPI_type 51 use m_build_info, only : abinit_version 52 use m_specialmsg, only : herald 53 use m_time, only : timein 54 use m_gwdefs, only : GW_TOLQ, GW_TOLQ0, GW_Q0_DEFAULT 55 use m_io_tools, only : prompt, flush_unit, open_file 56 use m_fstrings, only : int2char4, endswith, itoa, sjoin 57 use m_fft_mesh, only : g2ifft 58 use m_fftcore, only : get_cache_kb, getng 59 use m_fft, only : fourdp 60 use m_numeric_tools, only : iseven, cspint 61 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 62 use m_geometry, only : normv, metric 63 use m_gsphere, only : gsphere_t 64 use m_bz_mesh, only : kmesh_t, find_qmesh 65 use m_vcoul, only : vcoul_t 66 use m_ioarr, only : read_rhor 67 use m_io_screening, only : read_screening, hscr_t, ioscr_qmerge, ioscr_qrecover, ioscr_wmerge, ioscr_wremove 68 use m_ppmodel, only : ppm_init, ppm_free, setup_ppmodel, getem1_from_PPm_one_ggp, & 69 get_PPm_eigenvalues, ppmodel_t, cqratio 70 use m_model_screening, only : remove_phase 71 use m_screening, only : mkdump_er, em1results_free, em1results_print, decompose_epsm1, & 72 init_er_from_file, Epsilonm1_results 73 use m_wfd, only : test_charge 74 75 implicit none 76 77 !Local variables------------------------------- 78 !scalars 79 integer,parameter :: master=0,paral_kgb0=0,rdwr2=2,prtvol=0,cplex1=1 80 integer :: iomode,fform1,ifile,ierr,ii,ios,iqibz,iqf,nfiles 81 integer :: unt_dump,idx,ig1,ig2,iomega,ppmodel,npwe_asked,mqmem,io,unt_dump2 82 integer :: id_required,ikxc,approx_type,option_test,dim_kxcg,usexcnhat,usefinegrid 83 integer :: mgfft,nqlwl,nfft,igmax,comm,nq_selected,kptopt 84 integer :: choice,nfreq_tot,nfreqre,nfreqim,nfreqc,ifrq,imax 85 integer :: ig1_start,ig1_end,ig2_start,ig2_end,gmgp_idx,orig_npwe 86 real(dp) :: ucvol,boxcutmin,ecut,drude_plsmf,compch_fft,compch_sph 87 real(dp) :: nelectron_exp,freqremax,eps_diff,eps_norm,eps_ppm_norm 88 real(dp) :: value1,value2,factor,GN_drude_plsmf 89 real(dp) :: tcpu,tcpui,twall,twalli 90 real(gwp) :: phase 91 logical :: is_sus,is_scr,same_freqs,calc_epsilon,only_diag 92 character(len=1) :: ans 93 character(len=10) :: tagq 94 character(len=24) :: codename 95 character(len=500) :: msg 96 character(len=nctk_slen) :: varname 97 character(len=fnlen) :: fname_out,fname,fname_dump,fname_rho,prefix,fname_eigen,fname_dump2 98 type(hdr_type) :: hdr_rhor 99 type(abifile_t) :: abifile 100 type(hscr_t),pointer :: Hscr0 101 type(hscr_t),target :: Hscr_merge 102 type(MPI_type) :: MPI_enreg 103 type(kmesh_t) :: Kmesh,Qmesh 104 type(crystal_t) :: Cryst 105 type(gsphere_t) :: Gsphere 106 type(ppmodel_t) :: PPm 107 type(Epsilonm1_results) :: Er 108 type(vcoul_t), target :: Vcp 109 type(Dataset_type) :: Dtset 110 !arrays 111 integer :: ngfft(18) 112 integer,allocatable :: foundq(:),freq_indx(:,:) 113 real(dp),parameter :: k0(3) = [zero,zero,zero] 114 real(dp) :: gmet(3,3),gprimd(3,3),qdiff(3),rmet(3,3),qtmp(3),tsec(2) 115 real(dp),allocatable :: qlwl(:,:),real_omega(:),rhor(:,:),rhog(:,:),nhat(:,:) 116 real(dp),allocatable :: work(:),ftab(:),ysp(:,:),eint(:),qratio(:,:) 117 complex(gwpc),pointer :: vc_sqrt(:) 118 complex(gwpc),allocatable :: epsm1(:,:,:,:),kxcg(:,:) 119 complex(dpc),allocatable :: omega(:),em1_ppm(:),epsm1_eigen(:,:),ppm_eigen(:,:),rhoggp(:,:) 120 character(len=fnlen),allocatable :: filenames(:) 121 type(pawrhoij_type),allocatable :: pawrhoij(:) 122 type(hscr_t),target,allocatable :: Hscr_file(:) ! Cannot use allocatable as pathscale5 miscompiles the code 123 124 ! ************************************************************************* 125 126 ! Change communicator for I/O (mandatory!) 127 call abi_io_redirect(new_io_comm=xmpi_world) 128 129 ! Initialize MPI 130 call xmpi_init() 131 132 ! Initialize memory profiling if it is activated 133 ! if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0" 134 ! note that abimem.mocc files can easily be multiple GB in size so don't use this option normally 135 #ifdef HAVE_MEM_PROFILING 136 call abimem_init(0) 137 #endif 138 139 call timein(tcpui,twalli) 140 141 ! Default for sequential use 142 call initmpi_seq(MPI_enreg); comm = MPI_enreg%comm_world 143 144 is_sus=.FALSE.; is_scr=.FALSE. 145 146 ! Write greetings, and read the number of files === 147 codename='MRGSCR'//REPEAT(' ',18) 148 call herald(codename,abinit_version,std_out) 149 150 call prompt(' Enter the number of files to merge: ',nfiles) 151 ABI_CHECK(nfiles > 0, 'nfiles must be >0') 152 153 ABI_MALLOC(filenames,(nfiles)) 154 ABI_MALLOC(Hscr_file,(nfiles)) 155 156 if (nfiles == 1) then 157 call prompt(' Enter the name of the file to be analyzed: ',filenames(1)) 158 write(msg,'(7a)')ch10,& 159 ' Running single-file mode:',ch10,& 160 ' Checking the integrity of file: ',TRIM(filenames(1)),ch10,& 161 ' reporting the list of q-points that are missing. ' 162 call wrtout(std_out, msg) 163 164 if (nctk_try_fort_or_ncfile(filenames(1), msg) /= 0) then 165 ABI_ERROR(msg) 166 end if 167 168 else if (nfiles > 1) then 169 ! Read name of files to be merged and check for existence. 170 call prompt(' Enter the prefix for the final output file: ',fname_out) 171 172 do ifile=1,nfiles 173 write(msg,'(a,i4)')' Enter the name for the partial screening file no.',ifile 174 call prompt(msg,filenames(ifile)) 175 176 if (nctk_try_fort_or_ncfile(filenames(ifile), msg) /= 0) then 177 ABI_ERROR(msg) 178 end if 179 180 end do 181 end if 182 183 ! Read the header of each file. 184 do ifile=1,nfiles 185 iomode = IO_MODE_FORTRAN; if (endswith(filenames(ifile), ".nc")) iomode = IO_MODE_ETSF 186 187 call Hscr_file(ifile)%from_file(filenames(ifile), fform1, comm) 188 ABI_CHECK(fform1 /= 0, sjoin("fform == 0 in", filenames(ifile))) 189 190 abifile = abifile_from_fform(fform1) 191 if (abifile%fform == 0) then 192 ABI_ERROR(sjoin("Cannot find any abifile object associated to fform1:", itoa(fform1))) 193 end if 194 if (abifile%class /= "polariz" .and. abifile%class /= "epsm1") then 195 ABI_ERROR(sjoin('Error while reading header, fform= ',itoa(fform1))) 196 end if 197 is_scr = abifile%class == "epsm1" 198 is_sus = abifile%class == "polariz" 199 200 call Hscr_file(ifile)%print(unit=std_out,prtvol=1) 201 202 if (ifile == 1) call metric(gmet,gprimd,-1,rmet,Hscr_file(ifile)%Hdr%rprimd,ucvol) 203 end do !ifile 204 205 if (nfiles > 1) then 206 ! Put the correct ending on the output file 207 if (is_scr) fname_out=TRIM(fname_out)//'_SCR' 208 if (is_sus) fname_out=TRIM(fname_out)//'_SUS' 209 end if 210 211 ! Produce output file in netcdf format we are merging netcdf files. 212 if (iomode == IO_MODE_ETSF .and. .not. endswith(fname_out, ".nc")) fname_out = nctk_ncify(fname_out) 213 214 !============================ 215 !=== Merge multiple files === 216 !============================ 217 if (nfiles > 1) then 218 219 ! Check what kind of merging is to be performed 220 write(std_out,'(2(a))') ch10,' Do you want to merge q-points (= 1) ?' 221 write(std_out,'(a)') ' or do you want to merge frequencies (= 2) ?' 222 read(std_in,*)choice 223 224 select case(choice) 225 case(1) 226 write(std_out,'(3a)') ch10,' 1 => merging q-points',ch10 227 call ioscr_qmerge(nfiles, filenames, hscr_file, fname_out, hscr_merge) 228 229 case (2) 230 ! Merge frequencies 231 write(std_out,'(3a)') ch10,' 2 => merging frequency grids',ch10 232 !ABI_WARNING("Advanced user option, consistency in fform etc. will not be checked.") 233 234 write(std_out,'(2a)') ch10,' Enter freqremax [eV] for the merged file (Enter 0 to use all freq. found):' 235 read(std_in,*)freqremax 236 237 freqremax = freqremax/Ha_eV; if (freqremax<tol16) freqremax = HUGE(freqremax) 238 call ioscr_wmerge(nfiles, filenames, hscr_file, freqremax, fname_out, hscr_merge) 239 240 case default 241 ABI_ERROR("Invalid choice!") 242 end select 243 244 end if ! nfiles>1 245 246 ! Now check if the list of q-points is complete 247 ! Here we assume that the k-mesh reported in the header is the same as that used during the sigma calculation. 248 write(msg,'(3a)') ch10,' Checking if the list of q-points is complete. ',ch10 249 call wrtout(std_out, msg) 250 251 !call hscr_check_qpoints(hscr0) 252 253 Hscr0 => Hscr_file(1) 254 fname =filenames(1) 255 if (nfiles>1) then 256 Hscr0 => Hscr_merge 257 fname = fname_out 258 end if 259 260 cryst = HScr0%Hdr%get_crystal(remove_inv=.FALSE.) 261 262 kptopt=1 263 call Kmesh%init(Cryst,HScr0%Hdr%nkpt,Hscr0%Hdr%kptns,kptopt) 264 call Kmesh%print("K-mesh for the wavefunctions",prtvol=prtvol) 265 266 call find_qmesh(Qmesh,Cryst,Kmesh) 267 call qmesh%print("Q-mesh for the screening function",prtvol=prtvol) 268 269 ABI_MALLOC(foundq,(Qmesh%nibz)) 270 foundq(:)=0 271 do iqibz=1,Qmesh%nibz 272 do iqf=1,Hscr0%nqibz 273 qdiff(:)=Qmesh%ibz(:,iqibz)-Hscr0%qibz(:,iqf) 274 if (normv(qdiff,gmet,'G')<GW_TOLQ) foundq(iqibz)=foundq(iqibz)+1 275 end do 276 end do 277 278 if (ANY(foundq==0)) then 279 write(msg,'(6a)')ch10,& 280 ' File ',TRIM(fname),' is not complete ',ch10,& 281 ' The following q-points are missing:' 282 call wrtout(std_out, msg) 283 ii=0 284 do iqibz=1,Qmesh%nibz 285 if (foundq(iqibz)==0) then 286 ii=ii+1 287 write(msg,'(i3,a,3f12.6)')ii,') ',Qmesh%ibz(:,iqibz) 288 call wrtout(std_out, msg) 289 end if 290 end do 291 end if 292 293 if (ANY(foundq>1)) then 294 write(msg,'(6a)')ch10,& 295 ' File ',TRIM(fname),' is overcomplete ',ch10,& 296 ' The following q-points are present more than once:' 297 call wrtout(std_out, msg) 298 ii=0 299 do iqibz=1,Qmesh%nibz 300 if (foundq(iqibz)>1) then 301 ii=ii+1 302 write(msg,'(i3,a,3f12.6)')ii,') ',Qmesh%ibz(:,iqibz) 303 call wrtout(std_out, msg) 304 end if 305 end do 306 end if 307 308 if (ALL(foundq==1)) then 309 write(msg,'(5a)')ch10,& 310 '.File ',TRIM(fname),' contains a complete list of q-points ',ch10 311 call wrtout(std_out, msg) 312 end if 313 314 !===================== 315 !=== Recovery mode === 316 !===================== 317 if (nfiles==1) then 318 319 write(std_out,'(2(a))') ch10,' Do you want to recover a subset of q-points (= 1) ?' 320 write(std_out,'(a)') ' or extract the contents of the file (= 2) ?' 321 write(std_out,'(a)') ' or create dielectric function (SCR file)' 322 write(std_out,'(a)') ' and/or extract plasmon-pole parameters (= 3) ?' 323 write(std_out,'(a)') ' or remove real frequencies (= 4) ?' 324 write(std_out,'(a)') ' or remove imaginary frequencies (= 5) ?' 325 write(std_out,'(a)') ' or calculate a model screening (= 6) ?' 326 !write(std_out,'(a)') ' or interpolate the screening in k-space (= 8) ?' 327 !write(std_out,'(a)') ' or convert a netcd file to Fortran (= 9) ?' 328 read(std_in,*)choice 329 330 select case(choice) 331 case(1) 332 ! Recover subset of q-points -------------------------------------------------- 333 write(std_out,'(a)') ' 1 => Recovering subset of q-points' 334 call prompt(' Enter the number of q-points to be extracted: ',nq_selected) 335 call prompt(' Enter the name of the final output file: ',fname_out) 336 337 if (endswith(filenames(1), ".nc") .and. .not. endswith(fname_out, ".nc")) then 338 fname_out = nctk_ncify(fname_out) 339 call wrtout(std_out,"- Added .nc extension to output file as input data is in netcdf format.") 340 end if 341 342 call ioscr_qrecover(filenames(1), nq_selected, fname_out) 343 344 case(2) 345 ! Analyse file ---------------------------------------------------------------- 346 ABI_CHECK(iomode==IO_MODE_FORTRAN, "netcdf output not coded") 347 348 write(std_out,'(a)') ' 2 => Extraction of file contents' 349 350 ! Initialize the G-sphere. 351 call Gsphere%init(Cryst,Hscr0%npwe,gvec=Hscr0%gvec) 352 353 ABI_MALLOC_OR_DIE(epsm1,(Hscr0%npwe,Hscr0%npwe,Hscr0%nomega,1), ierr) 354 355 ! Give option to output epsilon instead of chi0 356 calc_epsilon = .FALSE. 357 if (is_sus) then 358 write(std_out,'(2a)') ch10,' You have provided a chi_0 file for analysis. Would you like to output' 359 write(std_out,'(2a)',advance='no') ' the dielectric function epsilon_GG'' ', '= delta_GG'' - v_G*chi0_GG''[Y/N] ? ' 360 read(std_in,*)ans 361 362 if (ans=='Y'.or.ans=='y') then 363 ! Initialise Coulomb terms 364 if (Er%Hscr%nqlwl==0) then 365 nqlwl=1 366 ABI_MALLOC(qlwl,(3,nqlwl)) 367 qlwl(:,1)= GW_Q0_DEFAULT 368 else 369 nqlwl=Er%Hscr%nqlwl 370 ABI_MALLOC(qlwl,(3,nqlwl)) 371 qlwl(:,:)=Er%Hscr%qlwl(:,1:nqlwl) 372 end if 373 374 Dtset%gw_icutcoul=3; Dtset%rcut=zero 375 Dtset%vcutgeo=(/zero,zero,zero/); 376 Dtset%boxcenter=(/zero,zero,zero/) 377 378 write(std_out,'(2a)',advance='no') ch10,' Was a Coulomb cutoff technique used [Y/N] ? ' 379 read(std_in,*)ans 380 if (ans=='Y'.or.ans=='y') then 381 write(std_out,'(2a)',advance='no') ' Enter gw_icutcoul: ' 382 read(std_in,*)Dtset%gw_icutcoul 383 write(std_out,'(2a)',advance='no') ' Enter vcutgeo: ' 384 read(std_in,*)Dtset%vcutgeo 385 write(std_out,'(2a)',advance='no') ' Enter boxcenter: ' 386 read(std_in,*)Dtset%boxcenter 387 end if 388 dtset%ecutsigx = -one 389 390 call Vcp%init(Gsphere,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,& 391 Dtset%vcutgeo,Dtset%ecutsigx,Hscr0%npwe,nqlwl,qlwl,comm) 392 ABI_FREE(qlwl) 393 394 calc_epsilon = .TRUE. 395 end if 396 end if 397 398 ig1 = 0; ig2 = 0 399 write(std_out,'(2(a),I0,a)',advance='NO') ch10,' Enter the starting index for G (1 - ',Hscr0%npwe,' ): ' 400 read(std_in,*)ig1_start 401 if (ig1_start<1.OR.ig1_start>Hscr0%npwe) then 402 ABI_ERROR(' Starting index out of bounds') 403 end if 404 write(std_out,'(a,I0,a,I0,a)',advance='NO') ' Enter the ending index for G ( ',ig1_start,' - ',Hscr0%npwe,' ): ' 405 read(std_in,*)ig1_end 406 if (ig1_end<ig1_start.OR.ig1_end>Hscr0%npwe) then 407 ABI_ERROR(' Ending index out of bounds') 408 end if 409 write(std_out,'(a,I0,a)',advance='NO') ' Enter the starting index for G'' (1 - ',Hscr0%npwe,' ): ' 410 read(std_in,*)ig2_start 411 if (ig2_start<1.OR.ig2_start>Hscr0%npwe) then 412 ABI_ERROR(' Starting index out of bounds') 413 end if 414 write(std_out,'(a,I0,a,I0,a)',advance='NO') ' Enter the ending index for G'' ( ',ig2_start,' - ',Hscr0%npwe,' ): ' 415 read(std_in,*)ig2_end 416 if (ig2_end<ig2_start.OR.ig2_end>Hscr0%npwe) then 417 ABI_ERROR(' Ending index out of bounds') 418 end if 419 420 only_diag = .FALSE. 421 write(std_out,'(a)',advance='no') ' Would you like to output only the diagonal [Y/N] ? ' 422 read(std_in,*)ans 423 if (ans=='Y'.or.ans=='y') only_diag = .TRUE. 424 425 do iqibz=1,Hscr0%nqibz 426 ! In the long wavelength limit we set q==0, because we still can use symmetries for the Body. 427 qtmp(:)=Hscr0%qibz(:,iqibz); if (normv(qtmp,Cryst%gmet,'G')<GW_TOLQ0) qtmp(:)=zero 428 429 ! FIXME 430 varname = "none" 431 call read_screening(varname,fname,Hscr0%npwe,1,Hscr0%nomega,epsm1,iomode,comm,iqiA=iqibz) 432 433 if (calc_epsilon) then ! Calculate epsilon 434 do iomega=1,Hscr0%nomega 435 if (iqibz==1) then 436 if (nqlwl>1) then 437 ABI_ERROR('nqlwl>1 not coded yet!') 438 end if 439 vc_sqrt => Vcp%vcqlwl_sqrt(:,iqibz) ! Use Coulomb term for q-->0 440 else 441 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 442 end if 443 do ig2=ig2_start,ig2_end 444 do ig1=ig1_start,ig1_end 445 epsm1(ig1,ig2,iomega,1) = -(vc_sqrt(ig1)**2)*epsm1(ig1,ig2,iomega,1) 446 end do ! ig1 447 epsm1(ig2,ig2,iomega,1) = one + epsm1(ig2,ig2,iomega,1) 448 end do ! ig2 449 end do ! iomega 450 end if ! Do we calculate epsilon 451 452 ! Find out the total number of frequencies along real/imaginary axes 453 ! and possibly in the z-plane 454 nfreqre=0; nfreqim=0; nfreqc=0; 455 do iomega=1,Hscr0%nomega 456 if (ABS(REAL(Hscr0%omega(iomega)))<tol8.AND. ABS(AIMAG(Hscr0%omega(iomega)))<tol8) nfreqre = nfreqre + 1 457 if (ABS(REAL(Hscr0%omega(iomega)))>tol8.AND. ABS(AIMAG(Hscr0%omega(iomega)))<tol8) nfreqre = nfreqre + 1 458 if (ABS(REAL(Hscr0%omega(iomega)))<tol8.AND. ABS(AIMAG(Hscr0%omega(iomega)))>tol8) nfreqim = nfreqim + 1 459 end do 460 if (Hscr0%nomega-nfreqre-nfreqim/=0) then 461 write(std_out,'(/,a)') ' WARNING: There are frequencies in the full complex plane.' 462 write(std_out,'(a)') ' The _SCR or _SUS file might not be suitable' 463 write(std_out,'(a,/)') ' for self-energy calculations.' 464 nfreqc = Hscr0%nomega-nfreqre-nfreqim 465 end if 466 write(std_out,'(2a,I0,a)') ch10,' Found ',Hscr0%nomega,' frequencies.' 467 write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10 468 if (nfreqc>0) then 469 write(std_out,'(a,I0)') ' There is a grid in the complex plane with ',nfreqc 470 write(std_out,'(2a)') ' extra frequencies in the list.',ch10 471 end if 472 473 ! Get Q index for name 474 call int2char4(iqibz,tagq) 475 ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!') 476 477 if (nfreqre>0) then ! Output real frequency axis 478 if (calc_epsilon) then 479 fname_dump=TRIM(fname)//'_EPS_Q'//TRIM(tagq) 480 else 481 fname_dump=TRIM(fname)//'_Q'//TRIM(tagq) 482 end if 483 484 if (open_file(fname_dump, msg, newunit=unt_dump, status='replace', form='formatted') /= 0) then 485 ABI_ERROR(msg) 486 end if 487 488 do ig1=ig1_start,ig1_end 489 do ig2=ig2_start,ig2_end 490 if (only_diag.AND.ig1/=ig2) CYCLE 491 write(unt_dump,'(2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,/)')& 492 '# ig1= ',ig1,' ig2= ',ig2,& 493 '# q = ',Hscr0%qibz(:,iqibz),& 494 '# G = ',Hscr0%gvec(:,ig1),' G''= ',Hscr0%gvec(:,ig2),& 495 '# omega [eV] Re Im ' 496 do iomega=1,nfreqre 497 write(unt_dump,'(f8.2,4x,2es16.8)') REAL(Hscr0%omega(iomega))*Ha_eV,& 498 REAL(epsm1(ig1,ig2,iomega,1)),AIMAG(epsm1(ig1,ig2,iomega,1)) 499 end do 500 write(unt_dump,*) 501 write(unt_dump,*) 502 end do !ig2 503 end do !ig1 504 close(unt_dump) 505 end if ! Output real frequency axis 506 507 if (nfreqim>0) then ! output imaginary frequency axis 508 if (calc_epsilon) then 509 fname_dump=TRIM(fname)//'_EPS_Imfrq_Q'//TRIM(tagq) 510 else 511 fname_dump=TRIM(fname)//'_Imfrq_Q'//TRIM(tagq) 512 end if 513 if (open_file(fname_dump,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then 514 ABI_ERROR(msg) 515 end if 516 do ig1=ig1_start,ig1_end 517 do ig2=ig2_start,ig2_end 518 if (only_diag.AND.ig1/=ig2) CYCLE 519 write(unt_dump,'(a,i4,2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,/)')& 520 '# index= ',idx,' ig1= ',ig1,' ig2= ',ig2,& 521 '# q = ',Hscr0%qibz(:,iqibz),& 522 '# G = ',Hscr0%gvec(:,ig1),' G''= ',Hscr0%gvec(:,ig2),& 523 '# omega [eV] Re Im ' 524 do iomega=nfreqre+1,nfreqre+nfreqim 525 write(unt_dump,'(f8.2,4x,2es16.8)') AIMAG(Hscr0%omega(iomega))*Ha_eV,epsm1(ig1,ig2,iomega,1) 526 end do 527 write(unt_dump,*) 528 write(unt_dump,*) 529 end do !ig2 530 end do !ig1 531 close(unt_dump) 532 end if ! Check for imaginary frequencies 533 534 ! Check for complex plane values 535 if (nfreqc>0) then 536 if (calc_epsilon) then 537 fname_dump=TRIM(fname)//'_EPS_ZPLANE_Q'//TRIM(tagq) 538 else 539 fname_dump=TRIM(fname)//'_ZPLANE_Q'//TRIM(tagq) 540 end if 541 542 if (open_file(fname_dump,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then 543 ABI_ERROR(msg) 544 end if 545 546 do ig1=ig1_start,ig1_end 547 do ig2=ig2_start,ig2_end 548 if (only_diag.AND.ig1/=ig2) CYCLE 549 write(unt_dump,'(a,i4,2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,/)')& 550 '# index= ',idx,' ig1= ',ig1,' ig2= ',ig2,& 551 '# q = ',Hscr0%qibz(:,iqibz),& 552 '# G = ',Hscr0%gvec(:,ig1),' G''= ',Hscr0%gvec(:,ig2),& 553 '# omega [eV] Re Im ' 554 do iomega=1,nfreqre 555 write(unt_dump,'(2(f8.2),4x,2es16.8)') REAL(Hscr0%omega(iomega))*Ha_eV,& 556 AIMAG(Hscr0%omega(iomega))*Ha_eV,epsm1(ig1,ig2,iomega,1) 557 end do 558 write(unt_dump,*) 559 do ios=1,nfreqim 560 do iomega=1,nfreqre 561 if (iomega==1) then 562 io = nfreqre + ios 563 else 564 io = nfreqre + nfreqim + (ios-1)*(nfreqre-1) + (iomega-1) 565 end if 566 write(unt_dump,'(2(f8.2),4x,2es16.8)') REAL(Hscr0%omega(io))*Ha_eV,& 567 AIMAG(Hscr0%omega(io))*Ha_eV,epsm1(ig1,ig2,io,1) 568 end do 569 write(unt_dump,*) 570 end do 571 write(unt_dump,*) 572 write(unt_dump,*) 573 end do !ig2 574 end do !ig1 575 close(unt_dump) 576 end if ! Check for complex plane freqs 577 578 end do !iqibz 579 580 ABI_FREE(epsm1) 581 call Gsphere%free() 582 583 case(3) 584 ! Extract dielectric function and plasmon-pole stuff -------------------------- 585 ABI_CHECK(iomode==IO_MODE_FORTRAN, "netcdf output not coded") 586 write(std_out,'(a)') ' 3 => Calculation of dielectric function and plasmon-pole model' 587 588 npwe_asked=Hscr0%npwe; mqmem=Hscr0%nqibz 589 call init_Er_from_file(Er,fname,mqmem,npwe_asked,comm) 590 591 ! Initialize the G-sphere === 592 call Gsphere%init(Cryst,Hscr0%npwe,gvec=Hscr0%gvec) 593 594 boxcutmin=two; igmax=Gsphere%shlim(Gsphere%nsh) 595 ecut=Er%Hscr%Hdr%ecutdg 596 597 call getng(boxcutmin,1,ecut,Gsphere%gmet,k0,MPI_enreg%me_fft,& 598 mgfft,nfft,ngfft,MPI_enreg%nproc_fft,Cryst%nsym,paral_kgb0,Cryst%symrel,Cryst%tnons) 599 600 ! I am using standard valued, it would be better to call indefo 601 ! ngfft(1:3)=Er%Hscr%Hdr%ngfft(1:3) 602 ngfft(7)=112 603 ngfft(8)=get_cache_kb() 604 nfft = PRODUCT(ngfft(1:3)) 605 606 Dtset%gw_icutcoul=3; Dtset%rcut=zero 607 Dtset%vcutgeo=(/zero,zero,zero/); Dtset%boxcenter=(/zero,zero,zero/) 608 Dtset%ecutsigx = -1 609 610 if (Er%Hscr%nqlwl==0) then 611 nqlwl=1 612 ABI_MALLOC(qlwl,(3,nqlwl)) 613 qlwl(:,1)= GW_Q0_DEFAULT 614 else 615 nqlwl=Er%Hscr%nqlwl 616 ABI_MALLOC(qlwl,(3,nqlwl)) 617 qlwl(:,:)=Er%Hscr%qlwl(:,1:nqlwl) 618 end if 619 620 call Vcp%init(Gsphere,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%gw_icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,& 621 Hscr0%npwe,nqlwl,qlwl,comm) 622 ABI_FREE(qlwl) 623 624 ! Get the density from an external file === 625 ! If meshes are not the same, do an FFT interpolation to have rhor on ngfft. 626 call prompt(' Enter name for external DEN (or PAWDEN) file: ', fname_rho) 627 628 ABI_MALLOC(rhor,(nfft,Hscr0%Hdr%nspden)) 629 ABI_MALLOC(pawrhoij,(Hscr0%Hdr%natom*Hscr0%Hdr%usepaw)) 630 631 call read_rhor(fname_rho, cplex1, nfft, Hscr0%Hdr%nspden, ngfft, 1, MPI_enreg, rhor, hdr_rhor, pawrhoij, comm) 632 633 call hdr_rhor%free() 634 call pawrhoij_free(pawrhoij) 635 ABI_FREE(pawrhoij) 636 637 ABI_MALLOC(rhog,(2,nfft)) 638 call fourdp(1,rhog,rhor(:,1),-1,MPI_enreg,nfft,1,ngfft,0) 639 640 ABI_MALLOC(nhat,(nfft,Hscr0%Hdr%nspden*Hscr0%Hdr%usepaw)) 641 compch_sph=greatest_real; compch_fft=greatest_real 642 usexcnhat=0; usefinegrid=0 643 644 nelectron_exp = Hscr0%Hdr%nelect 645 646 call test_charge(nfft,nelectron_exp,Hscr0%Hdr%nspden,rhor,Cryst%ucvol,& 647 Hscr0%Hdr%usepaw,usexcnhat,usefinegrid,compch_sph,compch_fft,drude_plsmf) 648 GN_drude_plsmf = drude_plsmf 649 650 ! Read and in case make Epsilon^{-1} according the the options specified 651 id_required=4; ikxc=0; approx_type=0; option_test=0; dim_kxcg=0 652 ABI_MALLOC(kxcg,(nfft,dim_kxcg)) 653 654 call prompt(' Enter prefix for output files: ',prefix) 655 fname_dump=TRIM(prefix)//'_SCR' 656 657 orig_npwe = Er%npwe 658 write(std_out,'(2a,I0)') ch10,' Number of plane waves is: ',Er%npwe 659 write(std_out,'(a)',advance='no') ' Would you like to change it [Y/N] ?' 660 read(std_in,*) ans 661 if (ans=='Y'.or.ans=='y') then 662 write(std_out,'(a)',advance='no') ' Enter new no. of plane waves (0 means use old value): ' 663 read(std_in,*) ii 664 if (ii>0.or.ii<=Er%npwe) Er%npwe = ii 665 if (ii<0.or.ii>Er%npwe) then 666 ABI_ERROR(' Wrong value for no. of plane waves!') 667 end if 668 end if 669 670 if (is_scr) Er%mqmem=1 671 if (is_sus) Er%mqmem=0 672 call mkdump_Er(Er,Vcp,Er%npwe,Gsphere%gvec,dim_kxcg,kxcg,id_required,approx_type,ikxc,option_test,& 673 fname_dump,iomode,nfft,ngfft,comm) 674 Er%mqmem=1 675 676 call em1results_print(Er) 677 678 write(std_out,'(2a)',advance='no') ch10,& 679 ' Would you like to calculate the eigenvalues of eps^{-1}_GG''(omega) [Y/N] ? ' 680 read(std_in,*) ans 681 682 if (ans=='Y'.or.ans=='y') then 683 ABI_MALLOC(epsm1_eigen,(Er%npwe,Er%nomega)) 684 imax = 10 685 if (Er%npwe < imax) imax = Er%npwe 686 do iqibz=1,Er%nqibz 687 call int2char4(iqibz,tagq) 688 ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!') 689 fname_eigen=TRIM(prefix)//'_EM1_EIG_Q'//TRIM(tagq) 690 if (open_file(fname_eigen,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then 691 ABI_ERROR(msg) 692 end if 693 call decompose_epsm1(Er,iqibz,epsm1_eigen) 694 write(unt_dump,'(a)') '# First (max 10) eigenvalues of eps^{-1}(omega)' 695 write(unt_dump,'(a,3f12.6)')'# q = ',Hscr0%qibz(:,iqibz) 696 write(unt_dump,'(a)') '# REAL omega [eV] REAL(eigen(esp^-1(1,w))) AIMAG(eigen(esp^-1(1,w)) ...' 697 do iomega=1,Er%nomega_r 698 write(unt_dump,'(21(es16.8))')REAL(Er%omega(iomega))*Ha_eV,& 699 (REAL(epsm1_eigen(ii,iomega)),ii=1,imax),(AIMAG(epsm1_eigen(ii,iomega)),ii=1,imax) 700 end do 701 close(unt_dump) 702 end do 703 ABI_FREE(epsm1_eigen) 704 ABI_FREE(kxcg) 705 end if ! Calculate eigenvalues 706 707 ! Analyze the PPmodel. 708 write(std_out,'(2a)') ch10,' Would you like to analyse plasmon-pole models [Y/N] ? ' 709 read(std_in,*)ans 710 711 if (ans=='Y'.or.ans=='y') then 712 write(std_out,'(2a,f6.2,a)') ch10,' Plasma frequency for GN PPM is: ',GN_drude_plsmf*Ha_eV, ' eV' 713 write(std_out,'(a)',advance='no') ' Would you like to change it [Y/N] ?' 714 read(std_in,*) ans 715 if (ans=='Y'.or.ans=='y') then 716 write(std_out,'(2a)',advance='no') ch10,' Enter plasma frequency [eV]: ' 717 read(std_in,*) GN_drude_plsmf 718 GN_drude_plsmf = GN_drude_plsmf/Ha_eV 719 end if 720 721 write(std_out,'(2a)') ch10,' Would you like to calculate the plasmon-pole model' 722 write(std_out,'(a)',advance='no') ' eigenvalues of eps^{-1}_GG''(omega) [Y/N] ? ' 723 read(std_in,*) ans 724 725 if (ans=='Y'.or.ans=='y') then 726 ABI_MALLOC(ppm_eigen,(PPm%npwc,Er%nomega)) 727 imax = 10; if (Er%npwe < imax) imax = Er%npwe 728 do iqibz=1,Er%nqibz 729 do ppmodel=1,2 730 731 call int2char4(iqibz,tagq) 732 ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!') 733 if (ppmodel==1) fname_dump=TRIM(prefix)//'_PPM_GN_EM1_EIG_Q'//TRIM(tagq) 734 if (ppmodel==2) fname_dump=TRIM(prefix)//'_PPM_HL_EM1_EIG_Q'//TRIM(tagq) 735 if (ppmodel==3) fname_dump=TRIM(prefix)//'_PPM_vdLH_EM1_EIG_Q'//TRIM(tagq) 736 if (ppmodel==4) fname_dump=TRIM(prefix)//'_PPM_EF_EM1_EIG_Q'//TRIM(tagq) 737 738 if (open_file(fname_eigen,msg,newunit=unt_dump,status='new',form='formatted') /= 0) then 739 ABI_ERROR(msg) 740 end if 741 742 call ppm_free(PPm) 743 if (ppmodel==1) then 744 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,GN_drude_plsmf,Dtset%gw_invalid_freq) 745 else 746 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,drude_plsmf,Dtset%gw_invalid_freq) 747 end if 748 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfft,Gsphere%gvec,ngfft,rhor(:,1),iqibz) 749 750 call get_PPm_eigenvalues(PPm,iqibz,Er%Hscr%zcut,Er%nomega,Er%omega,Vcp,ppm_eigen) 751 752 write(unt_dump,'(a)') '# First (max 10) eigenvalues of eps^{-1}(omega) from Plasmon-pole model' 753 write(unt_dump,'(a,3f12.6)')'# q = ',Hscr0%qibz(:,iqibz) 754 select case(ppmodel) 755 case(1) 756 write(unt_dump,'(a)') '# ppmodel = 1 : Godby - Needs' 757 case(2) 758 write(unt_dump,'(a)') '# ppmodel = 2 : Hybertsen - Louie' 759 case(3) 760 write(unt_dump,'(a)') '# ppmodel = 3 : von der Linden - Horsch' 761 case(4) 762 write(unt_dump,'(a)') '# ppmodel = 4 : Engel - Farid' 763 end select 764 write(unt_dump,'(a)') '# REAL omega [eV] REAL(eigen(ppm_eps^-1(1,w))) AIMAG(eigen(ppm_eps^-1(1,w)) ...' 765 do iomega=1,Er%nomega_r 766 write(unt_dump,'(21(es16.8))')REAL(Er%omega(iomega))*Ha_eV,& 767 (REAL(ppm_eigen(ii,iomega)),ii=1,imax),(AIMAG(ppm_eigen(ii,iomega)),ii=1,imax) 768 end do 769 close(unt_dump) 770 771 end do !ppmodel 772 end do ! iqibz 773 ABI_FREE(ppm_eigen) 774 end if ! Calculate PPM eigenvalues 775 776 ! Optionally output eps^{-1}_GG''(w) for a given set of GG' and gridpoints 777 write(std_out,'(2a)',advance='no') ch10,' Would you like to extract eps^{-1}_GG''(omega) for the PPM [Y/N] ?' 778 read(std_in,*) ans 779 780 if (ans=='Y'.or.ans=='y') then 781 ! Reconstruct e^{-1}_GG'(w) according to PPmodel for statistical analysis. 782 write(std_out,'(a)') ' Enter the number of frequency points in the' 783 write(std_out,'(a)') ' interval 0 - freqremax (0 means same as input file ): ' 784 read(std_in,*) nfreqre 785 if (nfreqre==0) then 786 nfreqre = Er%nomega_r 787 nfreqim = Er%nomega_i 788 nfreq_tot = Er%nomega 789 freqremax = REAL(Er%omega(Er%nomega_r)) 790 ABI_MALLOC(omega,(nfreq_tot)) 791 omega(:) = Er%omega(:) 792 same_freqs = .TRUE. 793 else 794 write(std_out,'(a)') ' Enter the value of freqremax (in eV): ' 795 read(std_in,*) freqremax 796 nfreqim = Er%nomega_i 797 nfreq_tot = nfreqre+nfreqim 798 ABI_MALLOC(omega,(nfreqre+Er%nomega_i)) 799 do iomega=1,nfreqre 800 omega(iomega) = CMPLX((freqremax/REAL((nfreqre-1)))*(iomega-1),zero) 801 end do 802 omega(nfreqre+1:nfreq_tot) = Er%omega(Er%nomega_r+1:Er%nomega) 803 same_freqs = .FALSE. 804 end if ! frequencies 805 806 do iqibz=1,Er%nqibz 807 qtmp(:)=Er%qibz(:,iqibz); if (normv(qtmp,Cryst%gmet,'G')<GW_TOLQ0) qtmp(:)=zero 808 call int2char4(iqibz,tagq) 809 ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!') 810 811 ! At this time only the Godby-Needs and Hybertsen-Louie models 812 ! TODO: Check the results from the others 813 do ppmodel=1,2 814 815 call ppm_free(PPm) 816 if (ppmodel==1) then 817 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,GN_drude_plsmf,Dtset%gw_invalid_freq) 818 else 819 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,drude_plsmf,Dtset%gw_invalid_freq) 820 end if 821 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,& 822 nfft,Gsphere%gvec,ngfft,rhor(:,1),iqibz) 823 824 ! Prepare file for data on real omega axis 825 if (ppmodel==1) fname_dump=TRIM(prefix)//'_PPM_w_GN_Q'//TRIM(tagq) 826 if (ppmodel==2) fname_dump=TRIM(prefix)//'_PPM_w_HL_Q'//TRIM(tagq) 827 if (ppmodel==3) fname_dump=TRIM(prefix)//'_PPM_w_vdLH_Q'//TRIM(tagq) 828 if (ppmodel==4) fname_dump=TRIM(prefix)//'_PPM_w_EF_Q'//TRIM(tagq) 829 830 if (open_file(fname_dump,msg,newunit=unt_dump,status='replace',form='formatted') /= 0) then 831 ABI_ERROR(msg) 832 end if 833 834 ! Prepare file for data on imaginary omega axis 835 if (ppmodel==1) fname_dump2=TRIM(prefix)//'_PPM_iw_GN_Q'//TRIM(tagq) 836 if (ppmodel==2) fname_dump2=TRIM(prefix)//'_PPM_iw_HL_Q'//TRIM(tagq) 837 if (ppmodel==3) fname_dump2=TRIM(prefix)//'_PPM_iw_vdLH_Q'//TRIM(tagq) 838 if (ppmodel==4) fname_dump2=TRIM(prefix)//'_PPM_iw_EF_Q'//TRIM(tagq) 839 840 if (open_file(fname_dump2,msg,newunit=unt_dump2,status='replace',form='formatted') /= 0) then 841 ABI_ERROR(msg) 842 end if 843 844 ABI_MALLOC(em1_ppm,(nfreq_tot)) 845 846 ig1 = 0; ig2 = 0 847 write(std_out,'(3a,I0,a,I0)') ch10,' Enter indices for G and G''.',& 848 'Entering 0 exits the loop. iqibz = ',iqibz,' ppmodel = ',ppmodel 849 850 do 851 write(std_out,'(2(a),I0,a)',advance='NO') ch10,' Enter index for G (1 - ',Er%npwe,' ): ' 852 read(std_in,*)ig1 853 if (ig1==0) EXIT 854 if (ig1<0.OR.ig1>Er%npwe) ABI_ERROR(' index out of bounds') 855 write(std_out,'(2(a),I0,a)',advance='NO') ch10,' Enter index for G'' (1 - ',Er%npwe,' ): ' 856 read(std_in,*)ig2 857 if (ig2==0) EXIT 858 if (ig2<0.OR.ig2>Er%npwe) ABI_ERROR(' index out of bounds') 859 860 ! Generate the PPM representation of epsilon^-1 861 call getem1_from_PPm_one_ggp(PPm,iqibz,Er%Hscr%zcut,nfreq_tot,omega,Vcp,em1_ppm,ig1,ig2) 862 863 write(unt_dump,'(a,I1)') '# epsilon^-1_GG''(omega) from ppmodel = ',ppmodel 864 write(unt_dump,'(2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,3F9.4,a,3F9.4,a,/a,f9.4,a,f9.4,a,/,a,/)')& 865 '# ig1= ',ig1,' ig2= ',ig2,& 866 '# q = ',Er%qibz(:,iqibz),& 867 '# G = ',Er%gvec(:,ig1),' G''= ',Er%gvec(:,ig2),& 868 '# G = (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig1)),& 869 ') G''= (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig2)),')',& 870 '# 1/2|G|^2 =',half*normv(Er%gvec(:,ig1),Cryst%gmet,'G')**2,& 871 ' Ha 1/2|G''|^2 =',half*normv(Er%gvec(:,ig2),Cryst%gmet,'G')**2,' Ha',& 872 '# omega [eV] Re Im ' 873 write(unt_dump2,'(a,I1)') '# epsilon^-1_GG''(iomega) from ppmodel = ',ppmodel 874 write(unt_dump2,'(2(a,i8),/,a,3f12.6,/,a,3i6,a,3i6,/,a,3F9.4,a,3F9.4,a,/a,f9.4,a,f9.4,a,/,a,/)')& 875 '# ig1= ',ig1,' ig2= ',ig2,& 876 '# q = ',Er%qibz(:,iqibz),& 877 '# G = ',Er%gvec(:,ig1),' G''= ',Er%gvec(:,ig2),& 878 '# G = (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig1)),& 879 ') G''= (',MATMUL(two_pi*Cryst%gmet,Er%gvec(:,ig2)),')',& 880 '# 1/2|G|^2 =',half*normv(Er%gvec(:,ig1),Cryst%gmet,'G')**2,& 881 ' Ha 1/2|G''|^2 =',half*normv(Er%gvec(:,ig2),Cryst%gmet,'G')**2,' Ha',& 882 '# iomega [eV] Re Im ' 883 884 do iomega=1,nfreqre 885 if (same_freqs) then 886 write(unt_dump,'(f8.2,4x,4es16.8)') REAL(omega(iomega))*Ha_eV,em1_ppm(iomega),& 887 Er%epsm1(ig1,ig2,iomega,iqibz) 888 else 889 write(unt_dump,'(f8.2,4x,2es16.8)') REAL(omega(iomega))*Ha_eV,em1_ppm(iomega) 890 end if 891 end do 892 ! First output the iomega = 0 point 893 write(unt_dump2,'(f8.2,4x,4es16.8)') AIMAG(omega(1))*Ha_eV,em1_ppm(1),& 894 Er%epsm1(ig1,ig2,1,iqibz) 895 ! Then the rest 896 do iomega=nfreqre+1,nfreq_tot 897 write(unt_dump2,'(f8.2,4x,4es16.8)') AIMAG(omega(iomega))*Ha_eV,em1_ppm(iomega),& 898 Er%epsm1(ig1,ig2,iomega,iqibz) 899 end do 900 write(unt_dump,*) 901 write(unt_dump,*) 902 write(unt_dump2,*) 903 write(unt_dump2,*) 904 end do ! Empty 905 ABI_FREE(em1_ppm) 906 close(unt_dump) 907 close(unt_dump2) 908 909 end do ! ppmodel 910 end do ! iqibz 911 ABI_FREE(omega) 912 end if ! Output epsilon for PPM 913 914 ! Optionally statistics for all PPMs 915 write(std_out,'(2a)',advance='no') ch10,' Would you like to output statistics for all PPMs [Y/N] ?' 916 read(std_in,*) ans 917 918 if (ans=='Y'.or.ans=='y') then 919 nfreqre = Er%nomega_r 920 nfreq_tot = Er%nomega 921 freqremax = REAL(Er%omega(Er%nomega_r)) 922 ABI_MALLOC(real_omega,(nfreqre)) 923 real_omega(:) = REAL(Er%omega(1:nfreqre)) 924 925 do iqibz=1,Er%nqibz 926 do ppmodel=1,2 927 928 qtmp(:)=Er%qibz(:,iqibz) 929 if (normv(qtmp,Cryst%gmet,'G')<GW_TOLQ0) qtmp(:)=zero 930 931 call ppm_free(PPm) 932 if (ppmodel==1) then 933 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,GN_drude_plsmf,Dtset%gw_invalid_freq) 934 else 935 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,ppmodel,drude_plsmf,Dtset%gw_invalid_freq) 936 end if 937 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,& 938 nfft,Gsphere%gvec,ngfft,rhor(:,1),iqibz) 939 940 ! Prepare ratios and density for the f-sum rule 941 ABI_MALLOC_OR_DIE(qratio,(orig_npwe,orig_npwe), ierr) 942 ABI_MALLOC_OR_DIE(rhoggp,(Er%npwe,Er%npwe), ierr) 943 944 call cqratio(orig_npwe,Gsphere%gvec,qtmp,Cryst%gmet,Cryst%gprimd,qratio) 945 ! Arrange n(G-G')->n(G,G') 946 ierr=0 947 do ig1=1,Er%npwe 948 do ig2=1,Er%npwe 949 gmgp_idx = g2ifft(Gsphere%gvec(:,ig1)-Gsphere%gvec(:,ig2),ngfft) 950 if (gmgp_idx/=0) then 951 rhoggp(ig1,ig2)=CMPLX(rhog(1,gmgp_idx),rhog(2,gmgp_idx)) 952 else 953 ierr=ierr+1 954 rhoggp(ig1,ig2)=czero 955 end if 956 end do 957 end do 958 if (ierr/=0) then 959 write(std_out,'(a,i0,a)')' Found ',ierr,' G1-G2 vectors falling outside the FFT box. ' 960 end if 961 962 ! Prepare files 963 call int2char4(iqibz,tagq) 964 ABI_CHECK((tagq(1:1)/='#'),'Bug: string length too short!') 965 if (ppmodel==1) fname_dump=TRIM(prefix)//'_norms_GN_Q'//TRIM(tagq) 966 if (ppmodel==2) fname_dump=TRIM(prefix)//'_norms_HL_Q'//TRIM(tagq) 967 if (open_file(fname_dump,msg, newunit=unt_dump, status='replace',form='formatted') /= 0) then 968 ABI_ERROR(msg) 969 end if 970 write(unt_dump,'(a)') '# Various norms integrated through spline interpolation' 971 write(unt_dump,'(a)') '# over all frequencies in the input file,' 972 write(unt_dump,'(a)') '# for all G and G'' vectors.' 973 write(unt_dump,'(a,I0)') '# ppmodel: ',ppmodel 974 write(unt_dump,'(a,I0)') '# Number of frequencies: ',nfreqre 975 write(unt_dump,'(a,f12.6)') '# Maximum frequency : ',freqremax 976 write(unt_dump,'(a)') '# Columns:' 977 write(unt_dump,'(2a)') '# ig1 ig2 |eps-eps_PPM|/|eps|',& 978 ' |eps-eps_PPM| |eps| |eps_PPM| G G''' 979 if (ppmodel==1) fname_dump2=TRIM(prefix)//'_f_sumrule_GN_Q'//TRIM(tagq) 980 if (ppmodel==2) fname_dump2=TRIM(prefix)//'_f_sumrule_HL_Q'//TRIM(tagq) 981 982 if (open_file(fname_dump2,msg,newunit=unt_dump2,status='replace',form='formatted') /= 0) then 983 ABI_ERROR(msg) 984 end if 985 986 write(unt_dump2,'(a)') '# The fulfillment of the f-sum rule: I(epsilon) =' 987 write(unt_dump2,'(a)') '# int_0^{inf}{omega*Im[epsilon_G,G''(omega)]}/C_qGG''' 988 write(unt_dump2,'(a)') '# C_qGG'' = ' 989 write(unt_dump2,'(a)') '# -Pi/2*omega_p^2*(q+G)*(q+G'')/|q+G|^2*n(G-G'')/n(0)' 990 write(unt_dump2,'(a)') '# for all G and G'' vectors.' 991 write(unt_dump2,'(a,I0)') '# ppmodel: ',ppmodel 992 write(unt_dump2,'(a,I0)') '# Number of frequencies: ',nfreqre 993 write(unt_dump2,'(a,f12.6)') '# Maximum frequency : ',freqremax 994 write(unt_dump2,'(a)') '# Columns:' 995 write(unt_dump2,'(3a)') '# ig1 ig2 I(epsilon)',& 996 ' I(eps_PPM) Re[n(G-G'')] Im[n(G-G'')] qratio I1*C_qGG''',& 997 ' Re[Omegatwsq] Im[Omegatwsq] Re[omegatw] Im[omegatw] |G| 1/2|G|^2' 998 999 ABI_MALLOC(em1_ppm,(nfreq_tot)) 1000 ABI_MALLOC(ftab,(nfreqre)) 1001 ABI_MALLOC(ysp,(3,nfreqre)) 1002 ABI_MALLOC(work,(nfreqre)) 1003 ABI_MALLOC(eint,(nfreqre)) 1004 1005 do ig1=1,Er%npwe 1006 write(std_out,'(2(a,I0))') ' ig1= ',ig1, ' of ',Er%npwe 1007 do ig2=1,Er%npwe 1008 !ig2 = ig1 1009 call getem1_from_PPm_one_ggp(PPm,iqibz,Er%Hscr%zcut,nfreq_tot,Er%omega,Vcp,em1_ppm,ig1,ig2) 1010 1011 ! Calculate norms in real 1012 eps_diff=0; eps_norm=0; eps_ppm_norm=0 1013 ftab(1:nfreqre) = ABS(Er%epsm1(ig1,ig2,1:nfreqre,iqibz)-em1_ppm(1:nfreqre)) 1014 call cspint(ftab,real_omega,nfreqre,real_omega(1),real_omega(nfreqre),ysp,eint,work,eps_diff) 1015 ftab(1:nfreqre) = ABS(Er%epsm1(ig1,ig2,1:nfreqre,iqibz)) 1016 call cspint(ftab,real_omega,nfreqre,real_omega(1),real_omega(nfreqre),ysp,eint,work,eps_norm) 1017 ftab(1:nfreqre) = ABS(em1_ppm(1:nfreqre)) 1018 call cspint(ftab,real_omega,nfreqre,real_omega(1),real_omega(nfreqre),ysp,eint,work,eps_ppm_norm) 1019 write(unt_dump,'(2i6,f12.4,3es14.4,6i4)') ig1,ig2,eps_diff/eps_norm,eps_diff,& 1020 eps_norm,eps_ppm_norm,Er%gvec(:,ig1),Er%gvec(:,ig2) 1021 1022 ! Evaluate the f-sum rule 1023 if (ig1==ig2) then 1024 ftab(1:nfreqre) = real_omega(1:nfreqre)*AIMAG(Er%epsm1(ig1,ig2,1:nfreqre,iqibz)) 1025 else 1026 ! Dephase first - HERE epsm1 is changed! 1027 call remove_phase(epsm1(ig1,ig2,:,1),Hscr_file(1)%nomega,phase) 1028 ftab(1:nfreqre) = real_omega(1:nfreqre)*AIMAG(Er%epsm1(ig1,ig2,1:nfreqre,iqibz)) 1029 end if 1030 1031 call cspint(ftab,real_omega,nfreqre,real_omega(1),& 1032 real_omega(nfreqre),ysp,eint,work,eps_diff) 1033 1034 if (ig1==ig2) then 1035 factor = -two*pi*pi*REAL(rhoggp(ig1,ig2))*qratio(ig1,ig2) 1036 else 1037 rhoggp(ig1,ig2) = CMPLX(COS(phase),-SIN(phase))*rhoggp(ig1,ig2) 1038 factor = -two*pi*pi*REAL(rhoggp(ig1,ig2))*qratio(ig1,ig2) 1039 end if 1040 1041 if (ABS(qratio(ig1,ig2))>zero) then 1042 value1 = eps_diff/factor 1043 if (ppmodel==1) then 1044 value2 = -pi*half*(REAL(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2))& 1045 /(REAL(PPm%omegatw(iqibz)%vals(ig1,ig2))))& 1046 /factor*(2*sqrt(pi*rhoggp(1,1))) 1047 else 1048 value2 = -pi*half*(SQRT(REAL(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2))))& 1049 /factor*(2*sqrt(pi*rhoggp(1,1))) 1050 end if 1051 else 1052 value1 = zero 1053 value2 = zero 1054 end if 1055 1056 write(unt_dump2,'(2i6,12es14.4)') ig1,ig2,value1,value2,& 1057 REAL(rhoggp(ig1,ig2)),AIMAG(rhoggp(ig1,ig2)),qratio(ig1,ig2),& 1058 eps_diff,REAL(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)),& 1059 AIMAG(PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)),& 1060 REAL(PPm%omegatw(iqibz)%vals(ig1,ig2)),& 1061 AIMAG(PPm%omegatw(iqibz)%vals(ig1,ig2)),& 1062 normv(Er%gvec(:,ig1),Cryst%gmet,'G'),& 1063 half*normv(Er%gvec(:,ig1),Cryst%gmet,'G')**2 1064 1065 end do !ig2 1066 end do !ig1 1067 1068 ABI_FREE(em1_ppm) 1069 ABI_FREE(ftab) 1070 ABI_FREE(ysp) 1071 ABI_FREE(work) 1072 ABI_FREE(eint) 1073 ABI_FREE(qratio) 1074 ABI_FREE(rhoggp) 1075 close(unt_dump); close(unt_dump2) 1076 1077 end do ! ppmodel 1078 end do ! iqibz 1079 1080 ABI_FREE(real_omega) 1081 end if ! Output statistics 1082 1083 call ppm_free(PPm) 1084 end if ! If ppmodel>0 1085 1086 ABI_FREE(rhor) 1087 ABI_FREE(rhog) 1088 ABI_FREE(nhat) 1089 1090 call Vcp%free() 1091 call em1results_free(Er) 1092 call Gsphere%free() 1093 1094 case(4) 1095 ! Remove real frequencies ---------------------------------------------------------- 1096 write(std_out,'(2(a))') ch10,' Do you want to remove every other real frequency (= 1) ?' 1097 write(std_out,'(a)') ' or specify for each real frequency individually (= 2) ?' 1098 write(std_out,'(a)') ' or remove ALL real frequencies (= 3) ?' 1099 read(std_in,*)choice 1100 1101 ! Calculate the total number of real freq 1102 nfreqre = 0; nfreqim = 0 1103 do ifrq=1,Hscr_file(1)%nomega 1104 ! If frequency is not imaginary, count. 1105 if (AIMAG(Hscr_file(1)%omega(ifrq)) < tol8) nfreqre = nfreqre + 1 1106 if (REAL(Hscr_file(1)%omega(ifrq)) < tol8 .and. AIMAG(Hscr_file(1)%omega(ifrq))>tol8) nfreqim = nfreqim + 1 1107 end do 1108 1109 nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq 1110 write(std_out,'(2a,I0,a)') ch10,' Found ',nfreq_tot,' frequencies.' 1111 write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10 1112 1113 ! Array with the index of frequencies to be kept. 1114 ABI_MALLOC(freq_indx,(nfreq_tot,1)) 1115 freq_indx = 0 1116 1117 select case(choice) 1118 case(1) 1119 ! Remove every other frequency 1120 write(std_out,'(2(a))') ch10,' Removing every other real frequency, i.e. every even one.' 1121 write(std_out,'(a)') ' If the total number of frequencies is odd, the first and last one will be kept.' 1122 write(std_out,'(a)') ' If the total number is even, the first one will still be in the final set.' 1123 1124 ! Test for no real frequencies 1125 ABI_CHECK(nfreqre /= 0, "No real frequencies in file!") 1126 1127 ii=nfreqre; nfreqre = 0 1128 do ifrq=1,ii 1129 if (.not. iseven(ifrq)) then 1130 nfreqre = nfreqre + 1; freq_indx(nfreqre,1) = ifrq 1131 end if 1132 end do 1133 write(std_out,'(2a,I0,a)') ch10,' ',nfreqre,' real frequencies will be kept.' 1134 1135 case(2) 1136 ! Specify freq. individually 1137 ii = nfreqre; nfreqre = 0 1138 do ifrq=1,ii 1139 write(std_out,'(a,f12.6,a)') ' Would you like to keep freq. at: ',REAL(Hscr_file(1)%omega(ifrq))*Ha_eV,' eV? [y/n]' 1140 read(std_in,*) ans 1141 if (ans=='Y'.or.ans=='y') then 1142 nfreqre = nfreqre + 1; freq_indx(nfreqre,1) = ifrq 1143 end if 1144 end do 1145 write(std_out,'(2a,I0,a)') ch10,' ',nfreqre,' real frequencies will be kept.' 1146 1147 case(3) 1148 ! Remove all real freq. 1149 nfreqre = 0 1150 1151 case default 1152 ABI_ERROR("Invalid choice!") 1153 end select 1154 1155 ! Add imaginary frequencies if any 1156 if (nfreqim > 0) then 1157 nfreqim = 0 1158 do ifrq=1,Hscr_file(1)%nomega 1159 if (AIMAG(Hscr_file(1)%omega(ifrq)) > tol8) then 1160 nfreqim = nfreqim + 1; freq_indx(nfreqre+nfreqim,1) = ifrq 1161 end if 1162 end do 1163 end if 1164 1165 nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq 1166 write(std_out,'(3(a,i0),a)')' Finally, we have ',nfreq_tot,' frequencies. ',nfreqre,' real, and ',nfreqim,' imaginary.' 1167 1168 call prompt(' Enter the full name of the final output file: ', fname_out) 1169 1170 if (endswith(filenames(1), ".nc") .and. .not. endswith(fname_out, ".nc")) then 1171 fname_out = nctk_ncify(fname_out) 1172 call wrtout(std_out,"- Added .nc extension to output file as input data is in netcdf format.") 1173 end if 1174 1175 call ioscr_wremove(filenames(1), hscr_file(1), fname_out, nfreq_tot, freq_indx, hscr_merge) 1176 ABI_FREE(freq_indx) 1177 1178 case(5) 1179 ! Remove imaginary frequencies ------------------------------------------------ 1180 1181 ! Calculate the total number of freq 1182 nfreqre = 0; nfreqim = 0 1183 do ifrq=1,Hscr_file(1)%nomega 1184 ! If frequency is not imaginary, count. 1185 if (AIMAG(Hscr_file(1)%omega(ifrq))<tol8) nfreqre = nfreqre + 1 1186 if (REAL(Hscr_file(1)%omega(ifrq))<tol8.and.AIMAG(Hscr_file(1)%omega(ifrq))>tol8) nfreqim = nfreqim + 1 1187 end do ! ifrq 1188 1189 ! Test for no real frequencies 1190 if (nfreqim == 0) then 1191 ABI_ERROR("No imaginary frequencies in file!") 1192 end if 1193 1194 nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq 1195 write(std_out,'(2a,I0,a)') ch10,' Found ',nfreq_tot,' frequencies.' 1196 write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10 1197 1198 ABI_MALLOC(freq_indx,(nfreq_tot,1)) 1199 freq_indx = 0 1200 1201 ! Specify freq. individually 1202 ii=nfreq_tot; nfreqim = 0 1203 do ifrq=nfreqre+1,ii 1204 write(std_out,'(a,f12.6,a)')& 1205 ' Would you like to keep imaginary freq. at: ',AIMAG(Hscr_file(1)%omega(ifrq))*Ha_eV,' eV? [y/n]' 1206 read(std_in,*) ans 1207 if (ans=='Y'.or.ans=='y') then 1208 nfreqim = nfreqim + 1; freq_indx(nfreqre+nfreqim,1) = ifrq 1209 end if 1210 end do ! ifrq 1211 write(std_out,'(2a,I0,a)') ch10,' ',nfreqim,' imaginary frequencies will be kept.' 1212 1213 ! Add real frequencies if any 1214 if (nfreqre > 0) then 1215 nfreqre = 0 1216 do ifrq=1,Hscr_file(1)%nomega 1217 if (AIMAG(Hscr_file(1)%omega(ifrq)) < tol8) then 1218 nfreqre = nfreqre + 1; freq_indx(nfreqre,1) = ifrq 1219 end if 1220 end do 1221 end if 1222 1223 nfreq_tot = nfreqre + nfreqim ! Here nfreq_tot becomes the *true* number of freq 1224 write(std_out,'(2a,I0,a)') ch10,' Finally, we have ',nfreq_tot,' frequencies.' 1225 write(std_out,'(2(a,I0),2a)') ' ',nfreqre,' real, and ',nfreqim,' imaginary.',ch10 1226 1227 call prompt(' Enter the full name of the final output file: ',fname_out) 1228 1229 if (endswith(filenames(1), ".nc") .and. .not. endswith(fname_out, ".nc")) then 1230 fname_out = nctk_ncify(fname_out) 1231 call wrtout(std_out,"- Added .nc extension to output file as input data is in netcdf format.") 1232 end if 1233 1234 call ioscr_wremove(filenames(1), hscr_file(1), fname_out, nfreq_tot, freq_indx, hscr_merge) 1235 1236 ABI_FREE(freq_indx) 1237 1238 case(6) 1239 ! Model screening ------------------------------------------------------------- 1240 ABI_ERROR("Model screening has been removed") 1241 1242 !case(9) 1243 ! TODO 1244 ! ! netcdf --> Fortran converter ------------------------------------------------------------- 1245 ! call prompt(' Enter the name of the final output Fortran file: ',fname_out) 1246 ! ! fname_out extension should be consistent with filenames(1) 1247 ! call ioscr_nc2fort(filenames(1), fname_out) 1248 1249 case default 1250 ! Bail if choice is wrong 1251 write(std_out,*) ' Invalid choice! Exiting...' 1252 goto 100 1253 end select 1254 1255 end if ! Single file mode 1256 1257 call timein(tcpu,twall) 1258 1259 tsec(1)=tcpu-tcpui 1260 tsec(2)=twall-twalli 1261 1262 write(std_out, '(a,a,a,f13.1,a,f13.1)' ) '-',ch10,'- Proc. 0 individual time (sec): cpu=',tsec(1),' wall=',tsec(2) 1263 1264 ! ===================== 1265 ! ==== Free memory ==== 1266 ! ===================== 1267 ABI_FREE(filenames) 1268 ABI_SFREE(kxcg) 1269 ABI_SFREE(foundq) 1270 1271 call cryst%free() 1272 call Kmesh%free() 1273 call Qmesh%free() 1274 call destroy_mpi_enreg(MPI_enreg) 1275 1276 nullify(Hscr0) 1277 call Hscr_merge%free() 1278 1279 do ifile=1,nfiles 1280 call Hscr_file(ifile)%free() 1281 end do 1282 ABI_FREE(Hscr_file) 1283 1284 call flush_unit(std_out) 1285 1286 call abinit_doctor("__mrgscr") 1287 1288 100 call xmpi_end() 1289 1290 end program mrgscr