Line data Source code
1 : !!****m* ABINIT/m_eph_driver
2 : !! NAME
3 : !! m_eph_driver
4 : !!
5 : !! FUNCTION
6 : !! Driver for EPH calculations
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2009-2026 ABINIT group (MG, MVer, GA)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_eph_driver
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 : use m_xomp
29 : use m_hdr
30 : use m_crystal
31 : use m_ebands
32 : use m_dtset
33 : use m_efmas_defs
34 : use m_dtfil
35 : use m_ddb
36 : use m_ddb_hdr
37 : use m_dvdb, only : dvdb_t
38 : use m_ifc
39 : use m_phonons
40 : use m_nctk
41 : use netcdf
42 :
43 : use defs_datatypes, only : pseudopotential_type
44 : use defs_abitypes, only : MPI_type
45 : use m_io_tools, only : file_exists, open_file
46 : use m_time, only : cwtime, cwtime_report
47 : use m_fstrings, only : strcat, sjoin, ftoa, itoa
48 : use m_fftcore, only : print_ngfft
49 : use m_wfk, only : wfk_read_ebands
50 : use m_rta, only : rta_driver, ibte_driver
51 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
52 : use m_pawang, only : pawang_type
53 : use m_pawrad, only : pawrad_type
54 : use m_pawtab, only : pawtab_type
55 : use m_paw_an, only : paw_an_type, paw_an_free !, paw_an_nullify, paw_an_init,
56 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
57 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
58 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, pawrhoij_symrhoij
59 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
60 : use m_phgamma, only : eph_phgamma
61 : use m_efmas, only : efmasdeg_free_array, efmasval_free_array, efmas_ncread
62 : use m_gkk, only : eph_gkk, ncwrite_v1qnu
63 : use m_phpi, only : eph_phpi
64 : use m_sigmaph, only : sigmaph
65 : use m_pspini, only : pspini
66 : use m_ephtk, only : ephtk_update_ebands
67 : use m_migdal_eliashberg, only : migdal_eliashberg_iso !, migdal_eliashberg_aniso
68 : use m_gstore, only : gstore_t, gstore_symmetrize
69 : use m_gstore_sigmaph, only : gstore_sigmaph
70 : use m_gstore_converters, only : gstore_convert
71 : use m_berry_curvature, only : berry_curvature
72 : use m_cumulant, only : cumulant_driver
73 : use m_frohlich, only : frohlich_t, frohlichmodel_zpr, frohlichmodel_polaronmass
74 : use m_gwpt, only : gwpt_run
75 : use m_varpeq, only : varpeq_run, varpeq_plot
76 : use m_eph_path, only : eph_path_run
77 : use m_wkk, only : wkk_run
78 :
79 : implicit none
80 :
81 : private
82 : !!***
83 :
84 : public :: eph
85 : !!***
86 :
87 : contains
88 : !!***
89 :
90 : !!****f* m_eph_driver/eph
91 : !! NAME
92 : !! eph
93 : !!
94 : !! FUNCTION
95 : !! Main routine to compute electron phonon coupling matrix elements and
96 : !! calculate related properties - superconducting Tc, phonon linewidths, electronic renormalization
97 : !! due to phonons and temperature effects...
98 : !!
99 : !! INPUTS
100 : !! acell(3)=Length scales of primitive translations (bohr)
101 : !! codvsn=Code version
102 : !! dtfil<datafiles_type>=Variables related to files.
103 : !! dtset<dataset_type>=All input variables for this dataset.
104 : !! pawang<pawang_type)>=PAW angular mesh and related data.
105 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
106 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
107 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
108 : !! Before entering the first time in the routine, a significant part of Psps has been initialized :
109 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm,
110 : !! and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in
111 : !! the call to pspini. The next time the code enters bethe_salpeter, Psps might be identical to the
112 : !! one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90.
113 : !! rprim(3,3)=Dimensionless real space primitive translations.
114 : !! xred(3,natom)=Reduced atomic coordinates.
115 : !!
116 : !! NOTES
117 : !!
118 : !! ON THE USE OF FFT GRIDS:
119 : !! =================
120 : !! In case of PAW:
121 : !! ---------------
122 : !! Two FFT grids are used:
123 : !! - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
124 : !! It is defined by nfft, ngfft, mgfft, ...
125 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
126 : !! - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
127 : !! It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
128 : !! In case of norm-conserving:
129 : !! ---------------------------
130 : !! - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
131 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
132 : !!
133 : !! SOURCE
134 :
135 109 : subroutine eph(acell, codvsn, dtfil, dtset, pawang, pawrad, pawtab, psps, rprim, xred)
136 :
137 : !Arguments ------------------------------------
138 : !scalars
139 : character(len=8),intent(in) :: codvsn
140 : type(datafiles_type),intent(in) :: dtfil
141 : type(dataset_type),intent(inout) :: dtset
142 : type(pawang_type),intent(inout) :: pawang
143 : type(pseudopotential_type),intent(inout) :: psps
144 : !arrays
145 : real(dp),intent(in) :: acell(3),rprim(3,3),xred(3,dtset%natom)
146 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
147 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
148 :
149 : !Local variables ------------------------------
150 : !scalars
151 : integer,parameter :: master = 0, selectz0 = 0, nsphere0 = 0, prtsrlr0 = 0
152 : integer,parameter :: with_cplex0 = 0, with_cplex1 = 1, with_cplex2 = 2
153 : integer :: ii,comm,nprocs,my_rank,psp_gencond,mgfftf,nfftf
154 : integer :: iblock_dielt_zeff, iblock_dielt, iblock_quadrupoles, ddb_nqshift, ierr, npert_miss
155 : integer :: omp_ncpus, work_size, nks_per_proc, lwsym, qptopt, ncid
156 : real(dp):: eff, mempercpu_mb, max_wfsmem_mb, nonscal_mem
157 : real(dp) :: ecore,ecut_eff,ecutdg_eff,gsqcutc_eff,gsqcutf_eff
158 : real(dp) :: cpu,wall,gflops
159 : logical :: use_wfk, use_wfq, use_dvdb, use_sigeph, use_drhodb, use_gstore, gstore_from_file
160 : character(len=500) :: msg
161 : character(len=fnlen) :: wfk0_path, wfq_path, ddb_filepath, dvdb_filepath, sigeph_filepath, path, drhodb_filepath, gstore_path
162 109 : type(hdr_type) :: wfk0_hdr, wfq_hdr
163 11227 : type(crystal_t) :: cryst, cryst_ddb
164 218 : type(ebands_t) :: ks_ebands, ks_ebands_kq, qp_ebands
165 109 : type(ddb_type) :: ddb, ddb_lw
166 5668 : type(ddb_hdr_type) :: ddb_hdr
167 17767 : type(dvdb_t) :: dvdb, drhodb
168 654 : type(ifc_type) :: ifc
169 : type(pawfgr_type) :: pawfgr
170 109 : type(mpi_type) :: mpi_enreg
171 1417 : type(phdos_t) :: phdos
172 436 : type(gstore_t) :: gstore
173 436 : type(frohlich_t) :: frohlich
174 : !arrays
175 : integer :: ngfftc(18), ngfftf(18), count_wminmax(2), units(2)
176 : real(dp),parameter :: k0(3)=zero
177 218 : real(dp) :: wminmax(2), dielt(3,3), zeff(3,3,dtset%natom), zeff_raw(3,3,dtset%natom)
178 218 : real(dp) :: qdrp_cart(3,3,3,dtset%natom)
179 218 : real(dp),allocatable :: ddb_qshifts(:,:), kpt_efmas(:,:)
180 109 : type(efmasdeg_type),allocatable :: efmasdeg(:)
181 109 : type(efmasval_type),allocatable :: efmasval(:,:)
182 : !type(pawfgrtab_type),allocatable :: pawfgrtab(:)
183 : !type(paw_ij_type),allocatable :: paw_ij(:)
184 : !type(paw_an_type),allocatable :: paw_an(:)
185 :
186 : !************************************************************************
187 :
188 : ! This part performs the initialization of the basic objects used to perform e-ph calculations:
189 : !
190 : ! 1) Crystal structure `cryst`
191 : ! 2) Ground state band energies: `ebands`
192 : ! 3) Interatomic force constants: `ifc`
193 : ! 4) DVDB database with the dvscf potentials
194 : ! 5) Pseudos and PAW basic objects.
195 : !
196 : ! Once we have these objects, we can call specialized routines for e-ph calculations.
197 : ! Notes:
198 : !
199 : ! * Any modification to the basic objects mentioned above should be done here (e.g. change of efermi)
200 : ! * This routines shall not allocate big chunks of memory. The CPU-demanding sections should be
201 : ! performed in the subdriver that will employ different MPI distribution schemes optimized for that particular task.
202 :
203 : DBG_ENTER('COLL')
204 :
205 109 : if (psps%usepaw == 1) then
206 0 : ABI_ERROR("PAW not implemented")
207 : ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
208 : end if
209 :
210 : ! abirules!
211 : if (.False.) write(std_out,*)acell,codvsn,rprim,xred
212 :
213 109 : comm = xmpi_world; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
214 327 : units = [std_out, ab_out]
215 :
216 : #ifndef HAVE_MPI_IBCAST
217 : do ii=1,5
218 : ABI_WARNING("Your MPI library does not provide MPI_IBCAST. Calculations parallelized over perturbations will be slow")
219 : end do
220 : #endif
221 :
222 : ! Initialize filenames
223 109 : wfk0_path = dtfil%fnamewffk
224 109 : wfq_path = dtfil%fnamewffq
225 109 : ddb_filepath = dtfil%filddbsin
226 :
227 : ! Use the ddb file as prefix if getdvdb or irddvb are not given in the input.
228 109 : dvdb_filepath = dtfil%fildvdbin
229 109 : if (dvdb_filepath == ABI_NOFILE) then
230 41 : dvdb_filepath = dtfil%filddbsin; ii = len_trim(dvdb_filepath); dvdb_filepath(ii-2:ii+1) = "DVDB"
231 : end if
232 :
233 109 : drhodb_filepath = dtfil%fildrhodbin
234 109 : if (drhodb_filepath == ABI_NOFILE) then
235 106 : drhodb_filepath = dtfil%filddbsin; ii = len_trim(drhodb_filepath); drhodb_filepath(ii-2:ii+1) = "DRHODB"
236 : end if
237 :
238 109 : sigeph_filepath = dtfil%filsigephin
239 :
240 909 : use_wfk = all(dtset%eph_task /= [0, 5, -5, 6, +15, -15, -16, 16])
241 109 : use_wfq = ((dtset%irdwfq /= 0 .or. dtset%getwfq /= 0 .or. dtset%getwfq_filepath /= ABI_NOFILE) .and. dtset%eph_frohlichm /= 1)
242 :
243 : ! If eph_task is needed and ird/get variables are not provided, assume WFQ == WFK
244 109 : if (any(dtset%eph_task == [2, -2, 3]) .and. .not. use_wfq) then
245 0 : wfq_path = wfk0_path
246 0 : use_wfq = .True.
247 : write(msg, "(4a)")&
248 0 : "eph_task requires WFQ but none among (irdwfq, getwfq, getwfq_filepath) is specified in input.", ch10, &
249 0 : "Will read WFQ wavefunctions from WFK file:", trim(wfk0_path)
250 0 : ABI_COMMENT(msg)
251 : end if
252 :
253 109 : use_dvdb = (dtset%eph_task /= 0 .and. dtset%eph_frohlichm /= 1 .and. abs(dtset%eph_task) /= 7 .and. dtset%eph_task /= 13)
254 109 : use_sigeph = (dtset%eph_task == 9)
255 109 : use_gstore = (dtset%eph_task == 13)
256 109 : use_drhodb = (dtset%eph_task == 17)
257 :
258 109 : if (my_rank == master) then
259 : ! GA: Let ddb object handle the error at reading time
260 : !if (.not. file_exists(ddb_filepath)) ABI_ERROR(sjoin("Cannot find DDB file:", ddb_filepath))
261 109 : if (use_dvdb .and. .not. file_exists(dvdb_filepath)) ABI_ERROR(sjoin("Cannot find DVDB file:", dvdb_filepath))
262 109 : if (use_sigeph .and. .not. file_exists(sigeph_filepath)) ABI_ERROR(sjoin("Cannot find SIGEPH file:", sigeph_filepath))
263 109 : if (use_drhodb .and. .not. file_exists(drhodb_filepath)) ABI_ERROR(sjoin("Cannot find DRHODB file:", drhodb_filepath))
264 :
265 : ! Accept WFK file in Fortran or netcdf format.
266 109 : if (use_wfk .and. nctk_try_fort_or_ncfile(wfk0_path, msg) /= 0) then
267 0 : ABI_ERROR(sjoin("Cannot find GS WFK file:", wfk0_path, ". Error:", msg))
268 : end if
269 : ! WFQ file
270 109 : if (use_wfq) then
271 2 : if (nctk_try_fort_or_ncfile(wfq_path, msg) /= 0) then
272 0 : ABI_ERROR(sjoin("Cannot find GS WFQ file:", wfq_path, ". Error:", msg))
273 : end if
274 : end if
275 : end if ! master
276 :
277 : ! Broadcast filenames (needed because they might have been changed if we are using netcdf files)
278 109 : if (use_wfk) then
279 97 : call xmpi_bcast(wfk0_path, master, comm, ierr)
280 97 : call wrtout(units, sjoin("- Reading GS states from WFK file:", wfk0_path))
281 : end if
282 109 : if (use_wfq) then
283 2 : call xmpi_bcast(wfq_path, master, comm, ierr)
284 2 : call wrtout(units, sjoin("- Reading GS states from WFQ file:", wfq_path) )
285 : end if
286 109 : call wrtout(units, sjoin("- Reading DDB from file:", ddb_filepath))
287 109 : if (use_dvdb) call wrtout(units, sjoin("- Reading DVDB from file:", dvdb_filepath))
288 109 : if (dtset%eph_frohlichm /= 0) call wrtout(units, sjoin("- Reading EFMAS information from file:", dtfil%fnameabi_efmas))
289 109 : call wrtout(units, ch10//ch10)
290 :
291 : ! autoparal section
292 : ! TODO: This just to activate autoparal in AbiPy. Lot of things should be improved.
293 109 : if (dtset%max_ncpus /= 0) then
294 0 : write(ab_out,'(a)')"--- !Autoparal"
295 0 : write(ab_out,"(a)")"# Autoparal section for EPH runs"
296 0 : write(ab_out,"(a)") "info:"
297 0 : write(ab_out,"(a,i0)")" autoparal: ",dtset%autoparal
298 0 : write(ab_out,"(a,i0)")" max_ncpus: ",dtset%max_ncpus
299 0 : write(ab_out,"(a,i0)")" nkpt: ",dtset%nkpt
300 0 : write(ab_out,"(a,i0)")" nsppol: ",dtset%nsppol
301 0 : write(ab_out,"(a,i0)")" nspinor: ",dtset%nspinor
302 0 : write(ab_out,"(a,i0)")" mband: ",dtset%mband
303 0 : write(ab_out,"(a,i0)")" eph_task: ",dtset%eph_task
304 :
305 0 : work_size = dtset%nkpt * dtset%nsppol
306 : ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI.
307 0 : nonscal_mem = zero
308 0 : max_wfsmem_mb = (two * dp * dtset%mpw * dtset%mband * dtset%nkpt * dtset%nsppol * dtset%nspinor * b2Mb) * 1.1_dp
309 :
310 : ! List of configurations.
311 : ! Assuming an OpenMP implementation with perfect speedup!
312 0 : write(ab_out,"(a)")"configurations:"
313 :
314 0 : do ii=1,dtset%max_ncpus
315 0 : nks_per_proc = work_size / ii
316 0 : nks_per_proc = nks_per_proc + mod(work_size, ii)
317 0 : eff = (one * work_size) / (ii * nks_per_proc)
318 : ! Add the non-scalable part and increase by 10% to account for other datastructures.
319 0 : mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp
320 :
321 0 : do omp_ncpus=1,1 !xomp_get_max_threads()
322 0 : write(ab_out,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus
323 0 : write(ab_out,"(a,i0)")" mpi_ncpus: ",ii
324 0 : write(ab_out,"(a,i0)")" omp_ncpus: ",omp_ncpus
325 0 : write(ab_out,"(a,f12.9)")" efficiency: ",eff
326 0 : write(ab_out,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb
327 : end do
328 : end do
329 0 : write(ab_out,'(a)')"..."
330 0 : call wrtout(std_out, "Stopping now!")
331 0 : call abi_abort('PERS', exit_status=0, print_config=.false.)
332 : end if
333 :
334 109 : call cwtime(cpu, wall, gflops, "start")
335 :
336 109 : if (use_wfk) then
337 : ! Construct crystal and ebands from the GS WFK file.
338 : ! TODO: If we want to update the eigenvalues to perform some kind of self-consistency,
339 : ! this operation should be done here! Either scissors operator or external file with QP energies.
340 97 : ks_ebands = wfk_read_ebands(wfk0_path, comm, out_hdr=wfk0_hdr)
341 97 : call wfk0_hdr%vs_dtset(dtset)
342 :
343 97 : cryst = wfk0_hdr%get_crystal()
344 97 : call cryst%print(header="crystal structure from WFK file")
345 :
346 : ! Here we change the GS bands (Fermi level, scissors operator ...)
347 : ! All the modifications to ks_ebands should be done here.
348 97 : call ephtk_update_ebands(dtset, ks_ebands, "Ground state energies")
349 :
350 97 : if (dtfil%filqpdatain /= ABI_NOFILE) then
351 : ! Read QP energies from an external file.
352 0 : call qp_ebands%read_qpdata(ks_ebands, dtfil%filqpdatain, comm)
353 0 : call ephtk_update_ebands(dtset, qp_ebands, "QP energies")
354 : else
355 97 : call ks_ebands%copy(qp_ebands)
356 : end if
357 :
358 : ! Need to update the WFK header to reflect the changes in ks_ebands.
359 : ! because we may need to write the header to ncfile
360 : ! NB: eigenvalues are not stored in the header.
361 97 : wfk0_hdr%occopt = ks_ebands%occopt
362 97 : call get_eneocc_vect(ks_ebands, "occ", wfk0_hdr%occ)
363 97 : wfk0_hdr%fermie = ks_ebands%fermie
364 97 : wfk0_hdr%nelect = ks_ebands%nelect
365 : end if
366 :
367 109 : if (use_wfq) then
368 : ! Read WFQ and construct ebands on the shifted grid.
369 2 : ks_ebands_kq = wfk_read_ebands(wfq_path, comm, out_hdr=wfq_hdr)
370 : ! GKA TODO: Have to construct a header with the proper set of q-shifted k-points then compare against dtset.
371 : !call wfq_hdr%vs_dtset(dtset)
372 2 : call wfq_hdr%free()
373 2 : call ephtk_update_ebands(dtset, ks_ebands_kq, "Ground state energies (K+Q)")
374 : end if
375 :
376 109 : call cwtime_report(" eph%init", cpu, wall, gflops)
377 :
378 : ! =======================================
379 : ! Output useful info on electronic bands
380 : ! =======================================
381 109 : call cwtime(cpu, wall, gflops, "start")
382 :
383 109 : if (my_rank == master) then
384 : ! Fermi Surface
385 109 : if (dtset%prtfsurf /= 0) then
386 1 : path = strcat(dtfil%filnam_ds(4), "_BXSF")
387 1 : call wrtout(units, sjoin("- Writing Fermi surface to file:", path))
388 1 : if (ks_ebands%write_bxsf(cryst, path) /= 0) then
389 0 : msg = "Cannot produce file for Fermi surface, check log file for more info"
390 0 : ABI_WARNING(msg)
391 0 : call wrtout(ab_out, msg)
392 : end if
393 : end if
394 :
395 : ! Nesting factor (requires qpath)
396 109 : if (dtset%prtnest /= 0 .and. dtset%ph_nqpath > 0) then
397 1 : path = strcat(dtfil%filnam_ds(4), "_NEST")
398 1 : call wrtout(ab_out, sjoin("- Writing nesting factor to file:", path))
399 1 : if (ks_ebands%write_nesting(cryst, path, dtset%prtnest, &
400 : dtset%tsmear, dtset%fermie_nest, dtset%ph_qpath(:,1:dtset%ph_nqpath), msg) /= 0) then
401 0 : ABI_WARNING(msg)
402 0 : call wrtout(ab_out,msg)
403 : end if
404 : end if
405 109 : if (use_wfk) call ks_ebands%write(dtset%prtebands, dtfil%filnam_ds(4))
406 : end if
407 :
408 109 : call cwtime_report(" eph%ebands_postprocess:", cpu, wall, gflops)
409 :
410 : ! Read the DDB file.
411 109 : if (use_wfk) then
412 97 : call ddb%from_file(ddb_filepath, ddb_hdr, cryst_ddb, comm, prtvol=dtset%prtvol)
413 :
414 : ! DDB cryst comes from DFPT --> no time-reversal if q /= 0
415 : ! Change the value so that we use the same as the GS part.
416 97 : cryst_ddb%timrev = cryst%timrev
417 97 : if (cryst%compare(cryst_ddb, header=" Comparing WFK crystal with DDB crystal") /= 0) then
418 0 : ABI_ERROR("Crystal structure from WFK and DDB do not agree! Check messages above!")
419 : end if
420 97 : call cryst_ddb%free()
421 : else
422 : ! Get crystal from DDB.
423 : ! Warning: We may loose precision in rprimd and xred because DDB in text format does not have enough significant digits.
424 12 : call ddb%from_file(ddb_filepath, ddb_hdr, cryst, comm, prtvol=dtset%prtvol)
425 12 : call ddb%set_brav(dtset%brav)
426 : end if
427 :
428 : ! Change the bravais lattice if needed
429 109 : call ddb%set_brav(dtset%brav)
430 :
431 : ! MR: a new ddb is necessary for the longwave quantities due to incompability of it with automatic reshapes
432 : ! that ddb%val and ddb%flg experience when passed as arguments of some routines
433 : ! GA: Should replace with ddb_hdr%with_d3E_lw
434 109 : iblock_quadrupoles = 0
435 8509 : qdrp_cart = zero
436 109 : if (ddb_hdr%has_d3E_lw) then
437 1 : lwsym = 1
438 1 : call ddb_lw_copy(ddb, ddb_lw, ddb_hdr)
439 : ! GA: FIXME Bad interface
440 1 : iblock_quadrupoles = ddb_lw%get_quadrupoles(ddb_hdr%ddb_version, lwsym, BLKTYP_d3E_lw, qdrp_cart)
441 1 : call ddb_lw%free()
442 : end if
443 :
444 : ! Set the q-shift for the DDB (well we mainly use gamma-centered q-meshes)
445 109 : ddb_nqshift = 1
446 545 : ABI_CALLOC(ddb_qshifts, (3, ddb_nqshift))
447 436 : ddb_qshifts(:,1) = dtset%ddb_shiftq(:)
448 :
449 : ! Get Dielectric Tensor
450 109 : iblock_dielt = ddb%get_dielt(dtset%rfmeth, dielt)
451 :
452 : ! Get Dielectric Tensor and Effective Charges
453 : ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
454 109 : iblock_dielt_zeff = ddb%get_dielt_zeff(cryst, dtset%rfmeth, dtset%chneut, selectz0, dielt, zeff, zeff_raw=zeff_raw)
455 109 : if (my_rank == master) then
456 109 : if (iblock_dielt_zeff == 0) then
457 49 : call wrtout(units, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", ddb_filepath))
458 49 : call wrtout(units, " Values initialized with zeros.")
459 : else
460 60 : call wrtout(units, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", ddb_filepath))
461 : end if
462 : end if
463 :
464 : ! The default value is 1. Here we set the flags to zero if Q* is not available.
465 109 : if (iblock_quadrupoles == 0) then
466 108 : dtset%dipquad = 0
467 108 : dtset%quadquad = 0
468 : end if
469 :
470 109 : if (my_rank == master) then
471 109 : if (iblock_quadrupoles == 0) then
472 108 : call wrtout(units, sjoin("- Cannot find quadrupole tensor in DDB file:", ddb_filepath))
473 108 : call wrtout(units, " Values initialized with zeros.")
474 : else
475 1 : call wrtout(units, sjoin("- Found quadrupole tensor in DDB file:", ddb_filepath))
476 : end if
477 : end if
478 :
479 : call ifc%init(cryst, ddb, &
480 : dtset%brav, dtset%asr, dtset%symdynmat, dtset%dipdip, dtset%rfmeth, &
481 : dtset%ddb_ngqpt, ddb_nqshift, ddb_qshifts, dielt, zeff, &
482 : qdrp_cart, nsphere0, dtset%rifcsph, prtsrlr0, dtset%enunit,1, comm, &
483 109 : dipquad=dtset%dipquad, quadquad=dtset%quadquad)
484 :
485 109 : ABI_FREE(ddb_qshifts)
486 218 : if (my_rank == master) call ifc%print([std_out])
487 :
488 : ! Output phonon band structure (requires qpath)
489 109 : if (dtset%prtphbands /= 0) call ifc_mkphbs(ifc, cryst, dtset, dtfil%filnam_ds(4), comm)
490 :
491 109 : if (dtset%prtphdos == 1) then
492 52 : call wrtout(std_out, " Computing Phonon DOS. Use prtphdos 0 to disable this part.")
493 52 : wminmax = zero
494 22 : do
495 : call phdos%init(cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, dtset%ph_smear, dtset%ph_ngqpt, &
496 74 : dtset%ph_nqshift, dtset%ph_qshift, "", wminmax, count_wminmax, comm)
497 193 : if (all(count_wminmax == 0)) exit
498 22 : wminmax(1) = wminmax(1) - abs(wminmax(1)) * 0.05; wminmax(2) = wminmax(2) + abs(wminmax(2)) * 0.05
499 22 : call phdos%free()
500 22 : write(msg, "(a, 2f8.5)") "Initial frequency mesh not large enough. Recomputing PHDOS with wmin, wmax: ",wminmax
501 22 : call wrtout(std_out, msg)
502 : end do
503 :
504 52 : if (my_rank == master) then
505 52 : if (dtset%prtvol > 0) then
506 : ! Disabled by default because it's slow and we use netcdf that is much better.
507 22 : path = strcat(dtfil%filnam_ds(4), "_PHDOS")
508 22 : call wrtout(units, sjoin("- Writing phonon DOS to file:", path))
509 22 : call phdos%print(path)
510 : end if
511 :
512 52 : path = strcat(dtfil%filnam_ds(4), "_PHDOS.nc")
513 52 : call wrtout(units, sjoin("- Writing phonon DOS to netcdf file:", path))
514 52 : NCF_CHECK_MSG(nctk_open_create(ncid, path, xmpi_comm_self), sjoin("Creating PHDOS.nc file:", path))
515 52 : NCF_CHECK(cryst%ncwrite(ncid))
516 52 : call phdos%ncwrite(ncid)
517 52 : NCF_CHECK(nf90_close(ncid))
518 : end if
519 52 : call phdos%free()
520 : end if ! prtphdos
521 :
522 109 : if (dtset%prtbltztrp == 1 .and. my_rank == master) then
523 1 : call ifc%outphbtrap(cryst, dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, dtfil%filnam_ds(4))
524 : ! BoltzTraP output files in GENEric format
525 1 : call ks_ebands%prtbltztrp(cryst, dtfil%filnam_ds(4))
526 : end if
527 :
528 : ! Output phonon isosurface in Xcrysden format.
529 109 : if (dtset%prtphsurf == 1) then
530 1 : path = strcat(dtfil%filnam_ds(4), "_PH.bxsf")
531 1 : call wrtout(units, sjoin("- Writing phonon frequencies in Xcrysden format to file:", path))
532 1 : call ifc%printbxsf(cryst, dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, path, comm)
533 : end if
534 :
535 109 : call cwtime_report(" eph%ifc:", cpu, wall, gflops)
536 :
537 : ! Initialize the object used to read DeltaVscf
538 109 : if (use_dvdb) then
539 96 : call dvdb%init(dvdb_filepath, dtset%gpu_option, comm)
540 96 : ABI_CHECK(dvdb%has_fields("pot1", msg), sjoin(dvdb_filepath, msg))
541 :
542 : ! DVDB cryst comes from DPPT --> no time-reversal if q /= 0
543 : ! Change the value so that we use the same as the GS part.
544 96 : dvdb%cryst%timrev = cryst%timrev
545 96 : if (cryst%compare(dvdb%cryst, header=" Comparing WFK crystal with DVDB crystal") /= 0) then
546 0 : ABI_ERROR("Crystal structure from WFK and DVDB do not agree! Check messages above!")
547 : end if
548 96 : dvdb%prtvol = dtset%prtvol
549 96 : if (dtset%prtvol > 10) dvdb%debug = .True.
550 :
551 : ! This to symmetrize the DFPT potentials.
552 96 : dvdb%symv1 = dtset%symv1scf
553 :
554 : ! Copy brav variable
555 96 : dvdb%brav = dtset%brav
556 :
557 : ! Select algorithm for generating the list of R-points and the weights used to compute W(r,R)
558 96 : dvdb%rspace_cell = dtset%dvdb_rspace_cell
559 :
560 : !call dvdb%load_ddb(dtset%prtvol, comm, ddb=ddb)
561 :
562 : ! Set qdamp from frohl_params
563 96 : dvdb%qdamp = dtset%dvdb_qdamp
564 :
565 : ! Set quadrupoles
566 7512 : dvdb%qstar = qdrp_cart; if (iblock_quadrupoles /= 0) dvdb%has_quadrupoles = .True.
567 :
568 : ! Set dielectric tensor, BECS and associated flags.
569 : ! This flag activates automatically the treatment of the long-range term in the Fourier interpolation
570 : ! of the DFPT potentials except when dvdb_add_lr == 0
571 96 : dvdb%add_lr = dtset%dvdb_add_lr
572 96 : if (iblock_dielt /= 0) then
573 611 : dvdb%has_dielt = .True.; dvdb%dielt = dielt
574 : end if
575 96 : if (iblock_dielt_zeff /= 0) then
576 2585 : dvdb%has_zeff = .True.; dvdb%zeff = zeff; dvdb%zeff_raw = zeff_raw
577 : end if
578 96 : if (.not. dvdb%has_dielt .or. .not. (dvdb%has_zeff .or. dvdb%has_quadrupoles)) then
579 49 : if (dvdb%add_lr /= 0) then
580 49 : dvdb%add_lr = 0
581 49 : ABI_WARNING("Setting dvdb_add_lr to 0. Long-range term won't be subtracted in Fourier interpolation.")
582 : end if
583 : end if
584 :
585 96 : if (dvdb%add_lr == 2) then
586 0 : if (dvdb%has_quadrupoles) then
587 0 : call wrtout(std_out, "dvdb_add_lr == 2 --> Quadrupoles are set to zero and won't be used in the interpolation")
588 : end if
589 0 : dvdb%has_quadrupoles = .False.
590 0 : dvdb%qstar = zero
591 : end if
592 :
593 96 : if (my_rank == master) then
594 192 : call dvdb%print([std_out], "DVDB FILE", dtset%prtvol)
595 96 : call dvdb%list_perts([-1, -1, -1], npert_miss)
596 96 : ABI_CHECK(npert_miss == 0, sjoin(itoa(npert_miss), "independent perturbation(s) are missing in the DVDB file!"))
597 : end if
598 : end if
599 :
600 109 : if (use_drhodb) then
601 : ! Store DRHODB as a DVDB object
602 3 : call drhodb%init(drhodb_filepath, dtset%gpu_option, comm)
603 3 : ABI_CHECK(drhodb%has_fields("den1", msg), sjoin(drhodb_filepath, msg))
604 :
605 : ! DVDB cryst comes from DPPT --> no time-reversal if q /= 0
606 : ! Change the value so that we use the same as the GS part.
607 3 : drhodb%cryst%timrev = cryst%timrev
608 3 : if (cryst%compare(drhodb%cryst, header=" Comparing WFK crystal with DRHODB crystal") /= 0) then
609 0 : ABI_ERROR("Crystal structure from WFK and DRHODB do not agree! Check messages above!")
610 : end if
611 3 : if (dtset%prtvol > 10) drhodb%debug = .True.
612 :
613 : ! This to symmetrize the DFPT densities.
614 3 : drhodb%symv1 = dtset%symv1scf
615 :
616 : ! Copy brav variable
617 3 : drhodb%brav = dtset%brav
618 :
619 : ! Select algorithm for generating the list of R-points and the weights used to compute W(r,R)
620 3 : drhodb%rspace_cell = dtset%dvdb_rspace_cell
621 :
622 : !call drhodb%load_ddb(dtset%prtvol, comm, ddb=ddb)
623 :
624 : ! Set qdamp, quadrupoles and all long-range terms to 0.
625 3 : drhodb%qdamp = 0
626 243 : drhodb%qstar = 0
627 3 : drhodb%has_quadrupoles = .False.
628 3 : drhodb%add_lr = 0
629 39 : drhodb%has_dielt = .False.; drhodb%dielt = 0
630 159 : drhodb%has_zeff = .False.; drhodb%zeff = 0; drhodb%zeff_raw = 0
631 :
632 3 : if (my_rank == master) then
633 6 : call drhodb%print([std_out], "DRHODB file", 0)
634 3 : call drhodb%list_perts([-1, -1, -1], npert_miss)
635 3 : ABI_CHECK(npert_miss == 0, sjoin(itoa(npert_miss), "independent perturbation(s) are missing in the DVDB file!"))
636 : end if
637 : end if
638 :
639 : call pawfgr_init(pawfgr, dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, &
640 109 : gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=cryst%gmet, k0=k0)
641 :
642 218 : call print_ngfft([std_out], ngfftc, header='Coarse FFT mesh used for the wavefunctions')
643 218 : call print_ngfft([std_out], ngfftf, header='Dense FFT mesh used for densities and potentials')
644 :
645 : ! Fake MPI_type for the sequential part.
646 109 : call initmpi_seq(mpi_enreg)
647 109 : call mpi_enreg%distribfft%init_seq('c', ngfftc(2), ngfftc(3), 'all')
648 109 : call mpi_enreg%distribfft%init_seq('f', ngfftf(2), ngfftf(3), 'all')
649 :
650 : ! I am not sure yet the EFMAS file will be needed as soon as eph_frohlichm/=0. To be decided later.
651 109 : if (dtset%eph_frohlichm /= 0) then
652 3 : NCF_CHECK(nctk_open_read(ncid, dtfil%fnameabi_efmas, xmpi_comm_self))
653 3 : call efmas_ncread(efmasdeg, efmasval, kpt_efmas, ncid)
654 3 : NCF_CHECK(nf90_close(ncid))
655 : end if
656 :
657 : ! ===========================================
658 : ! === Open and read pseudopotential files ===
659 : ! ===========================================
660 109 : call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm)
661 :
662 : ! Release nkpt-based arrays in dtset to decrease memory requirement if dense sampling.
663 : ! EPH routines should not access them after this point.
664 109 : if (all(dtset%eph_task /= [6, 10])) call dtset%free_nkpt_arrays()
665 :
666 : ! ====================================================
667 : ! === This is the real EPH stuff once all is ready ===
668 : ! ====================================================
669 :
670 109 : select case (dtset%eph_task)
671 : case (0)
672 : ! This is just to access the DDB post-processing tools for phonons.
673 5 : continue
674 :
675 : case (1)
676 : ! Compute phonon linewidths in metals.
677 : call eph_phgamma(wfk0_path, dtfil, ngfftc, ngfftf, dtset, cryst, ks_ebands, dvdb, ifc, &
678 5 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
679 :
680 : case (2, -2)
681 : ! Compute e-ph matrix elements (legacy version)
682 2 : ABI_CHECK(dtset%useylm == 0, "useylm != 0 not implemented/tested")
683 : call eph_gkk(wfk0_path, wfq_path, dtfil, ngfftc, ngfftf, dtset, cryst, ks_ebands, ks_ebands_kq, dvdb, ifc, &
684 2 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
685 :
686 : case (3)
687 : ! Compute phonon-electron self-energy.
688 0 : ABI_CHECK(dtset%useylm == 0, "useylm != 0 not implemented/tested")
689 : call eph_phpi(wfk0_path, wfq_path, dtfil, ngfftc, ngfftf, dtset, cryst, ks_ebands, ks_ebands_kq, dvdb, ifc, &
690 0 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
691 :
692 : case (4, -4)
693 : ! Compute electron-phonon self-energy (phonon contribution).
694 : call sigmaph(wfk0_path, dtfil, ngfftc, ngfftf, dtset, cryst, ks_ebands, dvdb, ifc, wfk0_hdr, &
695 58 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
696 :
697 : ! Compute transport properties in the RTA/IBTE only if sigma_erange has been used
698 155 : if (dtset%eph_task == -4 .and. any(abs(dtset%sigma_erange) > zero)) then
699 12 : if (dtset%ibte_prep > 0) then
700 1 : call ibte_driver(dtfil, ngfftc, dtset, ks_ebands, cryst, pawtab, psps, comm) ! Solve IBTE
701 : else
702 11 : call rta_driver(dtfil, ngfftc, dtset, ks_ebands, cryst, pawtab, psps, comm) ! Compute RTA
703 : end if
704 : end if
705 :
706 : case (24)
707 : ! Compute e-ph self-energy from GSTORE.nc file.
708 : call gstore_sigmaph(wfk0_path, ngfftc, ngfftf, dtset, dtfil, cryst, ks_ebands, qp_ebands, wfk0_hdr, dvdb, ifc, &
709 7 : pawfgr, pawtab, psps, mpi_enreg, comm)
710 :
711 : case (5, -5)
712 : ! Interpolate the DFPT potential.
713 : call dvdb%interpolate_and_write(dtset, dtfil%fnameabo_dvdb, ngfftc, ngfftf, cryst, &
714 8 : ifc%ngqpt, ifc%nqshft, ifc%qshft, comm)
715 :
716 : case (6)
717 : ! Estimate zero-point renormalization and temperature-dependent electronic structure using the Frohlich model.
718 2 : if (my_rank == master) call frohlichmodel_zpr(frohlich, cryst, dtset, efmasdeg, efmasval, ifc)
719 :
720 : case (7)
721 : ! Compute phonon-limited RTA from SIGEPH.nc file.
722 1 : call rta_driver(dtfil, ngfftc, dtset, qp_ebands, cryst, pawtab, psps, comm)
723 :
724 : case (8)
725 : ! Solve IBTE from SIGEPH.nc file.
726 0 : call ibte_driver(dtfil, ngfftc, dtset, qp_ebands, cryst, pawtab, psps, comm)
727 :
728 : case (9)
729 : ! Compute cumulant from SIGEPH.nc file.
730 2 : call cumulant_driver(dtfil, dtset, qp_ebands, cryst, comm)
731 :
732 : case (10)
733 : ! Estimate polaron effective mass in the triply-degenerate VB or CB cubic case
734 1 : if (my_rank == master) call frohlichmodel_polaronmass(frohlich, cryst, dtset, efmasdeg, efmasval, ifc)
735 :
736 : case (11)
737 : ! Compute and write e-ph matrix elements to GSTORE.nc file.
738 5 : if (dtfil%filgstorein /= ABI_NOFILE) then
739 : ! Init gstore from pre-existent file. gstore_gname and read_dw are not relevant here.
740 0 : call wrtout(units, sjoin(" Restarting GSTORE computation from:", dtfil%filgstorein))
741 : call gstore%from_ncpath(dtfil%filgstorein, with_cplex2, dtset, dtfil, cryst, qp_ebands, ifc, &
742 0 : "atom", dtset%gstore_gname, .False., comm)
743 : else
744 5 : gstore_path = strcat(dtfil%filnam_ds(4), "_GSTORE.nc")
745 5 : call gstore%init(gstore_path, dtset, dtfil, wfk0_hdr, cryst, qp_ebands, ifc, comm)
746 : end if
747 :
748 : call gstore%compute(wfk0_path, ngfftc, ngfftf, dtset, dtfil, cryst, qp_ebands, ifc, dvdb, &
749 5 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
750 :
751 5 : gstore_path = gstore%path
752 5 : call gstore%free()
753 :
754 5 : if (len(trim(dtset%gstore_convert)) /= 0) call gstore_convert(gstore_path, dtset, dtfil, cryst, qp_ebands, ifc, comm)
755 :
756 : ! Wannierize the e-ph matrix elements if the ABIWAN.nc file is provided.
757 5 : if (dtfil%filabiwanin /= ABI_NOFILE) then
758 : ! Init gstore from pre-existent file. gstore_gname and read_dw are not relevant here.
759 : call gstore%from_ncpath(gstore_path, with_cplex2, dtset, dtfil, cryst, qp_ebands, ifc, &
760 0 : "atom", dtset%gstore_gname, .False., comm)
761 0 : call gstore%wannierize_and_write_gwan(dvdb, dtfil)
762 0 : call gstore%free()
763 : end if
764 :
765 : !case (-11)
766 : ! Typical workflow for gstore with Wannierization:
767 : !
768 : ! 1) Wannierize with Abinit and wannier90 in library mode to get the ABIWAN.nc file.
769 : !
770 : ! 2) Pass ABIWAN.nc to the EPH code to compute GSTORE.nc only for the bands included in the wannierization step.
771 : !
772 : ! 3) Call gstore%wannierize_and_write_gwan to compute g(R_e, R_p) and save results to the GWAN.nc file.
773 : !
774 : ! 3) Start new job to compute properties with extra dense k/q-meshes (eph_ngkpt_fine and eph_ngqpt_fine)
775 : !
776 : ! - Init gstore object with extra dense meshes, possibly filtered and MPI-grid to distribute gvals.
777 : ! - Decide if gvals should be precomputed and stored or computed on the fly.
778 : ! - Read GWAN.nc file to build gstore%gqk(spin)%wan
779 : ! - Pass gstore object to the eph_task routines (what about ks_ebands)?
780 :
781 : ! call gstore%from_ncpath(gstore_path, with_cplex2, dtset, dtfil, cryst, ks_ebands, ifc, comm)
782 : ! call gstore%wannierize(dvdb, dtfil)
783 : ! call gstore%free()
784 : ! call gstore%init(gstore_path, dtset, dtfil, wfk0_hdr, cryst, ks_ebands, ifc, comm)
785 : ! call gstore%free()
786 :
787 : case (12, -12)
788 : ! Migdal-Eliashberg equations (isotropic or anisotropic case).
789 : ! Read |g(k,q)|^2 from GSTORE when present. With Wannier input, prepare
790 : ! the compact real-space vertex and interpolate it on demand in get_a2fw.
791 : call gstore%init_or_from_ncpath(merge(with_cplex1, with_cplex0, dtfil%filgstorein /= ABI_NOFILE), &
792 : dtset, dtfil, wfk0_hdr, cryst, qp_ebands, ifc, &
793 0 : "phonon", dtset%gstore_gname, .False., comm, gstore_from_file)
794 0 : if (gstore_from_file) then
795 0 : call wrtout(units, " Gstore built by reading a pre-existent GSTORE.nc file")
796 : else
797 0 : call wrtout(units, " Gstore built on the fly via Wannier interpolation (ABIWAN.nc + GWAN.nc)")
798 : end if
799 :
800 0 : if (dtset%eph_task == -12) call migdal_eliashberg_iso(gstore, dtset, dtfil)
801 : !if (dtset%eph_task == +12) call migdal_eliashberg_aniso(gstore, dtset, dtfil)
802 0 : call gstore%free()
803 :
804 : case (13)
805 : ! Variational polaron equations.
806 : call gstore%init_or_from_ncpath(with_cplex2, dtset, dtfil, wfk0_hdr, cryst, qp_ebands, ifc, &
807 9 : "phonon", dtset%gstore_gname, .False., comm, gstore_from_file)
808 9 : if (gstore_from_file) then
809 9 : call wrtout(units, sjoin(" Computing variational polaron equations from pre-existent GSTORE file:", dtfil%filgstorein))
810 : else
811 0 : call wrtout(units, " Gstore built on the fly via Wannier interpolation (ABIWAN.nc + GWAN.nc)")
812 : end if
813 9 : call varpeq_run(gstore, dtset, dtfil)
814 9 : call gstore%free()
815 :
816 : case (-13)
817 : ! Compute polaron wavefunctions and atomic displacements in the supercell and write results to files.
818 2 : call varpeq_plot(wfk0_path, ngfftc, dtset, dtfil, cryst, qp_ebands, pawtab, psps, comm)
819 :
820 : case (14)
821 : ! Molecular Berry Curvature.
822 : call gstore%init_or_from_ncpath(with_cplex2, dtset, dtfil, wfk0_hdr, cryst, qp_ebands, ifc, &
823 0 : "atom", dtset%gstore_gname, .False., comm, gstore_from_file)
824 0 : call wrtout(units, " Computing Berry curvature")
825 0 : if (gstore_from_file) then
826 0 : call wrtout(units, " Gstore built by reading a pre-existent GSTORE.nc file")
827 : else
828 0 : call wrtout(units, " Gstore built on the fly via Wannier interpolation (ABIWAN.nc + GWAN.nc)")
829 : end if
830 :
831 0 : call berry_curvature(gstore, dtset, dtfil)
832 0 : call gstore%free()
833 :
834 : case (15, -15)
835 : ! Write average of DFPT potentials to file.
836 2 : if (nprocs > 1) then
837 0 : ABI_WARNING("eph_task in [15, -15] (average of DFPT potentials) does not support nprocs > 1. Running in sequential.")
838 : end if
839 2 : dvdb%comm = xmpi_comm_self
840 2 : if (my_rank == master) then
841 2 : call dvdb%open_read(ngfftf, xmpi_comm_self)
842 2 : call dvdb%write_v1qavg(dtset, strcat(dtfil%filnam_ds(4), "_V1QAVG.nc"))
843 : end if
844 :
845 : case (-16, 16)
846 0 : if (nprocs > 1) then
847 0 : ABI_WARNING("eph_task in [16, -16] (test_phrotation) does not support nprocs > 1. Running in sequential.")
848 : end if
849 :
850 0 : qptopt = qp_ebands%kptopt; if (dtset%qptopt /= 0) qptopt = dtset%qptopt
851 0 : call test_phrotation(ifc, cryst, qptopt, dtset%ph_ngqpt, comm)
852 :
853 0 : dvdb%comm = xmpi_comm_self
854 0 : if (my_rank == master) then
855 0 : call dvdb%open_read(ngfftf, xmpi_comm_self)
856 : ! Compute \delta V_{q,nu)(r) and dump results to netcdf file.
857 0 : call ncwrite_v1qnu(dvdb, dtset, ifc, strcat(dtfil%filnam_ds(4), "_V1QNU.nc"))
858 : end if
859 :
860 : case (17)
861 : ! Compute e-ph matrix elements with the GWPT formalism.
862 : call gwpt_run(wfk0_path, dtfil, ngfftc, ngfftf, dtset, cryst, qp_ebands, dvdb, drhodb, ifc, wfk0_hdr, &
863 3 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
864 :
865 3 : if (len(trim(dtset%gstore_convert)) /= 0) then
866 0 : gstore_path = strcat(dtfil%filnam_ds(4), "_GSTORE.nc")
867 0 : call gstore_convert(gstore_path, dtset, dtfil, cryst, qp_ebands, ifc, comm)
868 : end if
869 :
870 : case (18)
871 : ! Compute e-ph matrix elements along path in the BZ.
872 2 : call eph_path_run(dtfil, dtset, cryst, ks_ebands, dvdb, ifc, pawfgr, pawang, pawrad, pawtab, psps, comm)
873 :
874 : case (19)
875 : ! Compute matrix elements of W_kk'.
876 :
877 0 : call wkk_run(wfk0_path, dtfil, ngfftc, ngfftf, dtset, cryst, qp_ebands, wfk0_hdr, pawtab, psps, mpi_enreg, comm)
878 :
879 : case (20)
880 0 : call gstore_convert(dtfil%filgstorein, dtset, dtfil, cryst, qp_ebands, ifc, comm)
881 :
882 : case default
883 109 : ABI_ERROR(sjoin("Unsupported value of eph_task:", itoa(dtset%eph_task)))
884 : end select
885 :
886 : !=====================
887 : !==== Free memory ====
888 : !=====================
889 109 : call cryst%free(); call dvdb%free(); call drhodb%free(); call ddb_hdr%free()
890 109 : call ddb%free(); call ifc%free(); call wfk0_hdr%free()
891 109 : call ks_ebands%free(); call ks_ebands_kq%free(); call qp_ebands%free()
892 109 : call pawfgr_destroy(pawfgr); call destroy_mpi_enreg(mpi_enreg)
893 :
894 109 : if (allocated(efmasdeg)) call efmasdeg_free_array(efmasdeg)
895 109 : if (allocated(efmasval)) call efmasval_free_array(efmasval)
896 109 : ABI_SFREE(kpt_efmas)
897 :
898 : ! Deallocation for PAW.
899 : if (dtset%usepaw == 1) then
900 : !call pawrhoij_free(pawrhoij)
901 : !ABI_FREE(pawrhoij)
902 : !call pawfgrtab_free(pawfgrtab)
903 : !ABI_FREE(pawfgrtab)
904 : !call paw_ij_free(paw_ij)
905 : !ABI_FREE(paw_ij)
906 : !call paw_an_free(paw_an)
907 : !ABI_FREE(paw_an)
908 : end if
909 :
910 : DBG_EXIT('COLL')
911 :
912 545 : end subroutine eph
913 : !!***
914 :
915 : end module m_eph_driver
916 : !!***
|
