Branch data Line data Source code
1 : : !!****m* ABINIT/m_phpi
2 : : !! NAME
3 : : !!
4 : : !! FUNCTION
5 : : !! Computation of phonon-electron self-energy.
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group (GKA)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : module m_phpi
22 : :
23 : : use, intrinsic :: iso_c_binding
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_xmpi
27 : : use m_errors
28 : : use m_ifc
29 : : use m_ebands
30 : : use m_nctk
31 : : use netcdf
32 : : use m_wfk
33 : : use m_ddb
34 : : use m_dvdb
35 : : use m_fft
36 : : use m_hamiltonian
37 : : use m_pawcprj
38 : : use m_dtset
39 : : use m_dtfil
40 : :
41 : : use defs_datatypes, only : pseudopotential_type
42 : : use defs_abitypes, only : mpi_type
43 : : use m_time, only : cwtime
44 : : use m_fstrings, only : sjoin, itoa, ftoa, ktoa, ltoa, strcat
45 : : use m_io_tools, only : iomode_from_fname
46 : : use m_cgtools, only : dotprod_g
47 : : use m_kg, only : getph
48 : : use m_fftcore, only : get_kg
49 : : use m_crystal, only : crystal_t
50 : : use m_bz_mesh, only : findqg0
51 : : use m_wfd, only : wfd_t
52 : : use m_pawang, only : pawang_type
53 : : use m_pawrad, only : pawrad_type
54 : : use m_pawtab, only : pawtab_type
55 : : use m_pawfgr, only : pawfgr_type
56 : : use m_getgh1c, only : getgh1c, rf_transgrid_and_pack, getgh1c_setup
57 : :
58 : : implicit none
59 : :
60 : : private
61 : : !!***
62 : :
63 : : public :: eph_phpi
64 : :
65 : :
66 : : contains !=================================================================================
67 : : !!***
68 : :
69 : : !!****f* m_phpi/eph_phpi
70 : : !! NAME
71 : : !! eph_phpi
72 : : !!
73 : : !! FUNCTION
74 : : !! Compute phonon-electron self-energy.
75 : : !!
76 : : !! INPUTS
77 : : !! wk0_path=String with the path to the GS unperturbed WFK file.
78 : : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
79 : : !! dtset<dataset_type>=All input variables for this dataset.
80 : : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
81 : : !! dvdb<dbdb_type>=Database with the DFPT SCF potentials.
82 : : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
83 : : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
84 : : !! pawang<pawang_type)>=PAW angular mesh and related data.
85 : : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
86 : : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
87 : : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
88 : : !! comm=MPI communicator.
89 : : !!
90 : : !! OUTPUT
91 : : !!
92 : : !! SOURCE
93 : :
94 : 0 : subroutine eph_phpi(wfk0_path,wfq_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands_k,ebands_kq,dvdb,ifc,&
95 : 0 : pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm)
96 : :
97 : : !Arguments ------------------------------------
98 : : !scalars
99 : : character(len=*),intent(in) :: wfk0_path, wfq_path
100 : : integer,intent(in) :: comm
101 : : type(datafiles_type),intent(in) :: dtfil
102 : : type(dataset_type),intent(in) :: dtset
103 : : type(crystal_t),intent(in) :: cryst
104 : : type(ebands_t),intent(in) :: ebands_k, ebands_kq
105 : : type(dvdb_t),intent(inout) :: dvdb
106 : : type(pawang_type),intent(in) :: pawang
107 : : type(pseudopotential_type),intent(in) :: psps
108 : : type(pawfgr_type),intent(in) :: pawfgr
109 : : type(ifc_type),intent(in) :: ifc
110 : : type(mpi_type),intent(in) :: mpi_enreg
111 : : !arrays
112 : : integer,intent(in) :: ngfft(18),ngfftf(18)
113 : : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
114 : : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
115 : :
116 : : !Local variables ------------------------------
117 : : !scalars
118 : : integer,parameter :: tim_getgh1c = 1,berryopt0 = 0, useylmgr1 = 0, master = 0
119 : : integer :: my_rank,nproc,iomode,mband,mband_kq,my_minb,my_maxb,nsppol,nkpt,nkpt_kq,idir,ipert
120 : : integer :: cplex,db_iqpt,natom,natom3,ipc,nspinor,onpw,imode
121 : : integer :: ib1,ib2,ik,ikq,spin,istwf_k,istwf_kq,npw_k,npw_kq
122 : : integer :: mpw,my_mpw,ierr,my_kstart,my_kstop,cnt
123 : : integer :: n1,n2,n3,n4,n5,n6,nspden
124 : : integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnlx1
125 : : integer :: nfft,nfftf,mgfft,mgfftf,nkpg,nkpg1
126 : : real(dp) :: cpu,wall,gflops
127 : : real(dp) :: ecut,eshift,eig0nk,eig0mkq,dotr,doti
128 : : real(dp) :: eta,f_nk,f_mkq,omega,wtk,gkk2,term1,term2
129 : : logical :: i_am_master,gen_eigenpb
130 : 0 : type(wfd_t) :: wfd_k,wfd_kq
131 : 0 : type(gs_hamiltonian_type) :: gs_hamkq
132 : 0 : type(rf_hamiltonian_type) :: rf_hamkq
133 : : character(len=500) :: msg
134 : : !arrays
135 : : integer :: g0_k(3), units(2)
136 : 0 : integer,allocatable :: kg_k(:,:),kg_kq(:,:),gtmp(:,:),nband(:,:),nband_kq(:,:),blkflg(:,:), wfd_istwfk(:)
137 : 0 : real(dp) :: kk(3),kq(3),qpt(3),phfrq(3*cryst%natom)
138 : 0 : real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom),displ_red(2,3,cryst%natom,3*cryst%natom)
139 : 0 : real(dp) :: Pi_ph(3*cryst%natom)
140 : 0 : real(dp),allocatable :: grad_berry(:,:),kinpw_kq(:),kpg_kq(:,:),kpg_k(:,:),dkinpw(:)
141 : 0 : real(dp),allocatable :: ffnl_k(:,:,:,:),ffnl_kq(:,:,:,:),ph3d_k(:,:,:),ph3d_kq(:,:,:)
142 : 0 : real(dp),allocatable :: v1scf(:,:,:,:),gkk(:,:,:,:,:), gkk_m(:,:,:)
143 : 0 : real(dp),allocatable :: bras_kq(:,:,:),kets_k(:,:,:),h1kets_kq(:,:,:)
144 : 0 : real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:)
145 : 0 : real(dp),allocatable :: ylm_kq(:,:),ylm_k(:,:),ylmgr_kq(:,:,:)
146 : 0 : real(dp),allocatable :: dummy_vtrial(:,:),gvnlx1(:,:)
147 : 0 : real(dp),allocatable :: gs1c(:,:)
148 : 0 : logical,allocatable :: bks_mask(:,:,:),bks_mask_kq(:,:,:),keep_ur(:,:,:),keep_ur_kq(:,:,:)
149 : 0 : type(pawcprj_type),allocatable :: cwaveprj0(:,:) !natom,nspinor*usecprj)
150 : : !************************************************************************
151 : :
152 [ # # ]: 0 : units = [std_out, ab_out]
153 : 0 : write(msg, '(3a)') ch10, "Computation of the real part of the phonon self-energy", ch10
154 : 0 : call wrtout(units, msg, do_flush=.True.)
155 : :
156 [ # # ]: 0 : if (psps%usepaw == 1) then
157 : 0 : ABI_ERROR("PAW not implemented")
158 : : ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
159 : : end if
160 : :
161 : 0 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
162 : 0 : i_am_master = (my_rank == master)
163 : :
164 : : ! Copy important dimensions
165 : 0 : natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands_k%nsppol; nspinor = ebands_k%nspinor; nspden = dtset%nspden
166 : 0 : nkpt = ebands_k%nkpt
167 : 0 : mband = ebands_k%mband
168 : 0 : nkpt_kq = ebands_kq%nkpt
169 : 0 : mband_kq = ebands_kq%mband
170 : 0 : ecut = dtset%ecut
171 : :
172 : : ! GKA TODO: Make sure there is a single q-point present.
173 [ # # ]: 0 : qpt = dtset%qptn(:)
174 : :
175 [ # # ]: 0 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
176 [ # # # # ]: 0 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
177 : 0 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
178 : 0 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
179 : :
180 : : ! Open the DVDB file
181 : 0 : call dvdb%open_read(ngfftf, xmpi_comm_self)
182 : :
183 : : ! Initialize the wave function descriptors.
184 : : ! For the time being, no memory distribution, each node has the full set of states.
185 : 0 : my_minb = 1; my_maxb = mband
186 : :
187 [ # # # # : 0 : ABI_MALLOC(nband, (nkpt, nsppol))
# # # # #
# # # # #
# # ]
188 [ # # # # : 0 : ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
# # # # #
# # # # #
# # # # #
# ]
189 [ # # # # : 0 : ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol))
# # # # #
# # # # #
# # ]
190 [ # # # # : 0 : nband=mband; bks_mask=.False.; keep_ur=.False.
# # # # #
# # # # #
# # ]
191 : :
192 [ # # # # : 0 : ABI_MALLOC(nband_kq, (nkpt_kq, nsppol))
# # # # #
# # # # #
# # ]
193 [ # # # # : 0 : ABI_MALLOC(bks_mask_kq,(mband_kq, nkpt_kq, nsppol))
# # # # #
# # # # #
# # # # #
# ]
194 [ # # # # : 0 : ABI_MALLOC(keep_ur_kq,(mband_kq, nkpt_kq ,nsppol))
# # # # #
# # # # #
# # ]
195 [ # # # # : 0 : nband_kq=mband_kq; bks_mask_kq=.False.; keep_ur_kq=.False.
# # # # #
# # # # #
# # ]
196 : :
197 : : ! Distribute the k-points over the processors
198 : 0 : call xmpi_split_work(nkpt,comm,my_kstart,my_kstop)
199 [ # # ]: 0 : do ik=1,nkpt
200 [ # # # # ]: 0 : if (.not. ((ik .ge. my_kstart) .and. (ik .le. my_kstop))) cycle
201 [ # # ]: 0 : kk = ebands_k%kptns(:,ik)
202 [ # # ]: 0 : kq = kk + qpt
203 : 0 : call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:),(/1,1,1/)) ! Find the index of the k+q point
204 [ # # # # ]: 0 : bks_mask(:,ik,:) = .True.
205 [ # # # # ]: 0 : bks_mask_kq(:,ikq,:) = .True.
206 : : end do
207 : :
208 : : ! Initialize the wavefunction descriptors
209 : :
210 : : ! Impose istwfk=1 for all k points. This is also done in respfn (see inkpts)
211 : : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
212 [ # # # # : 0 : ABI_MALLOC(wfd_istwfk, (nkpt))
# # # # #
# # # ]
213 [ # # ]: 0 : wfd_istwfk = 1
214 : :
215 : : call wfd_k%init(cryst,pawtab,psps,keep_ur,mband,nband,nkpt,nsppol,bks_mask,&
216 : : nspden,nspinor,ecut,dtset%ecutsm,dtset%dilatmx,wfd_istwfk,ebands_k%kptns,ngfft,&
217 : 0 : dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
218 [ # # ]: 0 : ABI_FREE(wfd_istwfk)
219 : :
220 [ # # ]: 0 : call wfd_k%print([std_out], header="Wavefunctions on the k-points grid")
221 : :
222 [ # # # # : 0 : ABI_MALLOC(wfd_istwfk, (nkpt_kq))
# # # # #
# # # ]
223 [ # # ]: 0 : wfd_istwfk = 1
224 : :
225 : : call wfd_kq%init(cryst,pawtab,psps,keep_ur_kq,mband_kq,nband_kq,nkpt_kq,nsppol,bks_mask_kq,&
226 : : nspden,nspinor,ecut,dtset%ecutsm,dtset%dilatmx,wfd_istwfk,ebands_kq%kptns,ngfft,&
227 : 0 : dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm)
228 : :
229 [ # # ]: 0 : ABI_FREE(wfd_istwfk)
230 : :
231 [ # # ]: 0 : call wfd_kq%print([std_out], header="Wavefunctions on the q-shifted k-points grid")
232 : :
233 [ # # ]: 0 : ABI_FREE(nband)
234 [ # # ]: 0 : ABI_FREE(bks_mask)
235 [ # # ]: 0 : ABI_FREE(keep_ur)
236 [ # # ]: 0 : ABI_FREE(nband_kq)
237 [ # # ]: 0 : ABI_FREE(bks_mask_kq)
238 [ # # ]: 0 : ABI_FREE(keep_ur_kq)
239 : :
240 : : ! Read wavefunctions on the k-points grid and q-shifted k-points grid.
241 : 0 : iomode = iomode_from_fname(wfk0_path)
242 : 0 : call wfd_k%read_wfk(wfk0_path,iomode)
243 : : if (.False.) call wfd_k%test_ortho(cryst,pawtab,unit=std_out,mode_paral="PERS")
244 : :
245 : 0 : iomode = iomode_from_fname(wfq_path)
246 : 0 : call wfd_kq%read_wfk(wfq_path,iomode)
247 : : if (.False.) call wfd_kq%test_ortho(cryst,pawtab,unit=std_out,mode_paral="PERS")
248 : :
249 : : ! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information on the coarse grid.
250 [ # # # # : 0 : ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
# # # # #
# # # #
# ]
251 : 0 : call getph(cryst%atindx,natom,n1,n2,n3,ph1d,cryst%xred)
252 : :
253 : : ! Find the appropriate value of mpw
254 : 0 : mpw = 0; cnt=0
255 [ # # ]: 0 : do spin=1,nsppol
256 [ # # ]: 0 : do ik=1,nkpt
257 [ # # ]: 0 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
258 [ # # ]: 0 : kk = ebands_k%kptns(:,ik)
259 [ # # ]: 0 : call get_kg(kk,1,ecut,cryst%gmet,onpw,gtmp)
260 [ # # ]: 0 : ABI_FREE(gtmp)
261 : 0 : mpw = max(mpw, onpw)
262 : : end do
263 : : end do
264 : 0 : cnt=0
265 [ # # ]: 0 : do spin=1,nsppol
266 [ # # ]: 0 : do ikq=1,nkpt_kq
267 [ # # ]: 0 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
268 [ # # ]: 0 : kq = ebands_kq%kptns(:,ikq)
269 [ # # ]: 0 : call get_kg(kq,1,ecut,cryst%gmet,onpw,gtmp)
270 [ # # ]: 0 : ABI_FREE(gtmp)
271 : 0 : mpw = max(mpw, onpw)
272 : : end do
273 : : end do
274 : 0 : my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
275 : :
276 : : ! Allow PW-arrays dimensioned with mpw
277 [ # # # # : 0 : ABI_MALLOC(kg_k, (3, mpw))
# # # # #
# # # ]
278 [ # # # # : 0 : ABI_MALLOC(kg_kq, (3, mpw))
# # # # ]
279 : :
280 : : ! Spherical Harmonics for useylm==1.
281 : : ! FIXME: These arrays should allocated with npw_k and npw_kq
282 [ # # # # : 0 : ABI_MALLOC(ylm_k,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
# # # # #
# # # # #
# # ]
283 [ # # # # : 0 : ABI_MALLOC(ylm_kq,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
# # # # #
# # # ]
284 [ # # # # : 0 : ABI_MALLOC(ylmgr_kq,(mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
# # # # ]
285 : :
286 : : ! TODO FOR PAW
287 : 0 : usecprj = 0
288 [ # # # # : 0 : ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
# # # # #
# ]
289 : :
290 : : ! Prepare call to getgh1c
291 : 0 : usevnl = 0
292 : 0 : optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
293 : 0 : optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
294 : 0 : opt_gvnlx1 = 0 ! gvnlx1 is output
295 [ # # ]: 0 : ABI_MALLOC(gvnlx1, (2,usevnl))
296 [ # # ]: 0 : ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4)))
297 : :
298 : : ! This part is taken from dfpt_vtorho
299 : : !==== Initialize most of the Hamiltonian (and derivative) ====
300 : : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
301 : : !2) Perform the setup needed for the non-local factors:
302 : : !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
303 : : !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
304 : :
305 : : call gs_hamkq%init(psps,pawtab,nspinor,NSPPOL,nspden,natom,&
306 : : dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,&
307 : : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
308 : 0 : usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option)
309 : :
310 : : ! Allocate vlocal. Note nvloc
311 : : ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data)
312 [ # # # # : 0 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
313 [ # # # # : 0 : vlocal = huge(one)
# # # # ]
314 : :
315 : : ! Allocate work space arrays.
316 [ # # # # : 0 : ABI_MALLOC(blkflg, (natom3,natom3))
# # # # #
# # # # #
# # ]
317 [ # # # # : 0 : ABI_CALLOC(dummy_vtrial, (nfftf,nspden))
# # # # #
# # # # #
# # # # #
# ]
318 : :
319 : 0 : call cwtime(cpu,wall,gflops,"start")
320 : :
321 : : ! Find the index of the q-point in the DVDB.
322 : 0 : db_iqpt = dvdb%findq(qpt)
323 : :
324 [ # # ]: 0 : if (db_iqpt /= -1) then
325 [ # # ]: 0 : if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found: ",ktoa(qpt)," in DVDB with index ",itoa(db_iqpt)))
326 : : ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
327 : : ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
328 : 0 : call dvdb%readsym_allv1(db_iqpt, cplex, nfftf, ngfftf, v1scf, comm)
329 : : else
330 : 0 : ABI_ERROR(sjoin("Could not find symmetric of q-point:", ktoa(qpt), "in DVDB"))
331 : : end if
332 : :
333 : : ! Allocate vlocal1 with correct cplex. Note nvloc
334 [ # # # # : 0 : ABI_MALLOC_OR_DIE(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
335 : :
336 : : ! Allocate el-ph coupling matrix elements
337 [ # # # # : 0 : ABI_MALLOC(gkk, (2, mband_kq, mband, natom, 3))
# # # # #
# # # # #
# # # # #
# # # ]
338 [ # # # # : 0 : ABI_MALLOC(gkk_m, (2, mband_kq, mband))
# # # # #
# # # # #
# # ]
339 : :
340 : : ! Compute displacement vectors and phonon frequencies
341 : 0 : call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
342 : :
343 : : ! Broadening parameter
344 [ # # ]: 0 : if (dtset%elph2_imagden .gt. tol12) then
345 : 0 : eta = dtset%elph2_imagden
346 : : else
347 : : eta = 0.0001_dp
348 : : end if
349 : :
350 : : ! Kpoints weights (not using symmetries at the moment)
351 : 0 : wtk = 1.0 / nkpt
352 : :
353 : : ! Initialize phonon self-energy
354 [ # # ]: 0 : Pi_ph = zero
355 : :
356 : : ! Examine the symmetries of the q wavevector
357 : : ! call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol)
358 : :
359 : : ! ----------------------------------------------------------------------------------------------- !
360 : : ! Begin loop over states
361 : : ! ----------------------------------------------------------------------------------------------- !
362 [ # # ]: 0 : do spin=1,nsppol
363 : :
364 : : ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
365 [ # # ]: 0 : do ipc=1,natom3
366 : : call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,&
367 : 0 : pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc))
368 : : end do
369 : :
370 : : ! Continue to initialize the Hamiltonian
371 : 0 : call gs_hamkq%load_spin(spin,vlocal=vlocal,with_nonlocal=.true.)
372 : :
373 [ # # ]: 0 : do ik=1,nkpt
374 : : ! Only do a subset a k-points
375 [ # # # # ]: 0 : if (.not. ((ik .ge. my_kstart) .and. (ik .le. my_kstop))) cycle
376 : :
377 : : ! Allocate workspace for wavefunctions. Make npw larger than expected.
378 [ # # # # : 0 : ABI_MALLOC(bras_kq, (2, mpw*nspinor, mband))
# # # # #
# # # # #
# # ]
379 [ # # # # : 0 : ABI_MALLOC(kets_k, (2, mpw*nspinor, mband))
# # # # #
# # # ]
380 [ # # # # : 0 : ABI_MALLOC(h1kets_kq, (2, mpw*nspinor, mband))
# # # # #
# # # ]
381 : :
382 [ # # ]: 0 : kk = ebands_k%kptns(:,ik)
383 [ # # ]: 0 : kq = kk + qpt
384 : 0 : call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:),(/1,1,1/)) ! Find the index of the k+q point
385 : :
386 : : ! Copy u_k(G)
387 : 0 : istwf_k = wfd_k%istwfk(ik); npw_k = wfd_k%npwarr(ik)
388 [ # # ]: 0 : ABI_CHECK(mpw >= npw_k, "mpw < npw_k")
389 [ # # # # ]: 0 : kg_k(:,1:npw_k) = wfd_k%kdata(ik)%kg_k
390 [ # # ]: 0 : do ib2=1,mband
391 : 0 : call wfd_k%copy_cg(ib2, ik, spin, kets_k(1,1,ib2))
392 : : end do
393 : :
394 : : ! Copy u_kq(G)
395 : 0 : istwf_kq = wfd_kq%istwfk(ikq); npw_kq = wfd_kq%npwarr(ikq)
396 [ # # ]: 0 : ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
397 [ # # # # ]: 0 : kg_kq(:,1:npw_kq) = wfd_kq%kdata(ikq)%kg_k
398 [ # # ]: 0 : do ib1=1,mband_kq
399 : 0 : call wfd_kq%copy_cg(ib1, ikq, spin, bras_kq(1,1,ib1))
400 : : end do
401 : :
402 : : ! if PAW, one has to solve a generalized eigenproblem
403 : : ! BE careful here because I will need sij_opt==-1
404 : 0 : gen_eigenpb = (psps%usepaw==1)
405 [ # # ]: 0 : sij_opt = 0; if (gen_eigenpb) sij_opt = 1
406 [ # # # # : 0 : ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2)))
# # # # #
# # # #
# ]
407 : :
408 [ # # # # : 0 : gkk = zero
# # # # #
# ]
409 : :
410 : : ! Loop over all 3*natom perturbations.
411 [ # # ]: 0 : do ipc=1,natom3
412 : 0 : idir = mod(ipc-1, 3) + 1
413 : 0 : ipert = (ipc - idir) / 3 + 1
414 : :
415 : : !write(msg, '(a,2i4)') "Treating ipert, idir = ", ipert, idir
416 : : !call wrtout(std_out, msg, do_flush=.True.)
417 : :
418 : : ! Prepare application of the NL part.
419 : 0 : call rf_hamkq%init(cplex,gs_hamkq,ipert,has_e1kbsc=.true.)
420 : 0 : call rf_hamkq%load_spin(spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.)
421 : :
422 : : ! This call is not optimal because there are quantities in out that do not depend on idir,ipert
423 : : call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,& ! In
424 : : cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,& ! In
425 : : npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq,& ! In
426 [ # # # # : 0 : dkinpw,nkpg,nkpg1,kpg_k,kpg_kq,kinpw_kq,ffnl_k,ffnl_kq,ph3d_k,ph3d_kq) ! Out
# # # # ]
427 : :
428 : : ! Calculate dvscf * psi_k, results stored in h1kets_kq on the k+q sphere.
429 : : ! Compute H(1) applied to GS wavefunction Psi(0)
430 [ # # ]: 0 : do ib2=1,mband
431 : 0 : eig0nk = ebands_k%eig(ib2,ik,spin)
432 : : ! Use scissor shift on 0-order eigenvalue
433 : 0 : eshift = eig0nk - dtset%dfpt_sciss
434 : :
435 : : call getgh1c(berryopt0,kets_k(:,:,ib2),cwaveprj0,h1kets_kq(:,:,ib2),&
436 : : & grad_berry,gs1c,gs_hamkq,gvnlx1,idir,ipert,(/eshift/),mpi_enreg,1,optlocal,&
437 [ # # ]: 0 : & optnl,opt_gvnlx1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
438 : : end do
439 : :
440 [ # # ]: 0 : ABI_FREE(kinpw_kq)
441 [ # # ]: 0 : ABI_FREE(kpg_kq)
442 [ # # ]: 0 : ABI_FREE(kpg_k)
443 [ # # ]: 0 : ABI_FREE(dkinpw)
444 [ # # ]: 0 : ABI_FREE(ffnl_k)
445 [ # # ]: 0 : ABI_FREE(ffnl_kq)
446 [ # # ]: 0 : ABI_FREE(ph3d_k)
447 [ # # ]: 0 : ABI_SFREE(ph3d_kq)
448 [ # # ]: 0 : ABI_SFREE(gs1c)
449 : :
450 : : ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation.
451 : : !The array eig1_k contains:
452 : : !
453 : : ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)> (NC psps)
454 : : ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW)
455 [ # # ]: 0 : do ib2=1,mband
456 [ # # ]: 0 : do ib1=1,mband_kq
457 : : call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bras_kq(1,1,ib1),h1kets_kq(1,1,ib2),&
458 : 0 : mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
459 [ # # ]: 0 : gkk(:,ib1,ib2,ipert,idir) = [dotr, doti]
460 : : end do
461 : : end do
462 : :
463 : : end do ! ipc
464 : :
465 : : ! Loop over 3*natom phonon branches.
466 [ # # ]: 0 : do imode=1,natom3
467 : 0 : omega = phfrq(imode)
468 : : ! Do not compute Pi for negative or too small frequencies
469 [ # # ]: 0 : if (omega .lt. tol6) cycle
470 [ # # # # : 0 : gkk_m = zero
# # ]
471 : :
472 : : ! Transform the gkk from atom,cart basis to mode basis
473 [ # # ]: 0 : do idir=1,3
474 [ # # ]: 0 : do ipert=1,natom
475 : : gkk_m(1,:,:) = gkk_m(1,:,:) &
476 : : + gkk(1,:,:,ipert,idir) * displ_red(1,idir,ipert,imode) &
477 [ # # # # ]: 0 : - gkk(2,:,:,ipert,idir) * displ_red(2,idir,ipert,imode)
478 : : gkk_m(2,:,:) = gkk_m(2,:,:) &
479 : : + gkk(1,:,:,ipert,idir) * displ_red(2,idir,ipert,imode) &
480 [ # # # # ]: 0 : + gkk(2,:,:,ipert,idir) * displ_red(1,idir,ipert,imode)
481 : : end do
482 : : end do
483 : :
484 [ # # # # : 0 : gkk_m = gkk_m / sqrt(two * omega)
# # ]
485 : :
486 : : ! sum contribution to phonon self-energy
487 [ # # ]: 0 : do ib2=1,mband
488 [ # # ]: 0 : do ib1=1,mband_kq
489 : 0 : f_nk = ebands_k%occ(ib2,ik,spin)
490 : 0 : f_mkq = ebands_kq%occ(ib1,ikq,spin)
491 [ # # ]: 0 : if (abs(f_mkq - f_nk) .le. tol12) cycle
492 : :
493 : 0 : eig0nk = ebands_k%eig(ib2,ik,spin)
494 : 0 : eig0mkq = ebands_kq%eig(ib1,ikq,spin)
495 : :
496 : 0 : gkk2 = gkk_m(1,ib1,ib2) ** 2 + gkk_m(2,ib1,ib2) ** 2
497 : :
498 : 0 : term1 = (f_mkq - f_nk) * (eig0mkq - eig0nk - omega) / ((eig0mkq - eig0nk - omega) ** 2 + eta ** 2)
499 : 0 : term2 = (f_mkq - f_nk) * (eig0mkq - eig0nk ) / ((eig0mkq - eig0nk ) ** 2 + eta ** 2)
500 : :
501 : 0 : Pi_ph(imode) = Pi_ph(imode) + wtk * gkk2 * (term1 - term2)
502 : : end do
503 : : end do
504 : :
505 : : end do ! imode
506 : :
507 [ # # ]: 0 : ABI_FREE(bras_kq)
508 [ # # ]: 0 : ABI_FREE(kets_k)
509 [ # # ]: 0 : ABI_FREE(h1kets_kq)
510 : : end do ! ikfs
511 : :
512 : 0 : call rf_hamkq%free()
513 : : end do ! spin
514 : :
515 : : ! Gather the k-points computed by all processes
516 : 0 : call xmpi_sum_master(Pi_ph,master,comm,ierr)
517 : :
518 : : ! Output the results
519 [ # # ]: 0 : if (i_am_master) then
520 : 0 : call out_phpi(ab_out, Pi_ph, phfrq, qpt, natom3)
521 : 0 : call out_phpi(std_out, Pi_ph, phfrq, qpt, natom3)
522 : : end if
523 : :
524 : 0 : if (i_am_master) call out_phpi_nc(dtfil, cryst, Pi_ph, phfrq, qpt, natom3)
525 : :
526 : : ! Free memory
527 : 0 : call cwtime(cpu,wall,gflops,"stop")
528 : :
529 : 0 : write(msg, '(3a)') "Computation of the real part of the phonon self-energy completed", ch10, &
530 : 0 : "--------------------------------------------------------------------------------"
531 : 0 : call wrtout(units, msg, do_flush=.True.)
532 : :
533 : : ! Free memory
534 [ # # ]: 0 : ABI_FREE(gkk)
535 [ # # ]: 0 : ABI_FREE(gkk_m)
536 [ # # ]: 0 : ABI_FREE(v1scf)
537 [ # # ]: 0 : ABI_FREE(vlocal1)
538 [ # # ]: 0 : ABI_FREE(gvnlx1)
539 [ # # ]: 0 : ABI_FREE(grad_berry)
540 [ # # ]: 0 : ABI_FREE(dummy_vtrial)
541 [ # # ]: 0 : ABI_FREE(ph1d)
542 [ # # ]: 0 : ABI_FREE(vlocal)
543 [ # # ]: 0 : ABI_FREE(kg_k)
544 [ # # ]: 0 : ABI_FREE(kg_kq)
545 [ # # ]: 0 : ABI_FREE(ylm_k)
546 [ # # ]: 0 : ABI_FREE(ylm_kq)
547 [ # # ]: 0 : ABI_FREE(ylmgr_kq)
548 [ # # ]: 0 : ABI_FREE(blkflg)
549 : :
550 : 0 : call gs_hamkq%free()
551 : 0 : call wfd_k%free()
552 : 0 : call wfd_kq%free()
553 : :
554 : 0 : call pawcprj_free(cwaveprj0)
555 [ # # # # : 0 : ABI_FREE(cwaveprj0)
# # # # #
# ]
556 : :
557 [ # # # # : 0 : end subroutine eph_phpi
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
558 : : !!***
559 : :
560 : : !----------------------------------------------------------------------
561 : :
562 : : !!****f* m_phpi/out_phpi
563 : : !! NAME
564 : : !! out_phpi
565 : : !!
566 : : !! FUNCTION
567 : : !! Output the phonon self-energy.
568 : : !!
569 : : !! INPUTS
570 : : !!
571 : : !! OUTPUT
572 : : !!
573 : : !! NOTES
574 : : !!
575 : : !! SOURCE
576 : :
577 : 0 : subroutine out_phpi(iout, Pi_ph, phfrq, qpt, natom3)
578 : :
579 : : !Arguments ------------------------------------
580 : : !scalars
581 : : integer,intent(in) :: iout
582 : : integer,intent(in) :: natom3
583 : : !arrays
584 : : real(dp),intent(in) :: Pi_ph(natom3),phfrq(natom3),qpt(3)
585 : :
586 : : !Local variables ------------------------------
587 : : !scalars
588 : : integer :: imode
589 : :
590 : 0 : write(iout,'(a)')' '
591 : : !write(iout,'(a)')' ----------------------------------------'
592 : : !write(iout,'(a)')' '
593 : 0 : write(iout,'(a)')' Phonon self-energy (Hartree)'
594 : 0 : write(iout,'(a)')' '
595 : 0 : write(iout,'(a,3f14.8)')' qpt =',qpt
596 : 0 : write(iout,'(a)')' '
597 : 0 : write(iout,'(1x,a,10x,a)')'omega','Pi(omega)'
598 : :
599 [ # # ]: 0 : do imode=1,natom3
600 : 0 : write(iout,'(1x,f12.8,1x,es14.6)') phfrq(imode), Pi_ph(imode)
601 : : end do
602 : :
603 : 0 : write(iout,'(a)')' '
604 : : !write(iout,'(a)')' ----------------------------------------'
605 : :
606 : 0 : end subroutine out_phpi
607 : : !!***
608 : :
609 : : !----------------------------------------------------------------------
610 : :
611 : : !!****f* m_phpi/out_phpi_nc
612 : : !! NAME
613 : : !! out_phpi_nc
614 : : !!
615 : : !! FUNCTION
616 : : !! Output the phonon self-energy in netCDF format.
617 : : !!
618 : : !! INPUTS
619 : : !!
620 : : !! OUTPUT
621 : : !!
622 : : !! NOTES
623 : : !!
624 : : !! SOURCE
625 : :
626 : 0 : subroutine out_phpi_nc(dtfil, cryst, Pi_ph, phfrq, qpt, natom3)
627 : :
628 : : !Arguments ------------------------------------
629 : : !scalars
630 : : integer,intent(in) :: natom3
631 : : type(datafiles_type), intent(in) :: dtfil
632 : : type(crystal_t),intent(in) :: cryst
633 : : !arrays
634 : : real(dp),intent(in) :: Pi_ph(natom3),phfrq(natom3),qpt(3)
635 : :
636 : : !Local variables ------------------------------
637 : : !scalars
638 : : integer :: natom,one_dim,cplex,cart_dir
639 : : integer :: ncid, ncerr
640 : : character(len=fnlen) :: fname
641 : :
642 : : ! Initialize NetCDF file.
643 [ # # ]: 0 : fname = strcat(dtfil%filnam_ds(4),"_Pi.nc")
644 [ # # ]: 0 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
645 : :
646 : : ! Write information of the crystal
647 [ # # ]: 0 : NCF_CHECK(cryst%ncwrite(ncid))
648 : :
649 : : ! Write the dimensions specified by ETSF
650 : 0 : one_dim = 1
651 : 0 : cplex = 2
652 : 0 : cart_dir = 3
653 : :
654 : 0 : natom = natom3 / 3
655 : :
656 : : ncerr = nctk_def_dims(ncid, [&
657 : : nctkdim_t('current_one_dim', one_dim), &
658 : : nctkdim_t('number_of_atoms', natom), &
659 : : nctkdim_t('number_of_cartesian_directions', cart_dir), &
660 : : nctkdim_t('number_of_perturbations', natom3), &
661 [ # # ]: 0 : nctkdim_t('cplex',cplex)], defmode=.True.)
662 [ # # ]: 0 : NCF_CHECK(ncerr)
663 : :
664 : : ! Create the arrays
665 : : ncerr = nctk_def_arrays(ncid, [&
666 : : nctkarr_t('q_point_reduced_coord', "dp", 'number_of_cartesian_directions'),&
667 : : nctkarr_t('phonon_frequencies', "dp", 'number_of_perturbations'), &
668 [ # # ]: 0 : nctkarr_t('phonon_self_energy_realpart', "dp", 'number_of_perturbations')])
669 [ # # ]: 0 : NCF_CHECK(ncerr)
670 : :
671 [ # # ]: 0 : NCF_CHECK(nctk_set_atomic_units(ncid, 'phonon_frequencies'))
672 [ # # ]: 0 : NCF_CHECK(nctk_set_atomic_units(ncid, 'phonon_self_energy_realpart'))
673 : :
674 : : ! Write data
675 [ # # ]: 0 : NCF_CHECK(nctk_set_datamode(ncid))
676 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'q_point_reduced_coord'), qpt))
677 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phonon_frequencies'), phfrq))
678 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phonon_self_energy_realpart'), Pi_ph))
679 : :
680 : : ! Close file
681 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
682 : :
683 : 0 : end subroutine out_phpi_nc
684 : : !!***
685 : :
686 : : !----------------------------------------------------------------------
687 : :
688 : : end module m_phpi
689 : : !!***
|
