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