TABLE OF CONTENTS


ABINIT/m_phpi [ Modules ]

[ Top ] [ Modules ]

NAME

FUNCTION

  Tools for the computation of phonon self-energy.

COPYRIGHT

  Copyright (C) 2008-2018 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 .

PARENTS

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 
25 module m_phpi
26 
27  use defs_basis
28  use defs_abitypes
29  use defs_datatypes
30  use m_abicore
31  use m_xmpi
32  use m_errors
33  use m_kptrank
34  use m_tetrahedron
35  use m_ifc
36  use m_ebands
37  use iso_c_binding
38  use m_nctk
39 #ifdef HAVE_NETCDF
40  use netcdf
41 #endif
42  use m_wfk
43  use m_ddb
44  use m_dvdb
45  use m_fft
46  use m_hamiltonian
47  use m_pawcprj
48 
49  use m_time,            only : cwtime
50  use m_fstrings,        only : sjoin, itoa, ftoa, ktoa, ltoa, strcat
51  use m_io_tools,        only : iomode_from_fname
52  use m_cgtools,         only : dotprod_g
53  use m_kg,              only : getph
54  use m_fftcore,         only : get_kg, kpgsph, sphere
55  use m_crystal,         only : crystal_t
56  use m_crystal_io,      only : crystal_ncwrite
57  use m_bz_mesh,         only : findqg0
58  use m_wfd,             only : wfd_init, wfd_free, wfd_print, wfd_t, wfd_test_ortho, wfd_copy_cg,&
59                                wfd_read_wfk, wfd_wave_free, wfd_rotate, wfd_reset_ur_cprj, wfd_get_ur
60  use m_pawang,          only : pawang_type
61  use m_pawrad,          only : pawrad_type
62  use m_pawtab,          only : pawtab_type
63  use m_pawfgr,          only : pawfgr_type
64  use m_eig2d,           only : gkk_t, gkk_init, gkk_ncwrite,gkk_free
65  use m_getgh1c,         only : getgh1c, rf_transgrid_and_pack, getgh1c_setup
66  use m_fourier_interpol, only : transgrid
67 ! use m_paw_an,          only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
68 ! use m_paw_ij,          only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
69 ! use m_pawfgrtab,       only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
70 ! use m_pawrhoij,        only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, symrhoij
71 ! use m_pawdij,          only : pawdij, symdij
72 ! use m_pawcprj,         only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
73 
74  implicit none
75 
76  private

m_phpi/eph_phpi [ Functions ]

[ Top ] [ m_phpi ] [ Functions ]

NAME

  eph_phpi

FUNCTION

  Compute phonon 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

PARENTS

      eph

NOTES

CHILDREN

SOURCE

117 subroutine eph_phpi(wfk0_path,wfq_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands_k,ebands_kq,dvdb,ifc,&
118                        pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm)
119 
120 
121 !This section has been created automatically by the script Abilint (TD).
122 !Do not modify the following lines by hand.
123 #undef ABI_FUNC
124 #define ABI_FUNC 'eph_phpi'
125 !End of the abilint section
126 
127  implicit none
128 
129 !Arguments ------------------------------------
130 !scalars
131  character(len=*),intent(in) :: wfk0_path, wfq_path
132  integer,intent(in) :: comm
133  type(datafiles_type),intent(in) :: dtfil
134  type(dataset_type),intent(in) :: dtset
135  type(crystal_t),intent(in) :: cryst
136  type(ebands_t),intent(in) :: ebands_k, ebands_kq
137  type(dvdb_t),intent(inout) :: dvdb
138  type(pawang_type),intent(in) :: pawang
139  type(pseudopotential_type),intent(in) :: psps
140  type(pawfgr_type),intent(in) :: pawfgr
141  type(ifc_type),intent(in) :: ifc
142  type(mpi_type),intent(in) :: mpi_enreg
143 !arrays
144  integer,intent(in) :: ngfft(18),ngfftf(18)
145  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
146  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
147 
148 !Local variables ------------------------------
149 !scalars
150  integer,parameter :: dummy_npw=1,nsig=1,tim_getgh1c=1,berryopt0=0,timrev1=1
151  integer,parameter :: useylmgr=0,useylmgr1=0,master=0
152  integer :: my_rank,nproc,iomode,mband,mband_kq,my_minb,my_maxb,nsppol,nkpt,nkpt_kq,idir,ipert
153  integer :: cplex,db_iqpt,natom,natom3,ipc,nspinor,onpw,imode
154  integer :: ib1,ib2
155  integer :: ik,ikq
156  integer :: spin,istwf_k,istwf_kq,npw_k,npw_kq
157  integer :: mpw,my_mpw,ierr,my_kstart,my_kstop,cnt
158  integer :: n1,n2,n3,n4,n5,n6,nspden
159  integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnl1
160  integer :: nfft,nfftf,mgfft,mgfftf,nkpg,nkpg1
161  real(dp) :: cpu,wall,gflops
162  real(dp) :: ecut,eshift,eig0nk,eig0mkq,dotr,doti
163  real(dp) :: eta,f_nk,f_mkq,omega,wtk,gkk2,term1,term2
164  logical :: i_am_master,gen_eigenpb
165  type(wfd_t) :: wfd_k,wfd_kq
166  type(gs_hamiltonian_type) :: gs_hamkq
167  type(rf_hamiltonian_type) :: rf_hamkq
168  character(len=500) :: msg
169 !arrays
170  integer :: g0_k(3),dummy_gvec(3,dummy_npw)
171  integer,allocatable :: kg_k(:,:),kg_kq(:,:),gtmp(:,:),nband(:,:),nband_kq(:,:),blkflg(:,:)
172  real(dp) :: kk(3),kq(3),qpt(3),phfrq(3*cryst%natom)
173  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom),displ_red(2,3,cryst%natom,3*cryst%natom)
174  real(dp) :: Pi_ph(3*cryst%natom)
175  real(dp),allocatable :: grad_berry(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),dkinpw(:)
176  real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:),ph3d(:,:,:),ph3d1(:,:,:)
177  real(dp),allocatable :: v1scf(:,:,:,:),gkk(:,:,:,:,:), gkk_m(:,:,:)
178  real(dp),allocatable :: bras_kq(:,:,:),kets_k(:,:,:),h1kets_kq(:,:,:)
179  real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:)
180  real(dp),allocatable :: ylm_kq(:,:),ylm_k(:,:),ylmgr_kq(:,:,:)
181  real(dp),allocatable :: dummy_vtrial(:,:),gvnl1(:,:)
182  real(dp),allocatable ::  gs1c(:,:)
183  logical,allocatable :: bks_mask(:,:,:),bks_mask_kq(:,:,:),keep_ur(:,:,:),keep_ur_kq(:,:,:)
184  type(pawcprj_type),allocatable  :: cwaveprj0(:,:) !natom,nspinor*usecprj)
185 
186 !************************************************************************
187 
188  write(msg, '(3a)') ch10, "Computation of the real part of the phonon self-energy", ch10
189  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
190  call wrtout(std_out, msg, "COLL", do_flush=.True.)
191 
192  if (psps%usepaw == 1) then
193    MSG_ERROR("PAW not implemented")
194    ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
195  end if
196 
197  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
198  i_am_master = (my_rank == master)
199 
200  ! Copy important dimensions
201  natom = cryst%natom
202  natom3 = 3 * natom
203  nsppol = ebands_k%nsppol
204  nspinor = ebands_k%nspinor
205  nspden = dtset%nspden
206  nkpt = ebands_k%nkpt
207  mband = ebands_k%mband
208  nkpt_kq = ebands_kq%nkpt
209  mband_kq = ebands_kq%mband
210  ecut = dtset%ecut
211 
212 ! GKA TODO: Make sure there is a single q-point present.
213  qpt = dtset%qptn(:)
214 
215  nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
216  nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
217  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
218  n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
219 
220 
221  ! Open the DVDB file
222  call dvdb_open_read(dvdb, ngfftf, xmpi_comm_self)
223 
224 
225  ! Initialize the wave function descriptors.
226  ! For the time being, no memory distribution, each node has the full set of states.
227  my_minb = 1; my_maxb = mband
228 
229  ABI_MALLOC(nband, (nkpt, nsppol))
230  ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
231  ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol))
232  nband=mband; bks_mask=.False.; keep_ur=.False.
233 
234  ABI_MALLOC(nband_kq, (nkpt_kq, nsppol))
235  ABI_MALLOC(bks_mask_kq,(mband_kq, nkpt_kq, nsppol))
236  ABI_MALLOC(keep_ur_kq,(mband_kq, nkpt_kq ,nsppol))
237  nband_kq=mband_kq; bks_mask_kq=.False.; keep_ur_kq=.False.
238 
239  ! Distribute the k-points over the processors
240  call xmpi_split_work(nkpt,comm,my_kstart,my_kstop,msg,ierr)
241  do ik=1,nkpt
242  if (.not. ((ik .ge. my_kstart) .and. (ik .le. my_kstop))) cycle
243 
244    kk = ebands_k%kptns(:,ik)
245    kq = kk + qpt
246    call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:),(/1,1,1/))  ! Find the index of the k+q point
247 
248    bks_mask(:,ik,:) = .True.
249    bks_mask_kq(:,ikq,:) = .True.
250  end do
251 
252  ! Initialize the wavefunction descriptors
253  call wfd_init(wfd_k,cryst,pawtab,psps,keep_ur,dtset%paral_kgb,dummy_npw,mband,nband,nkpt,nsppol,bks_mask,&
254    nspden,nspinor,dtset%ecutsm,dtset%dilatmx,ebands_k%istwfk,ebands_k%kptns,ngfft,&
255    dummy_gvec,dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm,opt_ecut=ecut)
256 
257  call wfd_print(wfd_k,header="Wavefunctions on the k-points grid",mode_paral='PERS')
258 
259  call wfd_init(wfd_kq,cryst,pawtab,psps,keep_ur_kq,dtset%paral_kgb,dummy_npw,mband_kq,nband_kq,nkpt_kq,nsppol,bks_mask_kq,&
260    nspden,nspinor,dtset%ecutsm,dtset%dilatmx,ebands_kq%istwfk,ebands_kq%kptns,ngfft,&
261    dummy_gvec,dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm,opt_ecut=ecut)
262 
263  call wfd_print(wfd_kq,header="Wavefunctions on the q-shifted k-points grid",mode_paral='PERS')
264 
265  ABI_FREE(nband)
266  ABI_FREE(bks_mask)
267  ABI_FREE(keep_ur)
268  ABI_FREE(nband_kq)
269  ABI_FREE(bks_mask_kq)
270  ABI_FREE(keep_ur_kq)
271 
272  ! Read wafefunctions on the k-points grid and q-shifted k-points grid.
273  iomode = iomode_from_fname(wfk0_path)
274  call wfd_read_wfk(wfd_k,wfk0_path,iomode)
275  if (.False.) call wfd_test_ortho(wfd_k,cryst,pawtab,unit=std_out,mode_paral="PERS")
276 
277  iomode = iomode_from_fname(wfq_path)
278  call wfd_read_wfk(wfd_kq,wfq_path,iomode)
279  if (.False.) call wfd_test_ortho(wfd_kq,cryst,pawtab,unit=std_out,mode_paral="PERS")
280 
281  ! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information on the coarse grid.
282  ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
283  call getph(cryst%atindx,natom,n1,n2,n3,ph1d,cryst%xred)
284 
285  ! Find the appropriate value of mpw
286  mpw = 0; cnt=0
287  do spin=1,nsppol
288    do ik=1,nkpt
289      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
290      kk = ebands_k%kptns(:,ik)
291      call get_kg(kk,1,ecut,cryst%gmet,onpw,gtmp)
292      ABI_FREE(gtmp)
293      mpw = max(mpw, onpw)
294    end do
295  end do
296  cnt=0
297  do spin=1,nsppol
298    do ikq=1,nkpt_kq
299      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
300      kq = ebands_kq%kptns(:,ikq)
301      call get_kg(kq,1,ecut,cryst%gmet,onpw,gtmp)
302      ABI_FREE(gtmp)
303      mpw = max(mpw, onpw)
304    end do
305  end do
306  my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
307 
308 
309  ! Allow PW-arrays dimensioned with mpw
310  ABI_MALLOC(kg_k, (3, mpw))
311  ABI_MALLOC(kg_kq, (3, mpw))
312 
313  ! Spherical Harmonics for useylm==1.
314  ABI_MALLOC(ylm_k,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
315  ABI_MALLOC(ylm_kq,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
316  ABI_MALLOC(ylmgr_kq,(mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
317 
318  ! TODO FOR PAW
319  usecprj = 0
320  ABI_DT_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
321 
322  ! Prepare call to getgh1c
323  usevnl = 0
324  optlocal = 1  ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
325  optnl = 2     ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
326  opt_gvnl1 = 0 ! gvnl1 is output
327  ABI_MALLOC(gvnl1, (2,usevnl))
328  ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4)))
329 
330  ! This part is taken from dfpt_vtorho
331  !==== Initialize most of the Hamiltonian (and derivative) ====
332  !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
333  !2) Perform the setup needed for the non-local factors:
334  !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
335  !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
336 
337  call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,NSPPOL,nspden,natom,&
338 &  dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,&
339 &  comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
340 &  usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
341 
342  ! Allocate vlocal. Note nvloc
343  ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data)
344  ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
345  vlocal = huge(one)
346 
347  ! Allocate work space arrays.
348  ABI_MALLOC(blkflg, (natom3,natom3))
349  ABI_CALLOC(dummy_vtrial, (nfftf,nspden))
350 
351 
352  call cwtime(cpu,wall,gflops,"start")
353 
354  ! Find the index of the q-point in the DVDB.
355  db_iqpt = dvdb_findq(dvdb, qpt)
356 
357  if (db_iqpt /= -1) then
358    if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found: ",ktoa(qpt)," in DVDB with index ",itoa(db_iqpt)))
359    ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
360    ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
361    call dvdb_readsym_allv1(dvdb, db_iqpt, cplex, nfftf, ngfftf, v1scf, comm)
362  else
363    if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Could not find: ",ktoa(qpt), "in DVDB - interpolating"))
364    ! Fourier interpolate of the potential
365    ABI_CHECK(any(abs(qpt) > tol12), "qpt cannot be zero if Fourier interpolation is used")
366    cplex = 2
367    ABI_MALLOC(v1scf, (cplex,nfftf,nspden,natom3))
368    call dvdb_ftinterp_qpt(dvdb, qpt, nfftf, ngfftf, v1scf, comm)
369  end if
370 
371  ! Allocate vlocal1 with correct cplex. Note nvloc
372  ABI_STAT_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr)
373  ABI_CHECK(ierr==0, "oom vlocal1")
374 
375  ! Allocate el-ph coupling matrix elements
376  ABI_MALLOC(gkk, (2, mband_kq, mband, natom, 3))
377  ABI_MALLOC(gkk_m, (2, mband_kq, mband))
378 
379  ! Compute displacement vectors and phonon frequencies
380  call ifc_fourq(ifc,cryst,qpt,phfrq,displ_cart, out_displ_red=displ_red)
381 
382  ! Broadening parameter
383  if (dtset%elph2_imagden .gt. tol12) then
384    eta = dtset%elph2_imagden
385  else
386    eta = 0.0001_dp
387  end if
388 
389  ! Kpoints weights (not using symmetries at the moment)
390  wtk = 1.0 / nkpt
391 
392  ! Initialize phonon self-energy
393  Pi_ph = zero
394 
395 
396  ! Examine the symmetries of the q wavevector
397  ! call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol)
398 
399  ! ----------------------------------------------------------------------------------------------- !
400  ! Begin loop over states
401  ! ----------------------------------------------------------------------------------------------- !
402  do spin=1,nsppol
403 
404    ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
405    do ipc=1,natom3
406      call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,&
407              pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc))
408    end do
409 
410    ! Continue to initialize the Hamiltonian
411    call load_spin_hamiltonian(gs_hamkq,spin,vlocal=vlocal,with_nonlocal=.true.)
412 
413    do ik=1,nkpt
414 
415      ! Only do a subset a k-points
416      if (.not. ((ik .ge. my_kstart) .and. (ik .le. my_kstop))) cycle
417 
418      ! Allocate workspace for wavefunctions. Make npw larger than expected.
419      ABI_MALLOC(bras_kq, (2, mpw*nspinor, mband))
420      ABI_MALLOC(kets_k, (2, mpw*nspinor, mband))
421      ABI_MALLOC(h1kets_kq, (2, mpw*nspinor, mband))
422 
423      kk = ebands_k%kptns(:,ik)
424 
425      kq = kk + qpt
426      call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:),(/1,1,1/))  ! Find the index of the k+q point
427 
428 
429      ! Copy u_k(G)
430      istwf_k = wfd_k%istwfk(ik); npw_k = wfd_k%npwarr(ik)
431      ABI_CHECK(mpw >= npw_k, "mpw < npw_k")
432      kg_k(:,1:npw_k) = wfd_k%kdata(ik)%kg_k
433      do ib2=1,mband
434        call wfd_copy_cg(wfd_k, ib2, ik, spin, kets_k(1,1,ib2))
435      end do
436 
437      ! Copy u_kq(G)
438      istwf_kq = wfd_kq%istwfk(ikq); npw_kq = wfd_kq%npwarr(ikq)
439      ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
440      kg_kq(:,1:npw_kq) = wfd_kq%kdata(ikq)%kg_k
441      do ib1=1,mband_kq
442        call wfd_copy_cg(wfd_kq, ib1, ikq, spin, bras_kq(1,1,ib1))
443      end do
444 
445      ! if PAW, one has to solve a generalized eigenproblem
446      ! BE careful here because I will need sij_opt==-1
447      gen_eigenpb = (psps%usepaw==1)
448      sij_opt = 0; if (gen_eigenpb) sij_opt = 1
449      ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2)))
450 
451      gkk = zero
452 
453      ! Loop over all 3*natom perturbations.
454      do ipc=1,natom3
455        idir = mod(ipc-1, 3) + 1
456        ipert = (ipc - idir) / 3 + 1
457 
458        !write(msg, '(a,2i4)')  "Treating ipert, idir = ", ipert, idir
459        !call wrtout(std_out, msg, "COLL", do_flush=.True.)
460 
461        ! Prepare application of the NL part.
462        call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.)
463            !&paw_ij1=paw_ij1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
464            !&mpi_spintab=mpi_enreg%my_isppoltab)
465        call load_spin_rf_hamiltonian(rf_hamkq,spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.)
466 
467        ! This call is not optimal because there are quantities in out that do not depend on idir,ipert
468        call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,&                   ! In
469          cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,&                           ! In
470          npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq,&                          ! In
471          dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1)                      ! Out
472 
473        ! Calculate dvscf * psi_k, results stored in h1kets_kq on the k+q sphere.
474        ! Compute H(1) applied to GS wavefunction Psi(0)
475        do ib2=1,mband
476          eig0nk = ebands_k%eig(ib2,ik,spin)
477          ! Use scissor shift on 0-order eigenvalue
478          eshift = eig0nk - dtset%dfpt_sciss
479 
480          call getgh1c(berryopt0,kets_k(:,:,ib2),cwaveprj0,h1kets_kq(:,:,ib2),&
481 &                     grad_berry,gs1c,gs_hamkq,gvnl1,idir,ipert,eshift,mpi_enreg,optlocal,&
482 &                     optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
483        end do
484 
485        ABI_FREE(kinpw1)
486        ABI_FREE(kpg1_k)
487        ABI_FREE(kpg_k)
488        ABI_FREE(dkinpw)
489        ABI_FREE(ffnlk)
490        ABI_FREE(ffnl1)
491        ABI_FREE(ph3d)
492 
493        if (allocated(gs1c)) then
494          ABI_FREE(gs1c)
495        end if
496 
497        if (allocated(ph3d1)) then
498          ABI_FREE(ph3d1)
499        end if
500 
501        ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation.
502        !The array eig1_k contains:
503        !
504        ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)>                           (NC psps)
505        ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW)
506        do ib2=1,mband
507          do ib1=1,mband_kq
508            call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bras_kq(1,1,ib1),h1kets_kq(1,1,ib2),&
509              mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
510 
511            gkk(1,ib1,ib2,ipert,idir) = dotr
512            gkk(2,ib1,ib2,ipert,idir) =  doti
513          end do
514        end do
515 
516      end do  ! ipc
517 
518 
519      ! Loop over 3*natom phonon branches.
520      do imode=1,natom3
521 
522        omega = phfrq(imode)
523 
524        ! Do not compute Pi for negative or too small frequencies
525        if (omega .lt. tol6) cycle
526 
527        gkk_m = zero
528 
529        ! Transform the gkk from atom,cart basis to mode basis
530        do idir=1,3
531          do ipert=1,natom
532             gkk_m(1,:,:) = gkk_m(1,:,:) &
533 &                          + gkk(1,:,:,ipert,idir) * displ_red(1,idir,ipert,imode) &
534 &                          - gkk(2,:,:,ipert,idir) * displ_red(2,idir,ipert,imode)
535             gkk_m(2,:,:) = gkk_m(2,:,:) &
536 &                          + gkk(1,:,:,ipert,idir) * displ_red(2,idir,ipert,imode) &
537 &                          + gkk(2,:,:,ipert,idir) * displ_red(1,idir,ipert,imode)
538          end do
539        end do
540 
541        gkk_m = gkk_m / sqrt(two * omega)
542 
543        ! sum contribution to phonon self-energy
544        do ib2=1,mband
545          do ib1=1,mband_kq
546 
547            f_nk = ebands_k%occ(ib2,ik,spin)
548            f_mkq = ebands_kq%occ(ib1,ikq,spin)
549 
550            if (abs(f_mkq - f_nk) .le. tol12) cycle
551 
552            eig0nk = ebands_k%eig(ib2,ik,spin)
553            eig0mkq = ebands_kq%eig(ib1,ikq,spin)
554 
555            gkk2 = gkk_m(1,ib1,ib2) ** 2 + gkk_m(2,ib1,ib2) ** 2
556 
557            term1 = (f_mkq - f_nk) * (eig0mkq - eig0nk - omega) / ((eig0mkq - eig0nk - omega) ** 2 + eta ** 2)
558            term2 = (f_mkq - f_nk) * (eig0mkq - eig0nk        ) / ((eig0mkq - eig0nk        ) ** 2 + eta ** 2)
559 
560            Pi_ph(imode) = Pi_ph(imode) + wtk * gkk2 * (term1 - term2)
561 
562          end do
563        end do
564 
565      end do  ! imode
566 
567      ABI_FREE(bras_kq)
568      ABI_FREE(kets_k)
569      ABI_FREE(h1kets_kq)
570 
571    end do ! ikfs
572 
573    call destroy_rf_hamiltonian(rf_hamkq)
574 
575  end do ! spin
576 
577  ! Gather the k-points computed by all processes
578  call xmpi_sum_master(Pi_ph,master,comm,ierr)
579 
580 
581  ! Output the results
582  if (i_am_master) then
583    call out_phpi(ab_out, Pi_ph, phfrq, qpt, natom3)
584    call out_phpi(std_out, Pi_ph, phfrq, qpt, natom3)
585  end if
586 
587 #ifdef HAVE_NETCDF
588  if (i_am_master) then
589    call out_phpi_nc(dtfil, cryst, Pi_ph, phfrq, qpt, natom3)
590  end if
591 #endif
592 
593  ! Free memory
594  call cwtime(cpu,wall,gflops,"stop")
595 
596  write(msg, '(3a)') "Computation of the real part of the phonon self-energy completed", ch10, &
597 &                   "--------------------------------------------------------------------------------"
598  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
599  call wrtout(std_out, msg, "COLL", do_flush=.True.)
600 
601  ! Free memory
602 
603  ABI_FREE(gkk)
604  ABI_FREE(gkk_m)
605  ABI_FREE(v1scf)
606  ABI_FREE(vlocal1)
607  ABI_FREE(gvnl1)
608  ABI_FREE(grad_berry)
609  ABI_FREE(dummy_vtrial)
610  ABI_FREE(ph1d)
611  ABI_FREE(vlocal)
612  ABI_FREE(kg_k)
613  ABI_FREE(kg_kq)
614  ABI_FREE(ylm_k)
615  ABI_FREE(ylm_kq)
616  ABI_FREE(ylmgr_kq)
617  ABI_FREE(blkflg)
618 
619  call destroy_hamiltonian(gs_hamkq)
620  call wfd_free(wfd_k)
621  call wfd_free(wfd_kq)
622 
623  call pawcprj_free(cwaveprj0)
624  ABI_DT_FREE(cwaveprj0)
625 
626 end subroutine eph_phpi

m_phpi/out_phpi [ Functions ]

[ Top ] [ m_phpi ] [ Functions ]

NAME

  out_phpi

FUNCTION

  Output the phonon self-energy.

INPUTS

OUTPUT

PARENTS

      m_phpi

NOTES

CHILDREN

SOURCE

651 subroutine out_phpi(iout, Pi_ph, phfrq, qpt, natom3)
652 
653 !Arguments ------------------------------------
654 !scalars
655 
656 !This section has been created automatically by the script Abilint (TD).
657 !Do not modify the following lines by hand.
658 #undef ABI_FUNC
659 #define ABI_FUNC 'out_phpi'
660 !End of the abilint section
661 
662  integer,intent(in) :: iout
663  integer,intent(in) :: natom3
664 !arrays
665  real(dp),intent(in) :: Pi_ph(natom3),phfrq(natom3),qpt(3)
666 
667 !Local variables ------------------------------
668 !scalars
669  integer :: imode
670 
671  write(iout,'(a)')' '
672  !write(iout,'(a)')' ----------------------------------------'
673  !write(iout,'(a)')' '
674  write(iout,'(a)')' Phonon self-energy (Hartree)'
675  write(iout,'(a)')' '
676  write(iout,'(a,3f14.8)')' qpt =',qpt
677  write(iout,'(a)')' '
678  write(iout,'(1x,a,10x,a)')'omega','Pi(omega)'
679 
680  do imode=1,natom3
681     write(iout,'(1x,f12.8,1x,es14.6)') phfrq(imode), Pi_ph(imode)
682  end do
683 
684  write(iout,'(a)')' '
685  !write(iout,'(a)')' ----------------------------------------'
686 
687 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

PARENTS

      m_phpi

NOTES

CHILDREN

SOURCE

712 subroutine out_phpi_nc(dtfil, cryst, Pi_ph, phfrq, qpt, natom3)
713 
714 #ifdef HAVE_NETCDF
715  use netcdf
716 #endif
717 
718 !Arguments ------------------------------------
719 !scalars
720 
721 !This section has been created automatically by the script Abilint (TD).
722 !Do not modify the following lines by hand.
723 #undef ABI_FUNC
724 #define ABI_FUNC 'out_phpi_nc'
725 !End of the abilint section
726 
727  integer,intent(in) :: natom3
728  type(datafiles_type), intent(in) :: dtfil
729  type(crystal_t),intent(in) :: cryst
730 !arrays
731  real(dp),intent(in) :: Pi_ph(natom3),phfrq(natom3),qpt(3)
732 
733 !Local variables ------------------------------
734 !scalars
735  integer :: natom,one_dim,cplex,cart_dir
736  integer :: ncid, ncerr
737  character(len=fnlen) :: fname
738 
739 #ifdef HAVE_NETCDF
740 
741  ! Initialize NetCDF file.
742  fname = strcat(dtfil%filnam_ds(4),"_Pi.nc")
743  NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
744 
745  ! Write information of the crystal
746  ncerr = crystal_ncwrite(cryst, ncid)
747  NCF_CHECK(ncerr)
748 
749  ! Write the dimensions specified by ETSF
750  one_dim = 1
751  cplex = 2
752  cart_dir = 3
753 
754  natom = natom3 / 3
755 
756  ncerr = nctk_def_dims(ncid, [&
757 & nctkdim_t('current_one_dim', one_dim), &
758 & nctkdim_t('number_of_atoms', natom), &
759 & nctkdim_t('number_of_cartesian_directions', cart_dir), &
760 & nctkdim_t('number_of_perturbations', natom3), &
761 & nctkdim_t('cplex',cplex)], defmode=.True.)
762  NCF_CHECK(ncerr)
763 
764  ! Create the arrays
765  ncerr = nctk_def_arrays(ncid, [&
766  nctkarr_t('q_point_reduced_coord', "dp", 'number_of_cartesian_directions'),&
767  nctkarr_t('phonon_frequencies', "dp", 'number_of_perturbations'), &
768  nctkarr_t('phonon_self_energy_realpart', "dp", 'number_of_perturbations')])
769  NCF_CHECK(ncerr)
770 
771  NCF_CHECK(nctk_set_atomic_units(ncid, 'phonon_frequencies'))
772  NCF_CHECK(nctk_set_atomic_units(ncid, 'phonon_self_energy_realpart'))
773 
774 ! Write data
775  NCF_CHECK(nctk_set_datamode(ncid))
776  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'q_point_reduced_coord'), qpt))
777  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phonon_frequencies'), phfrq))
778  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phonon_self_energy_realpart'), Pi_ph))
779 
780  ! Close file
781  NCF_CHECK(nf90_close(ncid))
782 
783 
784 #else
785  MSG_ERROR("NETCDF support required to write Pi.nc file.")
786 #endif
787 
788 end subroutine out_phpi_nc