TABLE OF CONTENTS


ABINIT/m_gkk [ Modules ]

[ Top ] [ Modules ]

NAME

FUNCTION

  Tools for the computation of electron-phonon coupling matrix elements (gkk)

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (GKA, MG)
  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 module m_gkk
25 
26  use defs_basis
27  use defs_abitypes
28  use m_abicore
29  use m_xmpi
30  use m_errors
31  use m_ifc
32  use m_ebands
33  use m_ddb
34  use m_dvdb
35  use m_fft
36  use m_hamiltonian
37  use m_pawcprj
38  use m_wfk
39  use m_nctk
40 #ifdef HAVE_NETCDF
41  use netcdf
42 #endif
43 
44  use m_time,           only : cwtime, sec2str
45  use m_io_tools,       only : iomode_from_fname
46  use m_fstrings,       only : itoa, sjoin, ktoa, ltoa, strcat
47  use m_symtk,          only : littlegroup_q
48  use m_fftcore,        only : ngfft_seq, get_kg, kpgsph, sphere
49  use defs_datatypes,   only : ebands_t, pseudopotential_type
50  use m_crystal,        only : crystal_t
51  use m_crystal_io,     only : crystal_ncwrite
52  use m_bz_mesh,        only : findqg0
53  use m_cgtools,        only : dotprod_g
54  use m_kg,             only : getph
55  use m_pawang,         only : pawang_type
56  use m_pawrad,         only : pawrad_type
57  use m_pawtab,         only : pawtab_type
58  use m_pawfgr,         only : pawfgr_type
59  use m_eig2d,          only : gkk_t, gkk_init, gkk_ncwrite, gkk_free
60  use m_wfd,            only : wfd_init, wfd_free, wfd_print, wfd_t, wfd_test_ortho, wfd_copy_cg,&
61                               wfd_read_wfk, wfd_wave_free, wfd_rotate, wfd_reset_ur_cprj, wfd_get_ur
62  use m_getgh1c,          only : getgh1c, rf_transgrid_and_pack, getgh1c_setup
63  use m_fourier_interpol, only : transgrid
64 ! use m_paw_an,          only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
65 ! use m_paw_ij,          only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
66 ! use m_pawfgrtab,       only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
67 ! use m_pawrhoij,        only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, symrhoij
68 ! use m_pawdij,          only : pawdij, symdij
69 ! use m_pawcprj,         only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
70 
71  implicit none
72 
73  private

m_gkk/eph_gkk [ Functions ]

[ Top ] [ m_gkk ] [ Functions ]

NAME

  eph_gkk

FUNCTION

  Compute electron-phonon coupling matrix elements.

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

      get_kg

SOURCE

115 subroutine eph_gkk(wfk0_path,wfq_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands_k,ebands_kq,dvdb,ifc,&
116                        pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm)
117 
118 
119 !This section has been created automatically by the script Abilint (TD).
120 !Do not modify the following lines by hand.
121 #undef ABI_FUNC
122 #define ABI_FUNC 'eph_gkk'
123 !End of the abilint section
124 
125  implicit none
126 
127 !Arguments ------------------------------------
128 !scalars
129  character(len=*),intent(in) :: wfk0_path, wfq_path
130  integer,intent(in) :: comm
131  type(datafiles_type),intent(in) :: dtfil
132  type(dataset_type),intent(in) :: dtset
133  type(crystal_t),intent(in) :: cryst
134  type(ebands_t),intent(in) :: ebands_k, ebands_kq
135  type(dvdb_t),target,intent(inout) :: dvdb
136  type(pawang_type),intent(in) :: pawang
137  type(pseudopotential_type),intent(in) :: psps
138  type(pawfgr_type),intent(in) :: pawfgr
139  type(ifc_type),intent(in) :: ifc
140  type(mpi_type),intent(inout) :: mpi_enreg
141 !arrays
142  integer,intent(in) :: ngfft(18),ngfftf(18)
143  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
144  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
145 
146 !Local variables ------------------------------
147 !scalars
148  integer,parameter :: dummy_npw=1,tim_getgh1c=1,berryopt0=0
149  integer,parameter :: useylmgr1=0,master=0
150  integer :: my_rank,nproc,iomode,mband,mband_kq,my_minb,my_maxb,nsppol,nkpt,nkpt_kq,idir,ipert
151  integer :: cplex,db_iqpt,natom,natom3,ipc,nspinor
152  integer :: ib1,ib2,band
153  integer :: ik,ikq,timerev_q
154  integer :: spin,istwf_k,istwf_kq,npw_k,npw_kq
155  integer :: mpw,mpw_k,mpw_kq,ierr,my_kstart,my_kstop,ncid
156  integer :: n1,n2,n3,n4,n5,n6,nspden
157  integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnl1
158  integer :: nfft,nfftf,mgfft,mgfftf,nkpg,nkpg1
159  real(dp) :: cpu,wall,gflops
160  real(dp) :: ecut,eshift,eig0nk,dotr,doti
161  logical :: i_am_master, gen_eigenpb
162  type(wfd_t) :: wfd_k, wfd_kq
163  type(gs_hamiltonian_type) :: gs_hamkq
164  type(rf_hamiltonian_type) :: rf_hamkq
165  character(len=500) :: msg
166  character(len=fnlen) :: fname, gkkfilnam
167 !arrays
168  integer :: g0_k(3),symq(4,2,cryst%nsym),dummy_gvec(3,dummy_npw)
169  integer,allocatable :: kg_k(:,:),kg_kq(:,:),nband(:,:),nband_kq(:,:),blkflg(:,:)
170  real(dp) :: kk(3),kq(3),qpt(3)
171  real(dp),allocatable :: grad_berry(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),dkinpw(:)
172  real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:),ph3d(:,:,:),ph3d1(:,:,:)
173  real(dp),allocatable :: v1scf(:,:,:,:),gkk(:,:,:,:,:)
174  real(dp),allocatable :: bras(:,:,:),kets(:,:,:),h1_kets(:,:,:)
175  real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:)
176  real(dp),allocatable :: ylm_kq(:,:),ylm_k(:,:),ylmgr_kq(:,:,:)
177  real(dp),allocatable :: dummy_vtrial(:,:),gvnl1(:,:)
178  real(dp),allocatable ::  gs1c(:,:)
179  logical,allocatable :: bks_mask(:,:,:),bks_mask_kq(:,:,:),keep_ur(:,:,:),keep_ur_kq(:,:,:)
180  type(pawcprj_type),allocatable  :: cwaveprj0(:,:) !natom,nspinor*usecprj)
181  type(gkk_t)     :: gkk2d
182 
183 !************************************************************************
184 
185  write(msg, '(2a)') "Computation of electron-phonon coupling matrix elements (gkk)", ch10
186  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
187  call wrtout(std_out, msg, "COLL", do_flush=.True.)
188 
189  if (psps%usepaw == 1) then
190    MSG_ERROR("PAW not implemented")
191    ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
192  end if
193 
194  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
195  i_am_master = (my_rank == master)
196 
197  ! Copy important dimensions
198  cplex = 2
199  natom = cryst%natom
200  natom3 = 3 * natom
201  nsppol = ebands_k%nsppol
202  nspinor = ebands_k%nspinor
203  nspden = dtset%nspden
204  nkpt = ebands_k%nkpt
205  mband = ebands_k%mband
206  nkpt_kq = ebands_kq%nkpt
207  mband_kq = ebands_kq%mband
208  ecut = dtset%ecut
209 
210 ! GKA TODO: Make sure there is a single q-point present.
211  qpt = dtset%qptn(:)
212 
213  nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
214  nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
215  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
216  n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
217 
218 
219  ! Open the DVDB file
220  call dvdb_open_read(dvdb, ngfftf, xmpi_comm_self)
221 
222 
223  ! Initialize the wave function descriptors.
224  ! For the time being, no memory distribution, each node has the full set of states.
225  my_minb = 1; my_maxb = mband
226 
227  ABI_MALLOC(nband, (nkpt, nsppol))
228  ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
229  ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol))
230  nband=mband; bks_mask=.False.; keep_ur=.False.
231 
232  ABI_MALLOC(nband_kq, (nkpt_kq, nsppol))
233  ABI_MALLOC(bks_mask_kq,(mband_kq, nkpt_kq, nsppol))
234  ABI_MALLOC(keep_ur_kq,(mband_kq, nkpt_kq ,nsppol))
235  nband_kq=mband_kq; bks_mask_kq=.False.; keep_ur_kq=.False.
236 
237 
238  ! Distribute the k-points over the processors
239  call xmpi_split_work(nkpt,comm,my_kstart,my_kstop,msg,ierr)
240  do ik=1,nkpt
241  if (.not. ((ik .ge. my_kstart) .and. (ik .le. my_kstop))) cycle
242 
243    kk = ebands_k%kptns(:,ik)
244    kq = kk + qpt
245    call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:),(/1,1,1/))  ! Find the index of the k+q point
246 
247    bks_mask(:,ik,:) = .True.
248    bks_mask_kq(:,ikq,:) = .True.
249  end do
250 
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  call find_mpw(mpw_k, ebands_k%kptns(:,:), nsppol, nkpt, cryst%gmet,ecut,comm)
287  call find_mpw(mpw_kq, ebands_kq%kptns(:,:), nsppol, nkpt_kq, cryst%gmet,ecut,comm)
288  mpw = max(mpw_k, mpw_kq)
289 
290  ! Allow PW-arrays dimensioned with mpw
291  ABI_MALLOC(kg_k, (3, mpw))
292  ABI_MALLOC(kg_kq, (3, mpw))
293 
294  ! Spherical Harmonics for useylm==1.
295  ABI_MALLOC(ylm_k,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
296  ABI_MALLOC(ylm_kq,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
297  ABI_MALLOC(ylmgr_kq,(mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
298 
299  ! TODO FOR PAW
300  usecprj = 0
301  ABI_DT_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
302 
303  ! Prepare call to getgh1c
304  usevnl = 0
305  optlocal = 1  ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
306  optnl = 2     ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
307  opt_gvnl1 = 0 ! gvnl1 is output
308  ABI_MALLOC(gvnl1, (2,usevnl))
309  ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4)))
310 
311  ! This part is taken from dfpt_vtorho
312  !==== Initialize most of the Hamiltonian (and derivative) ====
313  !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
314  !2) Perform the setup needed for the non-local factors:
315  !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
316  !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
317 
318  call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,natom,&
319 &  dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,&
320 &  usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda,&
321 &  comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
322 
323 ! Allocate vlocal. Note nvloc
324  ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
325 
326  ! Allocate work space arrays.
327  ABI_MALLOC(blkflg, (natom3,natom3))
328  ABI_CALLOC(dummy_vtrial, (nfftf,nspden))
329 
330 
331  call cwtime(cpu,wall,gflops,"start")
332 
333  ! Find the index of the q-point in the DVDB.
334  db_iqpt = dvdb_findq(dvdb, qpt)
335 
336  if (db_iqpt /= -1) then
337    if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found: ",ktoa(qpt)," in DVDB with index ",itoa(db_iqpt)))
338    ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
339    ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
340    call dvdb_readsym_allv1(dvdb, db_iqpt, cplex, nfftf, ngfftf, v1scf, comm)
341  else
342    if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Could not find: ",ktoa(qpt), "in DVDB - interpolating"))
343    ! Fourier interpolate of the potential
344    ABI_CHECK(any(abs(qpt) > tol12), "qpt cannot be zero if Fourier interpolation is used")
345 
346    ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
347    call dvdb_interpolate_v1scf(dvdb,cryst,qpt,ifc%ngqpt,ifc%nqshft,ifc%qshft, &
348    &                           nfft, ngfft, nfftf, ngfftf, v1scf, comm)
349 
350  end if
351 
352  ! Examine the symmetries of the q wavevector
353  call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol)
354 
355  ! Allocate vlocal1 with correct cplex. Note nvloc
356  ABI_STAT_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr)
357  ABI_CHECK(ierr==0, "oom vlocal1")
358 
359  ABI_MALLOC(gkk, (2*mband*nsppol,nkpt,1,1,mband_kq))
360 
361  ! ========================================================================== !
362  ! Begin loop on perturbations, spins and k-points
363  ! ========================================================================== !
364 
365  ! Loop over all 3*natom perturbations.
366  do ipc=1,natom3
367    idir = mod(ipc-1, 3) + 1
368    ipert = (ipc - idir) / 3 + 1
369 
370    write(msg, '(a,2i4)')  "Treating ipert, idir = ", ipert, idir
371    call wrtout(std_out, msg, "COLL", do_flush=.True.)
372 
373    gkk = zero
374 
375    do spin=1,nsppol
376 
377      ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
378      !do ipc=1,natom3
379      call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,&
380                pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc))
381      !end do
382 
383      ! Continue to initialize the Hamiltonian
384      call load_spin_hamiltonian(gs_hamkq,spin,vlocal=vlocal,with_nonlocal=.true.)
385 
386      ! Allocate workspace for wavefunctions. Make npw larger than expected.
387      ABI_MALLOC(bras, (2, mpw*nspinor, mband))
388      ABI_MALLOC(kets, (2, mpw*nspinor, mband))
389      ABI_MALLOC(h1_kets, (2, mpw*nspinor, mband))
390 
391 
392      ! GKA : This little block used to be right after the perturbation loop
393      ! Prepare application of the NL part.
394      call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.)
395                ! paw_ij1=paw_ij1,comm_atom=mpi_enreg%comm_atom,&
396                !&mpi_atmtab=mpi_enreg%my_atmtab,my_spintab=mpi_enreg%my_isppoltab)
397      call load_spin_rf_hamiltonian(rf_hamkq,spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.)
398 
399      do ik=1,nkpt
400 
401        ! Only do a subset a k-points
402        if (.not. ((ik .ge. my_kstart) .and. (ik .le. my_kstop))) cycle
403 
404        kk = ebands_k%kptns(:,ik)
405 
406        kq = kk + qpt
407        call findqg0(ikq,g0_k,kq,nkpt_kq,ebands_kq%kptns(:,:),(/1,1,1/))  ! Find the index of the k+q point
408 
409        ! Copy u_k(G)
410        istwf_k = wfd_k%istwfk(ik); npw_k = wfd_k%npwarr(ik)
411        ABI_CHECK(mpw >= npw_k, "mpw < npw_k")
412        kg_k(:,1:npw_k) = wfd_k%kdata(ik)%kg_k
413        do ib2=1,mband
414          call wfd_copy_cg(wfd_k, ib2, ik, spin, kets(1,1,ib2))
415        end do
416 
417        ! Copy u_kq(G)
418        istwf_kq = wfd_kq%istwfk(ikq); npw_kq = wfd_kq%npwarr(ikq)
419        ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
420        kg_kq(:,1:npw_kq) = wfd_kq%kdata(ikq)%kg_k
421        do ib1=1,mband_kq
422          call wfd_copy_cg(wfd_kq, ib1, ikq, spin, bras(1,1,ib1))
423        end do
424 
425        ! if PAW, one has to solve a generalized eigenproblem
426        ! BE careful here because I will need sij_opt==-1
427        gen_eigenpb = (psps%usepaw==1)
428        sij_opt = 0; if (gen_eigenpb) sij_opt = 1
429        ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2)))
430 
431        ! GKA : Previous loop on 3*natom perturbations used to start here
432 
433        ! This call is not optimal because there are quantities in out that do not depend on idir,ipert
434        call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,&                   ! In
435          cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,&                           ! In
436          npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq,&                          ! In
437          dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1)                      ! Out
438 
439        ! Calculate dvscf * psi_k, results stored in h1_kets on the k+q sphere.
440        ! Compute H(1) applied to GS wavefunction Psi(0)
441        do ib2=1,mband
442          eig0nk = ebands_k%eig(ib2,ik,spin)
443          ! Use scissor shift on 0-order eigenvalue
444          eshift = eig0nk - dtset%dfpt_sciss
445 
446          call getgh1c(berryopt0,kets(:,:,ib2),cwaveprj0,h1_kets(:,:,ib2),&
447 &                     grad_berry,gs1c,gs_hamkq,gvnl1,idir,ipert,eshift,mpi_enreg,optlocal,&
448 &                     optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
449        end do
450 
451        ABI_FREE(kinpw1)
452        ABI_FREE(kpg1_k)
453        ABI_FREE(kpg_k)
454        ABI_FREE(dkinpw)
455        ABI_FREE(ffnlk)
456        ABI_FREE(ffnl1)
457        ABI_FREE(ph3d)
458        ABI_FREE(gs1c)
459 
460        if (allocated(ph3d1)) then
461          ABI_FREE(ph3d1)
462        end if
463 
464        ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation.
465        !The array eig1_k contains:
466        !
467        ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)>                           (NC psps)
468        ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW)
469        do ib2=1,mband
470          do ib1=1,mband_kq
471            call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bras(1,1,ib1),h1_kets(1,1,ib2),&
472              mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
473            band = 2*ib2-1 + (spin-1) * 2 * mband
474            gkk(band,ik,1,1,ib1) = dotr
475            gkk(band+1,ik,1,1,ib1) = doti
476          end do
477        end do
478 
479      end do ! ikfs
480 
481      call destroy_rf_hamiltonian(rf_hamkq)
482 
483      ABI_FREE(bras)
484      ABI_FREE(kets)
485      ABI_FREE(h1_kets)
486 
487    end do ! spin
488 
489    ! Gather the k-points computed by all processes
490    call xmpi_sum_master(gkk,master,comm,ierr)
491 
492    ! Init a gkk_t object
493    call gkk_init(gkk,gkk2d,mband,nsppol,nkpt,1,1)
494 
495    ! Write the netCDF file.
496    call appdig(ipc,dtfil%fnameabo_gkk,gkkfilnam)
497    fname = strcat(gkkfilnam,".nc")
498 #ifdef HAVE_NETCDF
499    if (i_am_master) then
500      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file")
501      NCF_CHECK(crystal_ncwrite(cryst, ncid))
502      NCF_CHECK(ebands_ncwrite(ebands_k, ncid))
503      call gkk_ncwrite(gkk2d,qpt,1.0_dp, ncid)
504      NCF_CHECK(nf90_close(ncid))
505    end if
506 #endif
507 
508    ! Free memory
509    call gkk_free(gkk2d)
510 
511  end do ! ipc (loop over 3*natom atomic perturbations)
512 
513  call cwtime(cpu,wall,gflops,"stop")
514 
515  write(msg, '(2a)') "Computation of electron-phonon coupling matrix elements (gkk) completed", ch10
516  call wrtout(ab_out, msg, "COLL", do_flush=.True.)
517  call wrtout(std_out, msg, "COLL", do_flush=.True.)
518 
519  ! ========================================================================== !
520  ! Free memory
521  ! ========================================================================== !
522 
523  ABI_FREE(gkk)
524  ABI_FREE(v1scf)
525  ABI_FREE(vlocal1)
526  ABI_FREE(gvnl1)
527  ABI_FREE(grad_berry)
528  ABI_FREE(dummy_vtrial)
529  ABI_FREE(ph1d)
530  ABI_FREE(vlocal)
531  ABI_FREE(kg_k)
532  ABI_FREE(kg_kq)
533  ABI_FREE(ylm_k)
534  ABI_FREE(ylm_kq)
535  ABI_FREE(ylmgr_kq)
536  ABI_FREE(blkflg)
537 
538  call destroy_hamiltonian(gs_hamkq)
539  call wfd_free(wfd_k)
540  call wfd_free(wfd_kq)
541 
542  call pawcprj_free(cwaveprj0)
543  ABI_DT_FREE(cwaveprj0)
544 
545 end subroutine eph_gkk

m_gkk/find_mpw [ Functions ]

[ Top ] [ m_gkk ] [ Functions ]

NAME

  find_mpw

FUNCTION

  Look at all k-points and spins to find the maximum
  number of plane waves.

INPUTS

OUTPUT

PARENTS

      m_gkk

NOTES

CHILDREN

      get_kg

SOURCE

746 subroutine find_mpw(mpw, kpts, nsppol, nkpt, gmet, ecut, comm)
747 
748  use defs_basis
749  use m_abicore
750  use m_xmpi
751  use m_errors
752  use m_fftcore,         only : get_kg
753 
754 !This section has been created automatically by the script Abilint (TD).
755 !Do not modify the following lines by hand.
756 #undef ABI_FUNC
757 #define ABI_FUNC 'find_mpw'
758 !End of the abilint section
759 
760  implicit none
761 
762 !Arguments ------------------------------------
763 !scalars
764  integer,intent(out) :: mpw
765  integer,intent(in) :: nsppol, nkpt
766  integer,intent(in) :: comm
767  real(dp),intent(in) :: ecut
768 !arrays
769  real(dp),intent(in) :: kpts(3,nkpt)
770  real(dp),intent(in) :: gmet(3,3)
771 
772 !Local variables ------------------------------
773 !scalars
774  integer :: my_rank, cnt, nproc, ierr
775  integer :: ispin, ikpt
776  integer :: my_mpw, onpw
777  integer,allocatable :: gtmp(:,:)
778  real(dp) :: kpt(3)
779 
780  my_rank = xmpi_comm_rank(comm)
781  nproc = xmpi_comm_size(comm)
782 
783  mpw = 0; cnt=0
784  do ispin=1,nsppol
785    do ikpt=1,nkpt
786      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
787      kpt = kpts(:,ikpt)
788      call get_kg(kpt,1,ecut,gmet,onpw,gtmp)
789      ABI_FREE(gtmp)
790      mpw = max(mpw, onpw)
791    end do
792  end do
793  my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
794 
795 end subroutine find_mpw

m_gkk/ncwrite_v1qnu [ Functions ]

[ Top ] [ m_gkk ] [ Functions ]

NAME

  ncwrite_v1qnu

FUNCTION

  Compute \delta V_{q,nu)(r) and dump results to netcdf file.
  This routine should be called by a single processor.

 INPUT
  dvdb<dbdb_type>=Database with the DFPT SCF potentials.
  cryst(crystal_t)=Crystalline structure
  ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
  nqlist=Number of q-points
  qlist(3,nqlist)=List of q-points where \delta V_{q,nu)(r) is wanted.
    Potentials will be Fourier interpolated if qpt is not in DVDB.
  prtvol=Verbosity level
  path=Name of the netcdf file.

OUTPUT

  Only writing

PARENTS

CHILDREN

SOURCE

577 subroutine ncwrite_v1qnu(dvdb, cryst, ifc, nqlist, qlist, prtvol, path)
578 
579 
580 !This section has been created automatically by the script Abilint (TD).
581 !Do not modify the following lines by hand.
582 #undef ABI_FUNC
583 #define ABI_FUNC 'ncwrite_v1qnu'
584 !End of the abilint section
585 
586  implicit none
587 
588 !Arguments ------------------------------------
589 !scalars
590  integer,intent(in) :: nqlist,prtvol
591  type(dvdb_t),intent(inout) :: dvdb
592  type(crystal_t),intent(in) :: cryst
593  type(ifc_type),intent(in) :: ifc
594  character(len=*),intent(in) :: path
595 !arrays
596  real(dp),intent(in) :: qlist(3,nqlist)
597 
598 !Local variables-------------------------------
599 !scalars
600  integer :: ip,nu,iq,db_iqpt,do_ftv1q,cplex,ispden,nfftf,comm
601 #ifdef HAVE_NETCDF
602  integer :: ncid,ncerr
603 #endif
604  character(len=500) :: msg
605 !arrays
606  integer :: ngfftf(18)
607  real(dp) :: phfrq(3*cryst%natom),qpt(3)
608  real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom),displ_red(2,3*cryst%natom,3*cryst%natom)
609  real(dp),allocatable :: v1scf(:,:,:,:),v1qnu(:,:,:,:)
610 
611 !************************************************************************
612 
613  call wrtout(std_out, sjoin("Writing \delta V_{q,nu)(r) potentials to file:", path), do_flush=.True.)
614 
615  comm = xmpi_comm_self
616  call ngfft_seq(ngfftf, dvdb%ngfft3_v1(:,1)); nfftf = product(ngfftf(1:3))
617 
618  call dvdb_open_read(dvdb, ngfftf, xmpi_comm_self)
619 
620  ! Create netcdf file.
621 #ifdef HAVE_NETCDF
622  NCF_CHECK(nctk_open_create(ncid, path, comm))
623  NCF_CHECK(crystal_ncwrite(cryst, ncid))
624 
625  ! Add other dimensions.
626  ncerr = nctk_def_dims(ncid, [ &
627    nctkdim_t("nfftf", nfftf), nctkdim_t("nspden", dvdb%nspden), &
628    nctkdim_t("natom3", 3 * cryst%natom), nctkdim_t("nqlist", nqlist)], defmode=.True.)
629  NCF_CHECK(ncerr)
630 
631  ! Define arrays
632  ncerr = nctk_def_arrays(ncid, [ &
633    nctkarr_t("ngfftf", "int", "three"), &
634    nctkarr_t("qlist", "dp", "three, nqlist"), &
635    nctkarr_t("wqlist", "dp", "natom3, nqlist"), &
636    nctkarr_t("displ_cart", "dp", "two, natom3, natom3, nqlist"), &
637    nctkarr_t("v1qnu", "dp", "two, nfftf, nspden, natom3, nqlist")])
638  NCF_CHECK(ncerr)
639 
640  NCF_CHECK(nctk_set_datamode(ncid))
641  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngfftf"), ngfftf(1:3)))
642  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qlist"), qlist))
643 #endif
644 
645  ! Activate Fourier interpolation if q-points are not in the DVDB file.
646  ! TODO: handle q_bz = S q_ibz case by symmetrizing the potentials already available in the DVDB.
647  ! without performing FT interpolation.
648  do_ftv1q = 0
649  do iq=1,nqlist
650    if (dvdb_findq(dvdb, qlist(:,iq)) == -1) do_ftv1q = do_ftv1q + 1
651  end do
652  if (do_ftv1q /= 0) then
653    write(msg, "(2(a,i0),a)")"Will use Fourier interpolation of DFPT potentials [",do_ftv1q,"/",nqlist,"]"
654    call wrtout(std_out, msg)
655    call dvdb_ftinterp_setup(dvdb, ifc%ngqpt, 1, [zero,zero,zero], nfftf, ngfftf, comm)
656  end if
657 
658  ABI_MALLOC(v1qnu, (2, nfftf, dvdb%nspden, dvdb%natom3))
659 
660  do iq=1,nqlist
661    qpt = qlist(:,iq)
662    call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
663 
664    ! Find the index of the q-point in the DVDB.
665    db_iqpt = dvdb_findq(dvdb, qpt)
666 
667    ! TODO: handle q_bz = S q_ibz case by symmetrizing the potentials already available in the DVDB.
668    if (db_iqpt /= -1) then
669      if (prtvol > 0) call wrtout(std_out, sjoin("Found:", ktoa(qpt), "in DVDB with index", itoa(db_iqpt)))
670      ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
671      ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
672      call dvdb_readsym_allv1(dvdb, db_iqpt, cplex, nfftf, ngfftf, v1scf, comm)
673    else
674      if (prtvol > 0) call wrtout(std_out, sjoin("Could not find:", ktoa(qpt), "in DVDB - interpolating"))
675      ! Fourier interpolation of the potential
676      ABI_CHECK(any(abs(qpt) > tol12), "qpt cannot be zero if Fourier interpolation is used")
677      cplex = 2
678      ABI_MALLOC(v1scf, (cplex,nfftf, dvdb%nspden, dvdb%natom3))
679      call dvdb_ftinterp_qpt(dvdb, qpt, nfftf, ngfftf, v1scf, comm)
680    end if
681 
682    do nu=1,dvdb%natom3
683      ! v1qnu = \sum_{ka} phdispl{ka}(q,nu) D_{ka,q} V_scf(r)
684      ! NOTE: prefactor 1/sqrt(2 w(q,nu)) is not included in the potentials saved to file.
685      !v1qnu(2, nfftf, nspden, natom3), v1scf(cplex, nfftf, nspden, natom3)
686      v1qnu(:, :, :, nu) = zero
687      do ip=1,dvdb%natom3
688        do ispden=1,dvdb%nspden
689          if (cplex == 2) then
690            v1qnu(1, :, ispden, nu) = v1qnu(1, :, ispden, nu) + &
691              displ_red(1,ip,nu) * v1scf(1,:,ispden,ip) - displ_red(2,ip,nu) * v1scf(2,:,ispden,ip)
692            v1qnu(2, :, ispden, nu) = v1qnu(2, :, ispden, nu) + &
693              displ_red(2,ip,nu) * v1scf(1,:,ispden,ip) + displ_red(1,ip,nu) * v1scf(2,:,ispden,ip)
694          else
695            ! Gamma point. d(q) = d(-q)* --> d is real.
696            v1qnu(1, :, ispden, nu) = v1qnu(1, :, ispden, nu) + &
697              displ_red(1,ip,nu) * v1scf(1,:,ispden,ip)
698          end if
699        end do
700      end do
701    end do
702    ! Save results to file.
703 #ifdef HAVE_NETCDF
704    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wqlist"), phfrq, start=[1,iq]))
705    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "displ_cart"), displ_cart, start=[1,1,1,iq]))
706    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "v1qnu"), v1qnu, start=[1,1,1,1,iq]))
707 #endif
708    ABI_FREE(v1scf)
709  end do
710 
711  ABI_FREE(v1qnu)
712 #ifdef HAVE_NETCDF
713  NCF_CHECK(nf90_close(ncid))
714 #endif
715  call dvdb_close(dvdb)
716 
717  call wrtout(std_out, "dvqnu file written", do_flush=.True.)
718 
719 end subroutine ncwrite_v1qnu