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