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