TABLE OF CONTENTS


ABINIT/m_sigmaph [ Modules ]

[ Top ] [ Modules ]

NAME

  m_sigmaph

FUNCTION

  Compute the matrix elements of the Fan-Migdal Debye-Waller self-energy in the KS basis set.

COPYRIGHT

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

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 module m_sigmaph
26 
27  use defs_basis
28  use defs_abitypes
29  use iso_c_binding
30  use m_abicore
31  use m_xmpi
32  use m_errors
33  use m_hide_blas
34  use m_ifc
35  use m_ebands
36  use m_wfk
37  use m_ddb
38  use m_dvdb
39  use m_fft
40  use m_hamiltonian
41  use m_pawcprj
42  use m_wfd
43  use m_nctk
44 #ifdef HAVE_NETCDF
45  use netcdf
46 #endif
47 
48  use defs_datatypes,   only : ebands_t, pseudopotential_type
49  use m_time,           only : cwtime, sec2str
50  use m_fstrings,       only : itoa, ftoa, sjoin, ktoa, ltoa, strcat
51  use m_numeric_tools,  only : arth, c2r, get_diag, linfit
52  use m_io_tools,       only : iomode_from_fname
53  use m_special_funcs,  only : dirac_delta, gspline_t, gspline_new, gspline_eval, gspline_free
54  use m_fftcore,        only : ngfft_seq
55  use m_cgtools,        only : dotprod_g !set_istwfk
56  use m_cgtk,           only : cgtk_rotate
57  use m_crystal,        only : crystal_t
58  use m_crystal_io,     only : crystal_ncwrite
59  use m_kpts,           only : kpts_ibz_from_kptrlatt, kpts_timrev_from_kptopt, listkk
60  use m_lgroup,         only : lgroup_t, lgroup_new, lgroup_free
61  use m_fftcore,        only : get_kg
62  use m_kg,             only : getph
63  use m_getgh1c,        only : getgh1c, rf_transgrid_and_pack, getgh1c_setup
64  use m_pawang,         only : pawang_type
65  use m_pawrad,         only : pawrad_type
66  use m_pawtab,         only : pawtab_type
67  use m_pawfgr,         only : pawfgr_type
68  use m_fourier_interpol, only : transgrid
69  use m_symkpt,     only : symkpt
70 ! use m_paw_an,        only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
71 ! use m_paw_ij,        only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
72 ! use m_pawfgrtab,     only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
73 ! use m_pawrhoij,      only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, symrhoij
74 ! use m_pawdij,        only : pawdij, symdij
75 
76  implicit none
77 
78  private

m_sigmaph/gkknu_from_atm [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  gkknu_from_atm

FUNCTION

  Transform the gkk matrix elements from (atom, red_direction) basis to phonon-mode basis.

INPUTS

  nb1,nb2=Number of bands in gkk_atm matrix.
  nk=Number of k-points (usually 1)
  natom=Number of atoms.
  gkk_atm(2,nb1,nb2,3*natom)=EPH matrix elements in the atomic basis.
  phfrq(3*natom)=Phonon frequencies in Ha
  displ_red(2,3*natom,3*natom)=Phonon displacement in reduced coordinates.

OUTPUT

  gkk_nu(2,nb1,nb2,3*natom)=EPH matrix elements in the phonon-mode basis.
  num_smallw=Number of negative/too small frequencies that have been ignored
    by setting the corresponding gkk_nu to zero.

PARENTS

      m_sigmaph

CHILDREN

SOURCE

1241 subroutine gkknu_from_atm(nb1, nb2, nk, natom, gkk_atm, phfrq, displ_red, gkk_nu, num_smallw)
1242 
1243 
1244 !This section has been created automatically by the script Abilint (TD).
1245 !Do not modify the following lines by hand.
1246 #undef ABI_FUNC
1247 #define ABI_FUNC 'gkknu_from_atm'
1248 !End of the abilint section
1249 
1250  implicit none
1251 
1252 !Arguments ------------------------------------
1253 !scalars
1254  integer,intent(in) :: nb1, nb2, nk, natom
1255  integer,intent(out) :: num_smallw
1256 !arrays
1257  real(dp),intent(in) :: phfrq(3*natom),displ_red(2,3*natom,3*natom)
1258  real(dp),intent(in) :: gkk_atm(2,nb1,nb2,nk,3*natom)
1259  real(dp),intent(out) :: gkk_nu(2,nb1,nb2,nk,3*natom)
1260 
1261 !Local variables-------------------------
1262 !scalars
1263  integer :: nu,ipc
1264 
1265 ! *************************************************************************
1266 
1267  gkk_nu = zero; num_smallw = 0
1268 
1269  ! Loop over phonon branches.
1270  do nu=1,3*natom
1271    ! Ignore negative or too small frequencies
1272    if (phfrq(nu) < tol6) then
1273      num_smallw = num_smallw + 1; cycle
1274    end if
1275 
1276    ! Transform the gkk from atom, red_dir basis to phonon mode basis
1277    do ipc=1,3*natom
1278      gkk_nu(1,:,:,:,nu) = gkk_nu(1,:,:,:,nu) &
1279        + gkk_atm(1,:,:,:,ipc) * displ_red(1,ipc,nu) &
1280        - gkk_atm(2,:,:,:,ipc) * displ_red(2,ipc,nu)
1281      gkk_nu(2,:,:,:,nu) = gkk_nu(2,:,:,:,nu) &
1282        + gkk_atm(1,:,:,:,ipc) * displ_red(2,ipc,nu) &
1283        + gkk_atm(2,:,:,:,ipc) * displ_red(1,ipc,nu)
1284    end do
1285 
1286    gkk_nu(:,:,:,:,nu) = gkk_nu(:,:,:,:,nu) / sqrt(two * phfrq(nu))
1287  end do
1288 
1289 end subroutine gkknu_from_atm

m_sigmaph/nbe [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  nbe

FUNCTION

   Bose-Einstein statistic  1 / [(exp((e - mu)/ KT) - 1]

INPUTS

   ee=Single particle energy in Ha
   kT=Value of K_Boltzmann x T in Ha.
   mu=Chemical potential in Ha.

PARENTS

CHILDREN

SOURCE

1371 elemental real(dp) function nbe(ee, kT, mu)
1372 
1373 
1374 !This section has been created automatically by the script Abilint (TD).
1375 !Do not modify the following lines by hand.
1376 #undef ABI_FUNC
1377 #define ABI_FUNC 'nbe'
1378 !End of the abilint section
1379 
1380  implicit none
1381 
1382 !Arguments ------------------------------------
1383  real(dp),intent(in) :: ee, kT, mu
1384 
1385 !Local variables ------------------------------
1386  real(dp) :: ee_mu, arg
1387 ! *************************************************************************
1388 
1389  ee_mu = ee - mu
1390 
1391  !TODO: Find good tols.
1392  ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1393  if (kT > tol12) then
1394    arg = ee_mu / kT
1395    if (arg > tol12 .and. arg < 600._dp) then
1396      nbe = one / (exp(arg) - one)
1397    else
1398      nbe = zero
1399    end if
1400  else
1401    ! No condensate for T --> 0
1402    nbe = zero
1403  end if
1404 
1405 end function nbe

m_sigmaph/nfd [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  nfd

FUNCTION

  Fermi-Dirac statistic 1 / [(exp((e - mu)/ KT) + 1]

INPUTS

   ee=Single particle energy in Ha
   kT=Value of K_Boltzmann x T in Ha.
   mu=Chemical potential in Ha.

PARENTS

CHILDREN

SOURCE

1312 elemental real(dp) function nfd(ee, kT, mu)
1313 
1314 
1315 !This section has been created automatically by the script Abilint (TD).
1316 !Do not modify the following lines by hand.
1317 #undef ABI_FUNC
1318 #define ABI_FUNC 'nfd'
1319 !End of the abilint section
1320 
1321  implicit none
1322 
1323 !Arguments ------------------------------------
1324  real(dp),intent(in) :: ee, kT, mu
1325 
1326 !Local variables ------------------------------
1327  real(dp) :: ee_mu,arg
1328 ! *************************************************************************
1329 
1330  ee_mu = ee - mu
1331 
1332  !TODO: Find good tols.
1333  ! 1 kelvin [K] = 3.16680853419133E-06 Hartree
1334  if (kT > tol6) then
1335    arg = ee_mu / kT
1336    nfd = one / (exp(arg) + one)
1337  else
1338    ! Heaviside
1339    if (ee_mu > zero) then
1340      nfd = zero
1341    else if (ee < zero) then
1342      nfd = one
1343    else
1344      nfd = half
1345    end if
1346  end if
1347 
1348 end function nfd

m_sigmaph/sigmaph [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  sigmaph

FUNCTION

  Compute phonon-contribution to the 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

PARENTS

      eph

CHILDREN

SOURCE

 316 subroutine sigmaph(wfk0_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands,dvdb,ifc,&
 317                    pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm)
 318 
 319 
 320 !This section has been created automatically by the script Abilint (TD).
 321 !Do not modify the following lines by hand.
 322 #undef ABI_FUNC
 323 #define ABI_FUNC 'sigmaph'
 324 !End of the abilint section
 325 
 326  implicit none
 327 
 328 !Arguments ------------------------------------
 329 !scalars
 330  character(len=*),intent(in) :: wfk0_path
 331  integer,intent(in) :: comm
 332  type(datafiles_type),intent(in) :: dtfil
 333  type(dataset_type),intent(in) :: dtset
 334  type(crystal_t),intent(in) :: cryst
 335  type(ebands_t),intent(in) :: ebands
 336  type(dvdb_t),intent(inout) :: dvdb
 337  type(pawang_type),intent(in) :: pawang
 338  type(pseudopotential_type),intent(in) :: psps
 339  type(pawfgr_type),intent(in) :: pawfgr
 340  type(ifc_type),intent(in) :: ifc
 341  type(mpi_type),intent(in) :: mpi_enreg
 342 !arrays
 343  integer,intent(in) :: ngfft(18),ngfftf(18)
 344  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 345  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
 346 
 347 !Local variables ------------------------------
 348 !scalars
 349  integer,parameter :: dummy_npw=1,tim_getgh1c=1,berryopt0=0
 350  integer,parameter :: useylmgr=0,useylmgr1=0,master=0,ndat1=1
 351  integer :: my_rank,mband,my_minb,my_maxb,nsppol,nkpt,iq_ibz
 352  integer :: cplex,db_iqpt,natom,natom3,ipc,nspinor,nprocs
 353  integer :: ibsum_kq,ib_k,band,num_smallw,ibsum,ii,im,in,ndeg !,ib,nstates
 354  integer :: idir,ipert,ip1,ip2,idir1,ipert1,idir2,ipert2
 355  integer :: ik_ibz,ikq_ibz,isym_k,isym_kq,trev_k,trev_kq !,!timerev_q,
 356  integer :: spin,istwf_k,istwf_kq,istwf_kqirr,npw_k,npw_kq,npw_kqirr
 357  integer :: mpw,ierr,it !ipw
 358  integer :: n1,n2,n3,n4,n5,n6,nspden,do_ftv1q,nu
 359  integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnl1
 360  integer :: nfft,nfftf,mgfft,mgfftf,nkpg,nkpg1,nq
 361  integer :: nbcalc_ks,nbsum,bstart_ks,ikcalc
 362  real(dp),parameter :: tol_enediff=0.001_dp*eV_Ha
 363  real(dp) :: cpu,wall,gflops,cpu_all,wall_all,gflops_all
 364  real(dp) :: ecut,eshift,dotr,doti,dksqmax,weigth_q,rfact,alpha,beta,gmod2,hmod2
 365  complex(dpc) :: cfact,dka,dkap,dkpa,dkpap,my_ieta,cplx_ediff
 366  logical,parameter :: have_ktimerev=.True.
 367  logical :: isirr_k,isirr_kq,gen_eigenpb,isqzero
 368  type(wfd_t) :: wfd
 369  type(gs_hamiltonian_type) :: gs_hamkq
 370  type(rf_hamiltonian_type) :: rf_hamkq
 371  type(sigmaph_t) :: sigma
 372  type(gspline_t) :: gspl
 373  character(len=500) :: msg
 374 !arrays
 375  integer :: g0_k(3),g0_kq(3),dummy_gvec(3,dummy_npw)
 376  integer :: work_ngfft(18),gmax(3) !!g0ibz_kq(3),
 377  integer :: indkk_kq(1,6)
 378  integer,allocatable :: gtmp(:,:),kg_k(:,:),kg_kq(:,:),nband(:,:),distrib_bq(:,:),deg_ibk(:) !,degblock(:,:),
 379  real(dp) :: kk(3),kq(3),kk_ibz(3),kq_ibz(3),qpt(3),phfrq(3*cryst%natom),sqrt_phfrq0(3*cryst%natom)
 380  real(dp) :: lf(2),rg(2),res(2)
 381  real(dp) :: wqnu,nqnu,gkk2,eig0nk,eig0mk,eig0mkq,f_mkq
 382  !real(dp) :: kk(3),kq(3),kk_ibz(3),kq_ibz(3),qpt(3),phfrq(3*cryst%natom)
 383  !real(dp) :: wqnu,nqnu,gkk2,eig0nk,eig0mk,eig0mkq,ediff,f_mkq !,f_nk
 384  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom),displ_red(2,3,cryst%natom,3*cryst%natom)
 385  !real(dp) :: ucart(2,3,cryst%natom,3*cryst%natom)
 386  real(dp) :: d0mat(2,3*cryst%natom,3*cryst%natom)
 387  complex(dpc) :: cmat(3*cryst%natom,3*cryst%natom),cvec1(3*cryst%natom),cvec2(3*cryst%natom)
 388  real(dp),allocatable :: grad_berry(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),dkinpw(:)
 389  real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:),ph3d(:,:,:),ph3d1(:,:,:),v1scf(:,:,:,:)
 390  real(dp),allocatable :: gkk_atm(:,:,:),gkk_nu(:,:,:),dbwl_nu(:,:,:,:),gdw2_mn(:,:),gkk0_atm(:,:,:,:)
 391  complex(dpc),allocatable :: tpp(:,:),hka_mn(:,:,:),wmat1(:,:,:)
 392  real(dp),allocatable :: bra_kq(:,:),kets_k(:,:,:),h1kets_kq(:,:,:,:),cgwork(:,:)
 393  real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:)
 394  real(dp),allocatable :: ylm_kq(:,:),ylm_k(:,:),ylmgr_kq(:,:,:)
 395  real(dp),allocatable :: dummy_vtrial(:,:),gvnl1(:,:),work(:,:,:,:)
 396  real(dp),allocatable ::  gs1c(:,:),nqnu_tlist(:),dt_weights(:,:)
 397  complex(dpc),allocatable :: cfact_wr(:)
 398  logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
 399  type(pawcprj_type),allocatable  :: cwaveprj0(:,:)
 400 
 401 !************************************************************************
 402 
 403  if (psps%usepaw == 1) then
 404    MSG_ERROR("PAW not implemented")
 405    ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
 406  end if
 407 
 408  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
 409  call cwtime(cpu_all, wall_all, gflops_all, "start")
 410 
 411  ! Copy important dimensions
 412  natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor
 413  nspden = dtset%nspden; nkpt = ebands%nkpt
 414 
 415  nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
 416  nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
 417  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
 418  n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6)
 419 
 420  ! Get one-dimensional structure factor information on the coarse grid.
 421  ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
 422  call getph(cryst%atindx,natom,n1,n2,n3,ph1d,cryst%xred)
 423 
 424  ! Construct object to store final results.
 425  ecut = dtset%ecut ! dtset%dilatmx
 426  sigma = sigmaph_new(dtset, ecut, cryst, ebands, ifc, dtfil, comm)
 427  if (my_rank == master) call sigmaph_print(sigma, dtset, ab_out)
 428 
 429  ! This is the maximum number of PWs for all possible k+q treated.
 430  mpw = sigma%mpw; gmax = sigma%gmax
 431 
 432  ! Init work_ngfft
 433  gmax = gmax + 4 ! FIXME: this is to account for umklapp
 434  gmax = 2*gmax + 1
 435  call ngfft_seq(work_ngfft, gmax)
 436  !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
 437  ABI_MALLOC(work, (2, work_ngfft(4),work_ngfft(5),work_ngfft(6)))
 438 
 439  ! Initialize the wave function descriptor (read up to mband states where mband is defined by dtset%nband)
 440  ! For the time being, no memory distribution, each node has the full set of states.
 441  ! TODO: Distribute memory
 442  mband = dtset%mband
 443  my_minb = 1; my_maxb = mband
 444  ABI_MALLOC(nband, (nkpt, nsppol))
 445  ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
 446  ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol))
 447  nband=mband; bks_mask=.True.; keep_ur=.False.
 448 
 449  call wfd_init(wfd,cryst,pawtab,psps,keep_ur,dtset%paral_kgb,dummy_npw,mband,nband,nkpt,nsppol,bks_mask,&
 450    nspden,nspinor,dtset%ecutsm,dtset%dilatmx,ebands%istwfk,ebands%kptns,ngfft,&
 451    dummy_gvec,dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm,opt_ecut=ecut)
 452 
 453  call wfd_print(wfd,header="Wavefunctions for self-energy calculation.",mode_paral='PERS')
 454  ABI_FREE(nband)
 455  ABI_FREE(bks_mask)
 456  ABI_FREE(keep_ur)
 457 
 458  call wfd_read_wfk(wfd, wfk0_path, iomode_from_fname(wfk0_path))
 459  if (.False.) call wfd_test_ortho(wfd, cryst, pawtab, unit=std_out, mode_paral="PERS")
 460 
 461  ! TODO FOR PAW
 462  usecprj = 0
 463  ABI_DT_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
 464 
 465  ! Prepare call to getgh1c
 466  usevnl = 0
 467  optlocal = 1  ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
 468  optnl = 2     ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
 469  opt_gvnl1 = 0 ! gvnl1 is output
 470  ABI_MALLOC(gvnl1, (2,usevnl))
 471  ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4)))
 472 
 473  ! This part is taken from dfpt_vtorho
 474  !==== Initialize most of the Hamiltonian (and derivative) ====
 475  !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
 476  !2) Perform the setup needed for the non-local factors:
 477  !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
 478  !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
 479 
 480  call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,natom,&
 481 &  dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,&
 482 &  comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
 483 &  usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
 484 
 485  !PAW:allocate memory for non-symetrized 1st-order occupancies matrix (pawrhoij1)
 486  ! pawrhoij1_unsym => pawrhoij1
 487  ! if (psps%usepaw==1.and.iscf_mod>0) then
 488  !   if (paral_atom) then
 489  !     ABI_DATATYPE_ALLOCATE(pawrhoij1_unsym,(natom))
 490  !     cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=nspden
 491  !     call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,nspinor,&
 492  !       dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
 493  !   else
 494  !     pawrhoij1_unsym => pawrhoij1
 495  !     call pawrhoij_init_unpacked(pawrhoij1_unsym)
 496  !   end if
 497  ! end if
 498 
 499  ! Allocate work space arrays.
 500  ! Note nvloc in vlocal.
 501  ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data)
 502  ABI_CALLOC(dummy_vtrial, (nfftf,nspden))
 503  ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
 504  vlocal = huge(one)
 505  if (sigma%nwr > 0) then
 506    ABI_MALLOC(cfact_wr, (sigma%nwr))
 507  end if
 508  ABI_MALLOC(nqnu_tlist, (sigma%ntemp))
 509 
 510  ! Open the DVDB file
 511  call dvdb_open_read(dvdb, ngfftf, xmpi_comm_self)
 512  if (dtset%prtvol > 10) dvdb%debug = .True.
 513  ! This to symmetrize the DFPT potentials.
 514  if (dtset%symdynmat == 1) dvdb%symv1 = .True.
 515  call dvdb_print(dvdb, prtvol=dtset%prtvol)
 516 
 517  ! Compute gaussian spline.
 518  if (sigma%gfw_nomega > 0) then
 519    gspl = gspline_new(three * (sigma%gfw_mesh(2) - sigma%gfw_mesh(1)))
 520    ABI_MALLOC(dt_weights, (sigma%gfw_nomega, 2))
 521  end if
 522 
 523  ! Loop over k-points in Sigma_nk.
 524  do ikcalc=1,sigma%nkcalc
 525    kk = sigma%kcalc(:, ikcalc)
 526 
 527    ! Find IBZ(k) for q-point integration.
 528    call sigmaph_setup_kcalc(sigma, cryst, ikcalc)
 529    call wrtout(std_out, sjoin("Computing self-energy matrix elements for k-point:", ktoa(kk)))
 530    call wrtout(std_out, sjoin("Number of q-points in the IBZ(k):", itoa(sigma%nqibz_k)))
 531 
 532    ! Symmetry indices for kk.
 533    ik_ibz = sigma%kcalc2ibz(ikcalc,1); isym_k = sigma%kcalc2ibz(ikcalc,2)
 534    trev_k = sigma%kcalc2ibz(ikcalc,6); g0_k = sigma%kcalc2ibz(ikcalc,3:5)
 535    isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
 536    ABI_CHECK(isirr_k, "For the time being the k-point must be in the IBZ")
 537    kk_ibz = ebands%kptns(:,ik_ibz)
 538    npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz)
 539 
 540    ! Activate Fourier interpolation if q-points are not in the DVDB file.
 541    ! TODO: handle q_bz = S q_ibz case by symmetrizing the potentials already available in the DVDB.
 542    ! without performing FT interpolation.
 543    do_ftv1q = 0
 544    do iq_ibz=1,sigma%nqibz_k
 545      if (dvdb_findq(dvdb, sigma%qibz_k(:,iq_ibz)) == -1) do_ftv1q = do_ftv1q + 1
 546    end do
 547    if (do_ftv1q /= 0) then
 548      write(msg, "(2(a,i0),a)")"Will use Fourier interpolation of DFPT potentials [",do_ftv1q,"/",sigma%nqibz_k,"]"
 549      call wrtout(std_out, msg)
 550      call wrtout(std_out, sjoin("From ngqpt", ltoa(ifc%ngqpt), "to", ltoa(sigma%ngqpt)))
 551      call dvdb_ftinterp_setup(dvdb, ifc%ngqpt, 1, [zero,zero,zero], nfftf, ngfftf, comm)
 552    end if
 553 
 554    ! Allocate PW-arrays. Note mpw in kg_kq
 555    ABI_MALLOC(kg_k, (3, npw_k))
 556    kg_k = wfd%kdata(ik_ibz)%kg_k
 557    ABI_MALLOC(kg_kq, (3, mpw))
 558 
 559    ! Spherical Harmonics for useylm==1.
 560    ABI_MALLOC(ylm_k,(mpw, psps%mpsang**2 * psps%useylm))
 561    ABI_MALLOC(ylm_kq,(mpw, psps%mpsang**2 * psps%useylm))
 562    ABI_MALLOC(ylmgr_kq,(mpw, 3, psps%mpsang**2 * psps%useylm * useylmgr1))
 563 
 564    do spin=1,nsppol
 565      ! Bands in Sigma_nk to compute and number of bands in sum over states.
 566      bstart_ks = sigma%bstart_ks(ikcalc, spin)
 567      nbcalc_ks = sigma%nbcalc_ks(ikcalc, spin)
 568      nbsum = sigma%nbsum
 569 
 570      ! Zero self-energy matrix elements. Build frequency mesh for nk states.
 571      sigma%vals_e0ks = zero; sigma%dvals_de0ks = zero; sigma%dw_vals = zero
 572 
 573      if (sigma%nwr > 0) then
 574        sigma%vals_wr = zero
 575        do ib_k=1,nbcalc_ks
 576          band = ib_k + bstart_ks - 1
 577          eig0nk = ebands%eig(band,ik_ibz,spin) - sigma%wr_step * (sigma%nwr / 2)
 578          sigma%wrmesh_b(:,ib_k) = arth(eig0nk, sigma%wr_step, sigma%nwr)
 579        end do
 580      end if
 581      if (sigma%gfw_nomega > 0) sigma%gfw_vals = zero
 582      if (sigma%has_nuq_terms) sigma%vals_nuq = zero
 583 
 584      ! Allocate eph matrix elements.
 585      ABI_MALLOC(gkk_atm, (2, nbcalc_ks, natom3))
 586      ABI_MALLOC(gkk_nu, (2, nbcalc_ks, natom3))
 587      ABI_STAT_MALLOC(dbwl_nu, (2, nbcalc_ks, nbsum, natom3), ierr)
 588      ABI_CHECK(ierr == 0, "oom in dbwl_nu")
 589      dbwl_nu = zero
 590      ABI_MALLOC(gkk0_atm, (2, nbcalc_ks, nbsum, natom3))
 591      ABI_CHECK(ierr == 0, "oom in gkk0_atm")
 592      gkk0_atm = zero
 593 
 594      ! Load ground-state wavefunctions for which corrections are wanted (available on each node)
 595      ! TODO: symmetrize them if kk is not irred
 596      ABI_MALLOC(kets_k, (2, npw_k*nspinor, nbcalc_ks))
 597      do ib_k=1,nbcalc_ks
 598        band = ib_k + bstart_ks - 1
 599        call wfd_copy_cg(wfd, band, ik_ibz, spin, kets_k(1,1,ib_k))
 600      end do
 601 
 602      ! Continue to initialize the Hamiltonian
 603      call load_spin_hamiltonian(gs_hamkq,spin,vlocal=vlocal,with_nonlocal=.true.)
 604 
 605      ! Distribute q-points and bands.
 606      ABI_MALLOC(distrib_bq, (nbsum, sigma%nqibz_k))
 607      if (sigma%nqibz_k >= nprocs .and. mod(sigma%nqibz_k, nprocs) == 0) then
 608        do iq_ibz=1,sigma%nqibz_k
 609          distrib_bq(:, iq_ibz) = mod(iq_ibz, nprocs)
 610        end do
 611      else
 612        do it=1,nbsum*sigma%nqibz_k
 613          ibsum_kq = mod(it-1, nbsum) + 1; iq_ibz = (it - ibsum_kq) / nbsum + 1
 614          distrib_bq(ibsum_kq, iq_ibz) = mod(it, nprocs)
 615        end do
 616      end if
 617 
 618      ! Integrations over q-points in the IBZ(k)
 619      do iq_ibz=1,sigma%nqibz_k
 620        ! Quick-parallelization over q-points
 621        if (all(distrib_bq(1:nbsum, iq_ibz) /= my_rank)) cycle
 622 
 623        qpt = sigma%qibz_k(:,iq_ibz)
 624        isqzero = (sum(qpt**2) < tol14) !; if (isqzero) cycle
 625        call cwtime(cpu,wall,gflops,"start")
 626 
 627        ! Find the index of the q-point in the DVDB.
 628        db_iqpt = dvdb_findq(dvdb, qpt)
 629 
 630        ! TODO: handle q_bz = S q_ibz case by symmetrizing the potentials already available in the DVDB.
 631        if (db_iqpt /= -1) then
 632          if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found:", ktoa(qpt), "in DVDB with index", itoa(db_iqpt)))
 633          ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
 634          ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
 635          call dvdb_readsym_allv1(dvdb, db_iqpt, cplex, nfftf, ngfftf, v1scf, xmpi_comm_self)
 636        else
 637          if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Could not find:", ktoa(qpt), "in DVDB - interpolating"))
 638          ! Fourier interpolation of the potential
 639          ABI_CHECK(any(abs(qpt) > tol12), "qpt cannot be zero if Fourier interpolation is used")
 640          cplex = 2
 641          ABI_MALLOC(v1scf, (cplex,nfftf,nspden,natom3))
 642          call dvdb_ftinterp_qpt(dvdb, qpt, nfftf, ngfftf, v1scf, xmpi_comm_self)
 643          !v1scf = zero
 644        end if
 645 
 646        ! TODO: Make sure that symmetries in Q-space are preserved.
 647        ! Avoid fourq if qpt is in ddb
 648        ! Examine the symmetries of the q wavevector
 649        !call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol)
 650 
 651        ! Get phonon frequencies and displacements in reduced coordinates for this q-point
 652        call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
 653 
 654        ! Find k+q in the extended zone and extract symmetry info.
 655        ! Be careful here because there are two umklapp vectors to be considered:
 656        !
 657        !   k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz
 658        !
 659        kq = kk + qpt
 660        call listkk(dksqmax,cryst%gmet,indkk_kq,ebands%kptns,kq,ebands%nkpt,1,cryst%nsym,&
 661          1,cryst%symafm,cryst%symrel,sigma%timrev,use_symrec=.False.)
 662 
 663        if (dksqmax > tol12) then
 664          write(msg, '(3a,es16.6,7a)' )&
 665           "The WFK file cannot be used to compute self-energy corrections.",ch10,&
 666           "At least one of the k-points could not be generated from a symmetrical one. dksqmax: ",dksqmax, ch10,&
 667           "Q-mesh: ",ltoa(sigma%ngqpt),", K-mesh (from kptrlatt) ",ltoa(get_diag(dtset%kptrlatt)), ch10, &
 668           'Action: check your WFK file and (k,q) point input variables'
 669          MSG_ERROR(msg)
 670        end if
 671 
 672        ikq_ibz = indkk_kq(1,1); isym_kq = indkk_kq(1,2)
 673        trev_kq = indkk_kq(1, 6); g0_kq = indkk_kq(1, 3:5)
 674        isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0)) !; isirr_kq = .True.
 675        kq_ibz = ebands%kptns(:,ikq_ibz)
 676 
 677        ! Get npw_kq, kg_kq for k+q
 678        ! Be careful with time-reversal symmetry and istwf_kq
 679        if (isirr_kq) then
 680          ! Copy u_kq(G)
 681          istwf_kq = wfd%istwfk(ikq_ibz); npw_kq = wfd%npwarr(ikq_ibz)
 682          ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
 683          kg_kq(:,1:npw_kq) = wfd%kdata(ikq_ibz)%kg_k
 684        else
 685          ! Will Reconstruct u_kq(G) from the IBZ image.
 686          !istwf_kq = set_istwfk(kq); if (.not. have_ktimerev) istwf_kq = 1
 687          !call change_istwfk(from_npw,from_kg,from_istwfk,to_npw,to_kg,to_istwfk,n1,n2,n3,ndat1,from_cg,to_cg)
 688          istwf_kq = 1
 689          call get_kg(kq,istwf_kq,ecut,cryst%gmet,npw_kq,gtmp)
 690          ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
 691          kg_kq(:,1:npw_kq) = gtmp(:,:npw_kq)
 692          ABI_FREE(gtmp)
 693        end if
 694 
 695        ! Allocate array to store H1 |psi_nk> for all 3*natom perturbations
 696        ABI_STAT_MALLOC(h1kets_kq, (2, npw_kq*nspinor, nbcalc_ks, natom3), ierr)
 697        ABI_CHECK(ierr==0, "oom in h1kets_kq")
 698 
 699        ! Allocate vlocal1 with correct cplex. Note nvloc
 700        ABI_STAT_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr)
 701        ABI_CHECK(ierr==0, "oom in vlocal1")
 702 
 703        ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
 704        do ipc=1,natom3
 705          call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,&
 706            pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc))
 707        end do
 708        !vlocal1 = one
 709 
 710        ! if PAW, one has to solve a generalized eigenproblem
 711        ! BE careful here because I will need sij_opt==-1
 712        gen_eigenpb = (psps%usepaw==1)
 713        sij_opt = 0; if (gen_eigenpb) sij_opt = 1
 714        ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2)))
 715 
 716        ! Set up the spherical harmonics (Ylm) at k and k+q. See also dfpt_looppert
 717        !if (psps%useylm==1) then
 718        !   optder=0; if (useylmgr==1) optder=1
 719        !   call initylmg(cryst%gprimd,kg_k,kk,mkmem1,mpi_enreg,psps%mpsang,mpw,nband,mkmem1,&
 720        !     [npw_k],dtset%nsppol,optder,cryst%rprimd,ylm_k,ylmgr)
 721        !   call initylmg(cryst%gprimd,kg_kq,kq,mkmem1,mpi_enreg,psps%mpsang,mpw,nband,mkmem1,&
 722        !     [npw_kq],dtset%nsppol,optder,cryst%rprimd,ylm_kq,ylmgr_kq)
 723        !end if
 724 
 725        ! Loop over all 3*natom perturbations.
 726        ! In the inner loop, I calculate H1 * psi_k, stored in h1kets_kq on the k+q sphere.
 727        do ipc=1,natom3
 728          idir = mod(ipc-1, 3) + 1; ipert = (ipc - idir) / 3 + 1
 729 
 730          ! Prepare application of the NL part.
 731          call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.)
 732              !&paw_ij1=paw_ij1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
 733              !&mpi_spintab=mpi_enreg%my_isppoltab)
 734          call load_spin_rf_hamiltonian(rf_hamkq,spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.)
 735 
 736          ! This call is not optimal because there are quantities in out that do not depend on idir,ipert
 737          call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,&  ! In
 738            cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,&          ! In
 739            npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq,&         ! In
 740            dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1)     ! Out
 741 
 742          ! Compute H(1) applied to GS wavefunction Psi_nk(0)
 743          do ib_k=1,nbcalc_ks
 744            band = ib_k + bstart_ks - 1
 745            eig0nk = ebands%eig(band,ik_ibz,spin)
 746            ! Use scissor shift on 0-order eigenvalue
 747            eshift = eig0nk - dtset%dfpt_sciss
 748 
 749            call getgh1c(berryopt0,kets_k(:,:,ib_k),cwaveprj0,h1kets_kq(:,:,ib_k,ipc),&
 750              grad_berry,gs1c,gs_hamkq,gvnl1,idir,ipert,eshift,mpi_enreg,optlocal,&
 751              optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
 752          end do
 753 
 754          call destroy_rf_hamiltonian(rf_hamkq)
 755 
 756          ABI_FREE(kinpw1)
 757          ABI_FREE(kpg1_k)
 758          ABI_FREE(kpg_k)
 759          ABI_FREE(dkinpw)
 760          ABI_FREE(ffnlk)
 761          ABI_FREE(ffnl1)
 762          ABI_FREE(ph3d)
 763          if (allocated(ph3d1)) then
 764            ABI_FREE(ph3d1)
 765          end if
 766        end do ! ipc (loop over 3*natom atomic perturbations)
 767 
 768        ABI_FREE(gs1c)
 769 
 770        ! ==============
 771        ! Sum over bands
 772        ! ==============
 773        istwf_kqirr = wfd%istwfk(ikq_ibz); npw_kqirr = wfd%npwarr(ikq_ibz)
 774        ABI_MALLOC(bra_kq, (2, npw_kq*nspinor))
 775        ABI_MALLOC(cgwork, (2, npw_kqirr*nspinor))
 776 
 777        do ibsum_kq=1,nbsum
 778          if (distrib_bq(ibsum_kq, iq_ibz) /= my_rank) cycle
 779          ! This is to check whether the gkk elements in the degenerate subspace break symmetry
 780          !if (ibsum_kq >= bstart_ks .and. ibsum_kq <= bstart_ks + nbcalc_ks - 1) cycle
 781 
 782          ! symmetrize wavefunctions from IBZ (if needed).
 783          ! Be careful with time-reversal symmetry.
 784          if (isirr_kq) then
 785            ! Copy u_kq(G)
 786            call wfd_copy_cg(wfd, ibsum_kq, ikq_ibz, spin, bra_kq)
 787          else
 788            ! Reconstruct u_kq(G) from the IBZ image.
 789            !istwf_kq = set_istwfk(kq); if (.not. have_ktimerev) istwf_kq = 1
 790            !call change_istwfk(from_npw,from_kg,from_istwfk,to_npw,to_kg,to_istwfk,n1,n2,n3,ndat1,from_cg,to_cg)
 791 
 792            ! Use cgwork as workspace array, results stored in bra_kq
 793            !g0_kq =  g0ibz_kq + g0bz_kq
 794            call wfd_copy_cg(wfd, ibsum_kq, ikq_ibz, spin, cgwork)
 795            call cgtk_rotate(cryst, kq_ibz, isym_kq, trev_kq, g0_kq, nspinor, ndat1,&
 796                             npw_kqirr, wfd%kdata(ikq_ibz)%kg_k,&
 797                             npw_kq, kg_kq, istwf_kqirr, istwf_kq, cgwork, bra_kq, work_ngfft, work)
 798            !bra_kq = zero
 799          end if
 800 
 801          ! h1kets_kq is the main responsible for the symmetry breaking
 802          !h1kets_kq = one
 803          do ipc=1,natom3
 804            ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation.
 805            !The array eig1_k contains:
 806            !
 807            ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)>                           (NC psps)
 808            ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW)
 809            do ib_k=1,nbcalc_ks
 810              call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bra_kq,h1kets_kq(:,:,ib_k,ipc),&
 811                   mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 812              gkk_atm(:,ib_k,ipc) = [dotr, doti]
 813            end do
 814          end do
 815 
 816          ! Get gkk(kcalc, q, nu)
 817          call gkknu_from_atm(1, nbcalc_ks, 1, natom, gkk_atm, phfrq, displ_red, gkk_nu, num_smallw)
 818 
 819          ! Save data for Debye-Waller computation (performed outside the q-loop)
 820          ! dbwl_nu(2, nbcalc_ks, nbsum, natom3), gkk_nu(2, nbcalc_ks, natom3)
 821          if (isqzero) then
 822            gkk0_atm(:, :, ibsum_kq, :) = gkk_atm
 823            dbwl_nu(:, :, ibsum_kq, :) = gkk_nu
 824          end if
 825 
 826          ! Accumulate contribution to self-energy
 827          eig0mkq = ebands%eig(ibsum_kq,ikq_ibz,spin)
 828          f_mkq = ebands%occ(ibsum_kq,ikq_ibz,spin)
 829          if (nsppol == 1 .and. nspinor == 1 .and. nspden == 1) f_mkq = f_mkq * half
 830          weigth_q = sigma%wtq_k(iq_ibz)
 831 
 832          do nu=1,natom3
 833            ! Ignore acoustic or unstable modes.
 834            wqnu = phfrq(nu); if (wqnu < tol6) cycle
 835 
 836            ! Compute gaussian weights for Eliashberg function.
 837            if (sigma%gfw_nomega > 0) call gspline_eval(gspl, wqnu, sigma%gfw_nomega, sigma%gfw_mesh, dt_weights)
 838 
 839            do ib_k=1,nbcalc_ks
 840              band = ib_k + bstart_ks - 1
 841              eig0nk = ebands%eig(band,ik_ibz,spin)
 842              gkk2 = weigth_q * (gkk_nu(1,ib_k,nu) ** 2 + gkk_nu(2,ib_k,nu) ** 2)
 843 
 844              do it=1,sigma%ntemp
 845                !nqnu = one; f_mkq = one
 846                !write(std_out,*)wqnu,sigma%kTmesh(it)
 847                !write(std_out,*)eig0mkq,sigma%kTmesh(it),sigma%mu_e(it),eig0mkq-sigma%mu_e(it) / sigma%kTmesh(it)
 848                nqnu = nbe(wqnu, sigma%kTmesh(it), zero)
 849                ! TODO: Find good tolerances to treat limits
 850                !f_mkq = nfd(eig0mkq, sigma%kTmesh(it), sigma%mu_e(it))
 851                !if (nsppol == 1 .and. nspinor == 1 .and. nspden == 1) f_mkq = f_mkq * half
 852 
 853                ! Accumulate Sigma(w=eKS) for state ib_k
 854                !my_ieta = sigma%ieta * sign(one, ebands%fermie - eig0nk)
 855                my_ieta = sigma%ieta
 856                cfact =  (nqnu + f_mkq      ) / (eig0nk - eig0mkq + wqnu + my_ieta) + &
 857                         (nqnu - f_mkq + one) / (eig0nk - eig0mkq - wqnu + my_ieta)
 858 
 859                sigma%vals_e0ks(it, ib_k) = sigma%vals_e0ks(it, ib_k) + gkk2 * cfact
 860 
 861                ! Accumulate contribution to Eliashberg functions (Fan term)
 862                if (sigma%gfw_nomega > 0) then
 863                  sigma%gfw_vals(:, it, 1, ib_k) = sigma%gfw_vals(:, it, 1, ib_k) + &
 864                    gkk2 * cfact * dt_weights(:, 1)
 865                end if
 866                if (sigma%has_nuq_terms) then
 867                  !TODO: in principle iq_ibz --> iq_bz
 868                  sigma%vals_nuq(it, ib_k, nu, iq_ibz, 1) = sigma%vals_nuq(it, ib_k, nu, iq_ibz, 1) + &
 869                    gkk2 * cfact / weigth_q
 870                end if
 871 
 872                ! Accumulate dSigma(w)/dw(w=eKS) derivative for state ib_k
 873 #if 0
 874 !old derivative
 875                ! This to avoid numerical instability
 876                !my_ieta = five * sigma%ieta
 877                cfact = -(nqnu + f_mkq      ) / (eig0nk - eig0mkq + wqnu + my_ieta) ** 2 &
 878                        -(nqnu - f_mkq + one) / (eig0nk - eig0mkq - wqnu + my_ieta) ** 2
 879 
 880                ! Try to avoid numerical instability
 881                !dka = (eig0nk - eig0mkq + wqnu + my_ieta) ** 2 *  (eig0nk - eig0mkq - wqnu + my_ieta) ** 2
 882                !cfact = -(nqnu + f_mkq) * (eig0nk - eig0mkq - wqnu + my_ieta) ** 2 &
 883                !        -(nqnu - f_mkq + one) * (eig0nk - eig0mkq + wqnu + my_ieta) ** 2
 884                !cfact = cfact / dka
 885 
 886                sigma%dvals_de0ks(it, ib_k) = sigma%dvals_de0ks(it, ib_k) + gkk2 * cfact
 887 #else
 888                !cfact =  (nqnu + f_mkq      ) / (eig0nk - eig0mkq + wqnu + my_ieta) + &
 889                !         (nqnu - f_mkq + one) / (eig0nk - eig0mkq - wqnu + my_ieta)
 890                !sigma%vals_e0ks(it, ib_k) = sigma%vals_e0ks(it, ib_k) + gkk2 * cfact
 891                cfact = (eig0nk - eig0mkq + wqnu + my_ieta)
 892                gmod2 = cfact * dconjg(cfact)
 893                cfact = (eig0nk - eig0mkq - wqnu + my_ieta)
 894                hmod2 = cfact * dconjg(cfact)
 895 
 896                sigma%dvals_de0ks(it, ib_k) = sigma%dvals_de0ks(it, ib_k) + gkk2 * ( &
 897                  (nqnu + f_mkq)        * (gmod2 - two * (eig0nk - eig0mkq + wqnu) ** 2) / gmod2 ** 2 + &
 898                  (nqnu - f_mkq + one)  * (hmod2 - two * (eig0nk - eig0mkq - wqnu) ** 2) / hmod2 ** 2   &
 899                )
 900 #endif
 901 
 902                ! Accumulate Sigma(w) for state ib_k
 903                if (sigma%nwr > 0) then
 904                  my_ieta = sigma%ieta
 905                  cfact_wr(:) = (nqnu + f_mkq      ) / (sigma%wrmesh_b(:,ib_k) - eig0mkq + wqnu + my_ieta) + &
 906                                (nqnu - f_mkq + one) / (sigma%wrmesh_b(:,ib_k) - eig0mkq - wqnu + my_ieta)
 907                  sigma%vals_wr(:, it, ib_k) = sigma%vals_wr(:, it, ib_k) + gkk2 * cfact_wr(:)
 908                end if
 909              end do ! it
 910 
 911            end do ! ib_k
 912          end do ! nu
 913 
 914        end do ! ibsum_kq (sum over bands at k+q)
 915 
 916        ABI_FREE(bra_kq)
 917        ABI_FREE(cgwork)
 918        ABI_FREE(h1kets_kq)
 919        ABI_FREE(v1scf)
 920        ABI_FREE(vlocal1)
 921 
 922        call cwtime(cpu,wall,gflops,"stop")
 923        write(msg,'(2(a,i0),2(a,f8.2))')"q-point [",iq_ibz,"/",sigma%nqibz_k,"] completed. cpu:",cpu,", wall:",wall
 924        call wrtout(std_out, msg, do_flush=.True.)
 925      end do ! iq_ibz (sum over q-points)
 926 
 927      ABI_FREE(kets_k)
 928      ABI_FREE(gkk_atm)
 929      ABI_FREE(gkk_nu)
 930      ABI_FREE(distrib_bq)
 931 
 932      ! =========================
 933      ! Compute Debye-Waller term
 934      ! =========================
 935      call xmpi_sum(dbwl_nu, comm, ierr)
 936      call xmpi_sum(gkk0_atm, comm, ierr)
 937 
 938      ABI_MALLOC(deg_ibk, (nbcalc_ks))
 939      deg_ibk = 1
 940      ndeg = 1
 941 #if 0
 942      ! Count number of degeneracies.
 943      do ib_k=2,nbcalc_ks
 944        ib = ib_k + bstart_ks - 1
 945        if ( abs(ebands%eig(ib,ik_ibz,spin) - ebands%eig(ib-1,ik_ibz,spin) ) > tol_enediff) ndeg = ndeg + 1
 946      end do
 947 
 948      ! Build degblock table.
 949      ABI_MALLOC(degblock, (2, ndeg))
 950      ndeg = 1; degblock(1, 1) = 1
 951      do ib_k=2,nbcalc_ks
 952        ib = ib_k + bstart_ks - 1
 953        if ( abs(ebands%eig(ib,ik_ibz,spin) - ebands%eig(ib-1,ik_ibz,spin) ) > tol_enediff) then
 954          degblock(2, ndeg) = ib_k - 1
 955          ndeg = ndeg + 1
 956          degblock(1, ndeg) = ib_k
 957        end if
 958      end do
 959      degblock(2, ndeg) = nbcalc_ks
 960 
 961      !ABI_MALLOC(gkk0_atm, (2, nbcalc_ks, nbsum, natom3))
 962      do ii=1,ndeg
 963        !write(std_out,*)"ideg, ", ii, ", degblock", degblock(:,ii)
 964        nstates = degblock(2, ii) - degblock(1, ii) + 1
 965        do ib_k=degblock(1, ii), degblock(2, ii)
 966          deg_ibk(ib_k) = nstates
 967        end do
 968        if (nstates == 1) continue
 969        ! Use dbwl_nu as workspace
 970        dbwl_nu(:, 1, :, :) = zero
 971        do ib_k=degblock(1, ii), degblock(2, ii)
 972          dbwl_nu(:, 1, :, :) = dbwl_nu(:, 1, :, :) + gkk0_atm(:, ib_k, :, :)
 973        end do
 974        !dbwl_nu(:, 1, :, :) = dbwl_nu(:, 1, :, :) / nstates
 975        do ib_k=degblock(1, ii), degblock(2, ii)
 976          gkk0_atm(:, ib_k, :, :) = dbwl_nu(:, 1, :, :)
 977        end do
 978      end do
 979      ABI_FREE(degblock)
 980      deg_ibk(:) = 1
 981 #endif
 982 
 983      ABI_MALLOC(gdw2_mn, (nbsum, nbcalc_ks))
 984      ABI_MALLOC(tpp, (natom3, natom3))
 985      ABI_MALLOC(hka_mn, (natom3, nbsum, nbcalc_ks))
 986 
 987      qpt = zero
 988      call ifc_fourq(ifc, cryst, qpt, sqrt_phfrq0, displ_cart)
 989      where (sqrt_phfrq0 >= zero)
 990        sqrt_phfrq0 = sqrt(sqrt_phfrq0)
 991      else where
 992        sqrt_phfrq0 = zero
 993      end where
 994      d0mat = reshape(displ_cart, [2, natom3, natom3])
 995      ! cmat contains the displament vectors as complex array
 996      cmat = dcmplx(d0mat(1,:,:), d0mat(2,:,:))
 997      ! Multiply d by M to get e * M^{1/2}
 998      do nu=1,natom3
 999        do ip1=1,natom3
1000          idir1 = mod(ip1-1, 3) + 1; ipert1 = (ip1 - idir1) / 3 + 1
1001          rfact = ifc%amu(cryst%typat(ipert1)) * amu_emass !; rfact = one
1002          cmat(ip1, nu) = cmat(ip1, nu) * rfact
1003        end do
1004      end do
1005 
1006      ! Integral over IBZ. Note that here we can use IBZ(k=0).
1007      ! symsigma == 0 activates sum over full BZ for debugging purpose.
1008      ! For the time being, integrate over full BZ
1009      ! TODO Bug if symsigma 0 ?
1010      nq = sigma%nqibz
1011      if (sigma%symsigma == 0) nq = sigma%nqbz
1012      !if (sigma%symsigma == 0) nq = sigma%nqibz_k
1013      !do iq_ibz=1,sigma%nqibz
1014      !nq = sigma%nqbz
1015      do iq_ibz=1,nq
1016        if (mod(iq_ibz, nprocs) /= my_rank) cycle  ! MPI parallelism
1017 #if 0
1018        qpt = sigma%qbz(:,iq_ibz); weigth_q = one / sigma%nqbz
1019 #else
1020        if (abs(sigma%symsigma) == 1) then
1021          qpt = sigma%qibz(:,iq_ibz); weigth_q = sigma%wtq(iq_ibz)
1022          !qpt = sigma%qibz_k(:,iq_ibz); weigth_q = sigma%wtq_k(iq_ibz)
1023        else
1024          qpt = sigma%qbz(:,iq_ibz); weigth_q = one / sigma%nqbz
1025          !qpt = sigma%qibz_k(:,iq_ibz); weigth_q = sigma%wtq_k(iq_ibz)
1026        end if
1027 #endif
1028 
1029        ! Get phonons for this q-point.
1030        call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart)
1031 
1032        ! Compute hka_mn matrix with shape: (natom3, nbsum, nbcalc_ks))
1033        ! Needed for Giustino's equation.
1034        ! dbwl_nu(2, nbcalc_ks, nbsum, natom3)
1035        do in=1,nbcalc_ks
1036          do im=1,nbsum
1037            cvec1 = dcmplx(sqrt_phfrq0(:) * dbwl_nu(1, in, im, :), &
1038                           sqrt_phfrq0(:) * dbwl_nu(2, in, im, :))
1039            do ii=1,natom3
1040              cvec2 = cmat(ii,:)
1041              !hka_mn(ii, im, in) = xdotu(natom3, cvec2, 1, cvec1, 1)
1042              hka_mn(ii, im, in) = dot_product(dconjg(cvec2), cvec1)
1043              !hka_mn(ii, im, in) = dconjg(hka_mn(ii, im, in))
1044              !write(std_out,*)"hka_mn: ",hka_mn(ii, im, in)
1045            end do
1046          end do
1047        end do
1048 
1049        ! Sum over modes for this q-point.
1050        do nu=1,natom3
1051          ! Ignore acoustic or unstable modes.
1052          wqnu = phfrq(nu); if (wqnu < tol6) cycle
1053 
1054          ! Compute gaussian weights for Eliashberg function.
1055          if (sigma%gfw_nomega > 0) call gspline_eval(gspl, wqnu, sigma%gfw_nomega, sigma%gfw_mesh, dt_weights)
1056 
1057          ! Compute T_pp'(q,nu) matrix in cartesian coordinates.
1058          do ip2=1,natom3
1059            idir2 = mod(ip2-1, 3) + 1; ipert2 = (ip2 - idir2) / 3 + 1
1060            do ip1=1,natom3
1061              idir1 = mod(ip1-1, 3) + 1; ipert1 = (ip1 - idir1) / 3 + 1
1062              ! (k,a) (k,a')* + (k',a) (k',a')*
1063              dka   = dcmplx(displ_cart(1, idir1, ipert1, nu), displ_cart(2, idir1, ipert1, nu))
1064              dkap  = dcmplx(displ_cart(1, idir2, ipert1, nu), displ_cart(2, idir2, ipert1, nu))
1065              dkpa  = dcmplx(displ_cart(1, idir1, ipert2, nu), displ_cart(2, idir1, ipert2, nu))
1066              dkpap = dcmplx(displ_cart(1, idir2, ipert2, nu), displ_cart(2, idir2, ipert2, nu))
1067              tpp(ip1,ip2) = dka * dconjg(dkap) + dkpa * dconjg(dkpap)
1068              !write(std_out,*)"tpp: ",tpp(ip1, ip2)
1069            end do
1070          end do
1071 
1072          ! Giustino's equation in RMP
1073          ! gdw2_mn = diag(conjg(H) T H) / (2 w(qu,nu))
1074          ! tpp(natom3, natom3), hka_mn(natom3, nbsum, nbcalc_ks)
1075          ABI_MALLOC(wmat1, (natom3, nbsum, nbcalc_ks))
1076          call zgemm("N", "N", natom3, nbsum*nbcalc_ks, natom3, cone, tpp, natom3, hka_mn, natom3, czero, wmat1, natom3)
1077 
1078          do ibsum=1,nbsum
1079            do ib_k=1,nbcalc_ks
1080              !write(std_out,*)"maxval(wmat)",maxval(abs(wmat1(:, ibsum, ib_k)))
1081              cfact = dot_product(hka_mn(:, ibsum, ib_k), wmat1(:, ibsum, ib_k))
1082              !cfact = dot_product(wmat1(:, ibsum, ib_k), wmat1(:, ibsum, ib_k))
1083              !cfact = xdotc(natom3, hka_mn(:, ibsum, ib_k), 1, wmat1(:, ibsum, ib_k), 1)
1084              !cfact = xdotu(natom3, hka_mn(:, ibsum, ib_k), 1, wmat1(:, ibsum, ib_k), 1)
1085              gdw2_mn(ibsum, ib_k) = real(cfact) / (two * wqnu)
1086              !write(std_out, *)"gdw2_mn: ", gdw2_mn(ibsum, ib_k)
1087            end do
1088          end do
1089          ABI_FREE(wmat1)
1090          !gdw2_mn = zero
1091 
1092          ! Get phonon occupation for all temperatures.
1093          nqnu_tlist = nbe(wqnu, sigma%kTmesh(:), zero)
1094 
1095          ! Sum over bands and add (static) DW contribution for the different temperatures.
1096          do ibsum=1,nbsum
1097            eig0mk = ebands%eig(ibsum, ik_ibz, spin)
1098            do ib_k=1,nbcalc_ks
1099              band = ib_k + bstart_ks - 1
1100              eig0nk = ebands%eig(band, ik_ibz, spin)
1101              ! Handle n == m and degenerate states (either ignore or add broadening)
1102              cplx_ediff = (eig0nk - eig0mk)
1103              !if (abs(cplx_ediff) < tol6) cplx_ediff = cplx_ediff + sigma%ieta
1104              if (abs(cplx_ediff) < tol6) cycle
1105 
1106 #if 1
1107              ! Compute DW term following XG paper. Check prefactor.
1108              ! gkk0_atm(2, nbcalc_ks, nbsum, natom3)
1109              ! previous version
1110              gdw2_mn(ibsum, ib_k) = zero
1111              do ip2=1,natom3
1112                do ip1=1,natom3
1113                  gdw2_mn(ibsum, ib_k) = gdw2_mn(ibsum, ib_k) + tpp(ip1,ip2) * ( &
1114                    gkk0_atm(1, ib_k, ibsum, ip1) * gkk0_atm(1, ib_k, ibsum, ip2) + &
1115                    gkk0_atm(2, ib_k, ibsum, ip1) * gkk0_atm(2, ib_k, ibsum, ip2) + &
1116                    gkk0_atm(1, ib_k, ibsum, ip2) * gkk0_atm(1, ib_k, ibsum, ip1) + &
1117                    gkk0_atm(2, ib_k, ibsum, ip2) * gkk0_atm(2, ib_k, ibsum, ip1) &
1118                   )
1119                end do
1120 
1121                !do ip1=1,natom3
1122                !  gdw2_mn(ibsum, ib_k) = gdw2_mn(ibsum, ib_k) + tpp(ip1,ip2) * ( &
1123                !    gkk0_atm(1, ib_k, ibsum, ip1) * gkk0_atm(1, ib_k, ibsum, ip2) + &
1124                !    gkk0_atm(2, ib_k, ibsum, ip1) * gkk0_atm(2, ib_k, ibsum, ip2)  ) + &
1125                !    tpp(ip2, ip1) * ( &
1126                !    gkk0_atm(1, ib_k, ibsum, ip2) * gkk0_atm(1, ib_k, ibsum, ip1) + &
1127                !    gkk0_atm(2, ib_k, ibsum, ip2) * gkk0_atm(2, ib_k, ibsum, ip1) &
1128                !   )
1129                !end do
1130 
1131              end do
1132              !gdw2_mn(ibsum, ib_k) = gdw2_mn(ibsum, ib_k) * two
1133              gdw2_mn(ibsum, ib_k) = gdw2_mn(ibsum, ib_k) * two / deg_ibk(ib_k)
1134              !write(std_out,*)"gdw2_mn: ",gdw2_mn(ibsum, ib_k)
1135 #endif
1136              ! dbwl_nu(2, nbcalc_ks, nbsum, natom3), gkk_nu(2, nbcalc_ks, natom3)
1137 
1138              ! accumulate DW for each T, add it to Sigma(e0) and Sigma(w) as well
1139              do it=1,sigma%ntemp
1140                cfact = - weigth_q * gdw2_mn(ibsum, ib_k) * (two * nqnu_tlist(it) + one)  / cplx_ediff
1141                rfact = real(cfact)
1142                sigma%dw_vals(it, ib_k) = sigma%dw_vals(it, ib_k) + rfact
1143                sigma%vals_e0ks(it, ib_k) = sigma%vals_e0ks(it, ib_k) + rfact
1144                if (sigma%nwr > 0) sigma%vals_wr(:, it, ib_k) = sigma%vals_wr(:, it, ib_k) + rfact
1145 
1146                ! Optionally, accumulate contribution to Eliashberg functions (DW term)
1147                if (sigma%gfw_nomega > 0) then
1148                  sigma%gfw_vals(:, it, 2, ib_k) = sigma%gfw_vals(:, it, 2, ib_k) + &
1149                    rfact * dt_weights(:, 1)
1150                end if
1151 
1152                if (sigma%has_nuq_terms) then
1153                  ! TODO: in principle iq_ibz --> iq_bz if symsigma == 0
1154                  sigma%vals_nuq(it, ib_k, nu, iq_ibz, 2) = sigma%vals_nuq(it, ib_k, nu, iq_ibz, 2) + &
1155                    gkk2 * rfact / weigth_q
1156                end if
1157              end do
1158 
1159            end do
1160          end do
1161 
1162        end do ! nu
1163      end do ! iq_ibz
1164 
1165      ABI_FREE(deg_ibk)
1166      ABI_FREE(gdw2_mn)
1167      ABI_FREE(tpp)
1168      ABI_FREE(hka_mn)
1169      ABI_FREE(dbwl_nu)
1170      ABI_FREE(gkk0_atm)
1171 
1172      ! Collect results inside comm and write results for this (k-point, spin) to NETCDF file.
1173      call sigmaph_gather_and_write(sigma, ebands, ikcalc, spin, comm)
1174    end do ! spin
1175 
1176    ABI_FREE(kg_k)
1177    ABI_FREE(kg_kq)
1178    ABI_FREE(ylm_k)
1179    ABI_FREE(ylm_kq)
1180    ABI_FREE(ylmgr_kq)
1181  end do !ikcalc
1182 
1183  call cwtime(cpu_all, wall_all, gflops_all, "stop")
1184  call wrtout(std_out, "Computation of Sigma_eph completed", do_flush=.True.)
1185  call wrtout(std_out, sjoin("Total wall-time:", sec2str(cpu_all), ", Total cpu time:", sec2str(wall_all), ch10, ch10))
1186 
1187  ! Free memory
1188  ABI_FREE(gvnl1)
1189  ABI_FREE(grad_berry)
1190  ABI_FREE(dummy_vtrial)
1191  ABI_FREE(work)
1192  ABI_FREE(ph1d)
1193  ABI_FREE(vlocal)
1194  ABI_FREE(nqnu_tlist)
1195  if (sigma%nwr > 0) then
1196    ABI_FREE(cfact_wr)
1197  end if
1198  if (sigma%gfw_nomega > 0) then
1199    ABI_FREE(dt_weights)
1200  end if
1201 
1202  call gspline_free(gspl)
1203  call destroy_hamiltonian(gs_hamkq)
1204  call sigmaph_free(sigma)
1205  call wfd_free(wfd)
1206  call pawcprj_free(cwaveprj0)
1207  ABI_DT_FREE(cwaveprj0)
1208 
1209 end subroutine sigmaph

m_sigmaph/sigmaph_free [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  sigmaph_free

FUNCTION

  Deallocate dynamic memory

INPUTS

PARENTS

      m_sigmaph

CHILDREN

SOURCE

1942 subroutine sigmaph_free(self)
1943 
1944 
1945 !This section has been created automatically by the script Abilint (TD).
1946 !Do not modify the following lines by hand.
1947 #undef ABI_FUNC
1948 #define ABI_FUNC 'sigmaph_free'
1949 !End of the abilint section
1950 
1951  implicit none
1952 
1953 !Arguments ------------------------------------
1954  type(sigmaph_t),intent(inout) :: self
1955 
1956 !Local variables-------------------------------
1957  integer :: ii,jj,ideg
1958 
1959 ! *************************************************************************
1960 
1961  ! integer
1962  if (allocated(self%bstart_ks)) then
1963    ABI_FREE(self%bstart_ks)
1964  end if
1965  if (allocated(self%nbcalc_ks)) then
1966    ABI_FREE(self%nbcalc_ks)
1967  end if
1968  if (allocated(self%kcalc2ibz)) then
1969    ABI_FREE(self%kcalc2ibz)
1970  end if
1971 
1972  ! real
1973  if (allocated(self%kcalc)) then
1974    ABI_FREE(self%kcalc)
1975  end if
1976  if (allocated(self%kTmesh)) then
1977    ABI_FREE(self%kTmesh)
1978  end if
1979  if (allocated(self%mu_e)) then
1980    ABI_FREE(self%mu_e)
1981  end if
1982  if (allocated(self%wrmesh_b)) then
1983    ABI_FREE(self%wrmesh_b)
1984  end if
1985  if (allocated(self%qbz)) then
1986    ABI_FREE(self%qbz)
1987  end if
1988  if (allocated(self%qibz)) then
1989    ABI_FREE(self%qibz)
1990  end if
1991  if (allocated(self%wtq)) then
1992    ABI_FREE(self%wtq)
1993  end if
1994  if (allocated(self%qibz_k)) then
1995    ABI_FREE(self%qibz_k)
1996  end if
1997  if (allocated(self%wtq_k)) then
1998    ABI_FREE(self%wtq_k)
1999  end if
2000  if (allocated(self%gfw_mesh)) then
2001    ABI_FREE(self%gfw_mesh)
2002  end if
2003 
2004  ! complex
2005  if (allocated(self%vals_e0ks)) then
2006    ABI_FREE(self%vals_e0ks)
2007  end if
2008  if (allocated(self%dvals_de0ks)) then
2009    ABI_FREE(self%dvals_de0ks)
2010  end if
2011  if (allocated(self%dw_vals)) then
2012    ABI_FREE(self%dw_vals)
2013  end if
2014  if (allocated(self%vals_wr)) then
2015    ABI_FREE(self%vals_wr)
2016  end if
2017  if (allocated(self%gfw_vals)) then
2018    ABI_FREE(self%gfw_vals)
2019  end if
2020  if (allocated(self%vals_nuq)) then
2021    ABI_FREE(self%vals_nuq)
2022  end if
2023 
2024  ! types.
2025  if (allocated(self%degtab)) then
2026     do jj=1,size(self%degtab, dim=2)
2027       do ii=1,size(self%degtab, dim=1)
2028          do ideg=1,size(self%degtab(ii, jj)%bids)
2029            ABI_FREE(self%degtab(ii, jj)%bids(ideg)%vals)
2030          end do
2031          ABI_DT_FREE(self%degtab(ii, jj)%bids)
2032       end do
2033     end do
2034     ABI_DT_FREE(self%degtab)
2035  end if
2036 
2037  ! Close netcdf file.
2038 #ifdef HAVE_NETCDF
2039  if (self%ncid /= nctk_noid) then
2040    NCF_CHECK(nf90_close(self%ncid))
2041  end if
2042 #endif
2043 
2044 end subroutine sigmaph_free

m_sigmaph/sigmaph_gather_and_write [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  sigmaph_gather_and_write

FUNCTION

  Gather results from the MPI processes. Then master:

      1. Computes QP energies, Z factor and spectral function (if required).
      2. Saves results to file.

INPUTS

  ebands<ebands_t>=KS band energies.
  ikcalc=Index of the computed k-point
  spin=Spin index.
  comm=MPI communicator.

PARENTS

      m_sigmaph

CHILDREN

SOURCE

2139 subroutine sigmaph_gather_and_write(self, ebands, ikcalc, spin, comm)
2140 
2141 
2142 !This section has been created automatically by the script Abilint (TD).
2143 !Do not modify the following lines by hand.
2144 #undef ABI_FUNC
2145 #define ABI_FUNC 'sigmaph_gather_and_write'
2146 !End of the abilint section
2147 
2148  implicit none
2149 
2150 !Arguments ------------------------------------
2151  integer,intent(in) :: ikcalc, spin, comm
2152  type(sigmaph_t),target,intent(inout) :: self
2153  type(ebands_t),intent(in) :: ebands
2154 
2155 !Local variables-------------------------------
2156  integer,parameter :: master=0
2157  integer :: ideg,ib,it,ii,iw,nstates,ierr,my_rank,band,ik_ibz,ibc,ib_val,ib_cond
2158  real(dp) :: ravg,kse,kse_prev,dw,fan0,ks_gap,kse_val,kse_cond,qpe_adb,qpe_adb_val,qpe_adb_cond
2159  real(dp) :: smrt,alpha,beta,e0pde(9)
2160  complex(dpc) :: sig0c,zc,qpe,qpe_prev,qpe_val,qpe_cond,cavg1,cavg2
2161  character(len=500) :: msg
2162 !arrays
2163  integer :: shape3(3),shape4(4),shape5(5),shape6(6)
2164  real(dp), ABI_CONTIGUOUS pointer :: rdata3(:,:,:),rdata4(:,:,:,:),rdata5(:,:,:,:,:),rdata6(:,:,:,:,:,:)
2165  integer, ABI_CONTIGUOUS pointer :: bids(:)
2166  real(dp) :: qp_gaps(self%ntemp),qpadb_gaps(self%ntemp)
2167  real(dp),allocatable :: aw(:,:,:)
2168  complex(dpc),target :: qpadb_enes(self%ntemp, self%max_nbcalc),qp_enes(self%ntemp, self%max_nbcalc)
2169  real(dp) :: ks_enes(self%max_nbcalc),ze0_vals(self%ntemp, self%max_nbcalc)
2170 
2171 ! *************************************************************************
2172 
2173  my_rank = xmpi_comm_rank(comm)
2174 
2175  call xmpi_sum_master(self%vals_e0ks, master, comm, ierr)
2176  call xmpi_sum_master(self%dvals_de0ks, master, comm, ierr)
2177  call xmpi_sum_master(self%dw_vals, master, comm, ierr)
2178  if (self%nwr > 0) call xmpi_sum_master(self%vals_wr, master, comm, ierr)
2179  if (self%gfw_nomega > 0) call xmpi_sum_master(self%gfw_vals, master, comm, ierr)
2180  if (self%has_nuq_terms) call xmpi_sum_master(self%vals_nuq, master, comm, ierr)
2181 
2182  ! Only master writes
2183  if (my_rank /= master) return
2184 
2185  ik_ibz = self%kcalc2ibz(ikcalc, 1)
2186 
2187 #if 0
2188  ! This to compute derivative of Re Sigma with 9 points.
2189  write(ab_out, *)"Using finite derivative"
2190  iw = 1 + (self%nwr / 2) - 4
2191  e0pde = arth(zero, self%wr_step, 9)
2192  do ibc=1,self%nbcalc_ks(ikcalc,spin)
2193    do it=1,self%ntemp
2194      smrt = linfit(9, e0pde, real(self%vals_wr(iw:iw+8, it, ibc)), alpha, beta)
2195      self%dvals_de0ks(it, ibc) = alpha
2196      if (smrt > 0.1 / Ha_eV) then
2197        write(msg,'(3a,i0,a,i0,2a,2(f22.15,2a))')&
2198          'Values of Re Sig_c are not linear ',ch10,&
2199          'band index ibc = ',ibc,' spin component = ',spin,ch10,&
2200          'root mean square= ',smrt,ch10,&
2201          'estimated slope = ',alpha,ch10,&
2202          'Omega [eV] SigC [eV]'
2203        MSG_WARNING(msg)
2204      end if
2205    end do
2206  end do
2207 #endif
2208 
2209  if (self%symsigma == +1) then
2210    ! Average self-energy matrix elements in the degenerate subspace.
2211    do ideg=1,size(self%degtab(ikcalc, spin)%bids)
2212      bids => self%degtab(ikcalc, spin)%bids(ideg)%vals
2213      nstates = size(bids)
2214 
2215      do it=1,self%ntemp
2216        ! Average QP(T) and Z(T).
2217        cavg1 = sum(self%vals_e0ks(it, bids(:))) / nstates
2218        cavg2 = sum(self%dvals_de0ks(it, bids(:))) / nstates
2219        ravg = sum(self%dw_vals(it, bids(:))) / nstates
2220        do ii=1,nstates
2221          self%vals_e0ks(it, bids(ii)) = cavg1
2222          self%dvals_de0ks(it, bids(ii)) = cavg2
2223          self%dw_vals(it, bids(ii)) = ravg
2224        end do
2225 
2226        if (self%nwr > 0) then
2227          ! Average Sigma(omega, T)
2228          do iw=1,self%nwr
2229            cavg1 = sum(self%vals_wr(iw, it, bids(:))) / nstates
2230            do ii=1,nstates
2231              self%vals_wr(iw, it, bids(ii)) = cavg1
2232            end do
2233          end do
2234        end if
2235        !if (self%has_nuq_terms) then
2236        !  self%vals_nuq(it, bids(:), :, :, :) = sum(self%vals_nuq(it, bids(:), :, :, :)) / nstates
2237        !end if
2238      end do ! it
2239    end do ! ideg
2240  end if ! symsigma == +1
2241 
2242  ! Compute QP energies and Gaps.
2243  ib_val = nint(ebands%nelect / two); ib_cond = ib_val + 1
2244  kse_val = huge(one) * tol6; kse_cond = huge(one) * tol6
2245  qp_enes = huge(one) * tol6; qpadb_enes = huge(one) * tol6
2246  ks_enes = huge(one) * tol6; ze0_vals = huge(one) * tol6
2247  ks_gap = -one; qpadb_gaps = -one; qp_gaps = -one
2248 
2249  ! Write legend.
2250  if (ikcalc == 1 .and. spin == 1) then
2251    write(ab_out,"(a)")repeat("=", 80)
2252    write(ab_out,"(a)")"Final results in eV."
2253    write(ab_out,"(a)")"Notations:"
2254    write(ab_out,"(a)")"   eKS: Kohn-Sham energy. eQP: quasi-particle energy."
2255    write(ab_out,"(a)")"   eQP-eKS: Difference between the QP and the KS energy."
2256    write(ab_out,"(a)")"   SE1(eKS): Real part of the self-energy computed at the KS energy, SE2 for imaginary part."
2257    write(ab_out,"(a)")"   Z(eKS): Renormalization factor."
2258    write(ab_out,"(a)")"   FAN: Real part of the Fan term at eKS. DW: Debye-Waller term."
2259    write(ab_out,"(a)")"   DeKS: KS energy difference between this band and band-1, DeQP same meaning but for eQP."
2260    write(ab_out,"(a)")" "
2261    write(ab_out,"(a)")" "
2262  end if
2263 
2264  ! FIXME
2265  !do it=1,1
2266  do it=1,self%ntemp
2267    if (it == 1) then
2268      if (self%nsppol == 1) then
2269        write(ab_out,"(a)")sjoin("K-point:", ktoa(self%kcalc(:,ikcalc)))
2270      else
2271        write(ab_out,"(a)")sjoin("K-point:", ktoa(self%kcalc(:,ikcalc)), ", spin:", itoa(spin))
2272      end if
2273      write(ab_out,"(a)")"   B    eKS     eQP    eQP-eKS   SE1(eKS)  SE2(eKS)  Z(eKS)  FAN(eKS)   DW      DeKS     DeQP"
2274    end if
2275 
2276    do ibc=1,self%nbcalc_ks(ikcalc,spin)
2277      band = self%bstart_ks(ikcalc,spin) + ibc - 1
2278      kse = ebands%eig(band,ik_ibz,spin)
2279      ks_enes(ibc) = kse
2280      sig0c = self%vals_e0ks(it, ibc)
2281      dw = self%dw_vals(it, ibc)
2282      fan0 = real(sig0c) - dw
2283      ! TODO: Note that here I use the full Sigma including the imaginary part
2284      !zc = one / (one - self%dvals_de0ks(it, ibc))
2285      zc = one / (one - real(self%dvals_de0ks(it, ibc)))
2286      ze0_vals(it, ibc) = real(zc)
2287      qpe = kse + real(zc) * real(sig0c)
2288      qpe_adb = kse + real(sig0c)
2289      if (ibc == 1) then
2290        kse_prev = kse; qpe_prev = qpe
2291      end if
2292      if (band == ib_val) then
2293        kse_val = kse; qpe_val = qpe; qpe_adb_val = qpe_adb
2294      end if
2295      if (band == ib_cond) then
2296        kse_cond = kse; qpe_cond = qpe; qpe_adb_cond = qpe_adb
2297      end if
2298      if (it == 1) then
2299        !   B    eKS     eQP    eQP-eKS   SE1(eKS)  SE2(eKS)  Z(eKS)  FAN(eKS)   DW      DeKS     DeQP"
2300        write(ab_out, "(i4,10(f8.3,1x))") &
2301          band, kse * Ha_eV, real(qpe) * Ha_eV, (real(qpe) - kse) * Ha_eV, &
2302          real(sig0c) * Ha_eV, aimag(sig0c) * Ha_eV, real(zc), &
2303          fan0 * Ha_eV, dw * Ha_eV, (kse - kse_prev) * Ha_eV, real(qpe - qpe_prev) * Ha_eV
2304      end if
2305      if (ibc > 1) then
2306        kse_prev = kse; qpe_prev = qpe
2307      end if
2308      qpadb_enes(it, ibc) = qpe_adb
2309      qp_enes(it, ibc) = qpe
2310      if (kse_val /= huge(one) .and. kse_cond /= huge(one)) then
2311        ! We have enough states to compute the gap.
2312        if (it == 1) ks_gap = kse_cond - kse_val
2313        qpadb_gaps(it) = qpe_adb_cond - qpe_adb_val
2314        qp_gaps(it) = real(qpe_cond - qpe_val)
2315      end if
2316    end do ! ibc
2317 
2318    ! Print KS and QP gap
2319    if (it == 1 .and. kse_val /= huge(one) .and. kse_cond /= huge(one)) then
2320      write(ab_out, "(a)")" "
2321      write(ab_out, "(a,f8.3,1x,2(a,i0),a)")" KS gap: ",ks_gap * Ha_eV, &
2322        "(assuming bval:",ib_val," ==> bcond:",ib_cond,")"
2323      write(ab_out, "(2(a,f8.3),a)")" QP gap: ",qp_gaps(it) * Ha_eV," (adiabatic: ",qpadb_gaps(it) * Ha_eV, ")"
2324      write(ab_out, "(2(a,f8.3),a)")" QP_gap - KS_gap: ",(qp_gaps(it) - ks_gap) * Ha_eV,&
2325          " (adiabatic: ",(qpadb_gaps(it) - ks_gap) * Ha_eV, ")"
2326      write(ab_out, "(a)")" "
2327    end if
2328  end do ! it
2329 
2330  ! Write data to netcdf file **Only master writes**
2331  ! (use iso_c_binding to associate a real pointer to complex data because netcdf does not support complex types).
2332 
2333 #ifdef HAVE_NETCDF
2334  ! Write self-energy matrix elements for this (kpt, spin)
2335  ! Cannot use c_loc with gcc <= 4.8 due to internal compiler error so use c2r and stack memory.
2336  !shape3(1) = 2; shape4(1) = 2; shape5(1) = 2; shape6(1) = 2
2337 
2338  !shape3(2:) = shape(self%vals_e0ks); call c_f_pointer(c_loc(self%vals_e0ks), rdata3, shape3)
2339  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "vals_e0ks"), c2r(self%vals_e0ks), start=[1,1,1,ikcalc,spin]))
2340 
2341  !shape3(2:) = shape(self%dvals_de0ks); call c_f_pointer(c_loc(self%dvals_de0ks), rdata3, shape3)
2342  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "dvals_de0ks"), c2r(self%dvals_de0ks), start=[1,1,1,ikcalc,spin]))
2343 
2344  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "dw_vals"), self%dw_vals, start=[1,1,ikcalc,spin]))
2345 
2346  ! Dump QP energies and gaps for this (kpt, spin)
2347  !shape3(2:) = shape(qpadb_enes); call c_f_pointer(c_loc(qpadb_enes), rdata3, shape3)
2348  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "qpadb_enes"), c2r(qpadb_enes), start=[1,1,1,ikcalc,spin]))
2349 
2350  !shape3(2:) = shape(qp_enes); call c_f_pointer(c_loc(qp_enes), rdata3, shape3)
2351  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "qp_enes"), c2r(qp_enes), start=[1,1,1,ikcalc,spin]))
2352  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "ze0_vals"), ze0_vals, start=[1,1,ikcalc,spin]))
2353  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "ks_enes"), ks_enes, start=[1,ikcalc,spin]))
2354  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "ks_gaps"), ks_gap, start=[ikcalc,spin]))
2355  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "qpadb_gaps"), qpadb_gaps, start=[1,ikcalc,spin]))
2356  NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "qp_gaps"), qp_gaps, start=[1,ikcalc,spin]))
2357 
2358  ! Write frequency dependent data.
2359  if (self%nwr > 0) then
2360    NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "wrmesh_b"), self%wrmesh_b, start=[1,1,ikcalc,spin]))
2361 
2362    !shape4(2:) = shape(self%vals_wr); call c_f_pointer(c_loc(self%vals_wr), rdata4, shape4)
2363    NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "vals_wr"), c2r(self%vals_wr), start=[1,1,1,1,ikcalc,spin]))
2364 
2365    ! Compute spectral function.
2366    ! A = 1/pi [Im Sigma(ww)] / ([ww - ee - Re Sigma(ww)] ** 2 + Im Sigma(ww) ** 2])
2367    ABI_MALLOC(aw, (self%nwr, self%ntemp, self%max_nbcalc))
2368 
2369    do ib=1,self%nbcalc_ks(ikcalc,spin)
2370      band = self%bstart_ks(ikcalc, spin) + ib - 1
2371      kse = ebands%eig(band, ik_ibz, spin)
2372      do it=1,self%ntemp
2373        aw(:, it, ib) = piinv * abs(aimag(self%vals_wr(:, it, ib))) / &
2374          ((self%wrmesh_b(:, ib) - kse - real(self%vals_wr(:, it, ib))) ** 2 + aimag(self%vals_wr(:, it, ib)) ** 2)
2375      end do
2376    end do
2377    NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "spfunc_wr"), aw, start=[1,1,1,ikcalc,spin]))
2378    ABI_FREE(aw)
2379  end if
2380 
2381  ! Write Eliashberg functions
2382  if (self%gfw_nomega > 0) then
2383    !shape5(2:) = shape(self%gfw_vals); call c_f_pointer(c_loc(self%gfw_vals), rdata5, shape5)
2384    NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "gfw_vals"), c2r(self%gfw_vals), start=[1,1,1,1,1,ikcalc,spin]))
2385  end if
2386  if (self%has_nuq_terms) then
2387    !shape6(2:) = shape(self%vals_nuq); call c_f_pointer(c_loc(self%vals_nuq), rdata6, shape6)
2388    NCF_CHECK(nf90_put_var(self%ncid, nctk_idname(self%ncid, "vals_nuq"), c2r(self%vals_nuq), start=[1,1,1,1,1,1,ikcalc,spin]))
2389  end if
2390 #endif
2391 
2392 end subroutine sigmaph_gather_and_write

m_sigmaph/sigmaph_new [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  sigmaph_new

FUNCTION

  Creation method (allocates memory, initialize data from input vars).

INPUTS

  dtset<dataset_type>=All input variables for this dataset.
  ecut=Cutoff energy for wavefunctions.
  cryst<crystal_t>=Crystalline structure
  ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
  ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
  dtfil<datafiles_type>=variables related to files.
  comm=MPI communicator

PARENTS

CHILDREN

SOURCE

1432 type (sigmaph_t) function sigmaph_new(dtset, ecut, cryst, ebands, ifc, dtfil, comm) result(new)
1433 
1434 
1435 !This section has been created automatically by the script Abilint (TD).
1436 !Do not modify the following lines by hand.
1437 #undef ABI_FUNC
1438 #define ABI_FUNC 'sigmaph_new'
1439 !End of the abilint section
1440 
1441  implicit none
1442 
1443 !Arguments ------------------------------------
1444  integer,intent(in) :: comm
1445  real(dp),intent(in) :: ecut
1446  type(crystal_t),intent(in) :: cryst
1447  type(dataset_type),intent(in) :: dtset
1448  type(ebands_t),intent(in) :: ebands
1449  type(ifc_type),intent(in) :: ifc
1450  type(datafiles_type),intent(in) :: dtfil
1451 
1452 !Local variables ------------------------------
1453 !scalars
1454  integer,parameter :: master=0,brav1=1,occopt3=3,qptopt1=1
1455  integer :: my_rank,ik,my_nshiftq,my_mpw,cnt,nprocs,iq_ibz,ik_ibz,ndeg
1456  integer :: onpw,ii,ipw,ierr,it,spin,gap_err,ikcalc,gw_qprange,ibstop
1457  integer :: nk_found,ifo,jj
1458 #ifdef HAVE_NETCDF
1459  integer :: ncid,ncerr
1460 #endif
1461  real(dp),parameter :: spinmagntarget=-99.99_dp,tol_enediff=0.001_dp*eV_Ha
1462  real(dp) :: dksqmax,ph_wstep
1463  character(len=500) :: msg
1464  logical :: changed,found
1465  type(ebands_t) :: tmp_ebands
1466  type(gaps_t) :: gaps
1467 !arrays
1468  integer :: qptrlatt(3,3),indkk_k(1,6),my_gmax(3),kpos(6)
1469  integer :: val_indeces(ebands%nkpt, ebands%nsppol)
1470  integer,allocatable :: gtmp(:,:),degblock(:,:)
1471  real(dp) :: my_shiftq(3,1),kk(3),kq(3)
1472 
1473 ! *************************************************************************
1474 
1475  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1476 
1477  ! Copy important dimensions.
1478  new%nsppol = ebands%nsppol; new%nspinor = ebands%nspinor
1479 
1480  ! Broadening parameter from zcut
1481  new%ieta = + j_dpc * dtset%zcut
1482 
1483  ! Build (linear) mesh of K * temperatures. tsmesh(1:3) = [start, step, num]
1484  new%ntemp = nint(dtset%tmesh(3))
1485  ABI_CHECK(new%ntemp > 0, "ntemp <= 0")
1486  ABI_MALLOC(new%kTmesh, (new%ntemp))
1487  !write(std_out, *)"dtset%tmesh", dtset%tmesh
1488  new%kTmesh = arth(dtset%tmesh(1), dtset%tmesh(2), new%ntemp) * kb_HaK
1489 
1490  ! Number of bands included in self-energy summation.
1491  new%nbsum = dtset%mband
1492 
1493  gap_err = get_gaps(ebands, gaps)
1494  call gaps_print(gaps, unit=std_out)
1495  call ebands_report_gap(ebands, unit=std_out)
1496  val_indeces = get_valence_idx(ebands)
1497 
1498  ! Compute the chemical potential at the different physical temperatures with Fermi-Dirac.
1499  ABI_MALLOC(new%mu_e, (new%ntemp))
1500  new%mu_e = ebands%fermie
1501 
1502  ! TODO T == 0 >> SIGSEGV
1503 #if 0
1504  call ebands_copy(ebands, tmp_ebands)
1505  do it=1,new%ntemp
1506    call ebands_set_scheme(tmp_ebands, occopt3, new%kTmesh(it), spinmagntarget, dtset%prtvol)
1507    new%mu_e(it) = tmp_ebands%fermie
1508  end do
1509  call ebands_free(tmp_ebands)
1510 #endif
1511 
1512  ! Frequency mesh for sigma(w) and spectral function.
1513  ! TODO: Use GW variables but change default
1514  !dtset%freqspmin
1515  !dtset%freqspmax
1516  !dtset%nfreqsp
1517  new%nwr = 201
1518  ABI_CHECK(mod(new%nwr, 2) == 1, "nwr should be odd!")
1519  new%wr_step = 0.05 * eV_Ha
1520 
1521  ! Define q-mesh for integration.
1522  ! Either q-mesh from DDB (no interpolation) or eph_ngqpt_fine (Fourier interpolation if q not in DDB)
1523  new%ngqpt = dtset%ddb_ngqpt; my_nshiftq = 1; my_shiftq(:,1) = dtset%ddb_shiftq
1524  if (all(dtset%eph_ngqpt_fine /= 0)) then
1525    new%ngqpt = dtset%eph_ngqpt_fine; my_shiftq = 0
1526  end if
1527 
1528  ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
1529  qptrlatt = 0; qptrlatt(1,1) = new%ngqpt(1); qptrlatt(2,2) = new%ngqpt(2); qptrlatt(3,3) = new%ngqpt(3)
1530  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, my_nshiftq, my_shiftq, &
1531    new%nqibz, new%qibz, new%wtq, new%nqbz, new%qbz)
1532 
1533  ! Select k-point and bands where corrections are wanted
1534  ! if symsigma == 1, we have to include all degenerate states in the set
1535  ! because the final QP corrections will be obtained by averaging the results in the degenerate subspace.
1536  ! k-point and bands where corrections are wanted
1537  ! We initialize IBZ(k) here so that we have all the basic dimensions of the run and it's possible
1538  ! to distribuite the calculations among processors.
1539  new%symsigma = dtset%symsigma; new%timrev = kpts_timrev_from_kptopt(ebands%kptopt)
1540 
1541  ! TODO: debug gw_qprange > 0. Rename variable
1542  if (dtset%nkptgw /= 0) then
1543    ! Treat the k-points and bands specified in the input file.
1544    call wrtout(std_out, "Getting list of k-points for self-energy from kptgw and bdgw.")
1545    new%nkcalc = dtset%nkptgw
1546    ABI_MALLOC(new%kcalc, (3, new%nkcalc))
1547    ABI_MALLOC(new%bstart_ks, (new%nkcalc, new%nsppol))
1548    ABI_MALLOC(new%nbcalc_ks, (new%nkcalc, new%nsppol))
1549 
1550    new%kcalc = dtset%kptgw(:,1:new%nkcalc)
1551    do spin=1,new%nsppol
1552      new%bstart_ks(:,spin) = dtset%bdgw(1,1:new%nkcalc,spin)
1553      new%nbcalc_ks(:,spin) = dtset%bdgw(2,1:new%nkcalc,spin) - dtset%bdgw(1,1:new%nkcalc,spin) + 1
1554    end do
1555 
1556    ! Consistency check on bdgw and dtset%mband
1557    ierr = 0
1558    do spin=1,new%nsppol
1559      do ikcalc=1,new%nkcalc
1560        if (dtset%bdgw(2,ikcalc,spin) > dtset%mband) then
1561          ierr = ierr + 1
1562          write(msg,'(a,2i0,2(a,i0))')&
1563           "For (k, s) ",ikcalc,spin," bdgw= ",dtset%bdgw(2,ikcalc,spin), " > mband=",dtset%mband
1564          MSG_WARNING(msg)
1565        end if
1566      end do
1567    end do
1568    ABI_CHECK(ierr == 0, "Not enough bands in WFK file. See messages above. Aborting now.")
1569 
1570  else
1571    ! Use qp_range to select the interesting k-points and the corresponing bands.
1572    !
1573    !    0 --> Compute the QP corrections only for the fundamental and the optical gap.
1574    ! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num`
1575    !           bands above and below the Fermi level.
1576    ! -num --> Compute the QP corrections for all the k-points in the irreducible zone.
1577    !          Include all occupied states and `num` empty states.
1578 
1579    call wrtout(std_out, "nkptgw == 0 ==> Automatic selection of k-points and bands for the corrections.")
1580    gw_qprange = dtset%gw_qprange
1581 
1582    if (gap_err /=0 .and. gw_qprange == 0) then
1583      msg = "Problem while computing the fundamental and optical gap (likely metal). Will replace gw_qprange=0 with gw_qprange=1"
1584      MSG_WARNING(msg)
1585      gw_qprange = 1
1586    end if
1587 
1588    if (gw_qprange /= 0) then
1589      ! Include all the k-points in the IBZ.
1590      ! Note that kcalc == ebands%kptns so we can use a single ik index in the loop over k-points.
1591      ! No need to map kcalc onto ebands%kptns.
1592      new%nkcalc = ebands%nkpt
1593      ABI_MALLOC(new%kcalc, (3, new%nkcalc))
1594      ABI_MALLOC(new%bstart_ks, (new%nkcalc, new%nsppol))
1595      ABI_MALLOC(new%nbcalc_ks, (new%nkcalc, new%nsppol))
1596 
1597      new%kcalc = ebands%kptns
1598      new%bstart_ks = 1; new%nbcalc_ks = 6
1599 
1600      if (gw_qprange > 0) then
1601        ! All k-points: Add buffer of bands above and below the Fermi level.
1602        do spin=1,new%nsppol
1603          do ik=1,new%nkcalc
1604            new%bstart_ks(ik,spin) = max(val_indeces(ik,spin) - gw_qprange, 1)
1605            ii = max(val_indeces(ik,spin) + gw_qprange + 1, dtset%mband)
1606            new%nbcalc_ks(ik,spin) = ii - new%bstart_ks(ik,spin) + 1
1607          end do
1608        end do
1609 
1610      else
1611        ! All k-points: include all occupied states and -gw_qprange empty states.
1612        new%bstart_ks = 1
1613        do spin=1,new%nsppol
1614          do ik=1,new%nkcalc
1615            new%nbcalc_ks(ik,spin) = min(val_indeces(ik,spin) - gw_qprange, dtset%mband)
1616          end do
1617        end do
1618      end if
1619 
1620    else
1621      ! gw_qprange is not specified in the input.
1622      ! Include the optical and the fundamental KS gap.
1623      ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore
1624      ! we have compute the union of the k-points where the fundamental and the optical gaps are located.
1625      !
1626      ! Find the list of `interesting` kpoints.
1627      ABI_CHECK(gap_err == 0, "gw_qprange 0 cannot be used because I cannot find the gap (gap_err !=0)")
1628      nk_found = 1; kpos(1) = gaps%fo_kpos(1,1)
1629 
1630      do spin=1,new%nsppol
1631        do ifo=1,3
1632          ik = gaps%fo_kpos(ifo, spin)
1633          found = .False.; jj = 0
1634          do while (.not. found .and. jj < nk_found)
1635            jj = jj + 1; found = (kpos(jj) == ik)
1636          end do
1637          if (.not. found) then
1638            nk_found = nk_found + 1; kpos(nk_found) = ik
1639          end if
1640        end do
1641      end do
1642 
1643      ! Now we can define the list of k-points and the bands range.
1644      new%nkcalc = nk_found
1645      ABI_MALLOC(new%kcalc, (3, new%nkcalc))
1646      ABI_MALLOC(new%bstart_ks, (new%nkcalc, new%nsppol))
1647      ABI_MALLOC(new%nbcalc_ks, (new%nkcalc, new%nsppol))
1648 
1649      do ii=1,new%nkcalc
1650        ik = kpos(ii)
1651        new%kcalc(:,ii) = ebands%kptns(:,ik)
1652        do spin=1,new%nsppol
1653          new%bstart_ks(ii,spin) = val_indeces(ik,spin)
1654          new%nbcalc_ks(ii,spin) = 2
1655        end do
1656      end do
1657    end if ! gw_qprange
1658 
1659  end if ! nkptgw /= 0
1660 
1661  call gaps_free(gaps)
1662 
1663  ! The k-point and the symmetries relating the BZ k-point to the IBZ.
1664  ABI_MALLOC(new%kcalc2ibz, (new%nkcalc, 6))
1665  if (new%symsigma == 1) then
1666    ABI_DT_MALLOC(new%degtab, (new%nkcalc, new%nsppol))
1667  end if
1668 
1669  do ikcalc=1,new%nkcalc
1670    kk = new%kcalc(:,ikcalc)
1671    call listkk(dksqmax,cryst%gmet,indkk_k,ebands%kptns,kk,ebands%nkpt,1,cryst%nsym,&
1672       1,cryst%symafm,cryst%symrel,new%timrev,use_symrec=.False.)
1673 
1674    new%kcalc2ibz(ikcalc, :) = indkk_k(1, :)
1675 
1676    ik_ibz = indkk_k(1,1) !; isym_k = indkk_k(1,2)
1677    !trev_k = indkk_k(1, 6); g0_k = indkk_k(1, 3:5)
1678    !isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
1679    !kk_ibz = ebands%kptns(:,ik_ibz)
1680 
1681    if (dksqmax > tol12) then
1682       write(msg, '(4a,es16.6,7a)' )&
1683        "The WFK file cannot be used to compute self-energy corrections at kpoint: ",ktoa(kk),ch10,&
1684        "the k-point could not be generated from a symmetrical one. dksqmax: ",dksqmax, ch10,&
1685        "Q-mesh: ",ltoa(new%ngqpt),", K-mesh (from kptrlatt) ",ltoa(get_diag(dtset%kptrlatt)), ch10, &
1686        'Action: check your WFK file and (k,q) point input variables'
1687       MSG_ERROR(msg)
1688    end if
1689 
1690    ! Find the little-group of the k-point and initialize the weights for BZ integration.
1691    ! Examine the symmetries of the k-point
1692    ! TODO: See timerev_k
1693    !call littlegroup_q(cryst%nsym,kk,symk,cryst%symrec,cryst%symafm,timerev_k,prtvol=dtset%prtvol)
1694    !call wrtout(std_out, sjoin("Will compute", itoa(sigma%nqibz_k), "q-points in the IBZ(k)"))
1695 
1696    ! We will have to average the QP corrections over degenerate states if symsigma=1 is used.
1697    ! Here we make sure that all the degenerate states are included.
1698    ! Store also band indices of the degenerate sets, used to average final results.
1699    if (abs(new%symsigma) == 1) then
1700      do spin=1,new%nsppol
1701        ibstop = new%bstart_ks(ikcalc,spin) + new%nbcalc_ks(ikcalc,spin) - 1
1702        call enclose_degbands(ebands, ik_ibz, spin, new%bstart_ks(ikcalc,spin), ibstop, &
1703          changed, tol_enediff, degblock=degblock)
1704        if (changed) then
1705          new%nbcalc_ks(ikcalc,spin) = ibstop - new%bstart_ks(ikcalc,spin) + 1
1706          write(msg,'(2(a,i0),2a,2(1x,i0))')&
1707            "Not all the degenerate states at ikcalc= ",ikcalc,", spin= ",spin,ch10,&
1708            "were included in the bdgw set. bdgw has been changed to: ",new%bstart_ks(ikcalc,spin),ibstop
1709          MSG_COMMENT(msg)
1710        end if
1711        ! Store band indices used for averaging.
1712        ndeg = size(degblock, dim=2)
1713        ABI_DT_MALLOC(new%degtab(ikcalc, spin)%bids, (ndeg))
1714        do ii=1,ndeg
1715          cnt = degblock(2, ii) - degblock(1, ii) + 1
1716          ABI_DT_MALLOC(new%degtab(ikcalc, spin)%bids(ii)%vals, (cnt))
1717          new%degtab(ikcalc, spin)%bids(ii)%vals = [(jj, jj=degblock(1, ii), degblock(2, ii))]
1718        end do
1719        ABI_FREE(degblock)
1720      end do
1721    end if ! symsigma
1722 
1723  end do
1724 
1725  ! Now we can finally compute max_nbcalc
1726  new%max_nbcalc = maxval(new%nbcalc_ks)
1727 
1728  ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
1729  ! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
1730  ! used to symmetrize the wavefunctions in G-space.
1731  ! TODO: Should loop over IBZ(k)
1732  new%mpw = 0; new%gmax = 0; cnt = 0
1733  do ik=1,new%nkcalc
1734    kk = new%kcalc(:, ik)
1735    !call sigmaph_setup_kcalc(self, crystal, ik)
1736    do iq_ibz=1,new%nqbz
1737    !do iq_ibz=1,new%nqibz_k
1738      cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle
1739      !kq = kk + qibz_k(:,iq_ibz)
1740      kq = kk + new%qbz(:,iq_ibz)
1741 
1742      ! TODO: g0 umklapp here can enter into play!
1743      ! fstab should contains the max of the umlapp G-vectors.
1744      ! gmax could not be large enough!
1745      call get_kg(kq,1,ecut,cryst%gmet,onpw,gtmp)
1746      new%mpw = max(new%mpw, onpw)
1747      do ipw=1,onpw
1748        do ii=1,3
1749         new%gmax(ii) = max(new%gmax(ii), abs(gtmp(ii,ipw)))
1750        end do
1751      end do
1752      ABI_FREE(gtmp)
1753    end do
1754  end do
1755 
1756  my_mpw = new%mpw; call xmpi_max(my_mpw, new%mpw, comm, ierr)
1757  my_gmax = new%gmax; call xmpi_max(my_gmax, new%gmax, comm, ierr)
1758  call wrtout(std_out, sjoin('optimal value of mpw=', itoa(new%mpw)))
1759 
1760  ! ================================================================
1761  ! Allocate arrays used to store final results and set them to zero
1762  ! ================================================================
1763  ABI_CALLOC(new%vals_e0ks, (new%ntemp, new%max_nbcalc))
1764  ABI_CALLOC(new%dvals_de0ks, (new%ntemp, new%max_nbcalc))
1765  ABI_CALLOC(new%dw_vals, (new%ntemp, new%max_nbcalc))
1766 
1767  ! Frequency dependent stuff
1768  if (new%nwr > 0) then
1769    ABI_CALLOC(new%vals_wr, (new%nwr, new%ntemp, new%max_nbcalc))
1770    ABI_CALLOC(new%wrmesh_b, (new%nwr, new%max_nbcalc))
1771  end if
1772 
1773  ! Prepare calculation of generalized Eliashberg functions.
1774  new%gfw_nomega = 0; ph_wstep = one/Ha_cmm1
1775 
1776  if (new%gfw_nomega /= 0 .and. ph_wstep > zero) then
1777    new%gfw_nomega = nint((ifc%omega_minmax(2) - ifc%omega_minmax(1) ) / ph_wstep) + 1
1778    ABI_MALLOC(new%gfw_mesh, (new%gfw_nomega))
1779    new%gfw_mesh = arth(ifc%omega_minmax(1), ph_wstep, new%gfw_nomega)
1780    ABI_MALLOC(new%gfw_vals, (new%gfw_nomega, new%ntemp, 2, new%max_nbcalc))
1781  end if
1782 
1783  new%has_nuq_terms = .False.
1784  if (new%has_nuq_terms) then
1785    ABI_CALLOC(new%vals_nuq, (new%ntemp, new%max_nbcalc, 3*cryst%natom, new%nqbz, 2))
1786  end if
1787 
1788  ! Open netcdf file (only master work for the time being because cannot assume HDF5 + MPI-IO)
1789  ! This could create problems if MPI parallelism over (spin, nkptgw) ...
1790 #ifdef HAVE_NETCDF
1791  if (my_rank == master) then
1792    ! Master creates the netcdf file used to store the results of the calculation.
1793    NCF_CHECK(nctk_open_create(new%ncid, strcat(dtfil%filnam_ds(4), "_SIGEPH.nc"), xmpi_comm_self))
1794    ncid = new%ncid
1795 
1796    NCF_CHECK(crystal_ncwrite(cryst, ncid))
1797    NCF_CHECK(ebands_ncwrite(ebands, ncid))
1798 
1799    ! Add sigma_eph dimensions.
1800    ncerr = nctk_def_dims(ncid, [ &
1801      nctkdim_t("nkcalc", new%nkcalc), nctkdim_t("max_nbcalc", new%max_nbcalc), &
1802      nctkdim_t("nsppol", new%nsppol), nctkdim_t("ntemp", new%ntemp), nctkdim_t("natom3", 3 * cryst%natom), &
1803      nctkdim_t("nqibz", new%nqibz), nctkdim_t("nqbz", new%nqbz)], &
1804      defmode=.True.)
1805    NCF_CHECK(ncerr)
1806 
1807    if (new%nwr > 0) then
1808      NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("nwr", new%nwr)]))
1809    end if
1810    if (new%gfw_nomega > 0) then
1811      NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("gfw_nomega", new%gfw_nomega)]))
1812    end if
1813 
1814    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "symsigma", "nbsum"])
1815    NCF_CHECK(ncerr)
1816    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat"])
1817    NCF_CHECK(ncerr)
1818    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "ph_intmeth"], defmode=.True.)
1819    NCF_CHECK(ncerr)
1820    ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "eta", "wr_step"])
1821    NCF_CHECK(ncerr)
1822    ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie"])
1823    NCF_CHECK(ncerr)
1824    ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ph_wstep", "ph_smear"])
1825    NCF_CHECK(ncerr)
1826 
1827    ! Define arrays for results.
1828    ncerr = nctk_def_arrays(ncid, [ &
1829      nctkarr_t("ngqpt", "int", "three"), &
1830      nctkarr_t("eph_ngqpt_fine", "int", "three"), &
1831      nctkarr_t("ddb_ngqpt", "int", "three"), &
1832      nctkarr_t("ph_ngqpt", "int", "three"), &
1833      nctkarr_t("bstart_ks", "int", "nkcalc, nsppol"), &
1834      nctkarr_t("nbcalc_ks", "int", "nkcalc, nsppol"), &
1835      nctkarr_t("kcalc", "dp", "three, nkcalc"), &
1836      nctkarr_t("kcalc2ibz", "int", "nkcalc, six"), &
1837      nctkarr_t("kTmesh", "dp", "ntemp"), &
1838      nctkarr_t("mu_e", "dp", "ntemp"), &
1839      nctkarr_t("vals_e0ks", "dp", "two, ntemp, max_nbcalc, nkcalc, nsppol"), &
1840      nctkarr_t("dvals_de0ks", "dp", "two, ntemp, max_nbcalc, nkcalc, nsppol"), &
1841      nctkarr_t("dw_vals", "dp", "ntemp, max_nbcalc, nkcalc, nsppol"), &
1842      nctkarr_t("qpadb_enes", "dp", "two, ntemp, max_nbcalc, nkcalc, nsppol"), &
1843      nctkarr_t("qp_enes", "dp", "two, ntemp, max_nbcalc, nkcalc, nsppol"), &
1844      nctkarr_t("ze0_vals", "dp", "ntemp, max_nbcalc, nkcalc, nsppol"), &
1845      nctkarr_t("ks_enes", "dp", "max_nbcalc, nkcalc, nsppol"), &
1846      nctkarr_t("ks_gaps", "dp", "nkcalc, nsppol"), &
1847      nctkarr_t("qpadb_gaps", "dp", "ntemp, nkcalc, nsppol"), &
1848      nctkarr_t("qp_gaps", "dp", "ntemp, nkcalc, nsppol") &
1849    ])
1850    NCF_CHECK(ncerr)
1851 
1852    ! TODO: Check the mesh
1853    if (new%nwr > 0) then
1854      ! These arrays get two extra dimensions on file (nkcalc, nsppol).
1855      ncerr = nctk_def_arrays(ncid, [ &
1856        nctkarr_t("wrmesh_b", "dp", "nwr, max_nbcalc, nkcalc, nsppol"), &
1857        nctkarr_t("vals_wr", "dp", "two, nwr, ntemp, max_nbcalc, nkcalc, nsppol") &
1858      ])
1859      NCF_CHECK(ncerr)
1860    end if
1861    if (new%gfw_nomega > 0) then
1862      ncerr = nctk_def_arrays(ncid, [ &
1863        nctkarr_t("gfw_mesh", "dp", "gfw_nomega"), &
1864        nctkarr_t("gfw_vals", "dp", "two, gfw_nomega, ntemp, two, max_nbcalc, nkcalc, nsppol") &
1865      ])
1866      NCF_CHECK(ncerr)
1867    end if
1868    if (new%has_nuq_terms) then
1869      ncerr = nctk_def_arrays(ncid, [ &
1870        nctkarr_t("vals_nuq", "dp", "two, ntemp, max_nbcalc, natom3, nqbz, two, nkcalc, nsppol") &
1871      ])
1872      NCF_CHECK(ncerr)
1873    end if
1874 
1875    if (new%nwr > 1) then
1876      ! Make room for the spectral function.
1877      ncerr = nctk_def_arrays(ncid, [nctkarr_t("spfunc_wr", "dp", "nwr, ntemp, max_nbcalc, nkcalc, nsppol")])
1878      NCF_CHECK(ncerr)
1879    end if
1880 
1881    ! ======================================================
1882    ! Write data that do not depend on the (kpt, spin) loop.
1883    ! ======================================================
1884    NCF_CHECK(nctk_set_datamode(ncid))
1885    ncerr = nctk_write_iscalars(ncid, &
1886        [character(len=nctk_slen) :: "symsigma", "nbsum", "symdynmat", "ph_intmeth", "eph_intmeth", "eph_transport"], &
1887        [new%symsigma, new%nbsum, dtset%symdynmat, dtset%ph_intmeth, dtset%eph_intmeth, dtset%eph_transport])
1888    NCF_CHECK(ncerr)
1889    ncerr = nctk_write_dpscalars(ncid, &
1890      [character(len=nctk_slen) :: "eta", "wr_step"], &
1891      [aimag(new%ieta), new%wr_step])
1892    NCF_CHECK(ncerr)
1893 
1894    ncerr = nctk_write_dpscalars(ncid, &
1895      [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", "ph_wstep", "ph_smear"], &
1896      [dtset%eph_fsewin, dtset%eph_fsmear, dtset%eph_extrael, dtset%eph_fermie, dtset%ph_wstep, dtset%ph_smear])
1897    NCF_CHECK(ncerr)
1898 
1899    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngqpt"), new%ngqpt))
1900    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_ngqpt_fine"), dtset%eph_ngqpt_fine))
1901    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ddb_ngqpt"), dtset%ddb_ngqpt))
1902    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ph_ngqpt"), dtset%ph_ngqpt))
1903    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "bstart_ks"), new%bstart_ks))
1904    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nbcalc_ks"), new%nbcalc_ks))
1905    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc"), new%kcalc))
1906    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ibz"), new%kcalc2ibz))
1907    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kTmesh"), new%kTmesh))
1908    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mu_e"), new%mu_e))
1909    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eta"), aimag(new%ieta)))
1910    if (new%gfw_nomega > 0) then
1911      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gfw_mesh"), new%gfw_mesh))
1912    end if
1913    NCF_CHECK(nf90_close(ncid))
1914  end if ! master
1915 
1916  ! Now reopen the file in parallel.
1917  call xmpi_barrier(comm)
1918  NCF_CHECK(nctk_open_modify(new%ncid, strcat(dtfil%filnam_ds(4), "_SIGEPH.nc"), xmpi_comm_self))
1919  !NCF_CHECK(nctk_open_modify(new%ncid, strcat(dtfil%filnam_ds(4), "_SIGEPH.nc"), comm))
1920  NCF_CHECK(nctk_set_datamode(new%ncid))
1921 #endif
1922 
1923 end function sigmaph_new

m_sigmaph/sigmaph_print [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  sigmaph_print

FUNCTION

  Print self-energy and QP corrections for given (k-point, spin).

INPUTS

  dtset<dataset_type>=All input variables for this dataset.
  unt=Fortran unit number

PARENTS

      m_sigmaph

CHILDREN

SOURCE

2413 subroutine sigmaph_print(self, dtset, unt)
2414 
2415 
2416 !This section has been created automatically by the script Abilint (TD).
2417 !Do not modify the following lines by hand.
2418 #undef ABI_FUNC
2419 #define ABI_FUNC 'sigmaph_print'
2420 !End of the abilint section
2421 
2422  implicit none
2423 
2424 !Arguments ------------------------------------
2425  integer,intent(in) :: unt
2426  type(dataset_type),intent(in) :: dtset
2427  type(sigmaph_t),intent(in) :: self
2428 
2429 !Local variables-------------------------------
2430 !integer
2431  integer :: ikc,is
2432 
2433 ! *************************************************************************
2434 
2435  ! Write dimensions
2436  write(unt,"(a)")sjoin("Number of bands in e-ph self-energy:", itoa(self%nbsum))
2437  write(unt,"(a)")sjoin("Symsigma: ",itoa(self%symsigma), "Timrev:",itoa(self%timrev))
2438  write(unt,"(a)")sjoin("Imaginary shift in the denominator (zcut): ", ftoa(aimag(self%ieta) * Ha_eV, fmt="f5.3"), "[eV]")
2439  write(unt,"(a)")sjoin("Number of frequencies along the real axis:", itoa(self%nwr), ", Step:", ftoa(self%wr_step*Ha_eV), "[eV]")
2440  write(unt,"(a)")sjoin("Number of temperatures:", itoa(self%ntemp), &
2441    "From:", ftoa(self%kTmesh(1) / kb_HaK), "to", ftoa(self%kTmesh(self%ntemp) / kb_HaK), "[K]")
2442  write(unt,"(a)")sjoin("Ab-initio q-mesh from DDB file:", ltoa(dtset%ddb_ngqpt))
2443  write(unt,"(a)")sjoin("Q-mesh used for self-energy integration [ngqpt]:", ltoa(self%ngqpt))
2444  write(unt,"(a)")sjoin("Number of q-points in the IBZ:", itoa(self%nqibz))
2445  write(unt,"(a)")"List of K-points for self-energy corrections:"
2446  do ikc=1,self%nkcalc
2447    do is=1,self%nsppol
2448    if (self%nsppol == 2) write(unt,"(a,i1)")"... For spin: ",is
2449      write(unt, "(2(i4,2x),a,2(i4,1x))")&
2450        ikc, is, trim(ktoa(self%kcalc(:,ikc))), self%bstart_ks(ikc,is), self%bstart_ks(ikc,is) + self%nbcalc_ks(ikc,is) - 1
2451    end do
2452  end do
2453 
2454 end subroutine sigmaph_print

m_sigmaph/sigmaph_setup_kcalc [ Functions ]

[ Top ] [ m_sigmaph ] [ Functions ]

NAME

  sigmaph_setup_kcalc

FUNCTION

  Prepare calculations of self-energy matrix elements for ikcalc index.

INPUTS

  crystal<crystal_t> = Crystal structure.
  ikcalc=Index of the k-point to compute.

PARENTS

      m_sigmaph

CHILDREN

SOURCE

2065 subroutine sigmaph_setup_kcalc(self, cryst, ikcalc)
2066 
2067 
2068 !This section has been created automatically by the script Abilint (TD).
2069 !Do not modify the following lines by hand.
2070 #undef ABI_FUNC
2071 #define ABI_FUNC 'sigmaph_setup_kcalc'
2072 !End of the abilint section
2073 
2074  implicit none
2075 
2076 !Arguments ------------------------------------
2077  integer,intent(in) :: ikcalc
2078  type(crystal_t),intent(in) :: cryst
2079  type(sigmaph_t),intent(inout) :: self
2080 
2081 !Local variables-------------------------------
2082  type(lgroup_t) :: lgk
2083 
2084 ! *************************************************************************
2085 
2086  if (allocated(self%qibz_k)) then
2087    ABI_FREE(self%qibz_k)
2088  end if
2089  if (allocated(self%wtq_k)) then
2090    ABI_FREE(self%wtq_k)
2091  end if
2092 
2093  if (self%symsigma == 0) then
2094    ! Do not use symmetries in BZ sum_q --> nqibz_k == nqbz
2095    self%nqibz_k = self%nqbz
2096    ABI_MALLOC(self%qibz_k, (3, self%nqibz_k))
2097    ABI_MALLOC(self%wtq_k, (self%nqibz_k))
2098    self%qibz_k = self%qbz; self%wtq_k = one / self%nqbz
2099 
2100  else if (abs(self%symsigma) == 1) then
2101    ! Use the symmetries of the little group
2102    lgk = lgroup_new(cryst, self%kcalc(:, ikcalc), self%timrev, self%nqbz, self%qbz, self%nqibz, self%qibz)
2103 
2104    self%nqibz_k = lgk%nibz
2105    ABI_MALLOC(self%qibz_k, (3, self%nqibz_k))
2106    ABI_MALLOC(self%wtq_k, (self%nqibz_k))
2107    self%qibz_k = lgk%ibz; self%wtq_k = lgk%weights
2108    call lgroup_free(lgk)
2109  else
2110    MSG_ERROR(sjoin("Wrong symsigma:", itoa(self%symsigma)))
2111  end if
2112 
2113 end subroutine sigmaph_setup_kcalc

m_sigmaph/sigmaph_t [ Types ]

[ Top ] [ m_sigmaph ] [ Types ]

NAME

 sigmaph_t

FUNCTION

 Container for the (diagonal) matrix elements of the electronic self-energy (phonon contribution)
 computed in the KS representation i.e. Sigma_eph(omega, T, band, k, spin).
 Provides methods to compute the QP corrections, the spectral functions and save the results to file.

TODO

  1) Store (q, mu) or q-resolved contributions? More memory but could be useful for post-processing
  2) Use list of i.delta values instead of single shift? Useful for convergence studies.

SOURCE

110  type,private :: sigmaph_t
111 
112   integer :: nkcalc
113    ! Number of k-points computed.
114 
115   integer :: max_nbcalc
116    ! Maximum number of bands computed (max over nkcalc and spin).
117 
118   integer :: nsppol
119    ! Number of independent spin polarizations.
120 
121   integer :: nspinor
122    ! Number of spinorial components.
123 
124   !integer :: nsab
125    ! Number of spin components in self-energy matrix elements.
126    !   1 if (nspinor=1, nsppol=1),
127    !   2 if (nspinor=1, nsppol=2),
128    !   4 if (nspinor=2)
129 
130   integer :: nwr
131    ! Number of frequency points along the real axis for Sigma(w) and spectral function A(w)
132    ! Odd number so that the mesh is centered on the KS energy.
133    ! The spectral function is computed only if nwr > 1.
134 
135   integer :: ntemp
136    ! Number of temperatures.
137 
138   integer :: symsigma
139    ! 1 if matrix elements should be symmetrized.
140    ! Required when the sum over q in the BZ is replaced by IBZ(k).
141 
142   integer :: timrev
143    !  timrev=1 if the use of time-reversal is allowed; 0 otherwise
144 
145   integer :: nbsum
146    ! Number of bands used in sum over states.
147 
148   integer :: nqbz
149   ! Number of q-points in the BZ.
150 
151   integer :: nqibz
152   ! Number of q-points in the IBZ.
153 
154   integer :: nqibz_k
155   ! Number of q-points in the IBZ(k). Depends on ikcalc.
156 
157   integer :: ncid = nctk_noid
158   ! Netcdf file used to save results.
159 
160   integer :: mpw
161   ! Maximum number of PWs for all possible k+q
162 
163   complex(dpc) :: ieta
164    ! Used to shift the poles in the complex plane (Ha units)
165    ! Corresponds to `i eta` term in equations.
166 
167   real(dp) :: wr_step
168    ! Step of the linear mesh along the real axis (Ha units).
169 
170   integer :: gmax(3)
171 
172   integer :: ngqpt(3)
173    ! Number of divisions in the Q mesh in the BZ.
174 
175   integer,allocatable :: bstart_ks(:,:)
176    ! bstart_ks(nkcalc,nsppol)
177    ! Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
178    ! Depends on spin because all denerate states should be included when symmetries are used.
179 
180   integer,allocatable :: nbcalc_ks(:,:)
181    ! nbcalc_ks(nkcalc,nsppol)
182    ! Number of bands included in self-energy matrix elements for each k-point in kcalc.
183    ! Depends on spin because all denerate states should be included when symmetries are used.
184 
185   integer,allocatable :: kcalc2ibz(:,:)
186     !kcalc2ibz(nkcalc, 6))
187     ! Mapping kcalc --> ibz as reported by listkk.
188 
189   real(dp),allocatable :: kcalc(:,:)
190    ! kcalc(3, nkcalc)
191    ! List of k-points where the self-energy is computed.
192 
193   real(dp),allocatable :: qbz(:,:)
194   ! qbz(3,nqbz)
195   ! Reduced coordinates of the q-points in the full BZ.
196 
197   real(dp),allocatable :: qibz(:,:)
198   ! qibz(3,nqibz)
199   ! Reduced coordinates of the q-points in the IBZ.
200 
201   real(dp),allocatable :: wtq(:)
202   ! wtq(nqibz)
203   ! Weights of the q-points in the IBZ (normalized to one).
204 
205   real(dp),allocatable :: qibz_k(:,:)
206   ! qibz(3,nqibz_k)
207   ! Reduced coordinates of the q-points in the IBZ(k). Depends on ikcalc.
208 
209   real(dp),allocatable :: wtq_k(:)
210   ! wtq(nqibz_k)
211   ! Weights of the q-points in the IBZ(k) (normalized to one). Depends on ikcalc.
212 
213   real(dp),allocatable :: kTmesh(:)
214    ! kTmesh(ntemp)
215    ! List of temperatures (kT units).
216 
217   real(dp),allocatable :: mu_e(:)
218    ! mu_e(ntemp)
219    ! chemical potential of electrons for the different temperatures.
220 
221   real(dp),allocatable :: wrmesh_b(:,:)
222   ! wrmesh_b(nwr, max_nbcalc)
223   ! Frequency mesh along the real axis (Ha units) used for the different bands
224   ! This array depends on (ikcalc, spin)
225 
226   complex(dpc),allocatable :: vals_e0ks(:,:)
227    ! vals_e0ks(ntemp, max_nbcalc)
228    ! Sigma_eph(omega=eKS, kT, band) for fixed (kcalc, spin).
229 
230   complex(dpc),allocatable :: dvals_de0ks(:,:)
231    ! dvals_de0ks(ntemp, max_nbcalc) for fixed (kcalc, spin)
232    ! d Sigma_eph(omega, kT, band, kcalc, spin) / d omega (omega=eKS)
233 
234   real(dp),allocatable :: dw_vals(:,:)
235   !  dw_vals(ntemp, max_nbcalc) for fixed (kcalc, spin)
236   !  Debye-Waller term (static).
237 
238   !real(dp),allocatable :: qpadb_enes(:,:)
239   ! qp_adbenes(ntemp, max_nbcalc)
240   ! (Real) QP energies computed with the adiabatic formalism.
241 
242   !complex(dpc),allocatable :: qp_enes(:,:)
243   ! qp_enes(ntemp, max_nbcalc)
244   ! (Complex) QP energies computed with the non-adiabatic formalism.
245 
246   complex(dpc),allocatable :: vals_wr(:,:,:)
247    ! vals_wr(nwr, ntemp, max_nbcalc)
248    ! Sigma_eph(omega, kT, band) for given (k, spin).
249    ! enk_KS corresponds to nwr/2 + 1.
250    ! This array depends on (ikcalc, spin)
251 
252   integer :: gfw_nomega
253 
254   real(dp),allocatable :: gfw_mesh(:)
255    ! gfw_mesh(gfw_nomega)
256 
257   complex(dpc),allocatable :: gfw_vals(:,:,:,:)
258    !gfw_vals(gfw_nomega, ntemp, 2, max_nbcalc)
259 
260   logical :: has_nuq_terms
261 
262   complex(dpc),allocatable :: vals_nuq(:,:,:,:,:)
263    ! vals(ntemp, max_nbcalc, natom*3, nqbz, 2)
264    ! (nu, q)-resolved self-energy for given (k, spin)
265    ! First slice: Fan term at omega=e0_KS
266    ! Second slice: DW term
267 
268   type(degtab_t),allocatable :: degtab(:,:)
269    ! (nkcalc, nsppol)
270    ! Table used to average QP results in the degenerate subspace if symsigma == 1
271 
272  end type sigmaph_t