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