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 832 : 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 : !!***
|
