TABLE OF CONTENTS


ABINIT/m_phgamma [ Modules ]

[ Top ] [ Modules ]

NAME

FUNCTION

  Tools for the computation of phonon linewidths

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

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 module m_phgamma
25 
26  use defs_basis
27  use defs_abitypes
28  use m_abicore
29  use m_xmpi
30  use m_errors
31  use m_kptrank
32  use m_tetrahedron
33  use m_ifc
34  use m_ebands
35  use m_fstab
36  use iso_c_binding
37  use m_nctk
38 #ifdef HAVE_NETCDF
39  use netcdf
40 #endif
41 
42  use m_time,           only : cwtime
43  use m_fstrings,       only : toupper, itoa, sjoin, ktoa, ltoa, strcat
44  use m_numeric_tools,  only : arth, wrap2_pmhalf, simpson_int, simpson, bisect, mkherm, get_diag
45  use m_io_tools,       only : open_file
46  use m_symtk,          only : littlegroup_q
47  use m_geometry,       only : normv
48  use m_special_funcs,  only : dirac_delta
49  use m_fftcore,        only : ngfft_seq
50  use m_fft_mesh,       only : rotate_fft_mesh
51  !use m_cgtools,        only : set_istwfk
52  use m_cgtk,           only : cgtk_rotate
53  use m_kg,             only : getph
54  use m_dynmat,         only : d2sym3, symdyma, ftgam_init, ftgam, asrif9
55  use defs_datatypes,   only : ebands_t
56  use m_crystal,        only : crystal_t
57  use m_crystal_io,     only : crystal_ncwrite
58  use m_bz_mesh,        only : isamek, kpath_t, kpath_new, kpath_free, kpath_print
59  use m_special_funcs,  only : fermi_dirac
60  use m_kpts,           only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt, listkk
61  use defs_elphon,      only : gam_mult_displ, complete_gamma !, complete_gamma_tr
62  use m_getgh1c,          only : getgh1c, rf_transgrid_and_pack, getgh1c_setup
63  use m_fourier_interpol, only : transgrid
64 
65  implicit none
66 
67  private

m_phgamma/a2fw_ee_write [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_ee_write

FUNCTION

  Write alpha^2F(e,e',w) to an external file in text form

INPUTS

  a2f<a2fw_t>=Container storing the Eliashberg functions.
  basename=Filename for output.

OUTPUT

  Output is written to file. This routine should be called by one MPI proc.

PARENTS

      m_phgamma

CHILDREN

SOURCE

3173 subroutine a2fw_ee_write(a2f,basename)
3174 
3175 
3176 !This section has been created automatically by the script Abilint (TD).
3177 !Do not modify the following lines by hand.
3178 #undef ABI_FUNC
3179 #define ABI_FUNC 'a2fw_ee_write'
3180 !End of the abilint section
3181 
3182  implicit none
3183 
3184 !Arguments ------------------------------------
3185 !scalars
3186  character(len=*),intent(in) :: basename
3187  type(a2fw_t),intent(in) :: a2f
3188 
3189 !Local variables -------------------------
3190 !scalars
3191  integer :: iw,spin,unt,ii,mu
3192  integer :: iene, jene
3193  real(dp) :: ene1,ene2
3194  character(len=500) :: msg
3195  character(len=fnlen) :: path
3196 
3197 ! *********************************************************************
3198 
3199  ! Write spin-resolved a2F(e,e',w)
3200  path = strcat(basename, "_A2FEEW")
3201  if (open_file(path,msg,newunit=unt,form="formatted",status="unknown") /= 0) then
3202    MSG_ERROR(msg)
3203  end if
3204 
3205  call write_a2fw_header()
3206 
3207  write(unt,'(a)')"# en2, en1, Frequency, a2F_tot(w)"
3208  do iw=1,a2f%nomega
3209    do iene=1,a2f%nene
3210      ene1 = a2f%enemin + (iene-1)*a2f%deltaene
3211      do jene=1,a2f%nene
3212        ene2 = a2f%enemin + (jene-1)*a2f%deltaene
3213        write(unt,'(3E20.10,2x,E20.10)') ene2, ene1, a2f%omega(iw), sum(a2f%vals_ee(jene,iene,iw,:))   ! TOT
3214      end do
3215    end do
3216  end do
3217 
3218  close(unt)
3219 
3220  ! Write spin-resolved a2F(ef, ef, w)
3221  path = strcat(basename, "_A2FW_reference")
3222  if (open_file(path,msg,newunit=unt,form="formatted",status="unknown") /= 0) then
3223    MSG_ERROR(msg)
3224  end if
3225 
3226  call write_a2fw_header()
3227 
3228  write(unt,'(a)')"# Frequency, a2F_tot(ef,ef,w) for comparison with normal a2F(w)"
3229  iene = int(a2f%nene/2)
3230  jene = int(a2f%nene/2)
3231  do iw=1,a2f%nomega
3232    write(unt,'(E20.10,2x,E20.10)') a2f%omega(iw), sum(a2f%vals_ee(jene,iene,iw,:))   ! TOT
3233  end do
3234 
3235  close(unt)
3236 
3237 contains
3238 
3239 subroutine write_a2fw_header()
3240 
3241  ! Output the header.
3242 
3243 !This section has been created automatically by the script Abilint (TD).
3244 !Do not modify the following lines by hand.
3245 #undef ABI_FUNC
3246 #define ABI_FUNC 'write_a2fw_header'
3247 !End of the abilint section
3248 
3249  write(unt,'(a)')              '#'
3250  write(unt,'(a)')              '# ABINIT package: a2F(e,eprime,w) file'
3251  write(unt,'(a)')              '#'
3252  write(unt,'(a)')              '# a2F(e,eprime,w) function integrated over the FS. omega in a.u.'
3253  write(unt,'(2a)')             '# ngqpt: ',trim(ltoa(a2f%ngqpt))
3254  write(unt,'(a,i0,2(a,e16.6))')'# number of energies: ',a2f%nene," from enemin: ",a2f%enemin,&
3255 &    ' Ha with step ', a2f%deltaene
3256  write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f%nomega," between omega_min: ",a2f%omega_min,&
3257                                ' Ha and omega_max: ',a2f%omega_max,' Ha with step:',a2f%wstep
3258  write(unt,'(a,e16.6)')         '#  the smearing width for gaussians is ',a2f%smear
3259  write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f%n0)
3260  do spin=1,a2f%nsppol
3261    write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f%n0(spin)
3262  end do
3263  do ii=1,2
3264    write(unt,'(a)')     "# "
3265  end do
3266 
3267 end subroutine write_a2fw_header
3268 
3269 end subroutine a2fw_ee_write

m_phgamma/a2fw_free [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_free

FUNCTION

  Free the memory allocated in a2f

SIDE EFFECTS

  a2f<a2fw_t>=Structure storing the Eliashberg function a2F.

OUTPUT

PARENTS

      m_phgamma

CHILDREN

SOURCE

2181 subroutine a2fw_free(a2f)
2182 
2183 
2184 !This section has been created automatically by the script Abilint (TD).
2185 !Do not modify the following lines by hand.
2186 #undef ABI_FUNC
2187 #define ABI_FUNC 'a2fw_free'
2188 !End of the abilint section
2189 
2190  implicit none
2191 
2192 !Arguments ------------------------------------
2193  type(a2fw_t),intent(inout) :: a2f
2194 
2195 ! *********************************************************************
2196 
2197  ! @a2fw_t
2198 
2199  ! integer
2200  if (allocated(a2f%qshift)) then
2201    ABI_FREE(a2f%qshift)
2202  end if
2203 
2204  ! real
2205  if (allocated(a2f%n0)) then
2206    ABI_FREE(a2f%n0)
2207  end if
2208  if (allocated(a2f%omega)) then
2209    ABI_FREE(a2f%omega)
2210  end if
2211  if (allocated(a2f%vals)) then
2212    ABI_FREE(a2f%vals)
2213  end if
2214  if (allocated(a2f%vals_ee)) then
2215    ABI_FREE(a2f%vals_ee)
2216  end if
2217  if (allocated(a2f%lambdaw)) then
2218    ABI_FREE(a2f%lambdaw)
2219  end if
2220 
2221 end subroutine a2fw_free

m_phgamma/a2fw_init [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_init

FUNCTION

  Calculates the FS averaged alpha^2F(w) function

INPUTS

  cryst<crystal_t>=Info on the unit cell.
  ifc<ifc_type>=Interatomic force constants.
  gams<phgamma_t>=Structure storing the phonon linewidths.
  wstep=Step for linear frequency mesh in Ha.
  wminmax(2)=Minimum and maximum phonon frequency. Used to construct the linear mesh for A2F(w).
  intmeth=Integration method: 1 for gaussian, 2 for tetrahedra
  smear=Gaussian broadening used to approximate the Dirac delta.
  ngqpt(3)=Divisions of the Q-mesh used for interpolating the phonon linewidths (see also nqshift and qshift).
  nqshift=Number of shifts used to generated the Q-mesh.
  qshift(3,nqshift)=The shifts.
  comm=MPI communicator
  [qintp]=If set to False, ngqgpt, nqshift, qshift and qptop  are ignored and
     A2F(w) is computed from the IBZ values stored in gams. Default: True i.e use Fourier interpolation.
  [qptopt]=Controls the generation of the q-points. If not specified, the routine takes fully into account
    the symmetries of the system to generate the q points in the IBZone i.e. qptopt=1
    Other values of qptopt can be used for debugging purpose.

OUTPUT

  a2f<a2fw_t>=Structure storing the Eliashberg function a2F(w).

PARENTS

      m_phgamma

CHILDREN

SOURCE

2261 subroutine a2fw_init(a2f,gams,cryst,ifc,intmeth,wstep,wminmax,smear,ngqpt,nqshift,qshift,comm,&
2262   qintp,qptopt) ! optional
2263 
2264 
2265 !This section has been created automatically by the script Abilint (TD).
2266 !Do not modify the following lines by hand.
2267 #undef ABI_FUNC
2268 #define ABI_FUNC 'a2fw_init'
2269 !End of the abilint section
2270 
2271  implicit none
2272 
2273 !Arguments ------------------------------------
2274 !scalars
2275  integer,intent(in) :: intmeth,nqshift,comm
2276  integer,intent(in),optional :: qptopt
2277  real(dp),intent(in) :: wstep,smear
2278  logical,optional,intent(in) :: qintp
2279  type(phgamma_t),intent(inout) :: gams
2280  type(ifc_type),intent(in) :: ifc
2281  type(crystal_t),intent(in) :: cryst
2282  type(a2fw_t),target,intent(out) :: a2f
2283 !arrays
2284  integer,intent(in) :: ngqpt(3)
2285  real(dp),intent(in) :: wminmax(2),qshift(3,nqshift)
2286 
2287 !Local variables -------------------------
2288 !scalars
2289  integer,parameter :: bcorr0=0,master=0
2290  integer :: my_qptopt,iq_ibz,nqibz,ount,ii,my_rank,nproc,cnt
2291  integer :: mu,iw,natom3,nsppol,spin,ierr,nomega,nqbz
2292  integer :: iene, jene
2293  integer :: itemp, ntemp, jene_jump
2294 
2295  real(dp) :: cpu,wall,gflops
2296  real(dp) :: lambda_iso,omega,omega_log,xx,omega_min,omega_max,ww,mustar,tc_macmill
2297  real(dp) :: temp_el, min_temp, delta_temp, chempot, ene1, ene2, G0
2298  logical :: do_qintp
2299  character(len=500) :: msg
2300  type(t_tetrahedron) :: tetra
2301 !arrays
2302  integer :: qptrlatt(3,3),new_qptrlatt(3,3)
2303  real(dp),allocatable :: my_qshift(:,:)
2304  real(dp) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3)
2305  real(dp) :: invphfrq(gams%natom3)
2306  real(dp) :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3,gams%nsppol)
2307  real(dp) :: tmp_gam1(2,gams%natom3,gams%natom3)
2308  real(dp) :: tmp_gam2(2,gams%natom3,gams%natom3)
2309  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
2310  real(dp),allocatable :: tmp_a2f(:)
2311  real(dp), ABI_CONTIGUOUS pointer :: a2f_1d(:)
2312  real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:)
2313  real(dp),allocatable :: a2f_1mom(:),a2flogmom(:),a2flogmom_int(:),wdt(:,:)
2314  real(dp),allocatable :: lambda_tetra(:,:,:),phfreq_tetra(:,:)
2315  real(dp),allocatable :: tmp_gaussian(:,:)
2316  real(dp), allocatable :: a2feew_partial(:), a2feew_partial_int(:), a2feew_w(:), a2feew_w_int(:)
2317 
2318 ! *********************************************************************
2319 
2320  !@a2fw_t
2321  my_qptopt = 1; if (present(qptopt)) my_qptopt = qptopt
2322  do_qintp = .True.; if (present(qintp)) do_qintp = qintp
2323 
2324  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2325  nsppol = gams%nsppol; natom3 = gams%natom3
2326 
2327  call cwtime(cpu,wall,gflops,"start")
2328 
2329  if (do_qintp) then
2330    ! Generate the q-mesh by finding the IBZ and the corresponding weights.
2331    qptrlatt = 0
2332    do ii=1,3
2333      qptrlatt(ii,ii) = ngqpt(ii)
2334    end do
2335 
2336    call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, nqibz, qibz, wtq, nqbz, qbz, &
2337       new_kptrlatt=new_qptrlatt, new_shiftk=my_qshift)
2338    ABI_FREE(qbz)
2339 
2340    ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ.
2341    a2f%ngqpt = ngqpt; a2f%nqshift = size(my_qshift, dim=2)
2342    ABI_MALLOC(a2f%qshift, (3, a2f%nqshift))
2343    a2f%qshift = my_qshift
2344    ABI_FREE(my_qshift)
2345 
2346  else
2347    ! No interpolation. Use q-mesh parameters from gams%
2348    a2f%ngqpt = gams%ngqpt; a2f%nqshift = 1
2349    nqibz = gams%nqibz
2350    ABI_MALLOC(qibz,(3,nqibz))
2351    ABI_MALLOC(wtq,(nqibz))
2352    ABI_MALLOC(a2f%qshift,(3,a2f%nqshift))
2353    qibz = gams%qibz; wtq = gams%wtq
2354    a2f%qshift = zero ! Note: assuming q-mesh centered on Gamma.
2355  end if
2356 
2357  call cwtime(cpu,wall,gflops,"stop")
2358  write(msg,'(2(a,f8.2))')"a2fw_init, q-setup: ",cpu,", wall: ",wall
2359  call wrtout(std_out,msg,"COLL",do_flush=.True.)
2360 
2361  ! Define Min and max frequency for the mesh (enlarge it a bit)
2362  omega_min = wminmax(1); omega_max = wminmax(2)
2363  omega_min = omega_min - 0.1*abs(omega_min)
2364  if (omega_min >= zero) omega_min = one/Ha_meV
2365  omega_max = omega_max + 0.1*abs(omega_max)
2366 
2367  a2f%nsppol = nsppol; a2f%natom3 = gams%natom3; a2f%smear = smear
2368  a2f%omega_min = omega_min; a2f%omega_max = omega_max
2369  nomega = int((omega_max - omega_min) / wstep); a2f%nomega = nomega; a2f%wstep = wstep
2370  a2f%nene = gams%nene
2371  a2f%enemin = gams%enemin
2372  a2f%deltaene = gams%deltaene
2373 
2374  ABI_MALLOC(a2f%n0,(nsppol))
2375  a2f%n0=gams%n0
2376  ! Build linear mesh.
2377  ABI_MALLOC(a2f%omega, (nomega))
2378  a2f%omega = arth(omega_min,wstep,nomega)
2379  ABI_CALLOC(a2f%vals, (nomega,0:natom3, nsppol))
2380  ABI_CALLOC(a2f%lambdaw, (nomega,0:natom3, nsppol))
2381 
2382 #ifdef DEV_MJV
2383  ABI_CALLOC(a2f%vals_ee, (gams%nene,gams%nene,nomega,nsppol))
2384 #endif
2385 
2386  ABI_MALLOC(tmp_a2f, (nomega))
2387 
2388  if (intmeth == 2) then
2389    call cwtime(cpu,wall,gflops,"start")
2390 
2391    ! Prepare tetrahedron integration.
2392    qptrlatt = 0
2393    do ii=1,3
2394      qptrlatt(ii,ii) = a2f%ngqpt(ii)
2395    end do
2396 
2397    tetra = tetra_from_kptrlatt(cryst, my_qptopt, qptrlatt, a2f%nqshift, a2f%qshift, nqibz, qibz, msg, ierr)
2398    if (ierr/=0) MSG_ERROR(msg)
2399 
2400    ABI_STAT_MALLOC(lambda_tetra, (nqibz,natom3,nsppol), ierr)
2401    ABI_CHECK(ierr==0, "oom in lambda_tetra, use gaussians for A2F")
2402    lambda_tetra = zero
2403 
2404    ABI_STAT_MALLOC(phfreq_tetra, (nqibz,natom3), ierr)
2405    ABI_CHECK(ierr==0, "oom in phfreq_tetra, use gaussians for A2F")
2406    phfreq_tetra = zero
2407    cnt = 0
2408    do iq_ibz = 1, nqibz
2409      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2410      ! interpolated phonon freqs
2411      call ifc_fourq(ifc,cryst,qibz(:,iq_ibz),phfrq,displ_cart)
2412      ! save for tetrahedron interpolation
2413      phfreq_tetra(iq_ibz,:) = phfrq(:)
2414    end do
2415    call xmpi_sum(phfreq_tetra, comm, ierr)
2416 
2417    call cwtime(cpu,wall,gflops,"stop")
2418    write(msg,'(2(a,f8.2))')"a2fw_init%tetra, cpu",cpu,", wall: ",wall
2419    call wrtout(std_out,msg,"COLL",do_flush=.True.)
2420  end if
2421 
2422  call cwtime(cpu,wall,gflops,"start")
2423 
2424 #ifdef DEV_MJV
2425  open (unit=900, file="a2fvals_ee.dat")
2426  write (900,*) '# do_qintp ', do_qintp
2427 #endif
2428  ! Loop over spins and qpoints in the IBZ
2429  cnt = 0
2430  do spin=1,nsppol
2431    do iq_ibz=1,nqibz
2432      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2433 
2434      ! Interpolate or evaluate gamma directly.
2435 #ifdef DEV_MJV
2436      if (do_qintp) then
2437        call phgamma_interp(gams,cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee)
2438      else
2439        call phgamma_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee)
2440      end if
2441 #else
2442      if (do_qintp) then
2443        call phgamma_interp(gams,cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_ph,lambda_ph,displ_cart)
2444      else
2445        call phgamma_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph,displ_cart)
2446      end if
2447 #endif
2448 
2449      select case (intmeth)
2450      case (1)
2451        ! Gaussian: Add all contributions from the phonon modes at this qpoint to a2f
2452        ! (note that unstable modes are included).
2453        ABI_ALLOCATE(tmp_gaussian, (nomega, natom3))
2454        do mu=1,natom3
2455          tmp_a2f = zero
2456          do iw=1,nomega
2457            xx = a2f%omega(iw) - phfrq(mu)
2458            tmp_gaussian(iw,mu) = dirac_delta(xx, smear)
2459            tmp_a2f(iw) = tmp_a2f(iw) + tmp_gaussian(iw,mu) * lambda_ph(mu) * abs(phfrq(mu))
2460          end do
2461          a2f%vals(:,mu,spin) = a2f%vals(:,mu,spin) + tmp_a2f * wtq(iq_ibz)
2462        end do
2463 
2464 #ifdef DEV_MJV
2465        ! reset phfrq for low freq modes to
2466        invphfrq = zero
2467        do mu=1,natom3
2468          if (abs(phfrq(mu)) > EPH_WTOL) then
2469            invphfrq(mu) = one / abs(phfrq(mu))
2470          end if
2471        end do
2472 
2473        do iene= 1, gams%nene
2474          do jene= 1, gams%nene
2475            tmp_a2f = zero
2476 ! TODO: following block is just a GEMM
2477            do mu=1,natom3
2478              do iw=1,nomega
2479                tmp_a2f(iw) = tmp_a2f(iw) + tmp_gaussian(iw,mu) * gamma_ph_ee(jene,iene,mu,spin) * invphfrq(mu)
2480              end do
2481            end do
2482 
2483            a2f%vals_ee(jene, iene, :, spin) = a2f%vals_ee(jene, iene, :, spin) + tmp_a2f(:) * wtq(iq_ibz)
2484 if (iene == gams%nene/2 .and. jene == gams%nene/2) then
2485   write (900, '(a,E20.10,2x,2x,I6,3E20.10)') '#', wtq(iq_ibz), iq_ibz, invphfrq(1:3)
2486   do iw=1,nomega
2487     write (900, '(i6,2x,E20.10,2x,3E20.10,2x,3E20.10)') iw, a2f%vals_ee(jene, iene, iw, spin), &
2488 &        gamma_ph_ee(jene,iene,:,spin), tmp_gaussian(iw,1:3)
2489   end do
2490   write (900,*)
2491 end if
2492          end do
2493        end do
2494 #endif
2495        ABI_DEALLOCATE(tmp_gaussian)
2496 
2497 
2498      case (2)
2499        ! Tetra: store data.
2500        do mu=1,natom3
2501          lambda_tetra(iq_ibz, mu, spin) = lambda_ph(mu) * abs(phfrq(mu))
2502        end do
2503      end select
2504 
2505    end do ! iq_ibz
2506  end do ! spin
2507 #ifdef DEV_MJV
2508  close(900)
2509 #endif
2510 
2511  if (intmeth == 2) then
2512    ! Collect results on each node.
2513    call xmpi_sum(lambda_tetra, comm, ierr)
2514 
2515    ! workspace for tetra.
2516    ABI_MALLOC(wdt, (nomega, 2))
2517 
2518    ! For each mode get its contribution
2519    cnt = 0
2520    do spin=1,nsppol
2521      do mu=1,natom3
2522        do iq_ibz=1,nqibz
2523          cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! mpi-parallelism
2524 
2525          call tetra_get_onewk(tetra, iq_ibz, bcorr0, nomega, nqibz, phfreq_tetra(:,mu), &
2526            omega_min, omega_max, one, wdt)
2527 
2528          ! Accumulate (Integral of a2F is computed afterwards)
2529          a2f%vals(:,mu,spin) = a2f%vals(:,mu,spin) + wdt(:,1) * lambda_tetra(iq_ibz,mu,spin)
2530          !a2f%lambdaw(:,mu,spin) = a2f%lambdaw(:,mu,spin) + wdt(:,2) * lambda_tetra(iq_ibz, mu, spin)
2531       end do
2532      end do
2533    end do
2534 
2535    ! Free memory allocated for tetra.
2536    ABI_FREE(wdt)
2537    ABI_FREE(lambda_tetra)
2538    ABI_FREE(phfreq_tetra)
2539    call destroy_tetra(tetra)
2540  end if
2541 
2542  ! Collect final results on each node
2543  call xmpi_sum(a2f%vals, comm, ierr)
2544  do spin=1,nsppol
2545    a2f%vals(:,0,spin) = sum(a2f%vals(:,1:natom3,spin), dim=2)
2546 #ifdef DEV_MJV
2547    a2f%vals_ee(:,:,:,spin) = a2f%vals_ee(:,:,:,spin) / (two * pi * gams%n0(spin))
2548 #endif
2549 
2550    ! previously would divide by g(eF, spin)
2551    !a2f%vals(:,:,spin) = a2f%vals(:,:,spin) / (two_pi*a2f%n0(spin))
2552  end do
2553 
2554  !to avoid numerical noise uses a smoothing function
2555  ! TODO: Move smooth to m_numeric_tools and add ndat dimension.
2556  !do spin=1,nsppol
2557  !  do mu=0,natom3
2558  !    call smooth(a2f%vals(:,mu,spin),a2f%nomega,2)
2559  !  end do
2560  !end do
2561 
2562  call cwtime(cpu,wall,gflops,"stop")
2563  write(msg,'(2(a,f8.2))')"a2fw_init, a2f_eval: ",cpu,", wall: ",wall
2564  call wrtout(std_out,msg,"COLL",do_flush=.True.)
2565 
2566  ABI_MALLOC(a2f_1mom,(nomega))
2567  ABI_MALLOC(a2flogmom,(nomega))
2568  ABI_MALLOC(a2flogmom_int,(nomega))
2569 
2570  !call a2fw_print_info()
2571 
2572  ! Compute lambda(w) (resolved in spin and phonon mode).
2573  do spin=1,nsppol
2574    do mu=0,natom3
2575 
2576      do iw=1,nomega
2577        ww = a2f%omega(iw)
2578        if (abs(ww) > EPH_WTOL) then
2579           a2f_1mom(iw) = a2f%vals(iw,mu,spin) / abs(ww)
2580        else
2581           a2f_1mom(iw) = zero
2582        end if
2583        ! TODO: Strange that the first value of int(f) is not set to zero!
2584        ! FIXME: The frequency integration overestimates lambda if there are negative phonon frequencies!
2585      end do
2586 
2587      call simpson_int(nomega, wstep, a2f_1mom, a2f%lambdaw(:,mu,spin))
2588    end do
2589  end do
2590 
2591  ! print log moments of the alpha 2 F functions
2592  do spin=1,nsppol
2593    a2f_1d => a2f%vals(:,0,spin)
2594 
2595    lambda_iso = a2fw_moment(a2f,0,spin)
2596 
2597    ! Get log moment of alpha^2F.
2598    a2flogmom = zero
2599    do iw=1,nomega
2600      omega = a2f%omega(iw)
2601      if (abs(omega) > EPH_WTOL) then
2602        a2flogmom(iw) = (two/lambda_iso) * a2f_1d(iw)*log(abs(omega))/abs(omega)
2603      end if
2604    end do
2605    call simpson_int(nomega,wstep,a2flogmom,a2flogmom_int)
2606    omega_log = exp(a2flogmom_int(nomega))
2607 
2608    mustar = 0.12
2609    tc_macmill = omega_log/1.2_dp * exp((-1.04_dp*(one+lambda_iso)) / (lambda_iso-mustar*(one+0.62_dp*lambda_iso)))
2610 
2611    if (my_rank == master) then
2612      ount = std_out
2613      if (nsppol > 1) then
2614        write(msg,'(3a)') ch10,&
2615          'Warning: some of the following quantities should be integrated over spin', ch10
2616        call wrtout(ount,msg,'COLL')
2617      end if
2618 
2619      write(ount,'(a)')' Superconductivity: isotropic evaluation of parameters from electron-phonon coupling.'
2620      write(ount,'(a,es16.6)')' isotropic lambda = ',lambda_iso
2621      write(ount,'(a,es16.6,a,es16.6,a)' )' omegalog  = ',omega_log,' (Ha) ', omega_log*Ha_K, ' (Kelvin) '
2622      write(ount,'(a,es16.6,a,es16.6,a)')' MacMillan Tc = ',tc_macmill,' (Ha) ', tc_macmill*Ha_K, ' (Kelvin) '
2623      write(ount,"(a)")'    positive moments of alpha2F:'
2624      write(ount,'(a,es16.6)' )' lambda <omega^2> = ',a2fw_moment(a2f,2,spin)
2625      write(ount,'(a,es16.6)' )' lambda <omega^3> = ',a2fw_moment(a2f,3,spin)
2626      write(ount,'(a,es16.6)' )' lambda <omega^4> = ',a2fw_moment(a2f,4,spin)
2627      write(ount,'(a,es16.6)' )' lambda <omega^5> = ',a2fw_moment(a2f,5,spin)
2628    end if
2629  end do
2630 
2631 #ifdef DEV_MJV
2632  ! calculate the temperature dependence of the a2f(e,e',w) integrals (G_0(T_e) in PRL 110 016405 (2013) [[cite:Arnaud2013]]) 
2633  if (my_rank == master) then
2634    ntemp = 100
2635    min_temp = zero
2636    delta_temp = 40._dp ! Kelvin
2637    ABI_ALLOCATE (a2feew_partial, (a2f%nene))
2638    ABI_ALLOCATE (a2feew_partial_int, (a2f%nene))
2639    ABI_ALLOCATE (a2feew_w, (nomega))
2640    ABI_ALLOCATE (a2feew_w_int, (nomega))
2641 print *, "temp_el, G_0(T_e) in W/m^3/K, spin"
2642    do spin=1,nsppol
2643      do itemp = 1, ntemp
2644        temp_el = min_temp + (itemp-1)*delta_temp
2645 
2646        ! TODO: need to evolve the chemical potential with T, but I do not have access to the full DOS here!!
2647        ! possible fix using local information on DOS variation near E_F...
2648        chempot = a2f%enemin + half * a2f%nene * a2f%deltaene
2649 
2650        do iw=1,nomega
2651          omega = a2f%omega(iw)
2652          jene_jump = nint(omega/a2f%deltaene)
2653          a2feew_partial = zero
2654          do iene=1,a2f%nene
2655            ene1 = a2f%enemin + (iene-1)*a2f%deltaene
2656            ene2 = ene1 + omega
2657            a2feew_partial(iene) = a2f%vals_ee(min(a2f%nene,iene+jene_jump), iene, iw, spin) * &
2658 &             (fermi_dirac(ene1, chempot, temp_el/Ha_K) - fermi_dirac(ene2, chempot, temp_el/Ha_K))
2659          end do
2660          call simpson_int(a2f%nene, a2f%deltaene, a2feew_partial, a2feew_partial_int)
2661          a2feew_w(iw) = a2feew_partial_int(a2f%nene)
2662        end do
2663        call simpson_int(nomega,wstep,a2feew_w,a2feew_w_int)
2664        G0 = a2feew_w_int(nomega) * two_pi * a2f%n0(spin) / cryst%ucvol
2665        ! conversion factor for G0 to SI units =  Ha_J / Time_Sec / (Bohr_meter)**3 ~ 1.2163049915755545e+30
2666        print *, temp_el, G0  * kb_HaK / Time_Sec / (Bohr_meter)**3, spin !* Ha_J???
2667      end do
2668    end do
2669    ABI_DEALLOCATE (a2feew_partial)
2670    ABI_DEALLOCATE (a2feew_partial_int)
2671    ABI_DEALLOCATE (a2feew_w)
2672    ABI_DEALLOCATE (a2feew_w_int)
2673  end if
2674 #endif
2675 
2676  ABI_FREE(tmp_a2f)
2677  ABI_FREE(a2f_1mom)
2678  ABI_FREE(a2flogmom)
2679  ABI_FREE(a2flogmom_int)
2680  ABI_FREE(qibz)
2681  ABI_FREE(wtq)
2682 
2683 end subroutine a2fw_init

m_phgamma/a2fw_lambda_wij [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_lambda_wij

FUNCTION

  Compute \int dw [w x a2F(w)] (w_i-w_j)**2 + w**2)

INPUTS

  a2f<a2fw_t>=Structure storing the Eliashberg function.
  w1,wj=Matsubara frequencies (real) in Hartree
  spin=Spin index

PARENTS

CHILDREN

SOURCE

2934 real(dp) function a2fw_lambda_wij(a2f,wi,wj,spin) result(res)
2935 
2936 
2937 !This section has been created automatically by the script Abilint (TD).
2938 !Do not modify the following lines by hand.
2939 #undef ABI_FUNC
2940 #define ABI_FUNC 'a2fw_lambda_wij'
2941 !End of the abilint section
2942 
2943  implicit none
2944 
2945 !Arguments ------------------------------------
2946 !scalars
2947  integer,intent(in) :: spin
2948  real(dp),intent(in) :: wi,wj
2949  type(a2fw_t),intent(in) :: a2f
2950 
2951 !Local variables -------------------------
2952 !scalars
2953  integer :: iw
2954  real(dp) :: wd2,vv,omg
2955  logical :: iszero
2956 !arrays
2957  real(dp) :: values(a2f%nomega)
2958 
2959 ! *********************************************************************
2960 
2961  wd2 = (wi - wj)** 2
2962  iszero = (abs(wi - wj) < EPH_WTOL)
2963 
2964  do iw=1,a2f%nomega
2965    omg = a2f%omega(iw)
2966    vv = a2f%vals(iw,0,spin)
2967    if (abs(a2f%omega(iw)) > EPH_WTOL) then
2968      values(iw) = vv * omg / (wd2 + omg**2)
2969    else
2970      if (iszero) then ! TODO
2971        values(iw) = zero
2972      else
2973        values(iw) = vv * omg / (wd2 + omg**2)
2974      end if
2975    end if
2976  end do
2977 
2978  res = simpson(a2f%wstep, values)
2979  if (res < zero) res = zero
2980 
2981 end function a2fw_lambda_wij

m_phgamma/a2fw_logmoment [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_logmoment

FUNCTION

INPUTS

  a2f<a2fw_t>=Structure storing the Eliashberg function.
  nn
  spin=Spin index

OUTPUT

  a2fw_logmoment =

PARENTS

CHILDREN

SOURCE

2871 real(dp) function a2fw_logmoment(a2f,spin) result(res)
2872 
2873 
2874 !This section has been created automatically by the script Abilint (TD).
2875 !Do not modify the following lines by hand.
2876 #undef ABI_FUNC
2877 #define ABI_FUNC 'a2fw_logmoment'
2878 !End of the abilint section
2879 
2880  implicit none
2881 
2882 !Arguments ------------------------------------
2883 !scalars
2884  integer,intent(in) :: spin
2885  type(a2fw_t),intent(in) :: a2f
2886 
2887 !Local variables -------------------------
2888 !scalars
2889  integer :: iw
2890  real(dp) :: omg,lambda_iso
2891 !arrays
2892  real(dp) :: a2flogmom(a2f%nomega) !,a2flogmom_int(a2f%nomega)
2893 
2894 ! *********************************************************************
2895 
2896  ! Get log moment of alpha^2F.
2897  a2flogmom = zero
2898  lambda_iso = a2fw_moment(a2f,0,spin)
2899 
2900  do iw=1,a2f%nomega
2901    omg = a2f%omega(iw)
2902    if (abs(omg) > EPH_WTOL) then
2903      !a2flogmom(iw) = (two/lambda_iso) * a2f%vals(iw,spin) * log(abs(omg))/abs(omg)
2904    end if
2905  end do
2906 
2907  !call simpson_int(nomega,a2f%wstep,a2flogmom,a2flogmom_int)
2908  !res = exp(a2flogmom_int(nomega))
2909  res = zero
2910 
2911 end function a2fw_logmoment

m_phgamma/a2fw_moment [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_moment

FUNCTION

  Compute \int dw [a2F(w)/w] w^n
  From Allen PRL 59 1460 [[cite:Allen1987]] (See also [[cite:Grimvall1981]], Eq 6.72 page 175)

INPUTS

  a2f<a2fw_t>=Structure storing the Eliashberg function.
  nn=Value of n
  spin=The spin component

OUTPUT

  a2fw_moment = \int dw [a2F(w)/w] w^n
  [out_int(x)] = \int^{x} dw [a2F(w)/w] w^n

PARENTS

CHILDREN

SOURCE

2711 real(dp) function a2fw_moment(a2f,nn,spin,out_int)
2712 
2713 
2714 !This section has been created automatically by the script Abilint (TD).
2715 !Do not modify the following lines by hand.
2716 #undef ABI_FUNC
2717 #define ABI_FUNC 'a2fw_moment'
2718 !End of the abilint section
2719 
2720  implicit none
2721 
2722 !Arguments ------------------------------------
2723 !scalars
2724  integer,intent(in) :: spin,nn
2725  type(a2fw_t),intent(in) :: a2f
2726 !arrays
2727  real(dp),intent(out),optional :: out_int(a2f%nomega)
2728 
2729 !Local variables -------------------------
2730 !scalars
2731  integer :: iw !,nomega
2732  real(dp) :: omg,omg_nm1
2733 !arrays
2734  real(dp) :: ff(a2f%nomega),int_ff(a2f%nomega)
2735 
2736 ! *********************************************************************
2737 
2738  ! Construct the integrand function. [a2F(w)/w] w^n
2739  ff = zero; int_ff = zero
2740 
2741  if (nn-1 >= 0) then
2742    do iw=1,a2f%nomega
2743      omg = a2f%omega(iw)
2744      omg_nm1 = omg ** (nn-1)
2745      ff(iw) = a2f%vals(iw,0,spin) * omg_nm1
2746    end do
2747  else
2748    do iw=1,a2f%nomega
2749      omg = a2f%omega(iw)
2750      omg_nm1 = zero; if (abs(omg) > EPH_WTOL) omg_nm1 = omg**(nn-1)
2751      ff(iw) = a2f%vals(iw,0,spin) * omg_nm1
2752    end do
2753  end if
2754 
2755  ! Integration with simpson rule on a linear mesh.
2756  call simpson_int(a2f%nomega,a2f%wstep,ff,int_ff)
2757 
2758  a2fw_moment = int_ff(a2f%nomega)
2759  if (present(out_int)) out_int = int_ff
2760 
2761 end function a2fw_moment

m_phgamma/a2fw_solve_gap [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_solve_gap

FUNCTION

INPUTS

  a2f<a2fw_t>=Container storing the Eliashberg functions.
  ntemp=Number of temperatures
  temp_range = min and max temperatures
  wcut = frequency cutoff for Matsubara sums
  mustar= mustar parameter
  nstep=Max number of SCF steps
  reltol = relative tolerance accepted for exit of main Eliashberg loop
  comm=MPI communicator

TODO

  Implement python version in AbiPy

OUTPUT

PARENTS

CHILDREN

SOURCE

3301 subroutine a2fw_solve_gap(a2f,cryst,ntemp,temp_range,wcut,mustar,nstep,reltol,prefix,comm)
3302 
3303 
3304 !This section has been created automatically by the script Abilint (TD).
3305 !Do not modify the following lines by hand.
3306 #undef ABI_FUNC
3307 #define ABI_FUNC 'a2fw_solve_gap'
3308 !End of the abilint section
3309 
3310  implicit none
3311 
3312 !Arguments ------------------------------------
3313 !scalars
3314  integer,intent(in) :: ntemp,nstep,comm
3315  real(dp),intent(in) :: wcut,mustar,reltol
3316  character(len=*),intent(in) :: prefix
3317  type(a2fw_t),intent(inout) :: a2f
3318  type(crystal_t),intent(in) :: cryst
3319 !arrays
3320  real(dp),intent(in) :: temp_range(2)
3321 
3322 !Local variables -------------------------
3323 !scalars
3324  integer,parameter :: master=0
3325  integer :: istep,ii,jj,it,spin,iwn,nwm,conv !iw,
3326  integer :: my_rank,nproc,ncid,ncerr
3327  real(dp) :: summ,kT,tstep,abs_delta,rel_delta,alpha,gap
3328  !character(len=500) :: msg
3329 !arrays
3330  real(dp),allocatable :: wmts(:),lambda_ij(:,:),tmesh(:)
3331  real(dp),allocatable :: din(:),dout(:),zin(:),zout(:)
3332 
3333 ! *********************************************************************
3334 
3335  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3336 
3337  where (a2f%vals < zero) a2f%vals = zero
3338 
3339  ! Build linear mesh of temperatures.
3340  ABI_CHECK(ntemp > 1, "ntemp cannot be 1")
3341  tstep = (temp_range(2) - temp_range(1)) / (ntemp - 1)
3342  ABI_MALLOC(tmesh, (ntemp))
3343  tmesh = arth(temp_range(1), tstep, ntemp)
3344 
3345  ! Matsubara frequencies: i w_n = i (2n+1) pi T
3346  kT = kb_HaK * tmesh(1)
3347  nwm = 0
3348  do
3349    if ((2*nwm + 1) * pi * kT > wcut) exit
3350    nwm = nwm + 1
3351  end do
3352  ABI_CHECK(nwm /= 0, "Empy list of Matsubara frequencies, increase wcut")
3353 
3354 #ifdef HAVE_NETCDF
3355  ! Open the netcdf file used to store the results of the calculation.
3356  if (my_rank == master) then
3357    NCF_CHECK(nctk_open_create(ncid, strcat(prefix, "_ELIASHBERG.nc"), xmpi_comm_self))
3358    NCF_CHECK(crystal_ncwrite(cryst, ncid))
3359 
3360    ! Define dimensions.
3361    ncerr = nctk_def_dims(ncid, [&
3362      nctkdim_t("maxnum_matsubara_frequencies", nwm), nctkdim_t("num_temperatures", ntemp) &
3363      ], defmode=.True.)
3364    NCF_CHECK(ncerr)
3365 
3366    ncerr = nctk_def_arrays(ncid, [&
3367      nctkarr_t('temperatures', "dp", "num_temperatures"),&
3368      nctkarr_t('delta_imag_axis', "dp", "maxnum_matsubara_frequencies, num_temperatures"),&
3369      nctkarr_t('zeta_imag_axis', "dp", "maxnum_matsubara_frequencies, num_temperatures"), &
3370      nctkarr_t('delta_real_axis', "dp", "maxnum_matsubara_frequencies, num_temperatures"),&
3371      nctkarr_t('zeta_real_axis', "dp", "maxnum_matsubara_frequencies, num_temperatures")  &
3372      ])
3373    NCF_CHECK(ncerr)
3374 
3375    NCF_CHECK(nctk_set_datamode(ncid))
3376    NCF_CHECK(nf90_put_var(ncid, vid("temperatures"), tmesh))
3377  end if
3378 #endif
3379 
3380  do it=1,ntemp
3381    kT = kb_HaK * tmesh(it)
3382 
3383    ! Matsubara frequencies: i w_n = i (2n+1) pi T
3384    nwm = 0
3385    do
3386      if ((2*nwm + 1) * pi * kT > wcut) exit
3387      nwm = nwm + 1
3388    end do
3389    ABI_CHECK(nwm /= 0, "Empy list of Matsubara frequencies, increase wcut")
3390    write(std_out,*)"Number of matsubara frequencies",nwm
3391 
3392    ABI_MALLOC(wmts, (nwm))
3393    do ii=1,nwm
3394      wmts(ii) = (2*(ii-1) + 1) * pi * kT
3395    end do
3396 
3397    ! Compute lambda(n-n') kernel (symmetric)
3398    ABI_MALLOC(lambda_ij, (nwm, nwm))
3399    do spin=1,a2f%nsppol
3400      do jj=1,nwm
3401        do ii=1,jj
3402           lambda_ij(ii, jj) = a2fw_lambda_wij(a2f,  wmts(ii), wmts(jj), spin)
3403           if (ii /= jj) lambda_ij(jj, ii) = lambda_ij(ii, jj)
3404        end do
3405      end do
3406    end do
3407    !lambda_ij = lambda_ij * half
3408 
3409    ABI_MALLOC(din, (nwm))
3410    ABI_MALLOC(dout, (nwm))
3411    ABI_MALLOC(zin, (nwm))
3412    ABI_MALLOC(zout, (nwm))
3413 
3414    ! Initalize din
3415    if (it == 1) then
3416      din = 3.4*10-4 * eV_Ha
3417    else
3418      din = gap
3419    end if
3420    !dout = din
3421 
3422    conv = 0
3423    do istep=1,nstep
3424 
3425      !where (din < zero) din = zero
3426 
3427      do iwn=1,nwm
3428        summ = zero
3429        do jj=1,nwm
3430          summ = summ + wmts(jj) * lambda_ij(jj, iwn) / sqrt(wmts(jj)**2 + din(jj)**2)
3431        end do
3432        zout(iwn) = one + (pi * kT / wmts(iwn)) * summ
3433      end do
3434 
3435      do iwn=1,nwm
3436        summ = zero
3437        do jj=1,nwm
3438          summ = summ + din(jj) * (lambda_ij(jj, iwn) - mustar) / sqrt(wmts(jj)**2 + din(jj)**2)
3439        end do
3440        dout(iwn) = (pi * kT / zout(iwn)) * summ
3441      end do
3442 
3443      ! Test for convergence
3444      abs_delta = sum(abs(din))
3445      rel_delta = sum(abs(dout - din))
3446      !write(std_out,*)"rel_delta / abs_delta", (rel_delta / abs_delta)
3447 
3448      if ((rel_delta / abs_delta) < reltol) then
3449        conv = conv + 1
3450      else
3451        conv = 0
3452      end if
3453      if (conv == 2) exit
3454 
3455      ! TODO: Broyden mixing
3456      alpha = 0.2
3457      !alpha = 0.4
3458      !alpha = one
3459      din = alpha * dout + (one-alpha) * din
3460      zin = zout
3461    end do
3462 
3463    gap = dout(1)
3464    if (conv == 2) then
3465      write(std_out,*)"Converged at iteration: ",istep
3466    else
3467      write(std_out,*)"Not converged",rel_delta / abs_delta
3468    end if
3469    write(std_out,*)"T=",tmesh(it)," [K], gap ",gap*Ha_eV," [eV]"
3470 
3471    ! Write data to netcd file.
3472 #ifdef HAVE_NETCDF
3473    if (my_rank == master) then
3474      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "zeta_imag_axis"), zin, start=[1,it]))
3475      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "delta_imag_axis"), din, start=[1,it]))
3476    end if
3477 #endif
3478 
3479    !if (it == 1) then
3480    !  do iwn=1,nwm
3481    !    write(999,*)wmts(iwn),dout(iwn),zout(iwn)
3482    !  end do
3483    !end if
3484 
3485    ABI_FREE(din)
3486    ABI_FREE(dout)
3487    ABI_FREE(zin)
3488    ABI_FREE(zout)
3489 
3490    ABI_FREE(wmts)
3491    ABI_FREE(lambda_ij)
3492  end do ! it
3493 
3494  ABI_FREE(tmesh)
3495 
3496 #ifdef HAVE_NETCDF
3497  if (my_rank == master) then
3498    NCF_CHECK(nf90_close(ncid))
3499  end if
3500 #endif
3501 
3502 
3503 contains
3504  integer function vid(vname)
3505 
3506 
3507 !This section has been created automatically by the script Abilint (TD).
3508 !Do not modify the following lines by hand.
3509 #undef ABI_FUNC
3510 #define ABI_FUNC 'vid'
3511 !End of the abilint section
3512 
3513    character(len=*),intent(in) :: vname
3514    vid = nctk_idname(ncid, vname)
3515  end function vid
3516 
3517 end subroutine a2fw_solve_gap

m_phgamma/a2fw_t [ Types ]

[ Top ] [ m_phgamma ] [ Types ]

NAME

 a2fw_t

FUNCTION

 Store the Eliashberg function a2F(w).

SOURCE

218  type,public :: a2fw_t
219 
220   integer :: nomega
221   ! Number of frequency points in a2f(w).
222 
223   integer :: nsppol
224   ! Number of independent spin polarizations.
225 
226   integer :: natom3
227   ! Number of phonon modes.
228 
229   integer :: nene
230   ! Number of chemical potential values used for inelastic integration
231 
232   real(dp) :: enemin
233   ! Minimal chemical potential value used for inelastic integration Copied from fstab
234 
235   real(dp) :: deltaene
236   ! Chemical potential increment for inelastic integration Copied from fstab
237   !   for simplicity could be made equal to phonon frequency step
238 
239   real(dp) :: omega_min,omega_max
240   ! min and Max frequency (Ha) in the linear mesh.
241 
242   real(dp) :: wstep
243   ! Step of the linear mesh
244 
245   real(dp) :: smear
246   ! Gaussian broadening used to approximated the Dirac distribution.
247 
248   integer :: nqshift
249   ! Number of shifts in the q-mesh
250 
251   integer :: ngqpt(3)
252   ! The q-mesh used for calculating vals(w).
253 
254   real(dp),allocatable :: qshift(:,:)
255   ! qshift(3,nqshift)
256   ! The shifts used to generate the q-mesh.
257 
258   real(dp),allocatable :: n0(:)
259   ! n0(nsppol)
260   ! Electronic DOS at the Fermi level.
261 
262   real(dp),allocatable :: omega(:)
263   ! omega(nomega)
264   ! Frequency mesh in Hartree (linear).
265 
266   real(dp),allocatable :: vals(:,:,:)
267   ! vals(nomega,0:natom3,nsppol)
268   ! Eliashberg function
269   !   vals(w,1:natom3,1:nsppol): a2f(w) decomposed per phonon branch and spin
270   !   vals(w,0,1:nsppol): a2f(w) summed over phonons modes, decomposed in spin
271 
272   real(dp),allocatable :: vals_ee(:,:,:,:)
273   ! vals_ee(nene,nene,nomega,nsppol)
274   ! Eliashberg function
275   !   vals(e,e',w,0,1:nsppol): a2f(e,e',w) summed over phonons modes, decomposed in spin
276 
277   real(dp),allocatable :: lambdaw(:,:,:)
278   ! lambda(nomega,0:natom3,nsppol)
279 
280  end type a2fw_t
281 
282  public :: a2fw_free            ! Free the memory allocated in the structure.
283  public :: a2fw_init            ! Calculates the FS averaged alpha^2F(w) function.
284  public :: a2fw_write           ! Write alpha^2F(w) to an external file in text/netcdf format
285  public :: a2fw_moment          ! Compute moments of alpha^2F(w)/w .
286  !public :: a2fw_solve_gap       ! DEPRECATED

m_phgamma/a2fw_tr_free [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_tr_free

FUNCTION

  Free the memory allocated in a2f

SIDE EFFECTS

  a2f<a2fw_tr_t>=Structure storing the Eliashberg function a2F.

OUTPUT

PARENTS

      m_phgamma

CHILDREN

SOURCE

3541 subroutine a2fw_tr_free(a2f_tr)
3542 
3543 
3544 !This section has been created automatically by the script Abilint (TD).
3545 !Do not modify the following lines by hand.
3546 #undef ABI_FUNC
3547 #define ABI_FUNC 'a2fw_tr_free'
3548 !End of the abilint section
3549 
3550  implicit none
3551 
3552 !Arguments ------------------------------------
3553  type(a2fw_tr_t),intent(inout) :: a2f_tr
3554 
3555 ! *********************************************************************
3556 
3557  ! @a2fw_tr_t
3558 
3559  ! integer
3560  if (allocated(a2f_tr%qshift)) then
3561    ABI_FREE(a2f_tr%qshift)
3562  end if
3563 
3564  ! real
3565  if (allocated(a2f_tr%n0)) then
3566    ABI_FREE(a2f_tr%n0)
3567  end if
3568  if (allocated(a2f_tr%omega)) then
3569    ABI_FREE(a2f_tr%omega)
3570  end if
3571  if (allocated(a2f_tr%vals_in)) then
3572    ABI_FREE(a2f_tr%vals_in)
3573  end if
3574  if (allocated(a2f_tr%vals_out)) then
3575    ABI_FREE(a2f_tr%vals_out)
3576  end if
3577  if (allocated(a2f_tr%vals_tr)) then
3578    ABI_FREE(a2f_tr%vals_tr)
3579  end if
3580  if (allocated(a2f_tr%vals_tr_gen)) then
3581    ABI_FREE(a2f_tr%vals_tr_gen)
3582  end if
3583  if (allocated(a2f_tr%lambdaw_tr)) then
3584    ABI_FREE(a2f_tr%lambdaw_tr)
3585  end if
3586 
3587 end subroutine a2fw_tr_free

m_phgamma/a2fw_tr_init [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_tr_init

FUNCTION

  Calculates the FS averaged alpha^2F_tr,in,out(w) functions

INPUTS

  cryst<crystal_t>=Info on the unit cell.
  ifc<ifc_type>=Interatomic force constants.
  gams<phgamma_t>=Structure storing the phonon linewidths.
  wstep=Step for linear frequency mesh in Ha.
  wminmax(2)=Minimum and maximum phonon frequency. Used to construct the linear mesh for A2F(w).
  intmeth=Integration method: 1 for gaussian, 2 for tetrahedra
  smear=Gaussian broadening used to approximate the Dirac delta.
  ngqpt(3)=Divisions of the Q-mesh used for interpolating the phonon linewidths (see also nqshift and qshift).
  nqshift=Number of shifts used to generated the Q-mesh.
  qshift(3,nqshift)=The shifts.
  comm=MPI communicator
  [qintp]=If set to False, ngqgpt, nqshift, qshift and qptop  are ignored and
     A2F(w) is computed from the IBZ values stored in gams. Default: True i.e use Fourier interpolation.
  [qptopt]=Controls the generation of the q-points. If not specified, the routine takes fully into account
    the symmetries of the system to generate the q points in the IBZone i.e. qptopt=1
    Other values of qptopt can be used for debugging purpose.

OUTPUT

  a2f_tr<a2fw_tr_t>=Structure storing the Eliashberg transport function a2F_tr(w).

PARENTS

      m_phgamma

CHILDREN

SOURCE

3627 subroutine a2fw_tr_init(a2f_tr,gams,cryst,ifc,intmeth,wstep,wminmax,smear,ngqpt,nqshift,qshift,comm,&
3628   qintp,qptopt) ! optional
3629 
3630 
3631 !This section has been created automatically by the script Abilint (TD).
3632 !Do not modify the following lines by hand.
3633 #undef ABI_FUNC
3634 #define ABI_FUNC 'a2fw_tr_init'
3635 !End of the abilint section
3636 
3637  implicit none
3638 
3639 !Arguments ------------------------------------
3640 !scalars
3641  integer,intent(in) :: intmeth,nqshift,comm
3642  integer,intent(in),optional :: qptopt
3643  real(dp),intent(in) :: wstep,smear
3644  logical,optional,intent(in) :: qintp
3645  type(phgamma_t),intent(inout) :: gams
3646  type(ifc_type),intent(in) :: ifc
3647  type(crystal_t),intent(in) :: cryst
3648  type(a2fw_tr_t),target,intent(out) :: a2f_tr
3649 !arrays
3650  integer,intent(in) :: ngqpt(3)
3651  real(dp),intent(in) :: wminmax(2),qshift(3,nqshift)
3652 
3653 !Local variables -------------------------
3654 !scalars
3655  integer,parameter :: bcorr0=0,master=0
3656  integer :: my_qptopt,iq_ibz,nqibz,ount,ii,my_rank,nproc,cnt
3657  integer :: mu,iw,natom3,nsppol,spin,ierr,nomega,nqbz
3658  integer :: idir, jdir
3659  real(dp) :: cpu,wall,gflops
3660  real(dp) :: omega,xx,omega_min,omega_max,ww
3661  logical :: do_qintp
3662  character(len=500) :: msg
3663  type(t_tetrahedron) :: tetra
3664 !arrays
3665  integer :: qptrlatt(3,3),new_qptrlatt(3,3)
3666  real(dp),allocatable :: my_qshift(:,:)
3667  real(dp) :: lambda_iso(3,3), omega_log(3,3)
3668  real(dp) :: phfrq(gams%natom3)
3669  real(dp) :: gamma_in_ph(3,3,gams%natom3), gamma_out_ph(3,3,gams%natom3)
3670  real(dp) :: lambda_in_ph(3,3,gams%natom3), lambda_out_ph(3,3,gams%natom3)
3671  real(dp),allocatable :: tmp_a2f_in(:,:,:)
3672  real(dp),allocatable :: tmp_a2f_out(:,:,:)
3673  real(dp), ABI_CONTIGUOUS pointer :: a2f_tr_1d(:)
3674  real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:)
3675  real(dp),allocatable :: a2f_tr_1mom(:),a2f_tr_logmom(:),a2f_tr_logmom_int(:),wdt(:,:)
3676  real(dp),allocatable :: lambda_in_tetra(:,:,:,:,:),phfreq_tetra(:,:,:)
3677  real(dp),allocatable :: lambda_out_tetra(:,:,:,:,:)
3678 
3679 ! *********************************************************************
3680 
3681  !@a2fw_tr_t
3682  my_qptopt = 1; if (present(qptopt)) my_qptopt = qptopt
3683  do_qintp = .True.; if (present(qintp)) do_qintp = qintp
3684 
3685  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3686  nsppol = gams%nsppol; natom3 = gams%natom3
3687 
3688  call cwtime(cpu,wall,gflops,"start")
3689 
3690  if (do_qintp) then
3691    ! Generate the q-mesh by finding the IBZ and the corresponding weights.
3692    qptrlatt = 0
3693    do ii=1,3
3694      qptrlatt(ii,ii) = ngqpt(ii)
3695    end do
3696 
3697    call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, nqibz, qibz, wtq, nqbz, qbz, &
3698       new_kptrlatt=new_qptrlatt, new_shiftk=my_qshift)
3699    ABI_FREE(qbz)
3700 
3701    ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ.
3702    a2f_tr%ngqpt = ngqpt; a2f_tr%nqshift = size(my_qshift, dim=2)
3703    ABI_MALLOC(a2f_tr%qshift, (3, a2f_tr%nqshift))
3704    a2f_tr%qshift = my_qshift
3705    ABI_FREE(my_qshift)
3706 
3707  else
3708    ! No interpolation. Use q-mesh parameters from gams%
3709    a2f_tr%ngqpt = gams%ngqpt; a2f_tr%nqshift = 1
3710    nqibz = gams%nqibz
3711    ABI_MALLOC(qibz,(3,nqibz))
3712    ABI_MALLOC(wtq,(nqibz))
3713    ABI_MALLOC(a2f_tr%qshift,(3,a2f_tr%nqshift))
3714    qibz = gams%qibz; wtq = gams%wtq
3715    a2f_tr%qshift = zero ! Note: assuming q-mesh centered on Gamma.
3716  end if
3717 
3718  call cwtime(cpu,wall,gflops,"stop")
3719  write(msg,'(2(a,f8.2))')"a2fw_tr_init, q-setup: ",cpu,", wall: ",wall
3720  call wrtout(std_out,msg,"COLL",do_flush=.True.)
3721 
3722  ! Define Min and max frequency for the mesh (enlarge it a bit)
3723  omega_min = wminmax(1); omega_max = wminmax(2)
3724  omega_min = omega_min - 0.1*abs(omega_min)
3725  if (omega_min >= zero) omega_min = one/Ha_meV
3726  omega_max = omega_max + 0.1*abs(omega_max)
3727 
3728  a2f_tr%nsppol = nsppol; a2f_tr%natom3 = gams%natom3; a2f_tr%smear = smear
3729  a2f_tr%omega_min = omega_min; a2f_tr%omega_max = omega_max
3730  nomega = int((omega_max - omega_min) / wstep); a2f_tr%nomega = nomega; a2f_tr%wstep = wstep
3731 
3732  ABI_MALLOC(a2f_tr%n0,(nsppol))
3733  a2f_tr%n0=gams%n0
3734  ! Build linear mesh.
3735  ABI_MALLOC(a2f_tr%omega, (nomega))
3736  a2f_tr%omega = arth(omega_min,wstep,nomega)
3737  ABI_CALLOC(a2f_tr%vals_in, (nomega,3,3,0:natom3,nsppol))
3738  ABI_CALLOC(a2f_tr%vals_out, (nomega,3,3,0:natom3,nsppol))
3739  ABI_CALLOC(a2f_tr%vals_tr, (nomega,3,3,0:natom3,nsppol))
3740  ABI_CALLOC(a2f_tr%lambdaw_tr, (nomega,3,3,0:natom3, nsppol))
3741  ABI_MALLOC(tmp_a2f_in, (nomega,3,3))
3742  ABI_MALLOC(tmp_a2f_out, (nomega,3,3))
3743 
3744  if (intmeth == 2) then
3745    call cwtime(cpu,wall,gflops,"start")
3746 
3747    ! Prepare tetrahedron integration.
3748    qptrlatt = 0
3749    do ii=1,3
3750      qptrlatt(ii,ii) = a2f_tr%ngqpt(ii)
3751    end do
3752 
3753    tetra = tetra_from_kptrlatt(cryst, my_qptopt, qptrlatt, a2f_tr%nqshift, a2f_tr%qshift, nqibz, qibz, msg, ierr)
3754    if (ierr/=0) MSG_ERROR(msg)
3755 
3756    ABI_STAT_MALLOC(lambda_in_tetra, (nqibz,3,3,natom3,nsppol), ierr)
3757    ABI_CHECK(ierr==0, "oom in lambda_in_tetra, use gaussians for A2F_tr")
3758    ABI_STAT_MALLOC(lambda_out_tetra, (nqibz,3,3,natom3,nsppol), ierr)
3759    ABI_CHECK(ierr==0, "oom in lambda_out_tetra, use gaussians for A2F_tr")
3760    lambda_in_tetra = zero
3761    lambda_out_tetra = zero
3762 
3763    ABI_STAT_MALLOC(phfreq_tetra, (nqibz,natom3,nsppol), ierr)
3764    ABI_CHECK(ierr==0, "oom in phfreq_tetra, use gaussians for A2F_tr")
3765    phfreq_tetra = zero
3766 
3767    call cwtime(cpu,wall,gflops,"stop")
3768    write(msg,'(2(a,f8.2))')"a2fw_tr_init%tetra, cpu",cpu,", wall: ",wall
3769    call wrtout(std_out,msg,"COLL",do_flush=.True.)
3770  end if
3771 
3772  call cwtime(cpu,wall,gflops,"start")
3773 
3774  ! Loop over spins and qpoints in the IBZ
3775  cnt = 0
3776  do spin=1,nsppol
3777    do iq_ibz=1,nqibz
3778      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
3779 
3780      ! Interpolate or evaluate gamma directly.
3781      if (do_qintp) then
3782        call phgamma_vv_interp(gams,cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph)
3783      else
3784        call phgamma_vv_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph)
3785      end if
3786 
3787      select case (intmeth)
3788      case (1)
3789        ! Gaussian: Add all contributions from the phonon modes at this qpoint to a2f_tr
3790        ! (note that unstable modes are included).
3791        do mu=1,natom3
3792          tmp_a2f_in = zero
3793          tmp_a2f_out = zero
3794          do iw=1,nomega
3795            xx = a2f_tr%omega(iw) - phfrq(mu)
3796            tmp_a2f_in(iw,:,:)  = tmp_a2f_in(iw,:,:)  + dirac_delta(xx, smear) * lambda_in_ph(:,:,mu)  * abs(phfrq(mu))
3797            tmp_a2f_out(iw,:,:) = tmp_a2f_out(iw,:,:) + dirac_delta(xx, smear) * lambda_out_ph(:,:,mu) * abs(phfrq(mu))
3798          end do
3799          a2f_tr%vals_in(:,:,:,mu,spin)  = a2f_tr%vals_in(:,:,:,mu,spin)  + tmp_a2f_in(:,:,:) * wtq(iq_ibz)
3800          a2f_tr%vals_out(:,:,:,mu,spin) = a2f_tr%vals_out(:,:,:,mu,spin) + tmp_a2f_out(:,:,:) * wtq(iq_ibz)
3801        end do
3802 
3803      case (2)
3804        ! Tetra: store data.
3805        do mu=1,natom3
3806          lambda_in_tetra(iq_ibz, :,:, mu, spin) = lambda_in_ph(:,:,mu) * abs(phfrq(mu))
3807          lambda_out_tetra(iq_ibz, :,:, mu, spin) = lambda_out_ph(:,:,mu) * abs(phfrq(mu))
3808          phfreq_tetra(iq_ibz, mu, spin) = phfrq(mu)
3809        end do
3810      end select
3811 
3812    end do ! iq_ibz
3813  end do ! spin
3814 
3815  if (intmeth == 2) then
3816    ! Collect results on each node.
3817    call xmpi_sum(lambda_in_tetra, comm, ierr)
3818    call xmpi_sum(lambda_out_tetra, comm, ierr)
3819    call xmpi_sum(phfreq_tetra, comm, ierr)
3820 
3821    ! workspace for tetra.
3822    ABI_MALLOC(wdt, (nomega, 2))
3823 
3824 ! TODO: with the tetra_get_onewk call we can integrate this above and avoid allocating all of lambda_in_tetra and phfreq_tetra!!!
3825    ! For each mode get its contribution
3826    cnt = 0
3827    do spin=1,nsppol
3828      do mu=1,natom3
3829        do iq_ibz=1,nqibz
3830          cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! mpi-parallelism
3831 
3832          call tetra_get_onewk(tetra, iq_ibz, bcorr0, nomega, nqibz, phfreq_tetra(:,mu,spin), &
3833            omega_min, omega_max, one, wdt)
3834 
3835          ! Accumulate (Integral of a2F_tr is computed afterwards)
3836          do idir=1,3
3837            do jdir=1,3
3838              a2f_tr%vals_in(:,idir,jdir,mu,spin)  = a2f_tr%vals_in(:,idir,jdir,mu,spin)  &
3839 &               + wdt(:,1) * lambda_in_tetra(iq_ibz,idir,jdir,mu,spin)
3840              a2f_tr%vals_out(:,idir,jdir,mu,spin) = a2f_tr%vals_out(:,idir,jdir,mu,spin) &
3841 &               + wdt(:,1) * lambda_out_tetra(iq_ibz,idir,jdir,mu,spin)
3842              !a2f_tr%lambdaw(:,mu,spin) = a2f_tr%lambdaw(:,mu,spin) + wdt(:,2) * lambda_XX_tetra(iq_ibz, mu, spin)
3843            end do
3844          end do
3845        end do
3846      end do
3847    end do
3848 
3849    ! Free memory allocated for tetra.
3850    ABI_FREE(wdt)
3851    ABI_FREE(lambda_in_tetra)
3852    ABI_FREE(lambda_out_tetra)
3853    ABI_FREE(phfreq_tetra)
3854    call destroy_tetra(tetra)
3855  end if
3856 
3857  ! Collect final results on each node and divide by g(eF, spin)
3858 ! TODO: check whether this has to be reinstated, in particular for a2f_tr
3859  call xmpi_sum(a2f_tr%vals_in, comm, ierr)
3860  call xmpi_sum(a2f_tr%vals_out, comm, ierr)
3861 
3862 ! TODO: check normalization with N(0) and <v^2> from Savrasov eq 20
3863  do spin=1,nsppol
3864    a2f_tr%vals_in(:,:,:,0,spin)  = sum(a2f_tr%vals_in(:,:,:,1:natom3,spin), dim=4)
3865    a2f_tr%vals_out(:,:,:,0,spin) = sum(a2f_tr%vals_out(:,:,:,1:natom3,spin), dim=4)
3866    !a2f_tr%vals(:,:,spin) = a2f_tr%vals(:,:,spin) / (two_pi*a2f_tr%n0(spin))
3867  end do
3868 
3869  !to avoid numerical noise uses a smoothing function
3870  ! TODO: Move smooth to m_numeric_tools and add ndat dimension.
3871  !do spin=1,nsppol
3872  !  do mu=0,natom3
3873  !    call smooth(a2f_tr%vals(:,mu,spin),a2f_tr%nomega,2)
3874  !  end do
3875  !end do
3876 
3877  a2f_tr%vals_tr = a2f_tr%vals_out - a2f_tr%vals_in
3878 
3879  call cwtime(cpu,wall,gflops,"stop")
3880  write(msg,'(2(a,f8.2))')"a2fw_tr_init, a2f_eval: ",cpu,", wall: ",wall
3881  call wrtout(std_out,msg,"COLL",do_flush=.True.)
3882 
3883 ! NOTE: for the moment calculate the log moms etc on just the trace of the a2f functions
3884  ABI_MALLOC(a2f_tr_1mom,(nomega))
3885  ABI_MALLOC(a2f_tr_logmom,(nomega))
3886  ABI_MALLOC(a2f_tr_logmom_int,(nomega))
3887 
3888  !call a2fw_tr_print_info()
3889 
3890  ! Compute lambda(w) (resolved in spin mode and directions)
3891  do spin=1,nsppol
3892    do mu=0,natom3
3893      do idir=1,3
3894        do jdir=1,3
3895 
3896          do iw=1,nomega
3897            ww = a2f_tr%omega(iw)
3898            if (abs(ww) > EPH_WTOL) then
3899               a2f_tr_1mom(iw)  = a2f_tr%vals_tr(iw,idir,jdir,mu,spin) / abs(ww)
3900            else
3901               a2f_tr_1mom(iw) = zero
3902            end if
3903            ! TODO: Strange that the first value of int(f) is not set to zero!
3904            ! FIXME: The frequency integration overestimates lambda if there are negative phonon frequencies!
3905          end do
3906 
3907          call simpson_int(nomega, wstep, a2f_tr_1mom, a2f_tr%lambdaw_tr(:,idir,jdir,mu,spin))
3908        end do ! jdir
3909      end do ! idir
3910    end do ! mu
3911  end do ! spin
3912 
3913  if (my_rank == master) then
3914    ount = std_out
3915    if (nsppol > 1) then
3916      write(msg,'(3a)') ch10,&
3917        'Comment: some of the following quantities should be integrated over spin', ch10
3918      call wrtout(ount,msg,'COLL')
3919    end if
3920  end if
3921 
3922  do spin=1,nsppol
3923    lambda_iso = a2fw_tr_moment(a2f_tr,0,spin)
3924 
3925    do idir=1,3
3926      do jdir=1,3
3927        a2f_tr_1d => a2f_tr%vals_tr(:,idir,jdir,0,spin)
3928 
3929        ! Get log moment of alpha^2F.
3930        a2f_tr_logmom = zero
3931        do iw=1,nomega
3932          omega = a2f_tr%omega(iw)
3933          if (abs(omega) > EPH_WTOL .and. abs(lambda_iso(idir,jdir)) > EPH_WTOL) then
3934            a2f_tr_logmom(iw) = (two/lambda_iso(idir,jdir)) * a2f_tr_1d(iw)*log(abs(omega))/abs(omega)
3935          end if
3936        end do
3937        call simpson_int(nomega,wstep,a2f_tr_logmom,a2f_tr_logmom_int)
3938 !DEBUG
3939 !      write(std_out,*)' iw,nomega,greatest_real,a2f_tr_logmom_int(nomega)=',&
3940 !          & iw,nomega,greatest_real,a2f_tr_logmom_int(nomega)
3941 !ENDDEBUG
3942        if(abs(a2f_tr_logmom_int(nomega))<log(greatest_real*tol6))then
3943          omega_log(idir,jdir) = exp(a2f_tr_logmom_int(nomega))
3944        else
3945          omega_log(idir,jdir)=greatest_real*tol6
3946        endif
3947      end do
3948    end do
3949 
3950 
3951 ! TODO: make output only for irred values xx yy zz and top half of matrix
3952    if (my_rank == master) then
3953      write(ount,'(a)')' Evaluation of parameters analogous to electron-phonon coupling for 3x3 directions '
3954      write(ount,'(a,3(3es10.3,2x))') ' lambda = ',lambda_iso
3955      write(ount,'(a,3(3es10.3,2x),a)' )' omegalog  = ',omega_log,' (Ha) '
3956      write(ount,'(a,3(3es10.3,2x),a)' )'             ',omega_log*Ha_K, ' (Kelvin) '
3957      write(ount,"(a)")'    positive moments of alpha2F:'
3958      write(ount,'(a,3(3es10.3,2x))' )' lambda <omega^2> = ',a2fw_tr_moment(a2f_tr,2,spin)
3959      write(ount,'(a,3(3es10.3,2x))' )' lambda <omega^3> = ',a2fw_tr_moment(a2f_tr,3,spin)
3960      write(ount,'(a,3(3es10.3,2x))' )' lambda <omega^4> = ',a2fw_tr_moment(a2f_tr,4,spin)
3961      write(ount,'(a,3(3es10.3,2x))' )' lambda <omega^5> = ',a2fw_tr_moment(a2f_tr,5,spin)
3962    end if
3963  end do
3964 
3965  ABI_FREE(tmp_a2f_in)
3966  ABI_FREE(tmp_a2f_out)
3967  ABI_FREE(a2f_tr_1mom)
3968  ABI_FREE(a2f_tr_logmom)
3969  ABI_FREE(a2f_tr_logmom_int)
3970  ABI_FREE(qibz)
3971  ABI_FREE(wtq)
3972 
3973 end subroutine a2fw_tr_init

m_phgamma/a2fw_tr_moment [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_tr_moment

FUNCTION

  Compute \int dw [a2F_tr(w)/w] w^n
  From Allen PRL 59 1460 [[cite:Allen1987]] and later PRB papers (See also [[cite:Grimvall1981]] book)

INPUTS

  a2f_tr<a2fw_tr_t>=Structure storing the Eliashberg function.
  nn=Value of n
  spin=The spin component

OUTPUT

  a2fw_tr_moment = \int dw [a2F_tr(w)/w] w^n
  [out_int(x)] = \int^{x} dw [a2F_tr(w)/w] w^n

PARENTS

CHILDREN

SOURCE

2789 function a2fw_tr_moment(a2f_tr,nn,spin,out_int)
2790 
2791 
2792 !This section has been created automatically by the script Abilint (TD).
2793 !Do not modify the following lines by hand.
2794 #undef ABI_FUNC
2795 #define ABI_FUNC 'a2fw_tr_moment'
2796 !End of the abilint section
2797 
2798  implicit none
2799 
2800  real(dp), dimension(3,3) :: a2fw_tr_moment
2801 
2802 !Arguments ------------------------------------
2803 !scalars
2804  integer,intent(in) :: spin,nn
2805  type(a2fw_tr_t),intent(in) :: a2f_tr
2806 !arrays
2807  real(dp),intent(out),optional :: out_int(a2f_tr%nomega,3,3)
2808 
2809 !Local variables -------------------------
2810 !scalars
2811  integer :: iw !,nomega
2812  integer :: idir, jdir
2813  real(dp) :: omg,omg_nm1
2814 !arrays
2815  real(dp) :: ff(a2f_tr%nomega),int_ff(a2f_tr%nomega)
2816 
2817 ! *********************************************************************
2818 
2819  do jdir = 1, 3
2820    do idir = 1, 3
2821      ! Construct the integrand function. [a2F_tr(w)/w] w^n
2822      ff = zero; int_ff = zero
2823 
2824      if (nn-1 >= 0) then
2825        do iw=1,a2f_tr%nomega
2826          omg = a2f_tr%omega(iw)
2827          omg_nm1 = omg ** (nn-1)
2828          ff(iw) = a2f_tr%vals_tr(iw,idir,jdir,0,spin) * omg_nm1
2829        end do
2830      else
2831        do iw=1,a2f_tr%nomega
2832          omg = a2f_tr%omega(iw)
2833          omg_nm1 = zero; if (abs(omg) > EPH_WTOL) omg_nm1 = omg**(nn-1)
2834          ff(iw) = a2f_tr%vals_tr(iw,idir,jdir,0,spin) * omg_nm1
2835        end do
2836      end if
2837 
2838      ! Integration with simpson rule on a linear mesh.
2839      call simpson_int(a2f_tr%nomega,a2f_tr%wstep,ff,int_ff)
2840 
2841      a2fw_tr_moment(idir,jdir) = int_ff(a2f_tr%nomega)
2842      if (present(out_int)) out_int(:,idir,jdir) = int_ff(:)
2843    end do
2844  end do
2845 
2846 end function a2fw_tr_moment

m_phgamma/a2fw_tr_t [ Types ]

[ Top ] [ m_phgamma ] [ Types ]

NAME

 a2fw_tr_t

FUNCTION

 Store the Eliashberg transport spectral functions:
    a2F_trin(w, x, x')
    a2F_trout(w, x, x')
    a2F_tr(w, x, x') = in - out
    a2F_tr_gen(e, e', w, x, x')

SOURCE

302  type,public :: a2fw_tr_t
303 
304   integer :: nomega
305   ! Number of frequency points in a2f_tr(w).
306 
307   integer :: nene
308   ! Number of chemical potential values used for inelastic integration
309   ! Number of electron points in a2f_tr_gen(e,e',w).
310 
311   integer :: nsppol
312   ! Number of independent spin polarizations.
313 
314   integer :: natom3
315   ! Number of phonon modes.
316 
317   real(dp) :: enemin
318   ! Minimal chemical potential value used for inelastic integration
319 
320   real(dp) :: deltaene
321   ! Chemical potential increment for inelastic integration
322 
323   real(dp) :: omega_min,omega_max
324   ! min and Max frequency (Ha) in the linear mesh.
325 
326   real(dp) :: wstep
327   ! Step of the linear mesh
328 
329   real(dp) :: smear
330   ! Gaussian broadening used to approximated the Dirac distribution.
331 
332   integer :: nqshift
333   ! Number of shifts in the q-mesh
334 
335   integer :: ngqpt(3)
336   ! The q-mesh used for calculating vals(w).
337 
338   real(dp),allocatable :: qshift(:,:)
339   ! qshift(3,nqshift)
340   ! The shifts used to generate the q-mesh.
341 
342   real(dp),allocatable :: n0(:)
343   ! n0(nsppol)
344   ! Electronic DOS at the Fermi level.
345 
346   real(dp),allocatable :: omega(:)
347   ! omega(nomega)
348   ! Frequency mesh in Hartree (linear).
349 
350   real(dp),allocatable :: vals_in(:,:,:,:,:)
351   real(dp),allocatable :: vals_out(:,:,:,:,:)
352   ! vals_in(nomega,3,3,0:natom3,nsppol)
353   ! Eliashberg transport functions for in and out scattering
354   !   vals_in(w,3,3,1:natom3,1:nsppol): a2f_tr(w) decomposed per phonon branch and spin
355   !   vals_in(w,3,3,0,1:nsppol): a2f_tr(w) summed over phonons modes, decomposed in spin
356 
357   real(dp),allocatable :: vals_tr(:,:,:,:,:)
358   ! vals_tr(nomega,3,3,0:natom3,nsppol)
359   ! transport spectral function = in-out
360 
361   real(dp),allocatable :: vals_tr_gen(:,:,:,:,:,:,:)
362   ! vals(nene,nene,nomega,3,3,nsppol)
363   ! generalized transport spectral function from PB Allen Phys. Rev. Lett. 59, 1460 (1987) [[cite:Allen1987]]
364 
365   real(dp),allocatable :: lambdaw_tr(:,:,:,:,:)
366   ! lambda(nomega,3,3,0:natom3,nsppol)
367 
368  end type a2fw_tr_t
369 
370  public :: a2fw_tr_free            ! Free the memory allocated in the structure.
371  public :: a2fw_tr_init            ! Calculates the FS averaged alpha^2F_tr,in,out(w, x, x') functions.
372  public :: a2fw_tr_write           ! Write alpha^2F(w) to an external file in text/netcdf format

m_phgamma/a2fw_tr_write [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_tr_write

FUNCTION

  Write alpha^2F_tr(w) to an external file in text form

INPUTS

  a2f_tr<a2fw_tr_t>=Container storing the Eliashberg transport functions.
  basename=Filename for output.
  post=String appended to netcdf variables e.g. _qcoarse, _qintp
  ncid=Netcdf file handler. Set it to nctk_noid to disable output.

OUTPUT

  Output is written to file. This routine should be called by one MPI proc.

PARENTS

      m_phgamma

CHILDREN

SOURCE

4001 subroutine a2fw_tr_write(a2f_tr, basename, post, ncid)
4002 
4003 
4004 !This section has been created automatically by the script Abilint (TD).
4005 !Do not modify the following lines by hand.
4006 #undef ABI_FUNC
4007 #define ABI_FUNC 'a2fw_tr_write'
4008 !End of the abilint section
4009 
4010  implicit none
4011 
4012 !Arguments ------------------------------------
4013 !scalars
4014  integer,intent(in) :: ncid
4015  character(len=*),intent(in) :: basename, post
4016  type(a2fw_tr_t),intent(in) :: a2f_tr
4017 
4018 !Local variables -------------------------
4019 !scalars
4020  integer :: iw,spin,unt,ii,mu
4021  integer :: idir, jdir
4022 #ifdef HAVE_NETCDF
4023  integer :: ncerr
4024  character(len=500) :: dim1_name
4025 #endif
4026  character(len=500) :: msg
4027  character(len=fnlen) :: path
4028 
4029 ! *********************************************************************
4030 
4031  ! Write spin-resolved a2F_tr(w)
4032  path = strcat(basename, "_A2FW_tr")
4033  if (open_file(path,msg,newunit=unt,form="formatted",status="unknown") /= 0) then
4034    MSG_ERROR(msg)
4035  end if
4036 
4037  call write_a2fw_tr_header()
4038 
4039  if (a2f_tr%nsppol == 1) then
4040    write(unt,'(a)')"# Frequency, a2F_tr(w,i,j)"
4041    do iw=1,a2f_tr%nomega
4042      write(unt,'(e16.6, 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%vals_tr(iw,:,:,0,1)
4043    end do
4044    write(unt,'(3a)')ch10, ch10, "# Frequency, lambda(w,i,j)"
4045    do iw=1,a2f_tr%nomega
4046      write(unt,'(e16.6, 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%lambdaw_tr(iw,:,:,0,1)
4047    end do
4048 
4049  else
4050    write(unt,'(a)')"# Frequency, a2F_tr_tot(w,i,j), a2F_tr_spin1(w,i,j) ..."
4051    do iw=1,a2f_tr%nomega
4052      write(unt,'(e16.6, 2x, 3(3e16.6,1x), 2x,  3(3e16.6,1x), 2x,  3(3e16.6,1x))') a2f_tr%omega(iw), &
4053 &                 ((sum(a2f_tr%vals_tr(iw,idir,jdir,0,:)), idir=1,3), jdir=1,3) , &
4054 &                 a2f_tr%vals_tr(iw,:,:,0,1)      ,      &  ! UP
4055 &                 a2f_tr%vals_tr(iw,:,:,0,2)                ! DOWN
4056    end do
4057    write(unt,'(3a)') ch10, ch10, "# Frequency, lambda_tr_tot(w,i,j)dw, lambda_tr_spin1(w,i,j) ..."
4058    do iw=1,a2f_tr%nomega
4059      write(unt,'(e16.6, 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x,      3(3e16.6,1x), 2x,  3(3e16.6,1x))') a2f_tr%omega(iw), &
4060 &                 ((sum(a2f_tr%lambdaw_tr(iw,idir,jdir,0,:)), idir=1,3), jdir=1,3) , &  ! TOT
4061 &                 a2f_tr%lambdaw_tr(iw,:,:,0,1),      &  ! UP
4062 &                 a2f_tr%lambdaw_tr(iw,:,:,0,2)          ! DOWN
4063    end do
4064  end if
4065 
4066  close(unt)
4067 
4068  ! Write phonon mode contributions to a2F_tr(w,i,j)
4069  path = strcat(basename, "_PH_A2FW_tr")
4070  if (open_file(path,msg,newunit=unt,form="formatted",status="unknown") /= 0) then
4071    MSG_ERROR(msg)
4072  end if
4073  call write_a2fw_tr_header()
4074 
4075  if (a2f_tr%nsppol == 1) then
4076    do mu=0,a2f_tr%natom3
4077      write(unt,'(a,i0)')"# Phonon mode ",mu
4078      write(unt,'(a)')"# Frequency, a2F_tr(w)"
4079      do iw=1,a2f_tr%nomega
4080        write(unt,'(e16.6,2x,3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%vals_tr(iw,:,:,mu,1)
4081      end do
4082      do ii=1,2; write(unt,'(a)')""; end do
4083      write(unt,'(a)')"# Frequency, lambda_tr(w)"
4084      do iw=1,a2f_tr%nomega
4085        write(unt,'(e16.6,2x,3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%lambdaw_tr(iw,:,:,mu,1)
4086      end do
4087      do ii=1,2; write(unt,'(a)')""; end do
4088    end do
4089 
4090  else
4091    do mu=0,a2f_tr%natom3
4092      write(unt,'(a,i0)')"# Phonon mode ",mu
4093      write(unt,'(a)')"# Frequency, a2F_tot(w), a2F_spin1(w) ..."
4094      do iw=1,a2f_tr%nomega
4095        write(unt,*) a2f_tr%omega(iw), &
4096                     sum(a2f_tr%vals_tr(iw,:,:,mu,:)), &  ! TOT
4097                     a2f_tr%vals_tr(iw,:,:,mu,1)     , &  ! UP
4098                     a2f_tr%vals_tr(iw,:,:,mu,2)          ! DOWN
4099      end do
4100      write(unt,'(3a)')ch10,ch10,"# Frequency, lambda_tot(w)dw, lambda_spin1(w) ..."
4101      do iw=1,a2f_tr%nomega
4102        write(unt,*) a2f_tr%omega(iw), &
4103                      sum(a2f_tr%lambdaw_tr(iw,:,:,mu,:)), &  ! TOT
4104                      a2f_tr%lambdaw_tr(iw,:,:,mu,1),      &  ! UP
4105                      a2f_tr%lambdaw_tr(iw,:,:,mu,2)          ! DOWN
4106      end do
4107    end do
4108  end if
4109 
4110  close(unt)
4111 
4112  if (ncid /= nctk_noid) then
4113 #ifdef HAVE_NETCDF
4114    ! Define dimensions.
4115    dim1_name = strcat("a2ftr_nomega", post)
4116    ncerr = nctk_def_dims(ncid, [ &
4117      nctkdim_t(dim1_name, a2f_tr%nomega), &
4118      nctkdim_t("natom3p1", a2f_tr%natom3 + 1), nctkdim_t("number_of_spins", a2f_tr%nsppol) &
4119      ], defmode=.True.)
4120    NCF_CHECK(ncerr)
4121 
4122    ncerr = nctk_def_arrays(ncid, [ &
4123      nctkarr_t(strcat('a2ftr_mesh', post), "dp", dim1_name), &
4124      nctkarr_t(strcat('a2ftr_values', post), "dp", strcat(dim1_name, ", three, three, natom3p1, number_of_spins")), &
4125      nctkarr_t(strcat('a2ftr_lambdaw', post), "dp", strcat(dim1_name, ", three, three, natom3p1, number_of_spins")) &
4126      ])
4127    NCF_CHECK(ncerr)
4128 
4129    NCF_CHECK(nctk_set_datamode(ncid))
4130    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_mesh", post)), a2f_tr%omega))
4131    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_values", post)), a2f_tr%vals_tr))
4132    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_lambdaw", post)), a2f_tr%lambdaw_tr))
4133 #endif
4134  end if
4135 
4136 contains
4137 
4138 subroutine write_a2fw_tr_header()
4139 
4140  ! Output the header.
4141 
4142 !This section has been created automatically by the script Abilint (TD).
4143 !Do not modify the following lines by hand.
4144 #undef ABI_FUNC
4145 #define ABI_FUNC 'write_a2fw_tr_header'
4146 !End of the abilint section
4147 
4148  write(unt,'(a)')              '#'
4149  write(unt,'(a)')              '# ABINIT package: a2F_tr(w) file'
4150  write(unt,'(a)')              '#'
4151  write(unt,'(a)')              '# a2F_tr(w) function integrated over the FS. omega in a.u.'
4152  write(unt,'(2a)')             '# ngqpt: ',trim(ltoa(a2f_tr%ngqpt))
4153  write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f_tr%nomega," between omega_min: ",a2f_tr%omega_min,&
4154                                ' Ha and omega_max: ',a2f_tr%omega_max,' Ha with step:',a2f_tr%wstep
4155  write(unt,'(a,e16.6)')         '#  the smearing width for gaussians is ',a2f_tr%smear
4156  write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f_tr%n0)
4157  do spin=1,a2f_tr%nsppol
4158    write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f_tr%n0(spin)
4159  end do
4160  do ii=1,2; write(unt,'(a)')     "# "; end do
4161 
4162 end subroutine write_a2fw_tr_header
4163 
4164 end subroutine a2fw_tr_write

m_phgamma/a2fw_write [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 a2fw_write

FUNCTION

  Write alpha^2F(w) to an external file in text form

INPUTS

  a2f<a2fw_t>=Container storing the Eliashberg functions.
  basename=Filename for output.
  post=String appended to netcdf variables e.g. _qcoarse, _qintp
  ncid=Netcdf file handler. Set it to nctk_noid to disable output.

OUTPUT

  Output is written to file. This routine should be called by one MPI proc.

PARENTS

      m_phgamma

CHILDREN

SOURCE

3009 subroutine a2fw_write(a2f, basename, post, ncid)
3010 
3011 
3012 !This section has been created automatically by the script Abilint (TD).
3013 !Do not modify the following lines by hand.
3014 #undef ABI_FUNC
3015 #define ABI_FUNC 'a2fw_write'
3016 !End of the abilint section
3017 
3018  implicit none
3019 
3020 !Arguments ------------------------------------
3021 !scalars
3022  integer,intent(in) :: ncid
3023  character(len=*),intent(in) :: basename, post
3024  type(a2fw_t),intent(in) :: a2f
3025 
3026 !Local variables -------------------------
3027 !scalars
3028  integer :: iw,spin,unt,ii,mu
3029 #ifdef HAVE_NETCDF
3030  integer :: ncerr
3031  character(len=500) :: dim1_name
3032 #endif
3033  character(len=500) :: msg
3034  character(len=fnlen) :: path
3035 
3036 ! *********************************************************************
3037 
3038  ! Write spin-resolved a2F(w)
3039  path = strcat(basename, "_A2FW")
3040  if (open_file(path,msg,newunit=unt,form="formatted",status="unknown") /= 0) then
3041    MSG_ERROR(msg)
3042  end if
3043 
3044  call write_a2fw_header()
3045 
3046  if (a2f%nsppol == 1) then
3047    write(unt,'(a)')"# Frequency, a2F(w), lambda(w)"
3048    do iw=1,a2f%nomega
3049      write(unt,*) a2f%omega(iw), a2f%vals(iw,0,1), a2f%lambdaw(iw,0,1)
3050    end do
3051 
3052  else
3053    write(unt,'(a)')"# Frequency, a2F_tot(w), lambda_tot(w)dw, a2F_spin1(w), lambda_spin1(w) ..."
3054    do iw=1,a2f%nomega
3055      write(unt,*) a2f%omega(iw), &
3056                   sum(a2f%vals(iw,0,:)) , sum(a2f%lambdaw(iw,0,:)), &  ! TOT
3057                   a2f%vals(iw,0,1)      , a2f%lambdaw(iw,0,1),      &  ! UP
3058                   a2f%vals(iw,0,2)      , a2f%lambdaw(iw,0,2)          ! DOWN
3059    end do
3060  end if
3061 
3062  close(unt)
3063 
3064  ! Write phonon contributions to a2F(w)
3065  path = strcat(basename, "_PH_A2FW")
3066  if (open_file(path,msg,newunit=unt,form="formatted",status="unknown") /= 0) then
3067    MSG_ERROR(msg)
3068  end if
3069  call write_a2fw_header()
3070 
3071  if (a2f%nsppol == 1) then
3072    do mu=0,a2f%natom3
3073      write(unt,'(a,i0)')"# Phonon mode ",mu
3074      write(unt,'(a)')"# Frequency, a2F(w), lambda(w)"
3075      do iw=1,a2f%nomega
3076        write(unt,*) a2f%omega(iw), a2f%vals(iw,mu,1), a2f%lambdaw(iw,mu,1)
3077      end do
3078      do ii=1,2; write(unt,'(a)')""; end do
3079    end do
3080 
3081  else
3082    do mu=0,a2f%natom3
3083      write(unt,'(a,i0)')"# Phonon mode ",mu
3084      write(unt,'(a)')"# Frequency, a2F_tot(w), lambda_tot(w)dw, a2F_spin1(w), lambda_spin1(w) ..."
3085      do iw=1,a2f%nomega
3086        write(unt,*) a2f%omega(iw), &
3087                     sum(a2f%vals(iw,mu,:)), sum(a2f%lambdaw(iw,mu,:)), &  ! TOT
3088                     a2f%vals(iw,mu,1)     , a2f%lambdaw(iw,mu,1),      &  ! UP
3089                     a2f%vals(iw,mu,2)     , a2f%lambdaw(iw,mu,2)          ! DOWN
3090      end do
3091    end do
3092  end if
3093 
3094  close(unt)
3095 
3096  if (ncid /= nctk_noid) then
3097 #ifdef HAVE_NETCDF
3098    ! Define dimensions.
3099    dim1_name = strcat("a2f_nomega", post)
3100    ncerr = nctk_def_dims(ncid, [ &
3101      nctkdim_t(dim1_name, a2f%nomega), &
3102      nctkdim_t("natom3p1", a2f%natom3 + 1), nctkdim_t("number_of_spins", a2f%nsppol) &
3103      ], defmode=.True.)
3104    NCF_CHECK(ncerr)
3105 
3106    ncerr = nctk_def_arrays(ncid, [ &
3107      nctkarr_t(strcat('a2f_mesh', post), "dp", dim1_name), &
3108      nctkarr_t(strcat('a2f_values', post), "dp", strcat(dim1_name, ", natom3p1, number_of_spins")), &
3109      nctkarr_t(strcat('a2f_lambdaw', post), "dp", strcat(dim1_name, ", natom3p1, number_of_spins")) &
3110      ])
3111    NCF_CHECK(ncerr)
3112 
3113    NCF_CHECK(nctk_set_datamode(ncid))
3114    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_mesh", post)), a2f%omega))
3115    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_values", post)), a2f%vals))
3116    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_lambdaw", post)), a2f%lambdaw))
3117 #endif
3118  end if
3119 
3120 contains
3121 
3122 subroutine write_a2fw_header()
3123 
3124  ! Output the header.
3125 
3126 !This section has been created automatically by the script Abilint (TD).
3127 !Do not modify the following lines by hand.
3128 #undef ABI_FUNC
3129 #define ABI_FUNC 'write_a2fw_header'
3130 !End of the abilint section
3131 
3132  write(unt,'(a)')              '#'
3133  write(unt,'(a)')              '# ABINIT package: a2F(w) file'
3134  write(unt,'(a)')              '#'
3135  write(unt,'(a)')              '# a2F(w) function integrated over the FS. omega in a.u.'
3136  write(unt,'(2a)')             '# ngqpt: ',trim(ltoa(a2f%ngqpt))
3137  write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f%nomega," between omega_min: ",a2f%omega_min,&
3138                                ' Ha and omega_max: ',a2f%omega_max,' Ha with step:',a2f%wstep
3139  write(unt,'(a,e16.6)')         '#  the smearing width for gaussians is ',a2f%smear
3140  write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f%n0)
3141  do spin=1,a2f%nsppol
3142    write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f%n0(spin)
3143  end do
3144  do ii=1,2; write(unt,'(a)')     "# "; end do
3145 
3146 end subroutine write_a2fw_header
3147 
3148 end subroutine a2fw_write

m_phgamma/eph_phgamma [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

  eph_phgamma

FUNCTION

  Compute phonon linewidths in metals.

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

4199 subroutine eph_phgamma(wfk0_path,dtfil,ngfft,ngfftf,dtset,cryst,ebands,dvdb,ddk,ifc,&
4200                        pawfgr,pawang,pawrad,pawtab,psps,mpi_enreg,comm)
4201 
4202  use defs_basis
4203  use defs_datatypes
4204  use defs_abitypes
4205  use m_abicore
4206  use m_xmpi
4207  use m_errors
4208  use m_wfk
4209  use m_ddk
4210  use m_ddb
4211  use m_dvdb
4212  use m_ifc
4213  use m_fft
4214  use m_hamiltonian
4215  use m_pawcprj
4216 
4217  use m_time,            only : sec2str
4218  use m_fstrings,        only : sjoin, itoa, ftoa, ktoa
4219  use m_io_tools,        only : iomode_from_fname
4220  use m_cgtools,         only : dotprod_g
4221  use m_fftcore,         only : get_kg
4222  use m_crystal,         only : crystal_t
4223  use m_wfd,             only : wfd_init, wfd_free, wfd_print, wfd_t, wfd_test_ortho, wfd_copy_cg,&
4224                                wfd_read_wfk, wfd_wave_free, wfd_rotate, wfd_reset_ur_cprj, wfd_get_ur
4225  use m_pawang,          only : pawang_type
4226  use m_pawrad,          only : pawrad_type
4227  use m_pawtab,          only : pawtab_type
4228  use m_pawfgr,          only : pawfgr_type
4229 ! use m_paw_an,          only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
4230 ! use m_paw_ij,          only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
4231 ! use m_pawfgrtab,       only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_init
4232 ! use m_pawrhoij,        only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, symrhoij
4233 ! use m_pawdij,          only : pawdij, symdij
4234 ! use m_pawcprj,         only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
4235 
4236 !This section has been created automatically by the script Abilint (TD).
4237 !Do not modify the following lines by hand.
4238 #undef ABI_FUNC
4239 #define ABI_FUNC 'eph_phgamma'
4240 !End of the abilint section
4241 
4242  implicit none
4243 
4244 !Arguments ------------------------------------
4245 !scalars
4246  character(len=*),intent(in) :: wfk0_path
4247  integer,intent(in) :: comm
4248  type(datafiles_type),intent(in) :: dtfil
4249  type(dataset_type),intent(in) :: dtset
4250  type(crystal_t),intent(in) :: cryst
4251  type(ebands_t),intent(in) :: ebands
4252  type(dvdb_t),intent(inout) :: dvdb
4253  type(pawang_type),intent(in) :: pawang
4254  type(pseudopotential_type),intent(in) :: psps
4255  type(pawfgr_type),intent(in) :: pawfgr
4256  type(ifc_type),intent(in) :: ifc
4257  type(mpi_type),intent(in) :: mpi_enreg
4258 !arrays
4259  integer,intent(in) :: ngfft(18),ngfftf(18)
4260  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
4261  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
4262  type(ddk_t),intent(inout) :: ddk
4263 
4264 !Local variables ------------------------------
4265 !scalars
4266  integer,parameter :: dummy_npw=1,nsig=1,tim_getgh1c=1,berryopt0=0,timrev1=1
4267  integer,parameter :: useylmgr=0,useylmgr1=0,master=0,ndat1=1
4268  integer :: my_rank,nproc,iomode,mband,nsppol,nkpt,idir,ipert,iq_ibz
4269  integer :: cplex,db_iqpt,natom,natom3,ipc,ipc1,ipc2,nspinor,onpw
4270  integer :: bstart_k,bstart_kq,nband_k,nband_kq,ib1,ib2,band !band1,band2,
4271  integer :: ik_ibz,ik_bz,ikq_bz,ikq_ibz,isym_k,isym_kq,trev_k,trev_kq,timerev_q
4272  integer :: spin,istwf_k,istwf_kq,istwf_kirr,npw_k,npw_kq,npw_kirr
4273  integer :: ii,ipw,mpw,my_mpw,mnb,ierr,my_kstart,my_kstop,cnt,ncid
4274  integer :: isig,n1,n2,n3,n4,n5,n6,nspden,eph_scalprod,do_ftv1q
4275  integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnl1
4276  integer :: nfft,nfftf,mgfft,mgfftf,kqcount,nkpg,nkpg1,edos_intmeth
4277  integer :: iene, jene
4278 #ifdef HAVE_NETCDF
4279  integer :: ncerr
4280 #endif
4281  real(dp) :: cpu,wall,gflops
4282  real(dp) :: edos_step,edos_broad
4283  real(dp) :: ecut,eshift,eig0nk,dotr,doti,dksqmax
4284  logical,parameter :: have_ktimerev=.True.
4285  logical :: isirr_k,isirr_kq,gen_eigenpb
4286  type(wfd_t) :: wfd
4287  type(fstab_t),pointer :: fs
4288  type(gs_hamiltonian_type) :: gs_hamkq
4289  type(rf_hamiltonian_type) :: rf_hamkq
4290  type(edos_t) :: edos
4291  type(phgamma_t) :: gams
4292  type(a2fw_t) :: a2fw
4293  type(a2fw_tr_t) :: a2fw_tr
4294  character(len=500) :: msg
4295  character(len=fnlen) :: path
4296 !arrays
4297  integer :: g0_k(3),g0bz_kq(3),g0_kq(3),symq(4,2,cryst%nsym),dummy_gvec(3,dummy_npw)
4298  integer :: work_ngfft(18),gmax(3),my_gmax(3),gamma_ngqpt(3) !g0ibz_kq(3),
4299  integer,allocatable :: kg_k(:,:),kg_kq(:,:),gtmp(:,:),nband(:,:)
4300  integer :: indkk_kq(1,6)
4301  real(dp) :: kk(3),kq(3),kk_ibz(3),kq_ibz(3),qpt(3), phfrq(3*cryst%natom),lf(2),rg(2),res(2)
4302  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom),displ_red(2,3,cryst%natom,3*cryst%natom)
4303  real(dp) :: temp_tgam(2,3*cryst%natom,3*cryst%natom)
4304  real(dp) :: sigmas(2),wminmax(2)
4305  real(dp) :: n0(ebands%nsppol)
4306  real(dp),allocatable :: grad_berry(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),dkinpw(:)
4307  real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:),ph3d(:,:,:),ph3d1(:,:,:)
4308  real(dp),allocatable :: v1scf(:,:,:,:),tgam(:,:,:,:),gvals_qibz(:,:,:,:,:,:),gkk_atm(:,:,:,:)
4309  real(dp),allocatable :: bras_kq(:,:,:),kets_k(:,:,:),h1kets_kq(:,:,:)
4310  real(dp),allocatable :: ph1d(:,:),vlocal(:,:,:,:),vlocal1(:,:,:,:,:)
4311  real(dp),allocatable :: ylm_kq(:,:),ylm_k(:,:),ylmgr_kq(:,:,:)
4312  real(dp),allocatable :: dummy_vtrial(:,:),gvnl1(:,:),work(:,:,:,:)
4313  real(dp),allocatable ::  gs1c(:,:) !,gvnl_direc(:,:),pcon(:),sconjgr(:,:)
4314  !real(dp),allocatable :: eloc0_k(:),enl0_k(:),enl1_k(:),vlocal_tmp(:,:,:),vlocal1_tmp(:,:,:), rho1wfg(:,:),rho1wfr(:,:)
4315  real(dp),allocatable :: wt_k(:,:),wt_kq(:,:)
4316  real(dp),allocatable :: wt_k_en(:,:,:),wt_kq_en(:,:,:)
4317  logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
4318  type(fstab_t),target,allocatable :: fstab(:)
4319  type(pawcprj_type),allocatable  :: cwaveprj0(:,:) !natom,nspinor*usecprj)
4320  !real(dp),allocatable :: cwave0(:,:),gvnl1(:,:)
4321 
4322  real(dp), allocatable :: gvvvals_in_qibz(:,:,:,:,:,:,:)
4323  real(dp), allocatable :: gvvvals_out_qibz(:,:,:,:,:,:,:)
4324  real(dp), allocatable :: resvv_in(:,:)
4325  real(dp), allocatable :: resvv_out(:,:)
4326  real(dp), allocatable :: tgamvv_in(:,:,:,:,:),  vv_kk(:,:,:)
4327  real(dp), allocatable :: tgamvv_out(:,:,:,:,:), vv_kkq(:,:,:)
4328 
4329  ! for the Eliashberg solver
4330  integer :: ntemp
4331  real(dp) :: reltol,wcut
4332  real(dp) :: temp_range(2)
4333 
4334 !************************************************************************
4335 
4336  if (psps%usepaw == 1) then
4337    MSG_ERROR("PAW not implemented")
4338    ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
4339  end if
4340 
4341  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
4342 
4343  ! Copy important dimensions
4344  natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor; nspden = dtset%nspden
4345  nkpt = ebands%nkpt; mband=ebands%mband
4346 
4347  nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
4348  nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
4349  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
4350  n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
4351 
4352  ! Compute electron DOS.
4353  ! TODO: Optimize this part. Really slow if tetra and lots of points and/or bands
4354  ! Could just get DOS around efermi but then I cannot compute ef from IDOS.
4355  edos_intmeth = 2; if (dtset%prtdos == 1) edos_intmeth = 1
4356  !edos_intmeth = 1
4357  edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear
4358  edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha
4359  edos = ebands_get_edos(ebands,cryst,edos_intmeth,edos_step,edos_broad,comm)
4360 
4361  ! Store DOS per spin channels
4362  n0(:) = edos%gef(1:edos%nsppol)
4363  if (my_rank == master) then
4364    call edos_print(edos, unit=ab_out)
4365    path = strcat(dtfil%filnam_ds(4), "_EDOS")
4366    call wrtout(ab_out, sjoin("- Writing electron DOS to file:", path))
4367    call edos_write(edos, path)
4368  end if
4369 
4370  ! Find Fermi surface.
4371 
4372  ! array of widths: should this not imply nsig 2 in this whole routine?
4373  sigmas = [dtset%eph_fsmear, 2*dtset%eph_fsmear] !* eV_Ha
4374 
4375  ! Find Fermi surface k points
4376  ABI_DT_MALLOC(fstab, (nsppol))
4377  ! FIXME: kptopt, change setup of k-points if tetra: fist tetra weights then k-points on the Fermi surface.!
4378  call fstab_init(fstab, ebands, cryst, dtset%eph_fsewin, dtset%eph_intmeth, dtset%kptrlatt, &
4379    dtset%nshiftk, dtset%shiftk, comm)
4380  call fstab_print(fstab)
4381 
4382  ! now we can initialize the ddk velocities, on the FS grid only
4383  if (dtset%eph_transport > 0) then
4384    call ddk_read_fsvelocities(ddk, fstab, comm)
4385    call ddk_fs_average_veloc(ddk, ebands, fstab, sigmas)
4386  end if
4387 
4388  eph_scalprod = 0
4389  gamma_ngqpt = ifc%ngqpt
4390  ! TODO might this next condition not be any instead of all?
4391  if (all(dtset%eph_ngqpt_fine /= 0)) gamma_ngqpt = dtset%eph_ngqpt_fine
4392 
4393  ! TODO: Support nsig in phgamma_init
4394  call phgamma_init(gams,cryst,ifc,fstab(1),dtset%symdynmat,eph_scalprod,dtset%eph_transport,gamma_ngqpt,nsppol,nspinor,n0)
4395  call wrtout(std_out, sjoin("Will compute",itoa(gams%nqibz),"q-points in the IBZ"))
4396 
4397  ncid = nctk_noid
4398 #ifdef HAVE_NETCDF
4399  ! Open the netcdf file used to store the results of the calculation.
4400  if (my_rank == master) then
4401    NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), "_A2F.nc"), xmpi_comm_self))
4402    NCF_CHECK(crystal_ncwrite(cryst, ncid))
4403    NCF_CHECK(ebands_ncwrite(ebands, ncid))
4404    NCF_CHECK(edos_ncwrite(edos, ncid))
4405 
4406    ! Add eph dimensions.
4407    ncerr = nctk_def_dims(ncid, [nctkdim_t("nqibz", gams%nqibz), nctkdim_t("natom3", natom3)], defmode=.True.)
4408    NCF_CHECK(ncerr)
4409 
4410    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat"], defmode=.True.)
4411    NCF_CHECK(ncerr)
4412    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "ph_intmeth"], defmode=.True.)
4413    NCF_CHECK(ncerr)
4414    ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie"])
4415    NCF_CHECK(ncerr)
4416    ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ph_wstep", "ph_smear"])
4417    NCF_CHECK(ncerr)
4418 
4419    ! Define arrays for results.
4420    ncerr = nctk_def_arrays(ncid, [ &
4421      nctkarr_t("ngqpt", "int", "three"), &
4422      nctkarr_t("eph_ngqpt_fine", "int", "three"), &
4423      nctkarr_t("ddb_ngqpt", "int", "three"), &
4424      nctkarr_t("ph_ngqpt", "int", "three"), &
4425      ! linewidths in IBZ
4426      nctkarr_t('qibz', "dp", "number_of_reduced_dimensions, nqibz"), &
4427      nctkarr_t('wtq', "dp", "nqibz"), &
4428      nctkarr_t('phfreq_qibz', "dp", "natom3, nqibz"), &
4429      nctkarr_t('phdispl_cart_qibz', "dp", "two, natom3, natom3, nqibz"), &
4430      nctkarr_t('phgamma_qibz', "dp", "natom3, nqibz, number_of_spins"), &
4431      nctkarr_t('phlambda_qibz', "dp", "natom3, nqibz, number_of_spins") &
4432    ])
4433    NCF_CHECK(ncerr)
4434 
4435    ! ======================================================
4436    ! Write data that do not depend on the (kpt, spin) loop.
4437    ! ======================================================
4438    NCF_CHECK(nctk_set_datamode(ncid))
4439 
4440    ncerr = nctk_write_iscalars(ncid, &
4441        [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat", "ph_intmeth"], &
4442        [dtset%eph_intmeth, dtset%eph_transport, dtset%symdynmat, dtset%ph_intmeth])
4443    NCF_CHECK(ncerr)
4444    ncerr = nctk_write_dpscalars(ncid, &
4445      [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", "ph_wstep", "ph_smear"], &
4446      [dtset%eph_fsewin, dtset%eph_fsmear, dtset%eph_extrael, dtset%eph_fermie, dtset%ph_wstep, dtset%ph_smear])
4447    NCF_CHECK(ncerr)
4448 
4449    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qibz'), gams%qibz))
4450    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'wtq'), gams%wtq))
4451    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngqpt"), gamma_ngqpt))
4452    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_ngqpt_fine"), dtset%eph_ngqpt_fine))
4453    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ph_ngqpt"), dtset%ph_ngqpt))
4454    NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ddb_ngqpt"), dtset%ddb_ngqpt))
4455  end if
4456 #endif
4457  call edos_free(edos)
4458 
4459  ! Open the DVDB file
4460  call dvdb_open_read(dvdb, ngfftf, xmpi_comm_self)
4461 
4462  do_ftv1q = 0
4463  do iq_ibz=1,gams%nqibz
4464    qpt = gams%qibz(:,iq_ibz)
4465    if (dvdb_findq(dvdb, qpt) == -1) do_ftv1q = do_ftv1q + 1
4466 
4467    do spin=1,nsppol
4468      fs => fstab(spin)
4469      kqcount = 0
4470      do ik_bz=1,fs%nkfs
4471        kk = fs%kpts(:, ik_bz); kq = kk + qpt
4472        if (fstab_findkg0(fs, kq, g0bz_kq) == -1) cycle
4473        kqcount = kqcount + 1
4474      end do
4475      write(std_out,"((a,i0,2a,a,i0))")"For spin: ",spin,", qpt: ",trim(ktoa(qpt)),", number of (k,q) pairs: ",kqcount
4476    end do
4477  end do
4478  call wrtout(std_out, " ", do_flush=.True.)
4479 
4480  ! Activate Fourier interpolation if irred q-points are not in the DVDB file.
4481  if (do_ftv1q /= 0) then
4482    write(msg, "(2(a,i0),a)")"Will use Fourier interpolation of DFPT potentials [",do_ftv1q,"/",gams%nqibz,"]"
4483    call wrtout(std_out, msg)
4484    call wrtout(std_out, sjoin("From ngqpt", ltoa(ifc%ngqpt), "to", ltoa(gamma_ngqpt)))
4485    call dvdb_ftinterp_setup(dvdb,ifc%ngqpt,1,[zero,zero,zero],nfftf,ngfftf,comm)
4486  end if
4487 
4488  ! Initialize the wave function descriptor.
4489  ABI_MALLOC(nband, (nkpt, nsppol))
4490  ABI_MALLOC(bks_mask,(mband, nkpt, nsppol))
4491  ABI_MALLOC(keep_ur,(mband, nkpt ,nsppol))
4492  nband=mband; bks_mask=.False.; keep_ur=.False.
4493 
4494  ! Only wavefunctions on the FS are stored in wfd.
4495  ! Need all k-points on the FS because of k+q, spin is not distributed for the time being.
4496  ! One could reduce the memory allocated per MPI-rank via MPI-FFT or OpenMP...
4497  do spin=1,nsppol
4498    fs => fstab(spin)
4499    do ik_bz=1,fs%nkfs
4500      ik_ibz = fs%istg0(1, ik_bz)
4501      bstart_k = fs%bstcnt_ibz(1, ik_ibz); nband_k = fs%bstcnt_ibz(2, ik_ibz)
4502      bks_mask(bstart_k:bstart_k+nband_k-1, ik_ibz, spin) = .True.
4503    end do
4504  end do
4505 
4506  ! no memory distribution, each node has the full set of states.
4507  !bks_mask(1:mband,:,:) = .True.
4508 
4509  ecut = dtset%ecut ! dtset%dilatmx
4510  call wfd_init(wfd,cryst,pawtab,psps,keep_ur,dtset%paral_kgb,dummy_npw,mband,nband,nkpt,nsppol,bks_mask,&
4511    nspden,nspinor,dtset%ecutsm,dtset%dilatmx,ebands%istwfk,ebands%kptns,ngfft,&
4512    dummy_gvec,dtset%nloalg,dtset%prtvol,dtset%pawprtvol,comm,opt_ecut=ecut)
4513 
4514  call wfd_print(wfd,header="Wavefunctions on the Fermi Surface",mode_paral='PERS')
4515 
4516  ABI_FREE(nband)
4517  ABI_FREE(bks_mask)
4518  ABI_FREE(keep_ur)
4519 
4520  iomode = iomode_from_fname(wfk0_path)
4521 
4522  call wfd_read_wfk(wfd,wfk0_path,iomode)
4523 
4524  if (.False.) call wfd_test_ortho(wfd,cryst,pawtab,unit=std_out,mode_paral="PERS")
4525 
4526  ! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information on the coarse grid.
4527  ABI_MALLOC(ph1d, (2,3*(2*mgfft+1)*natom))
4528  call getph(cryst%atindx,natom,n1,n2,n3,ph1d,cryst%xred)
4529 
4530  ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
4531  ! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
4532  ! that will be used to symmetrize the wavefunctions in G-space.
4533  call cwtime(cpu,wall,gflops,"start")
4534  mpw = 0; gmax=0; cnt=0
4535  do iq_ibz=1,gams%nqibz
4536    qpt = gams%qibz(:,iq_ibz)
4537    do spin=1,nsppol
4538      fs => fstab(spin)
4539      do ik_bz=1,fs%nkfs
4540        cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
4541        kk = fs%kpts(:, ik_bz); kq = kk + qpt
4542        ! Do computation of G sphere, returning npw. Note istwfk==1.
4543        call get_kg(kk,1,ecut,cryst%gmet,onpw,gtmp)
4544        mpw = max(mpw, onpw)
4545        do ipw=1,onpw
4546          do ii=1,3
4547           gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw)))
4548          end do
4549        end do
4550        ABI_FREE(gtmp)
4551 
4552        ! TODO: g0 umklapp here can enter into play!
4553        ! fstab should contains the max of the umlapp G-vectors.
4554        ! gmax could not be large enough!
4555        call get_kg(kq,1,ecut,cryst%gmet,onpw,gtmp)
4556        mpw = max(mpw, onpw)
4557        do ipw=1,onpw
4558          do ii=1,3
4559           gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw)))
4560          end do
4561        end do
4562        ABI_FREE(gtmp)
4563      end do
4564    end do
4565  end do
4566  call cwtime(cpu,wall,gflops,"stop")
4567  write(msg,'(2(a,f8.2))')"gmax and mpw: ",cpu,", wall: ",wall
4568  call wrtout(std_out,msg,"COLL",do_flush=.True.)
4569 
4570  my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
4571  my_gmax = gmax; call xmpi_max(my_gmax, gmax, comm, ierr)
4572  call wrtout(std_out,sjoin('optimal value of mpw= ',itoa(mpw)),'COLL')
4573 
4574  ! Init work_ngfft
4575  gmax = gmax + 4 ! FIXME: this is to account for umklapp
4576  gmax = 2*gmax + 1
4577  call ngfft_seq(work_ngfft, gmax)
4578  write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
4579  ABI_STAT_MALLOC(work, (2, work_ngfft(4),work_ngfft(5),work_ngfft(6)), ierr)
4580  ABI_CHECK(ierr==0, 'out of memory in work')
4581 
4582  ! Allow PW-arrays dimensioned with mpw
4583  ABI_STAT_MALLOC(kg_k, (3, mpw), ierr)
4584  ABI_CHECK(ierr==0, 'out of memory in kg_k')
4585  ABI_STAT_MALLOC(kg_kq, (3, mpw), ierr)
4586  ABI_CHECK(ierr==0, 'out of memory in kg_kq')
4587 
4588 
4589  ! Spherical Harmonics for useylm==1.
4590  ABI_MALLOC(ylm_k,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
4591  ABI_MALLOC(ylm_kq,(mpw, psps%mpsang*psps%mpsang*psps%useylm))
4592  ABI_MALLOC(ylmgr_kq,(mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
4593 
4594  ! TODO FOR PAW
4595  usecprj = 0
4596  ABI_DT_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
4597 
4598 
4599  ! Prepare call to getgh1c
4600  usevnl = 0
4601  optlocal = 1  ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
4602  optnl = 2     ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
4603  opt_gvnl1 = 0 ! gvnl1 is output
4604  ABI_MALLOC(gvnl1, (2,usevnl))
4605  ABI_MALLOC(grad_berry, (2,nspinor*(berryopt0/4)))
4606 
4607  ! This part is taken from dfpt_vtorho
4608  !==== Initialize most of the Hamiltonian (and derivative) ====
4609  !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
4610  !2) Perform the setup needed for the non-local factors:
4611  !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
4612  !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
4613 
4614  call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,natom,&
4615 &  dtset%typat,cryst%xred,nfft,mgfft,ngfft,cryst%rprimd,dtset%nloalg,&
4616 &  comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
4617 &  usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
4618 
4619  !PAW:allocate memory for non-symetrized 1st-order occupancies matrix (pawrhoij1)
4620  ! pawrhoij1_unsym => pawrhoij1
4621  ! if (psps%usepaw==1.and.iscf_mod>0) then
4622  !   if (paral_atom) then
4623  !     ABI_DATATYPE_ALLOCATE(pawrhoij1_unsym,(natom))
4624  !     cplex_rhoij=max(cplex,dtset%pawcpxocc)
4625  !     nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
4626  !     call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,nspinor,&
4627  !       dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
4628  !   else
4629  !     pawrhoij1_unsym => pawrhoij1
4630  !     call pawrhoij_init_unpacked(pawrhoij1_unsym)
4631  !   end if
4632  ! end if
4633 
4634 ! Allocate vlocal. Note nvloc
4635  ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data)
4636  ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
4637  vlocal = huge(one)
4638 
4639  ! Allocate work space arrays.
4640  ABI_MALLOC(tgam, (2,natom3,natom3,nsig))
4641  ABI_CALLOC(dummy_vtrial, (nfftf,nspden))
4642 ! TODO: if we remove the nsig dependency we can remove this intermediate array
4643 ! and save a lot of memory
4644  ABI_STAT_MALLOC(gvals_qibz, (2,natom3,natom3,nsig,gams%nqibz,nsppol), ierr)
4645  ABI_CHECK(ierr==0, 'out of memory in gvals_qibz')
4646 
4647  if (dtset%eph_transport > 0) then
4648    ABI_MALLOC(tgamvv_in, (2,gams%ndir_transp**2,natom3,natom3,nsig))
4649    ABI_MALLOC(tgamvv_out, (2,gams%ndir_transp**2,natom3,natom3,nsig))
4650    ABI_MALLOC(resvv_in, (2,gams%ndir_transp**2))
4651    ABI_MALLOC(resvv_out, (2,gams%ndir_transp**2))
4652    ! TODO: if we remove the nsig dependency we can remove this intermediate array
4653    ! and save a lot of memory
4654    ABI_STAT_MALLOC(gvvvals_in_qibz, (2,gams%ndir_transp**2,natom3,natom3,nsig,gams%nqibz,nsppol), ierr)
4655    ABI_CHECK(ierr==0, 'out of memory in gvvvals_in_qibz')
4656    ABI_STAT_MALLOC(gvvvals_out_qibz, (2,gams%ndir_transp**2,natom3,natom3,nsig,gams%nqibz,nsppol), ierr)
4657    ABI_CHECK(ierr==0, 'out of memory in gvvvals_out_qibz')
4658  end if
4659 
4660 #ifdef DEV_MJV
4661  open (unit=800, file="wt_kq_en.dat")
4662  open (unit=801, file="wt_k_en.dat")
4663  open (unit=802, file="res_small.dat")
4664 #endif
4665  do iq_ibz=1,gams%nqibz
4666    qpt = gams%qibz(:,iq_ibz)
4667    tgam = zero
4668    if (dtset%eph_transport > 0) then
4669      tgamvv_in = zero
4670      tgamvv_out = zero
4671    end if
4672 
4673    call cwtime(cpu,wall,gflops,"start")
4674 
4675    ! Find the index of the q-point in the DVDB.
4676    db_iqpt = dvdb_findq(dvdb, qpt)
4677 
4678    if (db_iqpt /= -1) then
4679      if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Found: ",ktoa(qpt)," in DVDB with index ",itoa(db_iqpt)))
4680      ! Read or reconstruct the dvscf potentials for all 3*natom perturbations.
4681      ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
4682      call dvdb_readsym_allv1(dvdb, db_iqpt, cplex, nfftf, ngfftf, v1scf, comm)
4683    else
4684      if (dtset%prtvol > 0) call wrtout(std_out, sjoin("Could not find: ",ktoa(qpt), "in DVDB - interpolating"))
4685      ! Fourier interpolate of the potential
4686      ABI_CHECK(any(abs(qpt) > tol12), "qpt cannot be zero if Fourier interpolation is used")
4687      cplex = 2
4688      ABI_MALLOC(v1scf, (cplex,nfftf,nspden,natom3))
4689      call dvdb_ftinterp_qpt(dvdb, qpt, nfftf, ngfftf, v1scf, comm)
4690    end if
4691 
4692    ! Examine the symmetries of the q wavevector
4693    call littlegroup_q(cryst%nsym,qpt,symq,cryst%symrec,cryst%symafm,timerev_q,prtvol=dtset%prtvol)
4694 
4695    ! Allocate vlocal1 with correct cplex. Note nvloc
4696    ABI_STAT_MALLOC(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc,natom3), ierr)
4697    ABI_CHECK(ierr==0, "oom vlocal1")
4698 
4699 #ifdef DEV_MJV
4700    gams%vals_ee = zero
4701 #endif
4702    do spin=1,nsppol
4703      fs => fstab(spin)
4704 
4705      ! Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin.
4706      do ipc=1,natom3
4707        call rf_transgrid_and_pack(spin,nspden,psps%usepaw,cplex,nfftf,nfft,ngfft,gs_hamkq%nvloc,&
4708                  pawfgr,mpi_enreg,dummy_vtrial,v1scf(:,:,:,ipc),vlocal,vlocal1(:,:,:,:,ipc))
4709      end do
4710 
4711      ! Continue to initialize the Hamiltonian
4712      call load_spin_hamiltonian(gs_hamkq,spin,vlocal=vlocal,with_nonlocal=.true.)
4713 
4714      ! Allocate workspace for wavefunctions. Make npw larger than expected.
4715      ! maxnb is the maximum number of bands crossing the FS, used to dimension arrays.
4716      mnb = fs%maxnb
4717      ABI_CHECK(mnb <= ebands%mband, "mnb > ebands%mband")
4718      ABI_MALLOC(bras_kq, (2, mpw*nspinor, mnb))
4719      ABI_MALLOC(kets_k, (2, mpw*nspinor, mnb))
4720      ABI_MALLOC(h1kets_kq, (2, mpw*nspinor, mnb))
4721      ABI_MALLOC(gkk_atm, (2, mnb, mnb, natom3))
4722 
4723      ! The weights for FS integration.
4724      ABI_CALLOC(wt_k, (nsig, mnb))
4725      ABI_CALLOC(wt_kq, (nsig, mnb))
4726 
4727 !TODO: add flag around this
4728      ABI_CALLOC(wt_k_en, (nsig, mnb, gams%nene))
4729      ABI_CALLOC(wt_kq_en, (nsig, mnb, gams%nene))
4730 
4731      if (dtset%eph_transport > 0) then
4732        ABI_CALLOC(vv_kk, (gams%ndir_transp**2, mnb, mnb))
4733        ABI_CALLOC(vv_kkq, (gams%ndir_transp**2, mnb, mnb))
4734      end if
4735 
4736      ! =========================
4737      ! Integration over FS(spin)
4738      ! =========================
4739      call xmpi_split_work(fs%nkfs,comm,my_kstart,my_kstop,msg,ierr)
4740 
4741      do ik_bz=my_kstart,my_kstop
4742        ! The k-point and the symmetries relating the BZ points to the IBZ.
4743        kk = fs%kpts(:, ik_bz)
4744        ik_ibz = fs%istg0(1, ik_bz); isym_k = fs%istg0(2, ik_bz)
4745        trev_k = fs%istg0(3, ik_bz); g0_k = fs%istg0(4:6,ik_bz)
4746        isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
4747        kk_ibz = ebands%kptns(:,ik_ibz)
4748 
4749        ! Number of bands crossing the fermi level at k
4750        bstart_k = fs%bstcnt_ibz(1, ik_ibz); nband_k = fs%bstcnt_ibz(2, ik_ibz)
4751 
4752        ! Find k+q in the extended zone and extract symmetry info. cycle if k+q not in FS.
4753        ! Be careful here because there are two umklapp vectors to be considered:
4754        !
4755        !   k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz
4756        !
4757        kq = kk + qpt
4758        ikq_bz = fstab_findkg0(fs, kq, g0bz_kq); if (ikq_bz == -1) cycle
4759 
4760        !ikq_ibz = fs%istg0(1, ikq_bz); isym_kq = fs%istg0(2, ikq_bz)
4761        !trev_kq = fs%istg0(3, ikq_bz); g0ibz_kq = fs%istg0(4:6,ikq_bz)
4762        !isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0ibz_kq == 0))
4763        !g0_kq =  g0ibz_kq + g0bz_kq
4764 
4765        call listkk(dksqmax,cryst%gmet,indkk_kq,ebands%kptns,kq,ebands%nkpt,1,cryst%nsym,&
4766           1,cryst%symafm,cryst%symrel,timrev1,use_symrec=.False.)
4767 
4768        if (dksqmax > tol12) then
4769          write(msg, '(3a,es16.6,7a)' )&
4770           "The WFK file cannot be used to compute phonon linewidths.",ch10,&
4771           "At least one of the k-points could not be generated from a symmetrical one. dksqmax: ",dksqmax,ch10,&
4772           "Q-mesh: ",ltoa(gamma_ngqpt),", K-mesh (from kptrlatt) ",ltoa(get_diag(dtset%kptrlatt)), &
4773           'Action: check your WFK file and (k,q) point input variables'
4774           MSG_ERROR(msg)
4775        end if
4776 
4777        ikq_ibz = indkk_kq(1,1); isym_kq = indkk_kq(1,2)
4778        trev_kq = indkk_kq(1, 6); g0_kq = indkk_kq(1, 3:5)
4779        isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0))
4780        kq_ibz = ebands%kptns(:,ikq_ibz)
4781 
4782        ! Number of bands crossing the fermi level at k + q
4783        bstart_kq = fs%bstcnt_ibz(1, ikq_ibz); nband_kq = fs%bstcnt_ibz(2, ikq_ibz)
4784        ABI_CHECK(nband_k <= mnb .and. nband_kq <= mnb, "wrong nband")
4785 
4786        ! Get npw_k, kg_k and symmetrize wavefunctions from IBZ (if needed).
4787        ! Be careful with time-reversal symmetry.
4788        if (isirr_k) then
4789          ! Copy u_k(G)
4790          istwf_k = wfd%istwfk(ik_ibz); npw_k = wfd%npwarr(ik_ibz)
4791          ABI_CHECK(mpw >= npw_k, "mpw < npw_k")
4792          kg_k(:,1:npw_k) = wfd%kdata(ik_ibz)%kg_k
4793          do ib2=1,nband_k
4794            band = ib2 + bstart_k - 1
4795            call wfd_copy_cg(wfd, band, ik_ibz, spin, kets_k(1,1,ib2))
4796          end do
4797        else
4798          ! Reconstruct u_k(G) from the IBZ image.
4799          !istwf_k = set_istwfk(kk); if (.not. have_ktimerev) istwf_k = 1
4800          !call change_istwfk(from_npw,from_kg,from_istwfk,to_npw,to_kg,to_istwfk,n1,n2,n3,ndat,from_cg,to_cg)
4801          istwf_k = 1
4802          call get_kg(kk,istwf_k,ecut,cryst%gmet,npw_k,gtmp)
4803          ABI_CHECK(mpw >= npw_k, "mpw < npw_k")
4804          kg_k(:,1:npw_k) = gtmp(:,:npw_k)
4805          ABI_FREE(gtmp)
4806 
4807          ! Use h1kets_kq as workspace array, results stored in kets_k.
4808          istwf_kirr = wfd%istwfk(ik_ibz); npw_kirr = wfd%npwarr(ik_ibz)
4809          do ib2=1,nband_k
4810            band = ib2 + bstart_k - 1
4811            call wfd_copy_cg(wfd, band, ik_ibz, spin, h1kets_kq)
4812            call cgtk_rotate(cryst, kk_ibz, isym_k, trev_k, g0_k, nspinor, ndat1,&
4813                             npw_kirr, wfd%kdata(ik_ibz)%kg_k,&
4814                             npw_k, kg_k, istwf_kirr, istwf_k, h1kets_kq, kets_k(:,:,ib2), work_ngfft, work)
4815          end do
4816        end if
4817 
4818        ! Get npw_kq, kg_kq and symmetrize wavefunctions from IBZ (if needed).
4819        ! Be careful with time-reversal symmetry.
4820        if (isirr_kq) then
4821          ! Copy u_kq(G)
4822          istwf_kq = wfd%istwfk(ikq_ibz); npw_kq = wfd%npwarr(ikq_ibz)
4823          ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
4824          kg_kq(:,1:npw_kq) = wfd%kdata(ikq_ibz)%kg_k
4825          do ib1=1,nband_kq
4826            band = ib1 + bstart_kq - 1
4827            call wfd_copy_cg(wfd, band, ikq_ibz, spin, bras_kq(1,1,ib1))
4828          end do
4829        else
4830          ! Reconstruct u_kq(G) from the IBZ image.
4831          !istwf_kq = set_istwfk(kq); if (.not. have_ktimerev) istwf_kq = 1
4832          !call change_istwfk(from_npw,from_kg,from_istwfk,to_npw,to_kg,to_istwfk,n1,n2,n3,ndat,from_cg,to_cg)
4833          istwf_kq = 1
4834          call get_kg(kq,istwf_kq,ecut,cryst%gmet,npw_kq,gtmp)
4835          ABI_CHECK(mpw >= npw_kq, "mpw < npw_kq")
4836          kg_kq(:,1:npw_kq) = gtmp(:,:npw_kq)
4837          ABI_FREE(gtmp)
4838 
4839          ! Use h1kets_kq as workspace array, results stored in bras_kq
4840          istwf_kirr = wfd%istwfk(ikq_ibz); npw_kirr = wfd%npwarr(ikq_ibz)
4841          !g0_kq =  g0ibz_kq + g0bz_kq
4842          do ib1=1,nband_kq
4843            band = ib1 + bstart_kq - 1
4844            call wfd_copy_cg(wfd, band, ikq_ibz, spin, h1kets_kq)
4845            call cgtk_rotate(cryst, kq_ibz, isym_kq, trev_kq, g0_kq, nspinor, ndat1,&
4846                             npw_kirr, wfd%kdata(ikq_ibz)%kg_k,&
4847                             npw_kq, kg_kq, istwf_kirr, istwf_kq, h1kets_kq, bras_kq(:,:,ib1), work_ngfft, work)
4848          end do
4849        end if
4850 
4851        ! if PAW, one has to solve a generalized eigenproblem
4852        ! BE careful here because I will need sij_opt==-1
4853        gen_eigenpb = (psps%usepaw==1)
4854        sij_opt = 0; if (gen_eigenpb) sij_opt = 1
4855        ABI_MALLOC(gs1c, (2,npw_kq*nspinor*((sij_opt+1)/2)))
4856 
4857        ! Set up the spherical harmonics (Ylm) at k and k+q. See also dfpt_looppert
4858        !if (psps%useylm==1) then
4859        !   optder=0; if (useylmgr==1) optder=1
4860        !   call initylmg(cryst%gprimd,kg_k,kk,mkmem1,mpi_enreg,psps%mpsang,mpw,nband,mkmem1,&
4861        !     [npw_k],dtset%nsppol,optder,cryst%rprimd,ylm_k,ylmgr)
4862        !   call initylmg(cryst%gprimd,kg_kq,kq,mkmem1,mpi_enreg,psps%mpsang,mpw,nband,mkmem1,&
4863        !     [npw_kq],dtset%nsppol,optder,cryst%rprimd,ylm_kq,ylmgr_kq)
4864        !end if
4865 
4866        ! Loop over all 3*natom perturbations.
4867        do ipc=1,natom3
4868          idir = mod(ipc-1, 3) + 1; ipert = (ipc - idir) / 3 + 1
4869 
4870          ! Prepare application of the NL part.
4871          call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.)
4872              !&paw_ij1=paw_ij1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
4873              !&my_spintab=mpi_enreg%my_isppoltab)
4874          call load_spin_rf_hamiltonian(rf_hamkq,spin,vlocal1=vlocal1(:,:,:,:,ipc),with_nonlocal=.true.)
4875 
4876          ! This call is not optimal because there are quantities in out that do not depend on idir,ipert
4877          call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kk,kq,idir,ipert,&   ! In
4878            cryst%natom,cryst%rmet,cryst%gprimd,cryst%gmet,istwf_k,&           ! In
4879            npw_k,npw_kq,useylmgr1,kg_k,ylm_k,kg_kq,ylm_kq,ylmgr_kq,&          ! In
4880            dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1)      ! Out
4881 
4882          ! Calculate dvscf * psi_k, results stored in h1kets_kq on the k+q sphere.
4883          ! Compute H(1) applied to GS wavefunction Psi(0)
4884          do ib2=1,nband_k
4885            band = ib2 + bstart_k - 1
4886            eig0nk = ebands%eig(band,ik_ibz,spin)
4887            ! Use scissor shift on 0-order eigenvalue
4888            eshift = eig0nk - dtset%dfpt_sciss
4889 
4890            call getgh1c(berryopt0,kets_k(:,:,ib2),cwaveprj0,h1kets_kq(:,:,ib2),&
4891 &                       grad_berry,gs1c,gs_hamkq,gvnl1,idir,ipert,eshift,mpi_enreg,optlocal,&
4892 &                       optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
4893          end do
4894 
4895          call destroy_rf_hamiltonian(rf_hamkq)
4896 
4897          ABI_FREE(kinpw1)
4898          ABI_FREE(kpg1_k)
4899          ABI_FREE(kpg_k)
4900          ABI_FREE(dkinpw)
4901          ABI_FREE(ffnlk)
4902          ABI_FREE(ffnl1)
4903          ABI_FREE(ph3d)
4904          if (allocated(ph3d1)) then
4905            ABI_FREE(ph3d1)
4906          end if
4907 
4908          ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation.
4909          !The array eig1_k contains:
4910          !
4911          ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)>                           (NC psps)
4912          ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW)
4913          gkk_atm(:,:,:,ipc) = zero
4914          do ib2=1,nband_k
4915            do ib1=1,nband_kq
4916              call dotprod_g(dotr,doti,istwf_kq,npw_kq*nspinor,2,bras_kq(1,1,ib1),h1kets_kq(1,1,ib2),&
4917                mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
4918              gkk_atm(:,ib1,ib2,ipc) = [dotr, doti]
4919            end do
4920          end do
4921        end do ! ipc (loop over 3*natom atomic perturbations)
4922 
4923        ABI_FREE(gs1c)
4924 
4925        ! Sum over bands. Here weights come into play.
4926        ! Compute weights for FS integration.
4927        call fstab_weights_ibz(fs, ebands, ik_ibz, spin, sigmas, wt_k)
4928        call fstab_weights_ibz(fs, ebands, ikq_ibz, spin, sigmas, wt_kq)
4929 
4930 
4931        ! Accumulate results in tgam (sum over FS and bands).
4932        do ipc2=1,natom3
4933          do ipc1=1,natom3
4934            do ib2=1,nband_k
4935              do ib1=1,nband_kq
4936                lf = gkk_atm(:, ib1, ib2, ipc1)
4937                rg = gkk_atm(:, ib1, ib2, ipc2)
4938                res(1) = lf(1) * rg(1) + lf(2) * rg(2)
4939                res(2) = lf(1) * rg(2) - lf(2) * rg(1)
4940                ! Loop over smearing values.
4941                do isig=1,nsig
4942                  tgam(:,ipc1,ipc2,isig) = tgam(:,ipc1,ipc2,isig) &
4943 &                  + res(:)     * wt_kq(isig, ib1) * wt_k(isig, ib2)
4944                  !write(std_out,*)res, wt_kq(isig, ib,  wt_k(isig, ib2)
4945                end do
4946              end do
4947            end do
4948          end do
4949        end do
4950 
4951 !TODO: put a flag around this, in case it takes a lot of time
4952        do jene = 1, gams%nene
4953          call fstab_weights_ibz(fs, ebands, ik_ibz, spin, sigmas, wt_k_en(:,:,jene), iene=jene)
4954          call fstab_weights_ibz(fs, ebands, ikq_ibz, spin, sigmas, wt_kq_en(:,:,jene), iene=jene)
4955        end do
4956 
4957 #ifdef DEV_MJV
4958 write (800,*) wt_kq_en
4959 write (801,*) wt_k_en
4960        do ib2=1,nband_k
4961          do ib1=1,nband_kq
4962            do ipc2=1,natom3
4963              do ipc1=1,natom3
4964                lf = gkk_atm(:, ib1, ib2, ipc1)
4965                rg = gkk_atm(:, ib1, ib2, ipc2)
4966                res(1) = lf(1) * rg(1) + lf(2) * rg(2)
4967                res(2) = lf(1) * rg(2) - lf(2) * rg(1)
4968 if (sum(abs(res)) < tol6) then
4969   write (802,*) 'res small ib2,ib1,ipc2,ipc1, res ', ib2,ib1,ipc2,ipc1, res
4970 end if
4971                do jene = 1, gams%nene
4972                  do iene = 1, gams%nene
4973                    gams%vals_ee(:,iene,jene,ipc1,ipc2,iq_ibz,spin) = &
4974 &                    gams%vals_ee(:,iene,jene,ipc1,ipc2,iq_ibz,spin) + &
4975 &                    res(:) * wt_kq_en(1, ib1, iene) * wt_k_en(1, ib2, jene)
4976                  end do
4977                end do
4978              end do
4979            end do
4980          end do
4981        end do
4982 #endif
4983 
4984        if (dtset%eph_transport > 0) then
4985          ! TODO : could almost make this a BLAS call plus a reshape...
4986          do ib2 = 1,nband_k
4987            do ipc2 = 1,gams%ndir_transp
4988              do ib1 = 1,nband_kq
4989                do ipc1 = 1,gams%ndir_transp
4990                  vv_kk(ipc1+(ipc2-1)*gams%ndir_transp, ib1,ib2)  = ddk%velocity(ipc1,ib1,ik_bz,spin) &
4991 &                   * ddk%velocity(ipc2,ib2,ik_bz,spin) ! vk vk
4992                  vv_kkq(ipc1+(ipc2-1)*gams%ndir_transp, ib1,ib2) = ddk%velocity(ipc1,ib1,ikq_bz,spin) &
4993 &                   * ddk%velocity(ipc2,ib2,ik_bz,spin) ! vk vk+q
4994                end do
4995              end do
4996            end do
4997          end do
4998          ! Accumulate results in tgam (sum over FS and bands).
4999          do ipc2=1,natom3
5000            do ipc1=1,natom3
5001               do ib2=1,nband_k
5002                 do ib1=1,nband_kq
5003                   lf = gkk_atm(:, ib1, ib2, ipc1)
5004                   rg = gkk_atm(:, ib1, ib2, ipc2)
5005                   resvv_in(1,:) = res(1) * vv_kkq(:,ib1,ib2)
5006                   resvv_in(2,:) = res(2) * vv_kkq(:,ib1,ib2)
5007                   resvv_out(1,:) = res(1) * vv_kk(:,ib1,ib2)
5008                   resvv_out(2,:) = res(2) * vv_kk(:,ib1,ib2)
5009                   ! Loop over smearing values.
5010                   do isig=1,nsig
5011                     tgamvv_in(:,:,ipc1,ipc2,isig)  = tgamvv_in(:,:,ipc1,ipc2,isig)  &
5012 &                     + resvv_in(:,:)  * wt_kq(isig, ib1) * wt_k(isig, ib2)
5013                     tgamvv_out(:,:,ipc1,ipc2,isig) = tgamvv_out(:,:,ipc1,ipc2,isig) &
5014 &                     + resvv_out(:,:) * wt_kq(isig, ib1) * wt_k(isig, ib2)
5015                     !write(std_out,*)res, wt_kq(isig, ib,  wt_k(isig, ib2)
5016                   end do
5017                 end do
5018               end do
5019            end do
5020          end do
5021 
5022        end if ! add transport things
5023 
5024      end do ! ikfs
5025 
5026      ABI_FREE(wt_k)
5027      ABI_FREE(wt_kq)
5028      ABI_FREE(bras_kq)
5029      ABI_FREE(kets_k)
5030      ABI_FREE(h1kets_kq)
5031      ABI_FREE(gkk_atm)
5032 
5033 !TODO: add flag around this deallocation
5034      ABI_FREE(wt_k_en)
5035      ABI_FREE(wt_kq_en)
5036 
5037      call xmpi_sum(tgam, comm, ierr)
5038 #ifdef DEV_MJV
5039      call xmpi_sum(gams%vals_ee, comm, ierr)
5040 #endif
5041 
5042      if (dtset%eph_transport > 0) then
5043        ABI_FREE(vv_kk)
5044        ABI_FREE(vv_kkq)
5045 
5046        call xmpi_sum(tgamvv_in, comm, ierr)
5047        call xmpi_sum(tgamvv_out, comm, ierr)
5048      end if ! add transport things
5049 
5050      if (eph_scalprod == 1) then
5051        ! Get phonon frequencies and displacements for this q-point
5052        ! used in scalar product with H(1)_atom,idir  matrix elements
5053        call ifc_fourq(ifc, cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
5054 
5055        !call ifc_diagoq(ifc,cryst,qpt,phfrq,displ_cart,nanaqdir)
5056      end if
5057 
5058      do isig=1,nsig
5059        if (eph_scalprod == 1) then
5060          ! TODO: NotTested: moreover one should make sure that the phase of the eigenvectors is deterministic.
5061          ! Multiply by displacement matrices. Results are returned in temp_tgam.
5062          call gam_mult_displ(natom3, displ_red, tgam(:,:,:,isig), temp_tgam)
5063          tgam(:,:,:,isig) = temp_tgam
5064        end if
5065 
5066        ! Save results for this (q-point, spin)
5067        !write(std_out,*)tgam(:,:,:,isig)
5068        gvals_qibz(:,:,:,isig,iq_ibz,spin) = tgam(:,:,:,isig)
5069        if (dtset%eph_transport > 0) then
5070          gvvvals_in_qibz(:,:,:,:,isig,iq_ibz,spin)  = tgamvv_in(:,:,:,:,isig)
5071          gvvvals_out_qibz(:,:,:,:,isig,iq_ibz,spin) = tgamvv_out(:,:,:,:,isig)
5072        end if
5073      end do ! isig
5074 
5075    end do ! spin sppol
5076 
5077    ABI_FREE(v1scf)
5078    ABI_FREE(vlocal1)
5079 
5080    call cwtime(cpu,wall,gflops,"stop")
5081    write(msg,'(2(a,i0),2(a,f8.2))')"q-point [",iq_ibz,"/",gams%nqibz,"] completed. cpu:",cpu,", wall:",wall
5082    call wrtout(std_out, msg, do_flush=.True.)
5083  end do ! iq_ibz
5084 
5085 #ifdef DEV_MJV
5086 close(800)
5087 close(801)
5088 close(802)
5089 #endif
5090 
5091  ! Collect gvals_qibz on each node and divide by the total number of k-points in the full mesh.
5092  do spin=1,nsppol
5093    !call xmpi_sum(gvals_qibz, comm_qpts, ierr) ! for the moment only split over k??
5094    gvals_qibz(:,:,:,:,:,spin) = gvals_qibz(:,:,:,:,:,spin) / fstab(spin)%nktot
5095 #ifdef DEV_MJV
5096    gams%vals_ee(:,:,:,:,:,:,spin) = gams%vals_ee(:,:,:,:,:,:,spin) / fstab(spin)%nktot
5097 #endif
5098 
5099    if (dtset%eph_transport > 0) then
5100      gvvvals_in_qibz(:,:,:,:,:,:,spin) = gvvvals_in_qibz(:,:,:,:,:,:,spin) / fstab(spin)%nktot
5101      gvvvals_out_qibz(:,:,:,:,:,:,spin) = gvvvals_out_qibz(:,:,:,:,:,:,spin) / fstab(spin)%nktot
5102    end if
5103  end do
5104  call wrtout(std_out, "Computation of tgamma matrices completed", "COLL", do_flush=.True.)
5105 
5106  ! Free memory
5107  ABI_FREE(gvnl1)
5108  ABI_FREE(grad_berry)
5109  ABI_FREE(dummy_vtrial)
5110  ABI_FREE(work)
5111  ABI_FREE(ph1d)
5112  ABI_FREE(vlocal)
5113  ABI_FREE(kg_k)
5114  ABI_FREE(kg_kq)
5115  ABI_FREE(ylm_k)
5116  ABI_FREE(ylm_kq)
5117  ABI_FREE(ylmgr_kq)
5118  ABI_FREE(tgam)
5119  if (dtset%eph_transport > 0) then
5120    ABI_FREE(tgamvv_in)
5121    ABI_FREE(tgamvv_out)
5122    ABI_FREE(resvv_in)
5123    ABI_FREE(resvv_out)
5124  end if
5125 
5126  call destroy_hamiltonian(gs_hamkq)
5127  call wfd_free(wfd)
5128  do spin=1,ebands%nsppol
5129    call fstab_free(fstab(spin))
5130  end do
5131  ABI_DT_FREE(fstab)
5132 
5133  call pawcprj_free(cwaveprj0)
5134  ABI_DT_FREE(cwaveprj0)
5135 
5136  ! Initialize object used to interpolate linewidths, compute a2f, write results etc.
5137  ! TODO: also save values for isig > 1???
5138  gams%vals_qibz = gvals_qibz(:,:,:,1,:,:)
5139  !write(std_out,*)gvals_qibz
5140  ABI_FREE(gvals_qibz)
5141 
5142  if (dtset%eph_transport > 0) then
5143    gams%vals_in_qibz  = gvvvals_in_qibz(:,:,:,:,1,:,:)
5144    gams%vals_out_qibz = gvvvals_out_qibz(:,:,:,:,1,:,:)
5145    ABI_FREE(gvvvals_in_qibz)
5146    ABI_FREE(gvvvals_out_qibz)
5147  end if
5148 
5149  ! This call is not executed in elphon!
5150  if (gams%symgamma == 1) then
5151    do spin=1,gams%nsppol
5152      do iq_ibz=1,gams%nqibz
5153        call tgamma_symm(cryst,gams%qibz(:,iq_ibz),gams%vals_qibz(:,:,:,iq_ibz,spin))
5154      end do
5155    end do
5156  end if
5157 
5158  ! Print gamma(IBZ) to ab_out and ncid
5159  if (my_rank == master) call phgamma_print(gams, cryst, ifc, ncid)
5160 
5161  ! TODO: this should be the end of a subroutine to initialize the gams object, and the rest below called case by case,
5162  ! to modularize things further. The calculations of matrix elements above should also be encapsulated to avoid
5163  ! code duplication with the other eph sub-drivers
5164 
5165  ! Interpolate linewidths along the q-path.
5166  if (dtset%ph_nqpath <= 0) then
5167    write(msg, '(7a,es16.6,4a)' )&
5168     'You have not specified a path for the linewidth calculation - no interpolation or output will be done ',ch10,&
5169     'Action: check your input variables ph_nqpath and ph_qpath'
5170    MSG_ERROR(msg)
5171  end if
5172  call phgamma_linwid(gams,cryst,ifc,dtset%ph_ndivsm,dtset%ph_nqpath,dtset%ph_qpath,dtfil%filnam_ds(4),ncid,wminmax,comm)
5173 
5174  ! Compute a2Fw using ab-initio q-points (no interpolation)
5175  call a2fw_init(a2fw,gams,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,wminmax,dtset%ph_smear,&
5176    dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,comm,qintp=.False.,qptopt=1)
5177  if (my_rank == master) call a2fw_write(a2fw, strcat(dtfil%filnam_ds(4), "_NOINTP"), "_qcoarse", ncid)
5178 #ifdef DEV_MJV
5179  if (my_rank == master) call a2fw_ee_write(a2fw, strcat(dtfil%filnam_ds(4), "_NOINTP"))
5180 #endif
5181  call a2fw_free(a2fw)
5182 
5183  ! Compute a2Fw using Fourier interpolation.
5184  call a2fw_init(a2fw,gams,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,wminmax,dtset%ph_smear,&
5185    dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,comm,qptopt=1)
5186  if (my_rank == master) call a2fw_write(a2fw, dtfil%filnam_ds(4), "_qintp", ncid)
5187 #ifdef DEV_MJV
5188  if (my_rank == master) call a2fw_ee_write(a2fw, dtfil%filnam_ds(4))
5189 #endif
5190 
5191  ! TODO: Use KT mesh instead of T but read T from input.
5192  ntemp = 6
5193  temp_range = [0.6_dp, 1.2_dp]
5194  wcut = 10 * wminmax(2); reltol = 0.001
5195  !call a2fw_solve_gap(a2fw,cryst,dtset%tmesh,wcut,dtset%eph_mustar,dtset%nstep,reltol,dtfil%filnam_ds(4),comm)
5196  call a2fw_free(a2fw)
5197 
5198  ! Compute A2fw using Fourier interpolation and full BZ for debugging purposes.
5199  call a2fw_init(a2fw,gams,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,wminmax,dtset%ph_smear,&
5200    dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,comm,qptopt=3)
5201  if (my_rank == master) call a2fw_write(a2fw, strcat(dtfil%filnam_ds(4), "_A2FW_QPTOPT3"), "fake", nctk_noid)
5202  call a2fw_free(a2fw)
5203 
5204  if (dtset%eph_transport == 1) then
5205    ! Compute a2Fw_tr using ab-initio q-points (no interpolation)
5206    call a2fw_tr_init(a2fw_tr,gams,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,wminmax,dtset%ph_smear,&
5207      dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,comm,qintp=.False.,qptopt=1)
5208    if (my_rank == master) call a2fw_tr_write(a2fw_tr, strcat(dtfil%filnam_ds(4), "_NOINTP"), "_qcoarse", ncid)
5209    call a2fw_tr_free(a2fw_tr)
5210 
5211    ! Compute a2Fw_tr using Fourier interpolation.
5212    call a2fw_tr_init(a2fw_tr,gams,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,wminmax,dtset%ph_smear,&
5213      dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,comm,qptopt=1)
5214    if (my_rank == master) call a2fw_tr_write(a2fw_tr, dtfil%filnam_ds(4), "_qintp", ncid)
5215 
5216    ! calculate and output transport quantities
5217    call a2fw_tr_free(a2fw_tr)
5218 
5219    ! Compute A2fw_tr using Fourier interpolation and full BZ for debugging purposes.
5220    call a2fw_tr_init(a2fw_tr,gams,cryst,ifc,dtset%ph_intmeth,dtset%ph_wstep,wminmax,dtset%ph_smear,&
5221      dtset%ph_ngqpt,dtset%ph_nqshift,dtset%ph_qshift,comm,qptopt=3)
5222    if (my_rank == master) call a2fw_tr_write(a2fw_tr, strcat(dtfil%filnam_ds(4), "_A2FWTR_QPTOPT3"), "fake", nctk_noid)
5223 
5224    call a2fw_tr_free(a2fw_tr)
5225  end if
5226 
5227 #ifdef HAVE_NETCDF
5228  if (my_rank == master) then
5229    NCF_CHECK(nf90_close(ncid))
5230  end if
5231 #endif
5232 
5233  call phgamma_free(gams)
5234 
5235 end subroutine eph_phgamma

m_phgamma/phgamma_eval_qibz [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_eval_qibz

FUNCTION

 Compute the phonon linewidths for q-points in the IBZ without performing interpolation.

INPUTS

  gams<phgamma_t>
  cryst<crystal_t>=Crystal structure.
  ifc<ifc_type>=Interatomic force constants.
  iq_ibz=Index of the q-point in the IBZ array.
  spin=Spin index

OUTPUT

  phfrq(gams%natom3)=Phonon frequencies
  gamma_ph(gams%natom3)=Phonon linewidths.
  lambda_ph(gams%natom3)=Phonon linewidths.
  displ_cart(2,3,cry%natom,3*cryst%natom)=Phonon displacement in cartesian coordinates.

PARENTS

      m_phgamma

CHILDREN

SOURCE

810 subroutine phgamma_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee)
811 
812 
813 !This section has been created automatically by the script Abilint (TD).
814 !Do not modify the following lines by hand.
815 #undef ABI_FUNC
816 #define ABI_FUNC 'phgamma_eval_qibz'
817 !End of the abilint section
818 
819  implicit none
820 
821 !Arguments ------------------------------------
822 !scalars
823  integer,intent(in) :: iq_ibz,spin
824  type(phgamma_t),intent(inout) :: gams
825  type(crystal_t),intent(in) :: cryst
826  type(ifc_type),intent(in) :: ifc
827 !arrays
828  real(dp),intent(out) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3)
829  real(dp),intent(out),optional :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3)
830  real(dp),intent(out) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
831 
832 !Local variables-------------------------------
833 !scalars
834  integer,parameter :: qtor0=0
835  integer :: natom3,nu1,nu2
836  integer :: iene, jene
837  real(dp) :: diagerr,spinfact
838  character(len=500) :: msg
839  !arrays
840  real(dp) :: displ_red(2,gams%natom3,gams%natom3)
841  real(dp) :: img(gams%natom3),gam_now(2,gams%natom3**2)
842  real(dp) :: tmp_gam1(2,gams%natom3,gams%natom3),tmp_gam2(2,gams%natom3,gams%natom3)
843  real(dp) :: pheigvec(2*(3*cryst%natom)**2)
844 
845 ! *************************************************************************
846 
847  !@phgamma_t
848  natom3 = gams%natom3
849 
850  ! Get phonon frequencies and eigenvectors.
851  call ifc_fourq(ifc,cryst,gams%qibz(:,iq_ibz),phfrq,displ_cart,out_eigvec=pheigvec, out_displ_red=displ_red)
852 
853  select case (gams%eph_scalprod)
854  case (0)
855    ! If the matrices do not contain the scalar product with the displ_red vectors yet do it now.
856    tmp_gam2 = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2,natom3,natom3])
857    call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
858 
859    do nu1=1,natom3
860      gamma_ph(nu1) = tmp_gam1(1, nu1, nu1)
861      img(nu1) = tmp_gam1(2, nu1, nu1)
862      if (abs(img(nu1)) > tol8) then
863        write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
864        MSG_WARNING(msg)
865      end if
866    end do
867 
868 #ifdef DEV_MJV
869    if (present (gamma_ph_ee)) then
870      do iene = 1, gams%nene
871        do jene = 1, gams%nene
872          tmp_gam2 = reshape(gams%vals_ee(:,jene,iene,:,:,iq_ibz,spin), [2,natom3,natom3])
873          call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
874 
875          do nu1=1,natom3
876            gamma_ph_ee(jene,iene,nu1) = tmp_gam1(1, nu1, nu1)
877            img(nu1) = tmp_gam1(2, nu1, nu1)
878            if (abs(img(nu1)) > tol8) then
879              write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
880              MSG_WARNING(msg)
881            end if
882          end do
883        end do
884      end do
885    end if
886 #endif
887  case (1)
888    ! Diagonalize gamma matrix at qpoint (complex matrix).
889    ! MJV NOTE: gam_now is recast implicitly here to matrix
890    gam_now = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2,natom3**2])
891 
892    call ZGEMM('N','N',natom3,natom3,natom3,cone,gam_now, natom3,pheigvec,natom3,czero,tmp_gam1,natom3)
893    call ZGEMM('C','N',natom3,natom3,natom3,cone,pheigvec,natom3,tmp_gam1 ,natom3,czero,tmp_gam2,natom3)
894 
895    diagerr = zero
896    do nu2=1,natom3
897      gamma_ph(nu2) = tmp_gam2(1,nu2,nu2)
898 
899      do nu1=1,nu2-1
900        diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
901      end do
902      do nu1=nu2+1,natom3
903        diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
904      end do
905      diagerr = diagerr + abs(tmp_gam2(2,nu2,nu2))
906    end do
907 
908    if (diagerr > tol12) then
909      write (msg,'(a,es14.6)')'Numerical error in diagonalization of gamma with phon eigenvectors: ',diagerr
910      MSG_WARNING(msg)
911    end if
912 
913  case default
914    MSG_BUG(sjoin("Wrong value for eph_scalprod:",itoa(gams%eph_scalprod)))
915  end select
916 
917  ! Compute lambda
918  ! TODO : check this - looks like a factor of 2 wrt the inline documentation!
919  !spinfact should be 1 for a normal non sppol calculation without spinorbit
920  !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
921  !for nsppol 2 it should be 0.5 as we have 2 spin channels to sum
922  spinfact = two/(gams%nsppol*gams%nspinor)
923 
924  do nu1=1,gams%natom3
925    gamma_ph(nu1) =  gamma_ph(nu1) * pi * spinfact
926    lambda_ph(nu1) = zero
927    if (abs(phfrq(nu1)) > EPH_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
928 
929 #ifdef DEV_MJV
930    if (present(gamma_ph_ee)) then
931      gamma_ph_ee(:,:,nu1) =  gamma_ph_ee(:,:,nu1) * pi * spinfact
932    end if
933 #endif
934  end do
935 
936  ! This to avoid spurious results for the acoustic modes.
937  ! In principle, gamma(q) --> 0 and lambda(q) --> 0 for q --> 0
938  ! but the Fourier interpolated gammas do not fulfill this property so we set everything
939  ! to zero when we are inside a sphere or radius
940  if (normv(gams%qibz(:,iq_ibz), cryst%gmet, "G") < EPH_Q0TOL) then
941    gamma_ph(1:3) = zero
942    lambda_ph(1:3) = zero
943  end if
944 
945 end subroutine phgamma_eval_qibz

m_phgamma/phgamma_free [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_free

FUNCTION

  Free the dynamic memory in a <phgamma_t> datatype

PARENTS

      m_phgamma

CHILDREN

SOURCE

404 subroutine phgamma_free(gams)
405 
406 
407 !This section has been created automatically by the script Abilint (TD).
408 !Do not modify the following lines by hand.
409 #undef ABI_FUNC
410 #define ABI_FUNC 'phgamma_free'
411 !End of the abilint section
412 
413  implicit none
414 
415 !Arguments ------------------------------------
416  type(phgamma_t),intent(inout) :: gams
417 
418 ! *************************************************************************
419 
420  !@phgamma_t
421 
422  !integer
423 
424  !real
425  if (allocated(gams%n0)) then
426    ABI_FREE(gams%n0)
427  end if
428  if (allocated(gams%qibz)) then
429    ABI_FREE(gams%qibz)
430  end if
431  if (allocated(gams%wtq)) then
432    ABI_FREE(gams%wtq)
433  end if
434  if (allocated(gams%qbz)) then
435    ABI_FREE(gams%qbz)
436  end if
437  if (allocated(gams%rpt)) then
438    ABI_FREE(gams%rpt)
439  end if
440  if (allocated(gams%wghatm)) then
441    ABI_FREE(gams%wghatm)
442  end if
443  if (allocated(gams%vals_qibz)) then
444    ABI_FREE(gams%vals_qibz)
445  end if
446  if (allocated(gams%vals_bz)) then
447    ABI_FREE(gams%vals_bz)
448  end if
449  if (allocated(gams%vals_rpt)) then
450    ABI_FREE(gams%vals_rpt)
451  end if
452  if (allocated(gams%vals_in_qibz)) then
453    ABI_FREE(gams%vals_in_qibz)
454  end if
455  if (allocated(gams%vals_in_bz)) then
456    ABI_FREE(gams%vals_in_bz)
457  end if
458  if (allocated(gams%vals_in_rpt)) then
459    ABI_FREE(gams%vals_in_rpt)
460  end if
461  if (allocated(gams%vals_out_qibz)) then
462    ABI_FREE(gams%vals_out_qibz)
463  end if
464  if (allocated(gams%vals_out_bz)) then
465    ABI_FREE(gams%vals_out_bz)
466  end if
467  if (allocated(gams%vals_out_rpt)) then
468    ABI_FREE(gams%vals_out_rpt)
469  end if
470  if (allocated(gams%vals_ee)) then
471    ABI_FREE(gams%vals_ee)
472  end if
473 
474 end subroutine phgamma_free

m_phgamma/phgamma_init [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_init

FUNCTION

  Creation method for the phgamma_t datatype.

INPUTS

 cryst<crystal_t>
 ifc<ifc_type>=Interatomic force constants.
 symdynmat=1 to activa symmetrization of gamma matrices.
 eph_scalprod
 ngqpt(3)=Q-mesh divisions
 nsppol=Number of spin polarizations
 nspinor=Number of spinorial components.
 n0(nsppol)=Density of states at the Fermi level.

OUTPUT

 gams<phgamma_t>

PARENTS

      m_phgamma

CHILDREN

SOURCE

506 subroutine phgamma_init(gams,cryst,ifc,fstab,symdynmat,eph_scalprod,eph_transport,ngqpt,nsppol,nspinor,n0)
507 
508 
509 !This section has been created automatically by the script Abilint (TD).
510 !Do not modify the following lines by hand.
511 #undef ABI_FUNC
512 #define ABI_FUNC 'phgamma_init'
513 !End of the abilint section
514 
515  implicit none
516 
517 !Arguments ------------------------------------
518 !scalars
519  integer,intent(in) :: nsppol,nspinor,symdynmat,eph_scalprod
520  integer,intent(in) :: eph_transport
521  type(crystal_t),intent(in) :: cryst
522  type(ifc_type),intent(in) :: ifc
523  type(phgamma_t),intent(out) :: gams
524  type(fstab_t), intent(in) :: fstab
525 !arrays
526  integer,intent(in) :: ngqpt(3)
527  real(dp),intent(in) :: n0(nsppol)
528 
529 !Local variables-------------------------------
530 !scalars
531  integer,parameter :: qptopt1=1
532  integer :: ierr
533 !arrays
534  integer :: qptrlatt(3,3)
535 
536 ! *************************************************************************
537 
538  !@phgamma_t
539  ! Set basic dimensions.
540  gams%natom = cryst%natom; gams%natom3 = 3*cryst%natom; gams%nsppol = nsppol; gams%nspinor = nspinor
541  gams%symgamma = symdynmat; gams%eph_scalprod = eph_scalprod
542  !gams%asr = ifc%asr
543  gams%asr = 0
544 
545  gams%ndir_transp = 0; if (eph_transport > 0) gams%ndir_transp = 3
546 
547  ABI_MALLOC(gams%n0, (nsppol))
548  gams%n0 = n0
549 
550  gams%nene = fstab%nene
551  gams%enemin = fstab%enemin
552  gams%deltaene = fstab%deltaene
553 
554  ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
555  gams%ngqpt = ngqpt
556  qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
557  call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, 1, [zero, zero, zero], &
558    gams%nqibz, gams%qibz, gams%wtq, gams%nqbz, gams%qbz)
559 
560  ! Allocate matrices in the IBZ.
561  ABI_STAT_MALLOC(gams%vals_qibz, (2, gams%natom3, gams%natom3, gams%nqibz, nsppol), ierr)
562  ABI_CHECK(ierr==0, "out of memory in %vals_qibz")
563  gams%vals_qibz = zero
564 
565  ! TODO: if we remove the nsig dependency in the gvvvals_*_qibz we can remove
566  ! the intermediate array and save a lot of memory
567  ABI_STAT_MALLOC(gams%vals_in_qibz, (2,gams%ndir_transp**2,gams%natom3,gams%natom3,gams%nqibz,nsppol), ierr)
568  ABI_CHECK(ierr==0, "out of memory in %vals_in_qibz")
569  gams%vals_in_qibz = zero
570 
571  if (eph_transport > 0) then
572    ABI_STAT_MALLOC(gams%vals_out_qibz, (2,gams%ndir_transp**2,gams%natom3,gams%natom3,gams%nqibz,nsppol), ierr)
573    ABI_CHECK(ierr==0, "out of memory in %vals_out_qibz")
574    gams%vals_out_qibz = zero
575  end if
576 
577 #ifdef DEV_MJV
578  ABI_STAT_MALLOC(gams%vals_ee,(2,gams%nene,gams%nene,gams%natom3,gams%natom3,gams%nqibz,gams%nsppol), ierr)
579  ABI_CHECK(ierr==0, 'out of memory in gams%vals_ee')
580 #endif
581 
582  ! Prepare Fourier interpolation.
583  gams%gprim = ifc%gprim
584  gams%nrpt  = ifc%nrpt
585  ABI_MALLOC(gams%rpt,(3,gams%nrpt))
586  gams%rpt = ifc%rpt
587  ABI_MALLOC(gams%wghatm,(gams%natom,gams%natom,gams%nrpt))
588  gams%wghatm = ifc%wghatm
589 
590 end subroutine phgamma_init

m_phgamma/phgamma_interp [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_interp

FUNCTION

  Interpolate the phonon linewidths at a given q-point.

INPUTS

  gams<phgamma_t>
  cryst<crystal_t>=crystalline structure.
  ifc<ifc_type>=Interatomic force constants.
  spin=Spin index
  qpt(3)=q-point in reduced coordinates
  gamma_ph(3*natom)=Phonon linewidths

OUTPUT

  gamma_ph(gams%natom3)=Interpolated Phonon linewidths.
  lamda_ph(3*natom)=Lambda coefficients for the different phonon modes.
  phfrq(3*natom)=phonon frequencies at current q
  displ_cart(2,3,natom,3*natom) = Phonon displacement in Cartesian coordinates

PARENTS

      m_phgamma

CHILDREN

SOURCE

 978 subroutine phgamma_interp(gams,cryst,ifc,spin,qpt,phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee)
 979 
 980 
 981 !This section has been created automatically by the script Abilint (TD).
 982 !Do not modify the following lines by hand.
 983 #undef ABI_FUNC
 984 #define ABI_FUNC 'phgamma_interp'
 985 !End of the abilint section
 986 
 987  implicit none
 988 
 989 !Arguments ------------------------------------
 990 !scalars
 991  integer,intent(in) :: spin
 992  type(phgamma_t),intent(inout) :: gams
 993  type(crystal_t),intent(in) :: cryst
 994  type(ifc_type),intent(in) :: ifc
 995 !arrays
 996  real(dp),intent(in) :: qpt(3)
 997  real(dp),intent(out) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3)
 998  real(dp),intent(out),optional :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3)
 999  real(dp),intent(out) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
1000 
1001 !Local variables-------------------------------
1002 !scalars
1003  integer,parameter :: qtor0=0
1004  integer, save :: icall=0
1005  integer :: natom3,nu1,nu2
1006  real(dp) :: diagerr,spinfact
1007  character(len=500) :: msg
1008  !arrays
1009  real(dp) :: displ_red(2,gams%natom3,gams%natom3)
1010  real(dp) :: gam_now(2,gams%natom3**2),img(gams%natom3)
1011  real(dp) :: tmp_gam1(2,gams%natom3,gams%natom3),tmp_gam2(2,gams%natom3,gams%natom3)
1012  real(dp) :: pheigvec(2*gams%natom3**2)
1013  real(dp),allocatable :: coskr(:,:),sinkr(:,:)
1014 
1015 ! *************************************************************************
1016 
1017  ! Compute internal tables used for Fourier interpolation.
1018  if (.not.allocated(gams%vals_bz)) call phgamma_interp_setup(gams,cryst,"INIT")
1019 
1020 #ifdef DEV_MJV
1021  if (present(gamma_ph_ee) .and. icall == 0) then
1022    gamma_ph_ee = zero
1023    write (msg,'(2a)') "For the moment gams_ee matrix elements are not FT interpolated wrt q,',&
1024 &       ' only evaluated on the electron k grid. The resulting a2feew will be 0"
1025    MSG_WARNING(msg)
1026    icall = 1
1027  end if
1028 #endif
1029 
1030  !@phgamma_t
1031  natom3 = gams%natom3
1032 
1033  ! Taken from mkph_linwid
1034  ! This reduced version of ftgkk supposes the kpoints have been integrated
1035  ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt.
1036  ABI_MALLOC(coskr, (1,gams%nrpt))
1037  ABI_MALLOC(sinkr, (1,gams%nrpt))
1038  ! TODO: This is not optimal
1039  call ftgam_init(gams%gprim, 1, gams%nrpt, qpt, gams%rpt, coskr, sinkr)
1040 
1041  call ftgam(gams%wghatm,gam_now,gams%vals_rpt(:,:,:,spin),gams%natom,1,gams%nrpt,qtor0,coskr, sinkr)
1042 
1043  ! This call is not executed in elphon!
1044  if (gams%symgamma == 1) call tgamma_symm(cryst,qpt,gam_now)
1045 
1046  ABI_FREE(coskr)
1047  ABI_FREE(sinkr)
1048 
1049  ! Get phonon frequencies and eigenvectors.
1050  call ifc_fourq(ifc,cryst,qpt,phfrq,displ_cart,out_eigvec=pheigvec, out_displ_red=displ_red)
1051 
1052  select case (gams%eph_scalprod)
1053  case (0)
1054    ! If the matrices do not contain the scalar product with the displ_cart vectors yet do it now.
1055    tmp_gam2 = reshape (gam_now, [2,natom3,natom3])
1056    call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
1057 
1058    do nu1=1,natom3
1059      gamma_ph(nu1) = tmp_gam1(1, nu1, nu1)
1060      img(nu1) = tmp_gam1(2, nu1, nu1)
1061      if (abs(img(nu1)) > tol8) then
1062        write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
1063        MSG_WARNING(msg)
1064      end if
1065    end do
1066 
1067 
1068 !   do iene = 1, gams%nene
1069 !     do jene = 1, gams%nene
1070 !       tmp_gam2 = reshape (gam_now, [2,natom3,natom3])
1071 !       call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
1072 !       do nu1=1,natom3
1073 !         gamma_ph(nu1) = tmp_gam1(1, nu1, nu1)
1074 !         img(nu1) = tmp_gam1(2, nu1, nu1)
1075 !         if (abs(img(nu1)) > tol8) then
1076 !           write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
1077 !           MSG_WARNING(msg)
1078 !         end if
1079 !       end do
1080 !     end do
1081 !   end do
1082 
1083  case (1)
1084    ! Diagonalize gamma matrix at qpoint (complex matrix).
1085    ! MJV NOTE: gam_now is recast implicitly here to matrix
1086    call ZGEMM('N','N',natom3,natom3,natom3,cone,gam_now, natom3,pheigvec,natom3,czero,tmp_gam1,natom3)
1087    call ZGEMM('C','N',natom3,natom3,natom3,cone,pheigvec,natom3,tmp_gam1 ,natom3,czero,tmp_gam2,natom3)
1088 
1089    diagerr = zero
1090    do nu2=1,natom3
1091      gamma_ph(nu2) = tmp_gam2(1,nu2,nu2)
1092 
1093      do nu1=1,nu2-1
1094        diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1095      end do
1096      do nu1=nu2+1,natom3
1097        diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1098      end do
1099      diagerr = diagerr + abs(tmp_gam2(2,nu2,nu2))
1100    end do
1101 
1102    if (diagerr > tol12) then
1103      write (msg,'(a,es14.6)')'Numerical error in diagonalization of gamma with phon eigenvectors: ',diagerr
1104      MSG_WARNING(msg)
1105    end if
1106 
1107  case default
1108    MSG_BUG(sjoin("Wrong value for eph_scalprod:",itoa(gams%eph_scalprod)))
1109  end select
1110 
1111  ! Compute lambda
1112  !spinfact should be 1 for a normal non sppol calculation without spinorbit
1113  !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1114  !for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1115  spinfact = two/(gams%nsppol*gams%nspinor)
1116 
1117  ! Compute lambda
1118  do nu1=1,gams%natom3
1119    gamma_ph(nu1) = gamma_ph(nu1) * pi * spinfact
1120    lambda_ph(nu1) = zero
1121    if (abs(phfrq(nu1)) > EPH_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1122  end do
1123 
1124  ! This to avoid spurious results for the acoustic modes.
1125  ! In principle, gamma(q) --> 0 and lambda(q) --> 0 for q --> 0
1126  ! but the Fourier interpolated gammas do not fulfill this property so we set everything
1127  ! to zero when we are inside a sphere or radius
1128  if (normv(qpt, cryst%gmet, "G") < EPH_Q0TOL) then
1129    write(std_out,*)"Setting values to zero."
1130    gamma_ph(1:3) = zero
1131    lambda_ph(1:3) = zero
1132  end if
1133 
1134 end subroutine phgamma_interp

m_phgamma/phgamma_interp_setup [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_interp_setup

FUNCTION

  This routines performs the (allocation|deallocation) of the internal tables
  used to interpolate the linewidths in q-space

INPUTS

  action =
    "INIT" to allocate and compute the internal tables (default)
    "FREE" to deallocate the internal tables.

SIDE EFFECTS

  gams<phgamma_t>= gams%vals_bz and gams%vals_rpt, depending on action.

PARENTS

      m_phgamma

CHILDREN

SOURCE

1162 subroutine phgamma_interp_setup(gams,cryst,action)
1163 
1164 
1165 !This section has been created automatically by the script Abilint (TD).
1166 !Do not modify the following lines by hand.
1167 #undef ABI_FUNC
1168 #define ABI_FUNC 'phgamma_interp_setup'
1169 !End of the abilint section
1170 
1171  implicit none
1172 
1173 !Arguments ------------------------------------
1174 !scalars
1175  character(len=*),intent(in) :: action
1176  type(phgamma_t),intent(inout) :: gams
1177  type(crystal_t),intent(in) :: cryst
1178 
1179 !Local variables-------------------------------
1180 !scalars
1181  integer,parameter :: qtor1=1
1182  integer :: iq_bz,iq_ibz,isq_bz,spin,isym,ierr,ii
1183  !character(len=500) :: msg
1184  type(kptrank_type) :: qrank
1185 !arrays
1186  integer :: g0(3)
1187  integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:)
1188  real(dp) :: qirr(3),tmp_qpt(3)
1189  real(dp),allocatable :: coskr(:,:),sinkr(:,:)
1190  real(dp),allocatable :: gamma_qpt(:,:,:,:),atmfrc(:,:)
1191 
1192 ! *************************************************************************
1193 
1194  select case (toupper(action))
1195 
1196  case ("INIT")
1197    if (.not.allocated(gams%vals_bz)) then
1198      ABI_STAT_MALLOC(gams%vals_bz,(2,gams%natom3**2,gams%nqbz,gams%nsppol), ierr)
1199      ABI_CHECK(ierr==0, 'out of memory in gams%vals_bz')
1200      gams%vals_bz = zero
1201 
1202      ! Build tables needed by complete_gamma.
1203      call mkkptrank(gams%qbz,gams%nqbz,qrank)
1204 
1205      ! Compute index of IBZ q-point in the BZ array
1206      ABI_CALLOC(qirredtofull,(gams%nqibz))
1207 
1208      do iq_ibz=1,gams%nqibz
1209        qirr = gams%qibz(:,iq_ibz)
1210        iq_bz = kptrank_index(qrank, qirr)
1211        if (iq_bz /= -1) then
1212          ABI_CHECK(isamek(qirr,gams%qbz(:,iq_bz),g0), "isamek")
1213          qirredtofull(iq_ibz) = iq_bz
1214        else
1215          MSG_ERROR(sjoin("Full BZ does not contain IBZ q-point:", ktoa(qirr)))
1216        end if
1217      end do
1218 
1219      ! Build qpttoqpt table. See also mkqptequiv
1220      ABI_MALLOC(qpttoqpt,(2,cryst%nsym,gams%nqbz))
1221      qpttoqpt = -1
1222      do iq_bz=1,gams%nqbz
1223        do isym=1,cryst%nsym
1224          tmp_qpt = matmul(cryst%symrec(:,:,isym), gams%qbz(:,iq_bz))
1225 
1226          isq_bz = kptrank_index(qrank, tmp_qpt)
1227          if (isq_bz == -1) then
1228            MSG_ERROR("looks like no kpoint equiv to q by symmetry without time reversal!")
1229          end if
1230          qpttoqpt(1,isym,isq_bz) = iq_bz
1231 
1232          ! q --> -q
1233          tmp_qpt = -tmp_qpt
1234          isq_bz = kptrank_index(qrank, tmp_qpt)
1235          if (isq_bz == -1) then
1236            MSG_ERROR("looks like no kpoint equiv to q by symmetry with time reversal!")
1237          end if
1238          qpttoqpt(2,isym,isq_bz) = iq_bz
1239        end do
1240      end do
1241      call destroy_kptrank(qrank)
1242 
1243      ! Fill BZ array with IBZ data.
1244      do spin=1,gams%nsppol
1245        do iq_ibz=1,gams%nqibz
1246          iq_bz = qirredtofull(iq_ibz)
1247          gams%vals_bz(:,:,iq_bz,spin) = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2,gams%natom3**2])
1248        end do
1249      end do
1250 
1251      ! Complete vals_bz in the full BZ.
1252      ! FIXME: Change complete_gamma API to pass (..., nsppol)
1253      ABI_MALLOC(gamma_qpt, (2, gams%natom3**2, gams%nsppol, gams%nqbz))
1254      do spin=1,gams%nsppol
1255        gamma_qpt(:, :, spin, :) = gams%vals_bz(:, :, :, spin)
1256      end do
1257 
1258      call complete_gamma(cryst,gams%natom3,gams%nsppol,gams%nqibz,gams%nqbz,&
1259        gams%eph_scalprod,qirredtofull,qpttoqpt,gamma_qpt)
1260        !gams%eph_scalprod,qirredtofull,qpttoqpt,gams%vals_bz)
1261 
1262      do spin=1,gams%nsppol
1263        gams%vals_bz(:, :, :, spin) = gamma_qpt(:, :, spin, :)
1264      end do
1265      ABI_FREE(gamma_qpt)
1266 
1267      ! TODO: idem for vv_vals 3x3 matrices
1268 
1269      ABI_FREE(qirredtofull)
1270      ABI_FREE(qpttoqpt)
1271 
1272      ! This call is not executed in elphon!
1273      if (gams%symgamma == 1) then
1274        do spin=1,gams%nsppol
1275          do iq_bz=1,gams%nqbz
1276            call tgamma_symm(cryst,gams%qbz(:,iq_bz),gams%vals_bz(:,:,iq_bz,spin))
1277          end do
1278        end do
1279      end if
1280    end if ! first allocation and filling of arrays in qbz
1281 
1282    ! Now FT to real space too
1283    ! NOTE: gprim (not gprimd) is used for all FT interpolations,
1284    ! to be consistent with the dimensions of the rpt, which come from anaddb.
1285    ! TODO: this is needed only if FT is used, no when the linear interpolation is employed.
1286    if (.not.allocated(gams%vals_rpt)) then
1287      ABI_STAT_MALLOC(gams%vals_rpt,(2,gams%natom3**2,gams%nrpt,gams%nsppol), ierr)
1288      ABI_CHECK(ierr==0, 'out of memory in gams%vals_rpt')
1289      gams%vals_rpt = zero
1290 
1291      ! q --> r
1292      ABI_MALLOC(coskr, (gams%nqbz,gams%nrpt))
1293      ABI_MALLOC(sinkr, (gams%nqbz,gams%nrpt))
1294      call ftgam_init(gams%gprim, gams%nqbz, gams%nrpt, gams%qbz, gams%rpt, coskr, sinkr)
1295 
1296      do spin=1,gams%nsppol
1297        call ftgam(gams%wghatm,gams%vals_bz(:,:,:,spin),gams%vals_rpt(:,:,:,spin),gams%natom,gams%nqbz,&
1298           gams%nrpt,qtor1, coskr, sinkr)
1299 
1300        ! Enforce "acoustic" rule on vals_rpt
1301        ! This call is not executed in elphon!
1302        if (gams%asr /= 0) then
1303          ABI_MALLOC(atmfrc, (3*gams%natom*3*gams%natom,gams%nrpt))
1304          do ii=1,2
1305            atmfrc = gams%vals_rpt(ii,:,:,spin)
1306            !gals%vals_rpt(2,:,,spin) = zero
1307            call asrif9(gams%asr, atmfrc, gams%natom, gams%nrpt, gams%rpt, gams%wghatm)
1308            gams%vals_rpt(ii,:,:,spin) = atmfrc
1309          end do
1310          ABI_FREE(atmfrc)
1311        end if
1312      end do
1313 
1314      ABI_FREE(coskr)
1315      ABI_FREE(sinkr)
1316    end if ! allocation and filling of rpt as well
1317 
1318  case ("FREE")
1319    if (allocated(gams%vals_bz)) then
1320      ABI_FREE(gams%vals_bz)
1321    end if
1322    if (allocated(gams%vals_rpt)) then
1323      ABI_FREE(gams%vals_rpt)
1324    end if
1325    if (allocated(gams%vals_ee)) then
1326      ABI_FREE(gams%vals_ee)
1327    end if
1328 
1329  case default
1330    MSG_BUG(sjoin("Wrong action:",action))
1331  end select
1332 
1333 end subroutine phgamma_interp_setup

m_phgamma/phgamma_linwid [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_linwid

FUNCTION

  Interpolate the phonon linewidths along an arbitrary q-path (use Fourier interpolation).

INPUTS

  cryst<crystal_t>=Info on the unit cell and symmetries.
  ifc<ifc_type>=Interatomic force constants.
  ndivsm=Number of points used to sample the smallest segment
  nvert = Number of extrema in qverts
  qverts(3,nvert) = vertices of reciprocal space trajectory
  basename=name used to create the different output files (text format).
  ncid=Netcdf file handler (already open in the caller).
  comm=MPI communicator

OUTPUT

  wminmax=Minimum and max phonon frequency obtained on the path (Hartree units)

PARENTS

      m_phgamma

CHILDREN

SOURCE

1968 subroutine phgamma_linwid(gams,cryst,ifc,ndivsm,nvert,qverts,basename,ncid,wminmax,comm)
1969 
1970 
1971 !This section has been created automatically by the script Abilint (TD).
1972 !Do not modify the following lines by hand.
1973 #undef ABI_FUNC
1974 #define ABI_FUNC 'phgamma_linwid'
1975 !End of the abilint section
1976 
1977  implicit none
1978 
1979 !Arguments ------------------------------------
1980 !scalars
1981  integer,intent(in) :: nvert,ndivsm,comm,ncid
1982  type(crystal_t),intent(in) :: cryst
1983  type(ifc_type),intent(in) :: ifc
1984  type(phgamma_t),intent(inout) :: gams
1985  character(len=*),intent(in) :: basename
1986 !arrays
1987  real(dp),intent(in) :: qverts(3,nvert)
1988  real(dp),intent(out) :: wminmax(2)
1989 
1990 !Local variables-------------------------------
1991 !scalars
1992  integer,parameter :: master=0
1993  integer :: natom,ii,mu,iqpt,natom3,nsppol,ierr
1994  integer :: spin,unt,nqpt,nrpt,cnt,nproc,my_rank
1995 #ifdef HAVE_NETCDF
1996  integer :: ncerr
1997 #endif
1998  real(dp) :: omega_min,omega_max,wtmp,omega
1999  character(len=500) :: msg
2000  type(kpath_t) :: qpath
2001 !arrays
2002  real(dp) :: gamma_spin(gams%nsppol),lambda_spin(gams%nsppol)
2003  real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom)
2004  real(dp) :: phfrq(3*cryst%natom),gamma_ph(3*cryst%natom),lambda_ph(3*cryst%natom)
2005  real(dp) :: qpt(3),shift(3)
2006  real(dp),allocatable :: all_phfreq(:,:),all_gammaq(:,:,:),all_lambdaq(:,:,:),all_displ_cart(:,:,:,:)
2007 
2008 ! *********************************************************************
2009 
2010  DBG_ENTER("COLL")
2011 
2012  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2013  natom = cryst%natom; natom3 = gams%natom3; nsppol = gams%nsppol; nrpt = gams%nrpt
2014 
2015  ! Define the q-path along which phonon linwid will be interpolated.
2016  qpath = kpath_new(qverts, cryst%gprimd, ndivsm)
2017  nqpt = qpath%npts
2018 
2019  ! Allocate workspace arrays for MPI.
2020  ABI_CALLOC(all_phfreq, (natom3, nqpt))
2021  ABI_CALLOC(all_gammaq, (natom3, nqpt, nsppol))
2022  ABI_CALLOC(all_lambdaq, (natom3, nqpt, nsppol))
2023  ABI_CALLOC(all_displ_cart, (2, natom3, natom3, nqpt))
2024 
2025  ! initialize the minimum and maximum phonon frequency
2026  omega_min = huge(one); omega_max = -huge(one)
2027 
2028  ! Interpolation along specified path in q space (keep spin dep.)
2029  cnt = 0
2030  do spin=1,nsppol
2031    do iqpt=1,nqpt
2032      cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2033      call wrap2_pmhalf(qpath%points(:,iqpt), qpt, shift)
2034 
2035      ! Get phgamma
2036      call phgamma_interp(gams,cryst,ifc,spin,qpt,phfrq,gamma_ph,lambda_ph, displ_cart)
2037      all_gammaq(:, iqpt, spin) = gamma_ph
2038      all_lambdaq(:, iqpt, spin) = lambda_ph
2039      if (spin == 1) then
2040        all_phfreq(:, iqpt) = phfrq
2041        all_displ_cart(:, :, :, iqpt) = displ_cart
2042      end if
2043 
2044      ! Find max/min phonon frequency along path chosen
2045      ! presumed to be representative of full BZ to within 10 percent
2046      omega_min = min(omega_min, phfrq(1))
2047      omega_max = max(omega_max, phfrq(natom3))
2048    end do ! end iqpt do
2049  end do ! spin
2050 
2051  ! Collect results
2052  if (omega_min > tol12) omega_min = zero
2053  wtmp = omega_min; call xmpi_min(wtmp, omega_min, comm, ierr)
2054  wtmp = omega_max; call xmpi_max(wtmp, omega_max, comm, ierr)
2055  wminmax = [omega_min, omega_max]
2056 
2057  call xmpi_sum_master(all_gammaq, master, comm, ierr)
2058  call xmpi_sum_master(all_lambdaq, master, comm, ierr)
2059  call xmpi_sum_master(all_phfreq, master, comm, ierr)
2060  call xmpi_sum_master(all_displ_cart, master, comm, ierr)
2061 
2062  ! Master writes text file with final results.
2063  ! 3 * natom blocks, one block for each phonon mode.
2064  ! Each block contains:
2065  !   iqpt  omega_(q) gamma(q) lambda(q) nesting(q) ....
2066 
2067  if (xmpi_comm_rank(comm) == master) then
2068    if (open_file(strcat(basename, '_PHGAMMA'),msg,newunit=unt,form="formatted",status="unknown") /= 0) then
2069      MSG_ERROR(msg)
2070    end if
2071 
2072    write(unt,'(a)')     '#'
2073    write(unt,'(a)')     '# ABINIT package: E-PH band structure file. Hartree units'
2074    write(unt,'(a)')     '#'
2075    write(unt,'(a,i0,a)')'# Phonon frequencies, ph linewidths and lambda calculated on ',nqpt,' q-points'
2076    write(unt,"(a,i0)")  '# eph_scalprod = ',gams%eph_scalprod
2077    call kpath_print(qpath, header="Description of the q-path:", unit=unt, pre="#")
2078    do ii=1,2; write(unt,'(a)')     "# "; end do
2079 
2080    write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(gams%n0)
2081    do spin=1,nsppol
2082      write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",gams%n0(spin)
2083    end do
2084 
2085    do mu=1,natom3
2086      write(unt,'(a)')"#"
2087      if (nsppol == 1) write(unt,'(a,i0,a)')"# phonon mode [",mu,"] q-index omega gamma lambda"
2088      if (nsppol == 2) write(unt,'(a,i0,a)')&
2089        "# phonon mode [",mu,"] q-index omega gamma_tot lambda_tot gamma[spin=1] lambda[spin=1] gamma[2] lambda[2]"
2090      write(unt,'(a)')"#"
2091      do iqpt=1,nqpt
2092        omega = all_phfreq(mu, iqpt)
2093        gamma_spin = all_gammaq(mu, iqpt, :)
2094        lambda_spin = all_lambdaq(mu, iqpt, :)
2095        if (nsppol == 1) then
2096          write(unt,'(i8,3es16.6)' )iqpt,omega,gamma_spin(1),lambda_spin(1)
2097        else
2098          write(unt,'(i8,es20.6,6es16.6)' )iqpt,omega,&
2099             sum(gamma_spin),sum(lambda_spin),&
2100             gamma_spin(1),lambda_spin(1),&
2101             gamma_spin(2),lambda_spin(2)
2102        end if
2103      end do
2104    end do
2105 
2106    close(unt)
2107 
2108    ! Write data to netcdf file
2109    if (ncid /= nctk_noid) then
2110 #ifdef HAVE_NETCDF
2111      ncerr = nctk_def_dims(ncid, [&
2112        nctkdim_t("natom3", 3*natom), nctkdim_t("nqpath", nqpt), nctkdim_t("number_of_spins", nsppol) &
2113        ], defmode=.True.)
2114      NCF_CHECK(ncerr)
2115 
2116      ncerr = nctk_def_arrays(ncid, [&
2117        nctkarr_t('qpath', "dp", "number_of_reduced_dimensions, nqpath"), &
2118        nctkarr_t('phfreq_qpath', "dp", "natom3, nqpath"), &
2119        nctkarr_t('phdispl_cart_qpath', "dp", "two, natom3, natom3, nqpath"), &
2120        nctkarr_t('phgamma_qpath', "dp", "natom3, nqpath, number_of_spins"), &
2121        nctkarr_t('phlambda_qpath', "dp", "natom3, nqpath, number_of_spins")])
2122      NCF_CHECK(ncerr)
2123 
2124      NCF_CHECK(nctk_set_datamode(ncid))
2125      NCF_CHECK(nf90_put_var(ncid, vid("qpath"), qpath%points))
2126      NCF_CHECK(nf90_put_var(ncid, vid("phfreq_qpath"), all_phfreq))
2127      NCF_CHECK(nf90_put_var(ncid, vid("phdispl_cart_qpath"), all_displ_cart))
2128      NCF_CHECK(nf90_put_var(ncid, vid("phgamma_qpath"), all_gammaq))
2129      NCF_CHECK(nf90_put_var(ncid, vid("phlambda_qpath"), all_lambdaq))
2130 #endif
2131    end if
2132 
2133  end if ! master
2134 
2135  ABI_FREE(all_phfreq)
2136  ABI_FREE(all_gammaq)
2137  ABI_FREE(all_lambdaq)
2138  ABI_FREE(all_displ_cart)
2139 
2140  call kpath_free(qpath)
2141 
2142  DBG_EXIT("COLL")
2143 
2144 contains
2145  integer function vid(vname)
2146 
2147 !This section has been created automatically by the script Abilint (TD).
2148 !Do not modify the following lines by hand.
2149 #undef ABI_FUNC
2150 #define ABI_FUNC 'vid'
2151 !End of the abilint section
2152 
2153    character(len=*),intent(in) :: vname
2154    vid = nctk_idname(ncid, vname)
2155  end function vid
2156 
2157 end subroutine phgamma_linwid

m_phgamma/phgamma_print [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_print

FUNCTION

  Finalize the phgamma_t datatype.

INPUTS

 cryst<crystal_t>=Crystalline structure.
 ifc<ifc_type>=Interatomic force constants.
 ncid=Netcdf file handler (already open in the caller).

OUTPUT

 gams<phgamma_t>

PARENTS

      m_phgamma

CHILDREN

SOURCE

617 subroutine phgamma_print(gams,cryst,ifc,ncid)
618 
619 
620 !This section has been created automatically by the script Abilint (TD).
621 !Do not modify the following lines by hand.
622 #undef ABI_FUNC
623 #define ABI_FUNC 'phgamma_print'
624 !End of the abilint section
625 
626  implicit none
627 
628 !Arguments ------------------------------------
629 !scalars
630  type(phgamma_t),intent(inout) :: gams
631  type(crystal_t),intent(in) :: cryst
632  type(ifc_type),intent(in) :: ifc
633  integer,intent(in) :: ncid
634 
635 !Local variables-------------------------------
636 !scalars
637  integer :: iq_ibz,spin,mu
638  real(dp) :: lambda_tot
639  character(len=500) :: msg
640 !arrays
641  real(dp) :: phfrq(3*cryst%natom),gamma_ph(3*cryst%natom),lambda_ph(3*cryst%natom)
642  real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom)
643 
644 ! *************************************************************************
645 
646  ! ==========================================================
647  ! write data to files for each q point
648  ! ==========================================================
649  ! Compute total lambda
650  lambda_tot = zero
651  do spin=1,gams%nsppol
652    do iq_ibz=1,gams%nqibz
653 
654      ! Get phonon frequencies and eigenvectors.
655      call phgamma_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph,displ_cart)
656 
657      do mu=1,gams%natom3
658        lambda_tot = lambda_tot + lambda_ph(mu) * gams%wtq(iq_ibz)
659      end do
660 
661 #ifdef HAVE_NETCDF
662    ! Write data to netcdf file
663    if (ncid /= nctk_noid) then
664      !NCF_CHECK(nctk_set_datamode(ncid))
665      ! TODO: why does this depend on spin?????
666      if (spin == 1) then
667        NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qibz"), phfrq, start=[1, iq_ibz]))
668        NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_cart_qibz'), displ_cart, start=[1, 1, 1, iq_ibz]))
669      end if
670      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phgamma_qibz'), gamma_ph, start=[1, iq_ibz, spin]))
671      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phlambda_qibz'), lambda_ph, start=[1, iq_ibz, spin]))
672    end if
673 #endif
674 
675      ! Output to the main output file
676      if (gams%nsppol == 2) then
677        write(msg,'(2a,3es16.6,a,i1,a,a)')ch10,&
678          ' q-point =',gams%qibz(:, iq_ibz),'   spin = ',spin,ch10,&
679          ' Mode number    Frequency (Ha)  Linewidth (Ha)  Lambda(q,n)'
680      else
681        write(msg,'(2a,3es16.6,a,a)')ch10,&
682          ' q-point =',gams%qibz(:, iq_ibz),ch10,&
683          ' Mode number    Frequency (Ha)  Linewidth (Ha)  Lambda(q,n)'
684      end if
685      call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL')
686 
687      do mu=1,gams%natom3
688        write(msg,'(i5,es20.6,2es16.6)')mu,phfrq(mu),gamma_ph(mu),lambda_ph(mu)
689        call wrtout(std_out,msg,'COLL'); call wrtout(ab_out,msg,'COLL')
690      end do
691 
692    end do
693    ! Add blank lines to output files between spins
694    call wrtout(std_out,"",'COLL'); call wrtout(ab_out,"",'COLL')
695  end do
696 
697  write(ab_out,"(a,f8.4)")" lambda= ",lambda_tot
698 
699 end subroutine phgamma_print

m_phgamma/phgamma_t [ Types ]

[ Top ] [ m_phgamma ] [ Types ]

NAME

 phgamma_t

FUNCTION

 Provides methods for computing phonon linewidths, interpolating the results
 in q-space and evaluate superconducting properties.

SOURCE

 84  type,public :: phgamma_t
 85 
 86   integer :: natom
 87   ! Number of atoms per unit cell.
 88 
 89   integer :: natom3
 90   ! Number of phonon branches i.e. 3*natom.
 91 
 92   integer :: nsppol
 93   ! Number of independent spin polarizations.
 94 
 95   integer :: nspinor
 96   ! Number of spinorial components.
 97 
 98   integer :: nqibz
 99   ! Number of q-points in the IBZ.
100 
101   integer :: nqbz
102   ! Number of q-points in the BZ.
103 
104   integer :: eph_scalprod
105   !
106   integer :: nrpt
107   ! Number of points in the real space representation of the gamma matrices.
108 
109   integer :: symgamma
110   ! 1 if gamma matrices should be symmetrized by symdyma when using Fourier interpolation
111 
112   integer :: asr
113   ! If "Acoustic rule at gamma should be enforced.
114 
115   integer :: ndir_transp
116   ! 0 if no transport, otherwise 3
117 
118   integer :: ngqpt(3)
119   ! Number of divisions in the Q mesh.
120 
121   integer :: nene
122   ! Number of chemical potential values used for inelastic integration
123 
124   real(dp) :: enemin
125   ! Minimal chemical potential value used for inelastic integration Copied from fstab
126 
127   real(dp) :: deltaene
128   ! Chemical potential increment for inelastic integration Copied from fstab
129   !   for simplicity could be made equal to phonon frequency step
130 
131   real(dp) :: gprim(3,3)
132   ! Needed for Fourier interpolation.
133   ! NOTE: gprim (not gprimd) is used for all FT interpolations,
134   ! to be consistent with the dimensions of the rpt, which come from anaddb.
135 
136   real(dp),allocatable :: n0(:)
137    ! n0(gams%nsppol)=Density of states at the Fermi level.
138 
139   real(dp),allocatable :: qibz(:,:)
140   ! qibz(3,nqibz)
141   ! Reduced coordinates of the q-points in the IBZ.
142 
143   real(dp),allocatable :: wtq(:)
144   ! wtq(nqibz)
145   ! Weights of the q-points in the IBZ (normalized to one)
146 
147   real(dp),allocatable :: qbz(:,:)
148   ! qbz(3,nqbz)
149   ! Reduced coordinates of the q-points in the BZ.
150 
151   real(dp),allocatable :: rpt(:,:)
152   ! rpt(3,nrpt)
153   !  Reduced coordinates ***in terms of rprim*** of the lattice points used
154   !  for the Fourier transform of the phonon linewidths.
155 
156   real(dp),allocatable :: wghatm(:,:,:)
157   ! wghatm(natom,natom,nrpt)
158   ! Weights used in the FT of the phonon linewidths.
159 
160   real(dp),allocatable :: vals_qibz(:,:,:,:,:)
161   ! vals_qibz(2,natom3,natom3,nqibz,nsppol)) in reduced coordinates for each q-point in the IBZ.
162   ! matii_qibz {\tau'\alpha',\tau\alpha} = sum over k, ib1, ib2
163   !   <psi_{k+q,ib1} | H(1)_{\tau'\alpha'} | psi_{k,ib2}>*  \cdot
164   !   <psi_{k+q,ib1} | H(1)_{\tau \alpha } | psi_{k,ib2}>
165 
166   !NOTE: choice to put nsppol before or after nqbz is a bit arbitrary
167   !   abinit uses nband,nkpt,nsppol, but here for convenience nkpt_phon,nsppol,nqbz interpolation is on qpt
168   !MG: I think that nsppol should be the last dimension set call to ftgam.
169 
170   real(dp),allocatable :: vals_bz(:,:,:,:)
171   ! vals_bz(2,natom3**2,nqbz,nsppol)
172   ! tgamma matrices in reduced coordinates for each q in full zone
173   ! integrated over kpoints coeff and bands: still depends on qpt and spin.
174 
175   real(dp),allocatable :: vals_rpt(:,:,:,:)
176   ! vals_rpt(2,natom3**2,nrpt,nsppol)
177   ! tgamma matrices in real space in reduced coordinates.
178 
179   ! transport stuff with velocity factors
180   real(dp),allocatable :: vals_in_qibz(:,:,:,:,:,:)
181   real(dp),allocatable :: vals_out_qibz(:,:,:,:,:,:)
182   ! vals_XX_qibz(2,ndir_transp**2,natom3,natom3,nqibz,nsppol)) in reduced coordinates for each q-point in the IBZ.
183 
184   real(dp),allocatable :: vals_in_bz(:,:,:,:,:)
185   real(dp),allocatable :: vals_out_bz(:,:,:,:,:)
186   ! vals_XX_bz(2,ndir_transp**2,natom3**2,nqbz,nsppol)) in reduced coordinates for each q-point in the IBZ.
187 
188   real(dp),allocatable :: vals_in_rpt(:,:,:,:,:)
189   real(dp),allocatable :: vals_out_rpt(:,:,:,:,:)
190   ! vals_XX_rpt(2,ndir_transp**2,natom3**2,nrpt,nsppol)
191   ! tgamma matrices in real space in reduced coordinates.
192 
193   ! gamma matrices keeping full electron energy dependency
194   real(dp),allocatable :: vals_ee(:,:,:,:,:,:,:)
195   ! vals_eew(2, nene, nene, natom3, natom3, nqibz, nsppol)
196 
197  end type phgamma_t
198 
199  public :: phgamma_free          ! Free memory.
200  public :: phgamma_init          ! Creation method.
201  public :: phgamma_interp        ! Interpolates the phonon linewidths.
202  public :: phgamma_eval_qibz     ! Evaluate phonon linewidths without Fourier interpolation.
203  public :: phgamma_interp_setup  ! Compute the tables used for the interpolation in q-space.
204  public :: phgamma_linwid        ! Interpolate linewidths along an arbitrary q-path.

m_phgamma/phgamma_vv_eval_qibz [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_vv_eval_qibz

FUNCTION

 Compute the phonon linewidths times velocity squared, for q-points in the IBZ without performing interpolation.

INPUTS

  gams<phgamma_t>
  cryst<crystal_t>=Crystal structure.
  ifc<ifc_type>=Interatomic force constants.
  iq_ibz=Index of the q-point in the IBZ array.
  spin=Spin index

OUTPUT

  phfrq(gams%natom3)=Phonon frequencies
  gamma_in_ph(gams%natom3)=Phonon linewidths.
  gamma_out_ph(gams%natom3)=Phonon linewidths.
  lambda_in_ph(gams%natom3)=Phonon linewidths.
  lambda_out_ph(gams%natom3)=Phonon linewidths.

PARENTS

      m_phgamma

CHILDREN

SOURCE

1366 subroutine phgamma_vv_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph)
1367 
1368 
1369 !This section has been created automatically by the script Abilint (TD).
1370 !Do not modify the following lines by hand.
1371 #undef ABI_FUNC
1372 #define ABI_FUNC 'phgamma_vv_eval_qibz'
1373 !End of the abilint section
1374 
1375  implicit none
1376 
1377 !Arguments ------------------------------------
1378 !scalars
1379  integer,intent(in) :: iq_ibz,spin
1380  type(phgamma_t),intent(inout) :: gams
1381  type(crystal_t),intent(in) :: cryst
1382  type(ifc_type),intent(in) :: ifc
1383 !arrays
1384  real(dp),intent(out) :: phfrq(gams%natom3)
1385  real(dp),intent(out) :: gamma_in_ph(3,3,gams%natom3),lambda_in_ph(3,3,gams%natom3)
1386  real(dp),intent(out) :: gamma_out_ph(3,3,gams%natom3),lambda_out_ph(3,3,gams%natom3)
1387 
1388 !Local variables-------------------------------
1389 !scalars
1390  integer,parameter :: qtor0=0
1391  integer :: natom3,nu1,nu2
1392  integer :: idir, jdir, ii
1393  real(dp) :: diagerr,spinfact
1394  character(len=500) :: msg
1395  !arrays
1396  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
1397  real(dp) :: displ_red(2,gams%natom3,gams%natom3)
1398  real(dp) :: img(gams%natom3)
1399  real(dp) :: gam_now(2,gams%natom3,gams%natom3)
1400  real(dp) :: tmp_gam1(2,gams%natom3,gams%natom3),tmp_gam2(2,gams%natom3,gams%natom3)
1401  real(dp) :: pheigvec(2*(3*cryst%natom)**2)
1402 
1403 ! *************************************************************************
1404 
1405  !@phgamma_t
1406  natom3 = gams%natom3
1407 
1408  ! Get phonon frequencies and eigenvectors.
1409  call ifc_fourq(ifc,cryst,gams%qibz(:,iq_ibz),phfrq,displ_cart,out_eigvec=pheigvec, out_displ_red=displ_red)
1410 
1411  select case (gams%eph_scalprod)
1412  case (0)
1413    do jdir = 1,gams%ndir_transp
1414      do idir = 1,gams%ndir_transp
1415        ii = idir + gams%ndir_transp*(jdir-1)
1416        ! If the matrices do not contain the scalar product with the displ_red vectors yet do it now.
1417        tmp_gam2 = reshape(gams%vals_in_qibz(:,ii,:,:,iq_ibz,spin), [2,natom3,natom3])
1418        call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
1419 
1420        do nu1=1,natom3
1421          gamma_in_ph(idir,jdir,nu1) = tmp_gam1(1, nu1, nu1)
1422          img(nu1) = tmp_gam1(2, nu1, nu1)
1423          if (abs(img(nu1)) > tol8) then
1424            write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
1425            MSG_WARNING(msg)
1426          end if
1427        end do
1428 
1429        tmp_gam2 = reshape(gams%vals_out_qibz(:,ii,:,:,iq_ibz,spin), [2,natom3,natom3])
1430        call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
1431 
1432        do nu1=1,natom3
1433          gamma_out_ph(idir,jdir,nu1) = tmp_gam1(1, nu1, nu1)
1434          img(nu1) = tmp_gam1(2, nu1, nu1)
1435          if (abs(img(nu1)) > tol8) then
1436            write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
1437            MSG_WARNING(msg)
1438          end if
1439        end do
1440 
1441      end do ! idir
1442    end do ! jdir
1443 
1444  case (1)
1445    do jdir = 1,gams%ndir_transp
1446      do idir = 1,gams%ndir_transp
1447        ii = idir + gams%ndir_transp*(jdir-1)
1448        ! Diagonalize gamma matrix at qpoint (complex matrix).
1449        ! MJV NOTE: gam_now is recast implicitly here to matrix
1450        gam_now  = gams%vals_in_qibz(:,ii,:,:,iq_ibz,spin)
1451 
1452        call ZGEMM('N','N',natom3,natom3,natom3,cone,gam_now, natom3,pheigvec,natom3,czero,tmp_gam1,natom3)
1453        call ZGEMM('C','N',natom3,natom3,natom3,cone,pheigvec,natom3,tmp_gam1 ,natom3,czero,tmp_gam2,natom3)
1454 
1455        diagerr = zero
1456        do nu2=1,natom3
1457          gamma_in_ph(idir,jdir,nu2) = tmp_gam2(1,nu2,nu2)
1458 
1459          do nu1=1,nu2-1
1460            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1461          end do
1462          do nu1=nu2+1,natom3
1463            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1464          end do
1465          diagerr = diagerr + abs(tmp_gam2(2,nu2,nu2))
1466        end do
1467 
1468        if (diagerr > tol12) then
1469          write (msg,'(a,es14.6)')'Numerical error in diagonalization of gamma with phon eigenvectors: ',diagerr
1470          MSG_WARNING(msg)
1471        end if
1472 
1473        gam_now  = gams%vals_out_qibz(:,ii,:,:,iq_ibz,spin)
1474 
1475        call ZGEMM('N','N',natom3,natom3,natom3,cone,gam_now, natom3,pheigvec,natom3,czero,tmp_gam1,natom3)
1476        call ZGEMM('C','N',natom3,natom3,natom3,cone,pheigvec,natom3,tmp_gam1 ,natom3,czero,tmp_gam2,natom3)
1477 
1478        diagerr = zero
1479        do nu2=1,natom3
1480          gamma_out_ph(idir,jdir,nu2) = tmp_gam2(1,nu2,nu2)
1481 
1482          do nu1=1,nu2-1
1483            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1484          end do
1485          do nu1=nu2+1,natom3
1486            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1487          end do
1488          diagerr = diagerr + abs(tmp_gam2(2,nu2,nu2))
1489        end do
1490 
1491        if (diagerr > tol12) then
1492          write (msg,'(a,es14.6)')'Numerical error in diagonalization of gamma with phon eigenvectors: ',diagerr
1493          MSG_WARNING(msg)
1494        end if
1495      end do ! idir
1496    end do ! jdir
1497 
1498  case default
1499    MSG_BUG(sjoin("Wrong value for eph_scalprod:",itoa(gams%eph_scalprod)))
1500  end select
1501 
1502  ! Compute lambda
1503  do jdir = 1,gams%ndir_transp
1504    do idir = 1,gams%ndir_transp
1505      ! TODO : check this - looks like a factor of 2 wrt the inline documentation!
1506      !spinfact should be 1 for a normal non sppol calculation without spinorbit
1507      !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1508      !for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1509      spinfact = two/(gams%nsppol*gams%nspinor)
1510 
1511      do nu1=1,gams%natom3
1512        gamma_in_ph(idir,jdir,nu1)  =  gamma_in_ph(idir,jdir,nu1)  * pi * spinfact
1513        gamma_out_ph(idir,jdir,nu1) =  gamma_out_ph(idir,jdir,nu1) * pi * spinfact
1514        lambda_in_ph(idir,jdir,nu1)  = zero
1515        lambda_out_ph(idir,jdir,nu1) = zero
1516        if (abs(phfrq(nu1)) > EPH_WTOL) then
1517          lambda_in_ph(idir,jdir,nu1)  = gamma_in_ph(idir,jdir,nu1)  / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1518          lambda_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1519        end if
1520      end do
1521    end do ! idir
1522  end do ! jdir
1523 
1524 
1525 end subroutine phgamma_vv_eval_qibz

m_phgamma/phgamma_vv_interp [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_interp

FUNCTION

  Interpolate the linewidths at a given q-point.

INPUTS

  gams<phgamma_t>
  cryst<crystal_t>=crystalline structure.
  ifc<ifc_type>=Interatomic force constants.
  spin=Spin index
  qpt(3)=q-point in reduced coordinates
  gamma_ph(3*natom)=Phonon linewidths
  displ_cart(2,3,cryst%natom,3*cryst%natom)=Phonon displacement in cartesian coordinates.

OUTPUT

  gamma_ph(gams%natom3)=Interpolated Phonon linewidths.
  lamda_ph(3*natom)=Lambda coefficients for the different phonon modes.
  phfrq(3*natom)=phonon frequencies at current q

PARENTS

      m_phgamma

CHILDREN

SOURCE

1558 subroutine phgamma_vv_interp(gams,cryst,ifc,spin,qpt,phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph)
1559 
1560 
1561 !This section has been created automatically by the script Abilint (TD).
1562 !Do not modify the following lines by hand.
1563 #undef ABI_FUNC
1564 #define ABI_FUNC 'phgamma_vv_interp'
1565 !End of the abilint section
1566 
1567  implicit none
1568 
1569 !Arguments ------------------------------------
1570 !scalars
1571  integer,intent(in) :: spin
1572  type(phgamma_t),intent(inout) :: gams
1573  type(crystal_t),intent(in) :: cryst
1574  type(ifc_type),intent(in) :: ifc
1575 !arrays
1576  real(dp),intent(in) :: qpt(3)
1577  real(dp),intent(out) :: phfrq(gams%natom3)
1578  real(dp),intent(out) :: gamma_in_ph(3,3,gams%natom3),lambda_in_ph(3,3,gams%natom3)
1579  real(dp),intent(out) :: gamma_out_ph(3,3,gams%natom3),lambda_out_ph(3,3,gams%natom3)
1580 
1581 !Local variables-------------------------------
1582 !scalars
1583  integer,parameter :: qtor0=0
1584  integer :: natom3,nu1,nu2
1585  integer :: idir,jdir,ii
1586  real(dp) :: diagerr,spinfact
1587  character(len=500) :: msg
1588  !arrays
1589  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
1590  real(dp) :: displ_red(2,gams%natom3,gams%natom3),img(gams%natom3)
1591  real(dp) :: gam_in_now(2,3,3,gams%natom3**2)
1592  real(dp) :: gam_out_now(2,3,3,gams%natom3**2)
1593  real(dp) :: tmp_gam1(2,gams%natom3,gams%natom3),tmp_gam2(2,gams%natom3,gams%natom3)
1594  real(dp) :: pheigvec(2*gams%natom3**2)
1595  real(dp),allocatable :: coskr(:,:),sinkr(:,:)
1596 
1597 ! *************************************************************************
1598 
1599  ! Compute internal tables used for Fourier interpolation.
1600  if (.not.allocated(gams%vals_in_bz)) call phgamma_vv_interp_setup(gams,cryst,"INIT")
1601 
1602  !@phgamma_t
1603  natom3 = gams%natom3
1604 
1605  ! Get phonon frequencies and eigenvectors.
1606  call ifc_fourq(ifc,cryst,qpt,phfrq,displ_cart,out_eigvec=pheigvec, out_displ_red=displ_red)
1607 
1608  ! Taken from mkph_linwid
1609  ! This reduced version of ftgkk supposes the kpoints have been integrated
1610  ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt.
1611  ABI_MALLOC(coskr, (1,gams%nrpt))
1612  ABI_MALLOC(sinkr, (1,gams%nrpt))
1613  ! TODO: This is not optimal
1614  call ftgam_init(gams%gprim, 1, gams%nrpt, qpt, gams%rpt, coskr, sinkr)
1615 
1616  do idir=1,gams%ndir_transp
1617    do jdir=1,gams%ndir_transp
1618      ii = idir+gams%ndir_transp*(jdir-1)
1619 
1620      call ftgam(gams%wghatm,gam_in_now,gams%vals_in_rpt(:,ii,:,:,spin),gams%natom,1,gams%nrpt,qtor0,coskr, sinkr)
1621      call ftgam(gams%wghatm,gam_out_now,gams%vals_out_rpt(:,ii,:,:,spin),gams%natom,1,gams%nrpt,qtor0,coskr, sinkr)
1622 
1623      ! This call is not executed in elphon!
1624      ! TODO: needs to take into account the matrix nature of _in_ and _out_ gammas
1625      if (gams%symgamma == 1) then
1626        call tgamma_symm(cryst,qpt,gam_in_now)
1627        call tgamma_symm(cryst,qpt,gam_out_now)
1628      end if
1629 
1630      select case (gams%eph_scalprod)
1631      case (0)
1632        ! If the matrices do not contain the scalar product with the displ_cart vectors yet do it now.
1633        tmp_gam2 = reshape (gam_in_now, [2,natom3,natom3])
1634        call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
1635 
1636        do nu1=1,natom3
1637          gamma_in_ph(idir,jdir,nu1) = tmp_gam1(1, nu1, nu1)
1638          img(nu1) = tmp_gam1(2, nu1, nu1)
1639          if (abs(img(nu1)) > tol8) then
1640            write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
1641            MSG_WARNING(msg)
1642          end if
1643        end do
1644 
1645        tmp_gam2 = reshape (gam_out_now, [2,natom3,natom3])
1646        call gam_mult_displ(natom3, displ_red, tmp_gam2, tmp_gam1)
1647 
1648        do nu1=1,natom3
1649          gamma_out_ph(idir,jdir,nu1) = tmp_gam1(1, nu1, nu1)
1650          img(nu1) = tmp_gam1(2, nu1, nu1)
1651          if (abs(img(nu1)) > tol8) then
1652            write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch= ',nu1,', img= ',img(nu1)
1653            MSG_WARNING(msg)
1654          end if
1655        end do
1656 
1657      case (1)
1658        ! Diagonalize gamma matrix at qpoint (complex matrix).
1659        ! MJV NOTE: gam_now is recast implicitly here to matrix
1660        call ZGEMM('N','N',natom3,natom3,natom3,cone,gam_in_now,natom3,pheigvec,natom3,czero,tmp_gam1,natom3)
1661        call ZGEMM('C','N',natom3,natom3,natom3,cone,pheigvec,  natom3,tmp_gam1 ,natom3,czero,tmp_gam2,natom3)
1662 
1663        diagerr = zero
1664        do nu2=1,natom3
1665          gamma_in_ph(idir,jdir,nu2) = tmp_gam2(1,nu2,nu2)
1666 
1667          do nu1=1,nu2-1
1668            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1669          end do
1670          do nu1=nu2+1,natom3
1671            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1672          end do
1673          diagerr = diagerr + abs(tmp_gam2(2,nu2,nu2))
1674        end do
1675 
1676        if (diagerr > tol12) then
1677          write (msg,'(a,es14.6)')'Numerical error in diagonalization of gamma with phon eigenvectors: ',diagerr
1678          MSG_WARNING(msg)
1679        end if
1680 
1681        call ZGEMM('N','N',natom3,natom3,natom3,cone,gam_out_now,natom3,pheigvec,natom3,czero,tmp_gam1,natom3)
1682        call ZGEMM('C','N',natom3,natom3,natom3,cone,pheigvec,  natom3,tmp_gam1 ,natom3,czero,tmp_gam2,natom3)
1683 
1684        diagerr = zero
1685        do nu2=1,natom3
1686          gamma_out_ph(idir,jdir,nu2) = tmp_gam2(1,nu2,nu2)
1687 
1688          do nu1=1,nu2-1
1689            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1690          end do
1691          do nu1=nu2+1,natom3
1692            diagerr = diagerr + abs(tmp_gam2(1,nu1,nu2))+abs(tmp_gam2(2,nu1,nu2))
1693          end do
1694          diagerr = diagerr + abs(tmp_gam2(2,nu2,nu2))
1695        end do
1696 
1697        if (diagerr > tol12) then
1698          write (msg,'(a,es14.6)')'Numerical error in diagonalization of gamma with phon eigenvectors: ',diagerr
1699          MSG_WARNING(msg)
1700        end if
1701 
1702      case default
1703        MSG_BUG(sjoin("Wrong value for eph_scalprod:",itoa(gams%eph_scalprod)))
1704      end select
1705    end do
1706  end do
1707 
1708  ABI_FREE(coskr)
1709  ABI_FREE(sinkr)
1710 
1711  ! Compute lambda
1712  !spinfact should be 1 for a normal non sppol calculation without spinorbit
1713  !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1714  !for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1715  spinfact = two/(gams%nsppol*gams%nspinor)
1716 
1717  ! Compute lambda
1718  do nu1=1,gams%natom3
1719    do idir=1,gams%ndir_transp
1720      do jdir=1,gams%ndir_transp
1721        gamma_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) * pi * spinfact
1722        gamma_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) * pi * spinfact
1723        lambda_in_ph(idir,jdir,nu1) = zero
1724        lambda_out_ph(idir,jdir,nu1) = zero
1725        if (abs(phfrq(nu1)) > EPH_WTOL) then
1726          lambda_in_ph(idir,jdir,nu1)  = gamma_in_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1727          lambda_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1728        end if
1729      end do
1730    end do
1731  end do
1732 
1733 end subroutine phgamma_vv_interp

m_phgamma/phgamma_vv_interp_setup [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 phgamma_vv_interp_setup

FUNCTION

  This routines performs the (allocation|deallocation) of the internal tables
  used to interpolate the vv_linewidths in q-space

INPUTS

  action =
    "INIT" to allocate and compute the internal tables (default)
    "FREE" to deallocate the internal tables.

SIDE EFFECTS

  gams<phgamma_t>= gams%vals_in_bz and gams%vals_in_rpt, etc... depending on action.

PARENTS

      m_phgamma

CHILDREN

SOURCE

1761 subroutine phgamma_vv_interp_setup(gams,cryst,action)
1762 
1763 
1764 !This section has been created automatically by the script Abilint (TD).
1765 !Do not modify the following lines by hand.
1766 #undef ABI_FUNC
1767 #define ABI_FUNC 'phgamma_vv_interp_setup'
1768 !End of the abilint section
1769 
1770  implicit none
1771 
1772 !Arguments ------------------------------------
1773 !scalars
1774  character(len=*),intent(in) :: action
1775  type(phgamma_t),intent(inout) :: gams
1776  type(crystal_t),intent(in) :: cryst
1777 
1778 !Local variables-------------------------------
1779 !scalars
1780  integer,parameter :: qtor1=1
1781  integer :: iq_bz,iq_ibz,isq_bz,spin,isym,ierr
1782  integer :: ii, idir, jdir
1783  !character(len=500) :: msg
1784  type(kptrank_type) :: qrank
1785 !arrays
1786  integer :: g0(3)
1787  integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:)
1788  real(dp) :: qirr(3),tmp_qpt(3)
1789  real(dp),allocatable :: coskr(:,:),sinkr(:,:)
1790 
1791 ! *************************************************************************
1792 
1793  select case (toupper(action))
1794 
1795  case ("INIT")
1796    if (.not.allocated(gams%vals_in_bz)) then
1797      ABI_STAT_MALLOC(gams%vals_in_bz,(2,gams%ndir_transp**2,gams%natom3**2,gams%nqbz,gams%nsppol), ierr)
1798      ABI_CHECK(ierr==0, 'out of memory in gams%vals_in_bz')
1799      gams%vals_in_bz = zero
1800      ABI_STAT_MALLOC(gams%vals_out_bz,(2,gams%ndir_transp**2,gams%natom3**2,gams%nqbz,gams%nsppol), ierr)
1801      ABI_CHECK(ierr==0, 'out of memory in gams%vals_out_bz')
1802      gams%vals_out_bz = zero
1803 
1804      ! Build tables needed by complete_gamma.
1805      call mkkptrank(gams%qbz,gams%nqbz,qrank)
1806 
1807      ! Compute index of IBZ q-point in the BZ array
1808      ABI_CALLOC(qirredtofull,(gams%nqibz))
1809 
1810      do iq_ibz=1,gams%nqibz
1811        qirr = gams%qibz(:,iq_ibz)
1812        iq_bz = kptrank_index(qrank, qirr)
1813        if (iq_bz /= -1) then
1814          ABI_CHECK(isamek(qirr,gams%qbz(:,iq_bz),g0), "isamek")
1815          qirredtofull(iq_ibz) = iq_bz
1816        else
1817          MSG_ERROR(sjoin("Full BZ does not contain IBZ q-point:", ktoa(qirr)))
1818        end if
1819      end do
1820 
1821      ! Build qpttoqpt table. See also mkqptequiv
1822      ABI_MALLOC(qpttoqpt,(2,cryst%nsym,gams%nqbz))
1823      qpttoqpt = -1
1824      do iq_bz=1,gams%nqbz
1825        do isym=1,cryst%nsym
1826          tmp_qpt = matmul(cryst%symrec(:,:,isym), gams%qbz(:,iq_bz))
1827 
1828          isq_bz = kptrank_index(qrank, tmp_qpt)
1829          if (isq_bz == -1) then
1830            MSG_ERROR("looks like no kpoint equiv to q by symmetry without time reversal!")
1831          end if
1832          qpttoqpt(1,isym,isq_bz) = iq_bz
1833 
1834          ! q --> -q
1835          tmp_qpt = -tmp_qpt
1836          isq_bz = kptrank_index(qrank, tmp_qpt)
1837          if (isq_bz == -1) then
1838            MSG_ERROR("looks like no kpoint equiv to q by symmetry with time reversal!")
1839          end if
1840          qpttoqpt(2,isym,isq_bz) = iq_bz
1841        end do
1842      end do
1843      call destroy_kptrank(qrank)
1844 
1845      ! Fill BZ array with IBZ data.
1846      do spin=1,gams%nsppol
1847        do iq_ibz=1,gams%nqibz
1848          iq_bz = qirredtofull(iq_ibz)
1849          gams%vals_in_bz(:,:,:,iq_bz,spin) = reshape(gams%vals_in_qibz(:,:,:,:,iq_ibz,spin), &
1850 &             [2,gams%ndir_transp**2,gams%natom3**2])
1851          gams%vals_out_bz(:,:,:,iq_bz,spin) = reshape(gams%vals_out_qibz(:,:,:,:,iq_ibz,spin), &
1852 &             [2,gams%ndir_transp**2,gams%natom3**2])
1853        end do
1854      end do
1855 
1856      ! Complete vals_bz in the full BZ.
1857 !TODO!!! rotate the vv in and out matrices, according to the symmetry operation, instead of just copying them
1858 !     call complete_gamma_vv(cryst,gams%natom3,gams%nsppol,gams%nqibz,gams%nqbz,&
1859 !       gams%eph_scalprod,qirredtofull,qpttoqpt,gams%vals_in_bz(:,:,:,:,spin))
1860 !     call complete_gamma_vv(cryst,gams%natom3,gams%nsppol,gams%nqibz,gams%nqbz,&
1861 !       gams%eph_scalprod,qirredtofull,qpttoqpt,gams%vals_out_bz(:,:,:,:,spin))
1862 
1863 ! TODO: replace the above with these calls from anaddb
1864 !   call complete_gamma_tr(cryst,elph_ds%ep_scalprod,elph_ds%nbranch,elph_ds%nqptirred,&
1865 !&   elph_ds%nqpt_full,elph_ds%nsppol,elph_tr_ds%gamma_qpt_trout,elph_ds%qirredtofull,qpttoqpt)
1866 
1867 
1868      ! TODO: idem for vv_vals 3x3 matrices
1869 
1870      ABI_FREE(qirredtofull)
1871      ABI_FREE(qpttoqpt)
1872 
1873 !     ! This call is not executed in elphon!
1874 !     if (gams%symgamma == 1) then
1875 !       do spin=1,gams%nsppol
1876 !         do iq_bz=1,gams%nqbz
1877 !           call tgamma_symm_vv(cryst,gams%qbz(:,iq_bz),gams%vals_in_bz(:,:,:,iq_bz,spin))
1878 !           call tgamma_symm_vv(cryst,gams%qbz(:,iq_bz),gams%vals_out_bz(:,:,:,iq_bz,spin))
1879 !         end do
1880 !       end do
1881 !     end if
1882    end if
1883 
1884    ! Now FT to real space
1885    ! NOTE: gprim (not gprimd) is used for all FT interpolations,
1886    ! to be consistent with the dimensions of the rpt, which come from anaddb.
1887    ! TODO: this is needed only if FT is used, not when the linear interpolation is employed.
1888    if (.not.allocated(gams%vals_in_rpt)) then
1889      ABI_STAT_MALLOC(gams%vals_in_rpt,(2,gams%ndir_transp**2,gams%natom3**2,gams%nrpt,gams%nsppol), ierr)
1890      ABI_CHECK(ierr==0, 'out of memory in gams%vals_in_rpt')
1891      gams%vals_in_rpt = zero
1892      ABI_STAT_MALLOC(gams%vals_out_rpt,(2,gams%ndir_transp**2,gams%natom3**2,gams%nrpt,gams%nsppol), ierr)
1893      ABI_CHECK(ierr==0, 'out of memory in gams%vals_out_rpt')
1894      gams%vals_out_rpt = zero
1895 
1896      ! q --> r
1897      ABI_MALLOC(coskr, (gams%nqbz,gams%nrpt))
1898      ABI_MALLOC(sinkr, (gams%nqbz,gams%nrpt))
1899      call ftgam_init(gams%gprim, gams%nqbz, gams%nrpt, gams%qbz, gams%rpt, coskr, sinkr)
1900 
1901      do spin=1,gams%nsppol
1902        do idir=1,gams%ndir_transp
1903          do jdir=1,gams%ndir_transp
1904            ii = idir+gams%ndir_transp*(jdir-1)
1905 ! TODO: this is no contiguous in memory and will be slow. Make adapted ftgam?
1906            call ftgam(gams%wghatm,gams%vals_in_bz(:,ii,:,:,spin),gams%vals_in_rpt(:,ii,:,:,spin),gams%natom,gams%nqbz,&
1907              gams%nrpt,qtor1, coskr, sinkr)
1908            call ftgam(gams%wghatm,gams%vals_out_bz(:,ii,:,:,spin),gams%vals_out_rpt(:,ii,:,:,spin),gams%natom,gams%nqbz,&
1909              gams%nrpt,qtor1, coskr, sinkr)
1910          end do
1911        end do
1912      end do
1913 
1914 
1915      ABI_FREE(coskr)
1916      ABI_FREE(sinkr)
1917    end if
1918 
1919  case ("FREE")
1920    if (allocated(gams%vals_in_bz)) then
1921      ABI_FREE(gams%vals_in_bz)
1922    end if
1923    if (allocated(gams%vals_out_bz)) then
1924      ABI_FREE(gams%vals_out_bz)
1925    end if
1926    if (allocated(gams%vals_in_rpt)) then
1927      ABI_FREE(gams%vals_in_rpt)
1928    end if
1929    if (allocated(gams%vals_out_rpt)) then
1930      ABI_FREE(gams%vals_out_rpt)
1931    end if
1932 
1933  case default
1934    MSG_BUG(sjoin("Wrong action:",action))
1935  end select
1936 
1937 end subroutine phgamma_vv_interp_setup

m_phgamma/tgamma_symm [ Functions ]

[ Top ] [ m_phgamma ] [ Functions ]

NAME

 tgamma_symm

FUNCTION

  Symmetrize the tgamma matrix

INPUTS

 qpt(3)=phonon wavevector in reduced coordinates.
 cryst<crystal_t>=Crystalline structure.

PARENTS

      m_phgamma

CHILDREN

SOURCE

722 subroutine tgamma_symm(cryst,qpt,tgamma)
723 
724 
725 !This section has been created automatically by the script Abilint (TD).
726 !Do not modify the following lines by hand.
727 #undef ABI_FUNC
728 #define ABI_FUNC 'tgamma_symm'
729 !End of the abilint section
730 
731  implicit none
732 
733 !Arguments ------------------------------------
734 !scalars
735  type(crystal_t),intent(in) :: cryst
736 !arrays
737  real(dp),intent(in) :: qpt(3)
738  real(dp),intent(inout) :: tgamma(2,3*cryst%natom,3*cryst%natom)
739 
740 !Local variables-------------------------------
741 !scalars
742  integer :: ii,natom3,k
743 !arrays
744  real(dp) :: tgcart(2,3*cryst%natom,3*cryst%natom)
745  real(dp) :: umat(2,3*cryst%natom,3*cryst%natom),tmp_mat(2,3*cryst%natom,3*cryst%natom)
746 
747 ! *********************************************************************
748 
749  ! Build U matrix.
750  umat = zero; k = 1
751  do ii=1,cryst%natom
752    umat(1,k:k+2, k:k+2) = cryst%gprimd
753    k = k + 3
754  end do
755 
756  natom3 = 3 * cryst%natom
757 
758  ! Reduced --> Cartesian
759  call zgemm('N','N',natom3,natom3,natom3,cone,tgamma,natom3,umat,natom3,czero,tmp_mat,natom3)
760  call zgemm('T','N',natom3,natom3,natom3,cone,umat,natom3,tmp_mat,natom3,czero,tgcart,natom3)
761 
762  ! Make the matrix hermitian
763  call mkherm(tgcart,3*cryst%natom)
764 
765  ! Symmetrize tgamma matrix.
766  call symdyma(tgcart,cryst%indsym,cryst%natom,cryst%nsym,qpt,cryst%rprimd,cryst%symrel,cryst%symafm)
767 
768  umat = zero; k = 1
769  do ii=0,cryst%natom-1
770    umat(1,k:k+2, k:k+2) = cryst%rprimd
771    k = k + 3
772  end do
773 
774  ! Reduced --> Cartesian
775  call zgemm('N','N',natom3,natom3,natom3,cone,tgcart,natom3,umat,natom3,czero,tmp_mat,natom3)
776  call zgemm('T','N',natom3,natom3,natom3,cone,umat,natom3,tmp_mat,natom3,czero,tgamma,natom3)
777 
778 end subroutine tgamma_symm