Branch data 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 : : !!***
|
