diff --git a/src/libolle/mpi_wrappers.f90 b/src/libolle/mpi_wrappers.f90 index 671e06ca..33258a33 100644 --- a/src/libolle/mpi_wrappers.f90 +++ b/src/libolle/mpi_wrappers.f90 @@ -54,8 +54,8 @@ module mpi_wrappers !> mpi allreduce generic :: allreduce=>allreduce_int,allreduce_real,allreduce_1d_real,allreduce_2d_real,allreduce_3d_real,& allreduce_4d_real,allreduce_5d_real,allreduce_6d_real,allreduce_1d_int,allreduce_2d_int,allreduce_3d_int,& - allreduce_2d_complex,allreduce_3d_complex,allreduce_4d_complex,allreduce_5d_complex,allreduce_1d_logical,& - allreduce_int8 + allreduce_2d_complex,allreduce_3d_complex,allreduce_4d_complex,allreduce_5d_complex,allreduce_6d_complex,& + allreduce_1d_logical,allreduce_int8 procedure, private :: allreduce_int procedure, private :: allreduce_1d_int procedure, private :: allreduce_2d_int @@ -71,6 +71,7 @@ module mpi_wrappers procedure, private :: allreduce_3d_complex procedure, private :: allreduce_4d_complex procedure, private :: allreduce_5d_complex + procedure, private :: allreduce_6d_complex procedure, private :: allreduce_1d_logical procedure, private :: allreduce_int8 !> mpi broadcast @@ -1442,6 +1443,42 @@ subroutine allreduce_5d_complex(mw,operation,x,y,filename,linenumber) endif endif end subroutine +subroutine allreduce_6d_complex(mw,operation,x,y,filename,linenumber) + !> MPI helper + class(lo_mpi_helper), intent(inout) :: mw + !> what operation to do + character(len=*), intent(in) :: operation + !> array to allreduce + complex(r8), dimension(:,:,:,:,:,:), intent(inout) :: x + !> array to allreduce to, if omitted default to in-place + complex(r8), dimension(:,:,:,:,:,:), intent(out), optional :: y + !> filename we call from + character(len=*), intent(in), optional :: filename + !> line number we call from + integer, intent(in), optional :: linenumber + + integer :: mpiop + ! Fetch proper operation code + mpiop=mpi_operation_code(operation) + ! Do the actual communication + if ( present(y) ) then + if ( size(x) .ne. size(y) ) then + call lo_stop_gracefully(['mpi_allreduce inconsistent array sizes'],lo_exitcode_mpi,filename,linenumber,mw%comm) + endif + y=0.0_r8 + call mpi_allreduce(x,y,size(x),MPI_DOUBLE_COMPLEX,mpiop,mw%comm,mw%error) + else + call mpi_allreduce(MPI_IN_PLACE,x,size(x),MPI_DOUBLE_COMPLEX,mpiop,mw%comm,mw%error) + endif + ! Check that things went ok + if ( mw%error .ne. 0 ) then + if ( present(filename) .and. present(linenumber) ) then + call lo_stop_gracefully(['mpi_allreduce exit code '//tochar(mw%error)],lo_exitcode_mpi,filename,linenumber,mw%comm) + else + call lo_stop_gracefully(['mpi_allreduce exit code '//tochar(mw%error)],lo_exitcode_mpi,__FILE__,__LINE__,mw%comm) + endif + endif +end subroutine subroutine allreduce_1d_logical(mw,operation,i,j,filename,linenumber) !> MPI helper class(lo_mpi_helper), intent(inout) :: mw diff --git a/src/lineshape/dielscatter.f90 b/src/lineshape/dielscatter.f90 index dd284313..101a15f5 100644 --- a/src/lineshape/dielscatter.f90 +++ b/src/lineshape/dielscatter.f90 @@ -248,7 +248,7 @@ subroutine generate(dir, wp, di, p, qp, dr, fc, fct, se, isf, opts, tmr, mw, mem ! fetch probing energy and real and imaginary part of the self-energy Z = se%energy_axis(ie) dl = se%re_3ph(ie, imode) + se%re_4ph(ie, imode) - gm = se%im_3ph(ie, imode) + se%im_iso(ie, imode) + gm = se%im_3ph(ie, imode) + se%im_4ph(ie, imode) + se%im_iso(ie, imode) ! Zeroth order IR, selfenergy=0 c0 = diagonal_greensfunction(om, Z, 0.0_r8, smallsmearing) @@ -346,7 +346,7 @@ subroutine generate(dir, wp, di, p, qp, dr, fc, fct, se, isf, opts, tmr, mw, mem ! fetch probing energy and real and imaginary part of the self-energy Z = se%energy_axis(ie) dl = se%re_3ph(ie, b1) + se%re_4ph(ie, b1) - gm = se%im_3ph(ie, b1) + se%im_iso(ie, b1) + gm = se%im_3ph(ie, b1) + se%im_4ph(ie, b1) + se%im_iso(ie, b1) grfun(ie, b1) = aimag(diagonal_greensfunction(om1, Z, dl, gm)) end do end do diff --git a/src/lineshape/io.f90 b/src/lineshape/io.f90 index 73df05ee..17e3d4c9 100644 --- a/src/lineshape/io.f90 +++ b/src/lineshape/io.f90 @@ -161,6 +161,7 @@ subroutine write_lineshape_to_hdf5(se, p, qpt, wp, dir, di, qp, enhet, temperatu call h5%store_data(se%energy_axis*unitfactor, h5%group_id, 'frequency', enhet=trim(unitname)) call h5%store_data(se%im_iso*unitfactor, h5%group_id, 'imaginary_isotope_selfenergy', enhet=trim(unitname)) call h5%store_data(se%im_3ph*unitfactor, h5%group_id, 'imaginary_threephonon_selfenergy', enhet=trim(unitname)) + call h5%store_data(se%im_4ph*unitfactor, h5%group_id, 'imaginary_fourphonon_selfenergy', enhet=trim(unitname)) call h5%store_data(se%re_3ph*unitfactor, h5%group_id, 'real_threephonon_selfenergy', enhet=trim(unitname)) call h5%store_data(se%re_4ph*unitfactor, h5%group_id, 'real_fourphonon_selfenergy', enhet=trim(unitname)) @@ -179,7 +180,7 @@ subroutine write_lineshape_to_hdf5(se, p, qpt, wp, dir, di, qp, enhet, temperatu ! Imaginary part of Green's function do i = 1, se%n_mode if (wp%omega(i) .gt. lo_freqtol*10) then - call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, i) + se%im_iso(:, i), se%re_3ph(:, i) + se%re_4ph(:, i), wp%omega(i), dr0(:, i)) + call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, i) + se%im_4ph(:, i) + se%im_iso(:, i), se%re_3ph(:, i) + se%re_4ph(:, i), wp%omega(i), dr0(:, i)) else dr0(:, i) = 0.0_r8 end if diff --git a/src/lineshape/main.f90 b/src/lineshape/main.f90 index 268ae5c6..87a9b63b 100644 --- a/src/lineshape/main.f90 +++ b/src/lineshape/main.f90 @@ -72,7 +72,7 @@ program lineshape call fct%readfromfile(uc, 'infile.forceconstant_thirdorder') if (mw%talk) write (*, *) '... read third order forceconstant' end if - if (opts%fourthorder) then + if (opts%fourthorder_real .or. opts%fourthorder) then call fcf%readfromfile(uc, 'infile.forceconstant_fourthorder') if (mw%talk) write (*, *) '... read fourth order forceconstant' end if diff --git a/src/lineshape/options.f90 b/src/lineshape/options.f90 index 2bd48b32..a978d468 100644 --- a/src/lineshape/options.f90 +++ b/src/lineshape/options.f90 @@ -20,6 +20,7 @@ module options logical :: isotopescattering = .false. logical :: thirdorder = .false. logical :: fourthorder = .false. + logical :: fourthorder_real = .false. logical :: mct = .false. logical :: slightsmearing = .false. integer :: integrationtype = -lo_hugeint @@ -126,9 +127,12 @@ subroutine parse(opts) help='Switch of three-phonon scattering', & required=.false., act='store_true', def='.false.', error=lo_status) if (lo_status .ne. 0) stop - call cli%add(switch='--fourthorder', & + call cli%add(switch='--fourthorder_real', & help='Consider four-phonon contributions to the real part of the self-energy.', hidden=.true., & required=.false., act='store_true', def='.false.', error=lo_status) + call cli%add(switch='--fourthorder', & + help='Consider four-phonon contributions to the dynamic part of the self-energy.', hidden=.true., & + required=.false., act='store_true', def='.false.', error=lo_status) if (lo_status .ne. 0) stop call cli%add(switch='--remove_static_selfenergy', & help='Remove the static contribution to the self-energy as in the mode-coupling theory.', & @@ -237,6 +241,7 @@ subroutine parse(opts) opts%isotopescattering = .not. dumlog call cli%get(switch='--no_thirdorder_scattering', val=dumlog) opts%thirdorder = .not. dumlog + call cli%get(switch='--fourthorder_real', val=opts%fourthorder_real) call cli%get(switch='--fourthorder', val=opts%fourthorder) call cli%get(switch='--remove_static_selfenergy', val=opts%mct) call cli%get(switch='--nondiagonal', val=dumlog) @@ -336,7 +341,7 @@ subroutine parse(opts) end if ! If we are in the mode-coupling approach, the real part four phonon makes no sense - if (opts%mct .and. opts%fourthorder) then + if (opts%mct .and. opts%fourthorder_real) then write(*, *) 'There is no real part from the fourth order in the mode-coupling theory' stop end if diff --git a/src/lineshape/phonondamping.f90 b/src/lineshape/phonondamping.f90 index ceb2c88c..8f9a444a 100644 --- a/src/lineshape/phonondamping.f90 +++ b/src/lineshape/phonondamping.f90 @@ -41,7 +41,7 @@ module phonondamping ! energy axis for the self-energy real(r8), dimension(:), allocatable :: energy_axis ! Imaginary part of self energy (energy,mode) - real(r8), dimension(:, :), allocatable :: im_3ph, im_iso + real(r8), dimension(:, :), allocatable :: im_3ph, im_iso, im_4ph ! Real part of self energy (energy,mode) real(r8), dimension(:, :), allocatable :: re_3ph, re_4ph !> Direction of probe @@ -54,6 +54,7 @@ module phonondamping logical :: isotope_scattering = .false. logical :: thirdorder_scattering = .false. logical :: fourthorder_scattering = .false. + logical :: fourthorder_real = .false. logical :: diagonal = .false. logical :: skipsym = .false. ! auxiliary information about lineshapes: @@ -75,7 +76,8 @@ module phonondamping ! Prefactors, to make sure they are the same all the time. real(r8), parameter :: threephonon_prefactor = lo_pi/16.0_r8 -real(r8), parameter :: fourphonon_prefactor = 1.0_r8/8.0_r8 +real(r8), parameter :: fourphonon_real_prefactor = 1.0_r8/8.0_r8 +real(r8), parameter :: fourphonon_imag_prefactor = lo_pi/96.0_r8 real(r8), parameter :: isotope_prefactor = lo_pi/4.0_r8 ! dos interfaces @@ -199,6 +201,7 @@ function se_size_in_mem(se) result(mem) if (allocated(se%im_3ph)) mem = mem + storage_size(se%im_3ph)*size(se%im_3ph) if (allocated(se%im_iso)) mem = mem + storage_size(se%im_iso)*size(se%im_iso) if (allocated(se%re_3ph)) mem = mem + storage_size(se%re_3ph)*size(se%re_3ph) + if (allocated(se%im_4ph)) mem = mem + storage_size(se%im_4ph)*size(se%im_4ph) if (allocated(se%re_4ph)) mem = mem + storage_size(se%re_4ph)*size(se%re_4ph) if (allocated(se%xmid)) mem = mem + storage_size(se%xmid)*size(se%xmid) if (allocated(se%xlo)) mem = mem + storage_size(se%xlo)*size(se%xlo) diff --git a/src/lineshape/phonondamping_dos.f90 b/src/lineshape/phonondamping_dos.f90 index cb822a2f..4ef3e943 100644 --- a/src/lineshape/phonondamping_dos.f90 +++ b/src/lineshape/phonondamping_dos.f90 @@ -146,7 +146,7 @@ module subroutine get_intensity_as_dos(pd, qpd, drd, uc, fc, fct, fcf, ise, sf, if (drd%iq(iq)%omega(imode) .lt. lo_freqtol) cycle ! Get the spectral function buf1 = 0.0_r8 - call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, imode) + se%im_iso(:, imode), se%re_3ph(:, imode) + se%re_4ph(:, imode), drd%iq(iq)%omega(imode), buf1) + call evaluate_spectral_function(se%energy_axis, se%im_3ph(:, imode) + se%im_4ph(:, imode) + se%im_iso(:, imode), se%re_3ph(:, imode) + se%re_4ph(:, imode), drd%iq(iq)%omega(imode), buf1) buf1 = buf1/lo_trapezoid_integration(se%energy_axis, buf1) ! Stash the unsmeared spectral function @@ -173,7 +173,9 @@ module subroutine get_intensity_as_dos(pd, qpd, drd, uc, fc, fct, fcf, ise, sf, ! Store shifts and linewidth buf_shift_3rd(imode, iq) = lo_linear_interpolation(se%energy_axis, se%re_3ph(:, imode), drd%iq(iq)%omega(imode)) buf_shift_4th(imode, iq) = lo_linear_interpolation(se%energy_axis, se%re_4ph(:, imode), drd%iq(iq)%omega(imode)) - buf_linewidth(imode, iq) = lo_linear_interpolation(se%energy_axis, se%im_3ph(:, imode), drd%iq(iq)%omega(imode)) + lo_linear_interpolation(se%energy_axis, se%im_iso(:, imode), drd%iq(iq)%omega(imode)) + buf_linewidth(imode, iq) = lo_linear_interpolation(se%energy_axis, se%im_3ph(:, imode), drd%iq(iq)%omega(imode)) + & + lo_linear_interpolation(se%energy_axis, se%im_iso(:, imode), drd%iq(iq)%omega(imode)) + & + lo_linear_interpolation(se%energy_axis, se%im_4ph(:, imode), drd%iq(iq)%omega(imode)) end do ! Store spectral function diff --git a/src/lineshape/phonondamping_fourthorder.f90 b/src/lineshape/phonondamping_fourthorder.f90 index 8f3e9a30..7b4a393a 100644 --- a/src/lineshape/phonondamping_fourthorder.f90 +++ b/src/lineshape/phonondamping_fourthorder.f90 @@ -13,7 +13,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d type(lo_crystalstructure), intent(in) :: uc !> harmonic properties for q',q'',q'' type(lo_phonon_dispersions), intent(in) :: dr - !> third order force constant + !> fourth order force constant type(lo_forceconstant_fourthorder), intent(in) :: fcf !> temperature real(r8), intent(in) :: temperature @@ -238,7 +238,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d ctr = ctr + 1 if (mod(ctr, mw%n) .ne. mw%r) cycle - prefactor = qp%ip(iq)%integration_weight*fourphonon_prefactor + prefactor = qp%ip(iq)%integration_weight*fourphonon_real_prefactor ! pre-transform the matrix element call pretransform_phi4(fcf, qpoint%r, qp%ip(iq)%r, ptf) @@ -278,7 +278,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d ctr = ctr + 1 if (mod(ctr, mw%n) .ne. mw%r) cycle - prefactor = qp%ap(iq)%integration_weight*fourphonon_prefactor + prefactor = qp%ap(iq)%integration_weight*fourphonon_real_prefactor ! pre-transform the matrix element call pretransform_phi4(fcf, qpoint%r, qp%ap(iq)%r, ptf) @@ -432,7 +432,7 @@ subroutine fourphonon_selfenergy(qpoint, wp, gp, qp, uc, temperature, dr, fcf, d ! if ( minval(omega) .gt. omegathres ) then ! f0=real(fcf%scatteringamplitude(omega,egv,-qv2,-qv3,-qv4),r8) ! f0=f0*(1.0_r8+lo_planck(temperature,omega(3))+lo_planck(temperature,omega(4)) ) - ! f0=f0*qp%ap(q)%integration_weight*fourphonon_prefactor + ! f0=f0*qp%ap(q)%integration_weight*fourphonon_real_prefactor ! else ! f0=0.0_r8 ! endif diff --git a/src/lineshape/phonondamping_gaussian.f90 b/src/lineshape/phonondamping_gaussian.f90 index f9e2bc97..2a48a49f 100644 --- a/src/lineshape/phonondamping_gaussian.f90 +++ b/src/lineshape/phonondamping_gaussian.f90 @@ -372,3 +372,143 @@ subroutine threephonon_imaginary_selfenergy_gaussian(wp, se, sr, qp, dr, tempera call mem%deallocate(buf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) if (verbosity .gt. 0) call lo_progressbar(' ... threephonon imaginary selfenergy', dr%n_mode*qp%n_full_point, dr%n_mode*qp%n_full_point, walltime() - t0) end subroutine + + +!> fourphonon self energy gaussian +subroutine fourphonon_imaginary_selfenergy_gaussian(wp, se, sr, qp, dr, temperature, mw, mem, verbosity) + !> harmonic properties at this q-point + type(lo_phonon_dispersions_qpoint), intent(in) :: wp + !> self-energy + type(lo_phonon_selfenergy), intent(inout) :: se + !> scattering rates + type(lo_listofscatteringrates), intent(inout) :: sr + !> qpoint mesh + class(lo_qpoint_mesh), intent(in) :: qp + !> harmonic dispersions + type(lo_phonon_dispersions), intent(in) :: dr + !> temperature + real(r8), intent(in) :: temperature + !> MPI helper + type(lo_mpi_helper), intent(inout) :: mw + !> memory tracker + type(lo_mem_helper), intent(inout) :: mem + !> talk a lot? + integer, intent(in) :: verbosity + + real(r8), dimension(:), allocatable :: buf + real(r8) :: psisquare + real(r8) :: sigma, dom, om2, om3, om4, invf, s2, s3, s4 + real(r8) :: plf1, plf2, t0, pref + real(r8) :: n2, n3, n4 + integer :: q2, q3, ctr, i, ii, jj, b1, b2, b3, b4, ilo, ihi + + ! set the unit + t0 = walltime() + + invf = se%n_energy/se%energy_axis(se%n_energy) + + ctr = 0 + se%im_4ph = 0.0_r8 + call mem%allocate(buf, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + buf = 0.0_r8 + + qloopfull: do q2 = 1, qp%n_full_point + do q3=1, qp%n_full_point + ! Prefactor for this q-point + pref = fourphonon_imag_prefactor * qp%ap(q2)%integration_weight * qp%ap(q3)%integration_weight + do b1 = 1, dr%n_mode + do b2 = 1, dr%n_mode + ! Make it parallel + ctr = ctr + 1 + if (mod(ctr, mw%n) .ne. mw%r) cycle + + do b3 = 1, dr%n_mode + do b4 = 1, dr%n_mode + ! reset buffer + buf = 0.0_r8 + ilo = lo_hugeint + ihi = -lo_hugeint + ! Get the smearing parameter, in case it is adaptive + select case (se%integrationtype) + case (1) + s2 = dr%default_smearing(b2) + s3 = dr%default_smearing(b3) + s4 = dr%default_smearing(b4) + sigma = sqrt(s2**2 + s3**2 + s4**2) + case (2) + s2 = qp%adaptive_sigma(qp%ap(q2)%radius, dr%aq(q2)%vel(:, b2), dr%default_smearing(b2), se%smearing_prefactor) + s3 = qp%adaptive_sigma(qp%ap(q3)%radius, dr%aq(q3)%vel(:, b3), dr%default_smearing(b3), se%smearing_prefactor) + s4 = qp%adaptive_sigma(qp%ap(q3)%radius, sr%vel4(:, b4, q2, q3), dr%default_smearing(b4), se%smearing_prefactor) + sigma = sqrt(s2**2 + s3**2 + s4**2) + end select + ! Fetch frequencies + om2 = dr%aq(q2)%omega(b2) + om3 = dr%aq(q3)%omega(b3) + om4 = sr%omega4(b4, q2, q3) + + ! Define occupation numbers + n2 = lo_planck(temperature, om2) + n3 = lo_planck(temperature, om3) + n4 = lo_planck(temperature, om4) + + plf1 = (n2 + 1) * (n3 + 1) * (n4 + 1) + plf2 = n2 * n3 * n4 + dom = om2 + om3 + om4 + ! Delta(BigOmega - om2 - om3 - om4) + ii = max(floor((dom - 4 * sigma)*invf), 1) + jj = min(ceiling((dom + 4 * sigma) * invf), se%n_energy) + ilo = min(ilo, ii) + ihi = max(ihi, jj) + do i = ii, jj + buf(i) = buf(i) + (plf1 - plf2) * lo_gauss(se%energy_axis(i), dom, sigma) + end do + + plf1 = n2 * (n3 + 1) * (n4 + 1) + plf2 = (n2 + 1) * n3 * n4 + dom = -om2 + om3 + om4 + ! Delta(BigOmega - om2 + om3 + om4) + ii = max(floor((dom - 4 * sigma)*invf), 1) + jj = min(ceiling((dom + 4 * sigma)*invf), se%n_energy) + ilo = min(ilo, ii) + ihi = max(ihi, jj) + do i = ii, jj + buf(i) = buf(i) + 3 * (plf1 - plf2) * lo_gauss(se%energy_axis(i), dom, sigma) + end do + + ! Increment the self-energy + if (ilo .lt. ihi) then + psisquare = abs(sr%psi_4ph(b1, b2, b3, b4, q2, q3) * conjg(sr%psi_4ph(b1, b2, b3, b4, q2, q3))) * pref + se%im_4ph(ilo:ihi, b1) = se%im_4ph(ilo:ihi, b1) + buf(ilo:ihi)*psisquare + end if + end do + end do + end do + end do + if (verbosity .gt. 0) then + if (lo_trueNtimes(ctr, 127, qp%n_full_point*qp%n_full_point*dr%n_mode*dr%n_mode)) then + call lo_progressbar(' ... fourphonon imaginary selfenergy', ctr, dr%n_mode*dr%n_mode*qp%n_full_point*qp%n_full_point) + end if + end if + end do + end do qloopfull + + ! Sum it up + call mw%allreduce('sum', se%im_4ph) + + ! Fix degeneracies ? + do b1 = 1, dr%n_mode + buf = 0.0_r8 + do i = 1, wp%degeneracy(b1) + b2 = wp%degenmode(i, b1) + buf = buf + se%im_4ph(:, b2) + end do + buf = buf / real(wp%degeneracy(b1), r8) + do i = 1, wp%degeneracy(b1) + b2 = wp%degenmode(i, b1) + se%im_4ph(:, b2) = buf + end do + end do + call mem%deallocate(buf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + if (verbosity .gt. 0) call lo_progressbar(' ... fourphonon imaginary selfenergy', dr%n_mode*qp%n_full_point,& + dr%n_mode*qp%n_full_point, walltime() - t0) +end subroutine diff --git a/src/lineshape/phonondamping_generation.f90 b/src/lineshape/phonondamping_generation.f90 index 870fd9f4..7c3ca7aa 100644 --- a/src/lineshape/phonondamping_generation.f90 +++ b/src/lineshape/phonondamping_generation.f90 @@ -59,6 +59,7 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op se%isotope_scattering = opts%isotopescattering ! include isotope scattering? se%thirdorder_scattering = opts%thirdorder ! include threephonon term? se%fourthorder_scattering = opts%fourthorder ! include fourphonon term? + se%fourthorder_real = opts%fourthorder_real ! include fourphonon term? se%diagonal = opts%diagonal ! calculate only the diagonal part of the self-energy se%qdir = qdir ! Make note on the probe direction ! Make sure we know exactly when to skip symmetry things @@ -70,6 +71,7 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op allocate (se%im_3ph(se%n_energy, se%n_mode)) allocate (se%im_iso(se%n_energy, se%n_mode)) allocate (se%re_3ph(se%n_energy, se%n_mode)) + allocate (se%im_4ph(se%n_energy, se%n_mode)) allocate (se%re_4ph(se%n_energy, se%n_mode)) ! Get the energy range if (se%integrationtype .eq. 4) then @@ -82,16 +84,18 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op se%im_3ph = 0.0_r8 se%im_iso = 0.0_r8 se%re_3ph = 0.0_r8 + se%im_4ph = 0.0_r8 se%re_4ph = 0.0_r8 ! Now things should be clean and nice. Maybe say what we are about to do? if (verbosity .gt. 1) then - write (*, *) ' liso:', se%isotope_scattering - write (*, *) ' lthirdorder:', se%thirdorder_scattering - write (*, *) ' lfourthorder:', se%fourthorder_scattering - write (*, *) ' sigma:', se%smearing_prefactor*lo_mean(dr%default_smearing)*lo_frequency_hartree_to_thz, ' Thz' - write (*, *) ' temperature: ', tochar(opts%temperature), ' K' - write (*, *) ' frequencies:' + write (*, *) ' liso:', se%isotope_scattering + write (*, *) ' lthirdorder:', se%thirdorder_scattering + write (*, *) ' lfourthorder_real:', se%fourthorder_real + write (*, *) ' lfourthorder:', se%fourthorder_scattering + write (*, *) ' sigma:', se%smearing_prefactor*lo_mean(dr%default_smearing)*lo_frequency_hartree_to_thz, ' Thz' + write (*, *) ' temperature: ', tochar(opts%temperature), ' K' + write (*, *) ' frequencies:' do i = 1, dr%n_mode write (*, *) ' mode ', tochar(i, -3), ', omega: ', tochar(wp%omega(i)*lo_frequency_Hartree_to_THz, 6), ' THz' end do @@ -169,34 +173,35 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op call tmr%tock('three-phonon integrals') - ! Maybe fourthorder things if (se%fourthorder_scattering) then - fourthorder: block - real(r8), dimension(:), allocatable :: delta - integer :: j - - allocate (delta(dr%n_mode)) - call fourphonon_selfenergy(qpoint, wp, gp, qp, uc, opts%temperature, dr, fcf, delta, se%skipsym, mw, mem, verbosity) - do j = 1, se%n_mode - se%re_4ph(:, j) = delta(j) - end do - deallocate (delta) - end block fourthorder + call sr%generate_fourthorder(qpoint, wp, gp, qp, dr, uc, fc, fcf, & + se%skipsym, .false., opts%grid, tmr, mw, mem, verbosity) + select case (se%integrationtype) + case (1:2) + call fourphonon_imaginary_selfenergy_gaussian(wp, se, sr, qp, dr, opts%temperature, mw, mem, verbosity) + case (3:5) + call lo_stop_gracefully(['Integration type not implemented for fourth order scattering'], lo_exitcode_param, __FILE__, __LINE__, mw%comm) + case default + call lo_stop_gracefully(['Unknown integration type'], lo_exitcode_param, __FILE__, __LINE__, mw%comm) + end select end if - call tmr%tock('four-phonon self-energy') + call tmr%tock('four-phonon integrals') ! finalize to ensure that it's reasonable. sanity: block ! Make sure the selfenergy is zero where it's supposed to be. First ! make sure it's at least not negative. se%im_3ph = max(se%im_3ph, 0.0_r8) + se%im_4ph = max(se%im_4ph, 0.0_r8) se%im_iso = max(se%im_iso, 0.0_r8) ! zero at zero se%im_3ph(1, :) = 0.0_r8 + se%im_4ph(1, :) = 0.0_r8 se%im_iso(1, :) = 0.0_r8 ! zero at the end se%im_3ph(se%n_energy, :) = 0.0_r8 + se%im_4ph(se%n_energy, :) = 0.0_r8 se%im_iso(se%n_energy, :) = 0.0_r8 end block sanity @@ -255,7 +260,74 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op end do end if end block kktransform - call tmr%tock('Kramers-Kronig transformation') + call tmr%tock('Kramers-Kronig transformation 3ph') + end if + + + ! Kramers-Kronig-transform the imaginary part to get the real for the four phonon + if (se%fourthorder_scattering) then + kktransform_4ph: block + real(r8) :: pref = 2.0_r8/lo_pi + complex(r8), dimension(:), allocatable :: z0 + complex(r8) :: eta + real(r8), dimension(:), allocatable :: x, xs, y0 + real(r8) :: xp, dlx + integer :: imode, ie, ctr + + call mem%allocate(x, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(xs, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(z0, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(y0, se%n_energy, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + + x = se%energy_axis ! copy of x-axis + xs = x**2 ! x^2, precalculated + dlx = x(2) - x(1) ! prefactor for riemann integral + eta = lo_imag*(x(2) - x(1))*1E-8_r8 ! small imaginary thing to avoid divergence + + ctr = 0 + do imode = 1, se%n_mode + do ie = 1, se%n_energy + ctr = ctr + 1 + if (mod(ctr, mw%n) .ne. mw%r) cycle + xp = x(ie) + y0 = se%im_4ph(:, imode)*x + ! To anyone reading this code: + ! The correct way to compute the principal value would be + ! z0 = (xp + eta)**2 - xs + ! However, it doesn't matter for the Im -> Re transform because we + ! don't really hit 0. + z0 = xp**2 - xs + eta + y0 = real(y0/z0, r8) + y0(1) = y0(1)*0.5_r8 + y0(se%n_energy) = y0(se%n_energy)*0.5_r8 + se%re_4ph(ie, imode) = se%re_4ph(ie, imode) + sum(y0)*dlx + end do + end do + call mw%allreduce('sum', se%re_4ph) + se%re_4ph = se%re_4ph*pref + + call mem%deallocate(x, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(xs, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(z0, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(y0, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + end block kktransform_4ph + call tmr%tock('Kramers-Kronig transformation 4ph') + end if + + ! Do the fourth order real part after the Kramer-Kronig to not mess with the (maybe) already computed real part + if (se%fourthorder_real) then + fourthorder_real: block + real(r8), dimension(:), allocatable :: delta + integer :: j + + allocate (delta(dr%n_mode)) + call fourphonon_selfenergy(qpoint, wp, gp, qp, uc, opts%temperature, dr, fcf, delta, se%skipsym, mw, mem, verbosity) + do j = 1, se%n_mode + se%re_4ph(:, j) = se%re_4ph(:, j) + delta(j) + end do + deallocate (delta) + call tmr%tock('four-phonon real part') + end block fourthorder_real end if ! Correct the memory matrix to get back to the orthogonal time one @@ -342,7 +414,7 @@ subroutine generate(se, qpoint, qdir, wp, uc, fc, fct, fcf, ise, isf, qp, dr, op if (mod(imode, mw%n) .ne. mw%r) cycle ! Evaluate rough spectral function - yim = se%im_3ph(:, imode) + se%im_iso(:, imode) + yim = se%im_3ph(:, imode) + se%im_4ph(:, imode) + se%im_iso(:, imode) yre = se%re_3ph(:, imode) + se%re_4ph(:, imode) call evaluate_spectral_function(se%energy_axis, yim, yre, wp%omega(imode), ysf) ! Find peaks in the spectral function @@ -400,6 +472,7 @@ subroutine generate_interp(se, qpoint, qdir, wp, uc, ise, opts, mw, mem) se%isotope_scattering = opts%isotopescattering ! include isotope scattering? se%thirdorder_scattering = opts%thirdorder ! include threephonon term? se%fourthorder_scattering = opts%fourthorder ! include fourphonon term? + se%fourthorder_real = opts%fourthorder_real ! include fourphonon term? se%diagonal = opts%diagonal ! calculate only the diagonal part of the self-energy se%skipsym = opts%qpointpath ! If on a path, skip using symmetry. For now, should reintroduce it later. se%qdir = qdir ! Make note on the probe direction diff --git a/src/lineshape/phonondamping_grid.f90 b/src/lineshape/phonondamping_grid.f90 index 5d7b8bf9..16b99293 100644 --- a/src/lineshape/phonondamping_grid.f90 +++ b/src/lineshape/phonondamping_grid.f90 @@ -142,7 +142,7 @@ module subroutine get_selfenergy_on_closed_grid(sf, tc, pd, qp, dr, uc, fc, fct, if (mod(imode, mw%n) .ne. mw%r) cycle ! Get the spectral functions if (dr%iq(iq)%omega(imode) .gt. lo_freqtol) then - buf_sigmaIm(:, imode) = se%im_3ph(:, imode) + se%im_iso(:, imode) + buf_sigmaIm(:, imode) = se%im_3ph(:, imode) + se%im_4ph(:, imode) + se%im_iso(:, imode) buf_sigmaRe(:, imode) = se%re_3ph(:, imode) + se%re_4ph(:, imode) call taperfn_im(se%energy_axis, opts%maxf, dr%iq(iq)%omega(imode), buf_taper) buf_sigmaIm(:, imode) = buf_sigmaIm(:, imode)*buf_taper @@ -414,7 +414,7 @@ module subroutine get_selfenergy_on_points(qvec, wp, qp, dr, uc, fc, fct, fcf, i call se%generate(qpoint, qdir, wp(iq), uc, fc, fct, fcf, ise, isf, qp, dr, opts, tmr, mw, mem, verbosity=-1) ! Store it sigmaRe(:, :, iq) = se%re_3ph + se%re_4ph - sigmaIm(:, :, iq) = se%im_3ph + se%im_iso + sigmaIm(:, :, iq) = se%im_3ph + se%im_4ph + se%im_iso ! Report if (verbosity .gt. 0) then write (*, *) 'did', iq, 'out of', size(qvec, 2) diff --git a/src/lineshape/phonondamping_path.f90 b/src/lineshape/phonondamping_path.f90 index 323ece5a..4a46788f 100644 --- a/src/lineshape/phonondamping_path.f90 +++ b/src/lineshape/phonondamping_path.f90 @@ -159,7 +159,7 @@ module subroutine spectral_function_along_path(bs, uc, fc, fct, fcf, ise, qp, dr if (mw%r .eq. solrnk) then ! Store self-energy in buffers buf_re(:, :, jpt) = se%re_3ph + se%re_4ph - buf_im(:, :, jpt) = se%im_3ph + se%im_iso + buf_im(:, :, jpt) = se%im_3ph + se%im_4ph + se%im_iso ! Store the shifts from third and fourth order while they are available do imode = 1, dr%n_mode buf_shift_3rd(imode, jpt) = lo_linear_interpolation(se%energy_axis, se%re_3ph(:, imode), bs%p(iqp)%omega(imode)) diff --git a/src/lineshape/scatteringrates.f90 b/src/lineshape/scatteringrates.f90 index fb4b0a25..3f4eb0d5 100644 --- a/src/lineshape/scatteringrates.f90 +++ b/src/lineshape/scatteringrates.f90 @@ -17,25 +17,43 @@ module scatteringrates public :: lo_listofscatteringrates type lo_listofscatteringrates + !> For third order !> frequency at q'' real(r8), dimension(:, :), allocatable :: omega3 !> group velocity at q'' real(r8), dimension(:, :, :), allocatable :: vel3 !> three-phonon matrix elements, not squared (mode1,mode2,mode3,q) complex(r8), dimension(:, :, :, :), allocatable :: psi_3ph + !> q-vectors are needed for evaluation of spectral functions + real(r8), dimension(:, :), allocatable :: qvec3 + !> eigenvectors are needed for spectral functions + complex(r8), dimension(:, :, :), allocatable :: egv3 + + !> For isotope !> isotope matrix elements (mode1,mode2,q) real(r8), dimension(:, :, :), allocatable :: psi_iso + !> make a note if q1 is gamma logical :: atgamma = .false. + !> For fourth order + !> frequency at q4 + real(r8), dimension(:, :, :), allocatable :: omega4 + !> group velocity at q4 + real(r8), dimension(:, :, :, :), allocatable :: vel4 + !> four-phonon matrix elements, not squared (mode1, mode2, mode3, mode4, q3, q4) + complex(r8), dimension(:, :, :, :, :, :), allocatable :: psi_4ph !> q-vectors are needed for evaluation of spectral functions - real(r8), dimension(:, :), allocatable :: qvec3 + real(r8), dimension(:, :, :, :), allocatable :: qvec4 !> eigenvectors are needed for spectral functions - complex(r8), dimension(:, :, :), allocatable :: egv3 + complex(r8), dimension(:, :, :, :, :), allocatable :: egv4 + contains !> calculate matrix elements procedure :: generate + !> calculate matrix elements only for four phonons + procedure :: generate_fourthorder end type !> build a flattened helper to make things fast @@ -43,10 +61,12 @@ module scatteringrates !> eigenvectors? complex(r8), dimension(:, :), allocatable :: ugv1 complex(r8), dimension(:, :, :), allocatable :: ugv2, ugv3 + complex(r8), dimension(:, :, :, :), allocatable :: ugv4 !> pre-transformed matrix element complex(r8), dimension(:), allocatable :: ptf_phi complex(r8), dimension(:), allocatable :: evp1 complex(r8), dimension(:), allocatable :: evp2 + complex(r8), dimension(:), allocatable :: evp3 end type contains @@ -863,4 +883,373 @@ subroutine pretransform_phi(fct, q2, q3, ptf) ! end function ! end function + +subroutine generate_fourthorder(sr, qpoint, ompoint, gpoint, qp, dr, uc, fc, fcf, skipsym, eigenvectors, & + closedgrid, tmr, mw, mem, verbosity) + !> scattering rates + class(lo_listofscatteringrates), intent(out) :: sr + !> q-point in question + type(lo_qpoint) :: qpoint + !> harmonic properties at the relevant q-point + type(lo_phonon_dispersions_qpoint), intent(in) :: ompoint + !> harmonic properties at Gamma with the right direction + type(lo_phonon_dispersions_qpoint), intent(in) :: gpoint + !> qpoint mesh + class(lo_qpoint_mesh), intent(in) :: qp + !> phonon dispersions + type(lo_phonon_dispersions), intent(in) :: dr + !> crystal structure + type(lo_crystalstructure), intent(in) :: uc + !> second order forceconstant + type(lo_forceconstant_secondorder), intent(inout) :: fc + !> third order force constant + type(lo_forceconstant_fourthorder), intent(in) :: fcf + !> calculate everything? + logical, intent(in) :: skipsym, eigenvectors, closedgrid + !> timer + type(lo_timer), intent(inout) :: tmr + !> MPI communicator + type(lo_mpi_helper), intent(inout) :: mw + !> Memory tracker + type(lo_mem_helper), intent(inout) :: mem + !> how much to talk + integer, intent(in) :: verbosity + + ! Just keep a list of q4, temporarily + type(flathelper) :: fh + real(r8), dimension(:, :, :), allocatable :: qvec4 + complex(r8), dimension(:, :, :, :), allocatable :: egv4 + real(r8) :: timer, t0, t1 + integer :: ind_gamma_full + + ! Start timer + timer = walltime() + t0 = timer + t1 = timer + + call tmr%tick() + + init: block + integer :: i + + ! Locate the index for Gamma on the q-grid + do i=1, qp%n_full_point + if (lo_sqnorm(qp%ap(i)%r) .lt. lo_sqtol) then + ind_gamma_full = i + exit + end if + end do + end block init + + fourthfreq: block + type(lo_phonon_dispersions_qpoint) :: op4 + real(r8), dimension(3) :: qv1, qv2, qv3, qv4, v0 + integer :: q2, q3, idxq, ctr + + ! Make some space + allocate (sr%omega4(dr%n_mode, qp%n_full_point, qp%n_full_point)) + allocate (sr%vel4(3, dr%n_mode, qp%n_full_point, qp%n_full_point)) + call mem%allocate(egv4, [dr%n_mode, dr%n_mode, qp%n_full_point, qp%n_full_point], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%allocate(qvec4, [3, qp%n_full_point, qp%n_full_point], persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + sr%omega4 = 0.0_r8 + sr%vel4 = 0.0_r8 + egv4 = 0.0_r8 + qvec4 = 0.0_r8 + + ctr = 0 + do q2 = 1, qp%n_full_point + ! make it parallel + ctr = ctr + 1 + if (mod(ctr, mw%n) .ne. mw%r) cycle + + do q3 = 1, qp%n_full_point + ! Get the q-vectors + qv1 = qpoint%r + qv2 = qp%ap(q2)%r + qv3 = qp%ap(q3)%r + qv4 = -qv1 - qv2 - qv3 + if (closedgrid) then + ! A closed grid, fetch data from the grid to the fourth point + qv4 = -qv1 - qv2 - qv3 + qv4 = matmul(uc%inv_reciprocal_latticevectors, qv4) + qv4 = lo_clean_fractional_coordinates(qv4) + idxq = index_on_grid(qp, qv4) + sr%omega4(:, q2, q3) = dr%aq(idxq)%omega + sr%vel4(:, :, q2, q3) = dr%aq(idxq)%vel + egv4(:, :, q2, q3) = dr%aq(idxq)%egv + qvec4(:, q2, q3) = qv4 + else + ! Not a closed grid, just fetch everything + ! Get harmonic things + call op4%generate(fc, uc, mem, qvec=qv4) + sr%omega4(:, q2, q3) = op4%omega + sr%vel4(:, :, q2, q3) = op4%vel + egv4(:, :, q2, q3) = op4%egv + qvec4(:, q2, q3) = qv4 + + v0 = matmul(uc%inv_reciprocal_latticevectors, qv4) + if (sum(abs(v0 - anint(v0))) .lt. lo_sqtol) then + ! we are at Gamma, Replace with proper Gamma + qvec4(:, q2, q3) = 0.0_r8 + egv4(:, :, q2, q3) = gpoint%egv + sr%vel4(:, :, q2, q3) = gpoint%vel + sr%omega4(:, q2, q3) = gpoint%omega + end if + end if + end do + end do + call mw%allreduce('sum', sr%omega4) + call mw%allreduce('sum', sr%vel4) + call mw%allreduce('sum', egv4) + call mw%allreduce('sum', qvec4) + if (verbosity .gt. 0) then + t1 = walltime() + write (*, *) '' + write (*, *) "... got q'' and q''' frequencies (", tochar(t1 - t0), "s)" + t0 = t1 + end if + end block fourthfreq + + call tmr%tock("q''' harmonic properties") + + pretransform: block + real(r8) :: f0, f1 + integer :: iq, jq, imode, jmode, iatom, ix, ialpha, ctr + + ! Some temporary space for scaled eigenvectors + allocate (fh%ugv1(dr%n_mode, dr%n_mode)) + allocate (fh%ugv2(dr%n_mode, dr%n_mode, qp%n_full_point)) + allocate (fh%ugv4(dr%n_mode, dr%n_mode, qp%n_full_point, qp%n_full_point)) + fh%ugv1 = 0.0_r8 + fh%ugv2 = 0.0_r8 + fh%ugv4 = 0.0_r8 + + ! First fix the one at the main q-point, always the same + ! just pre-multiply with masses and frequencies + do imode = 1, dr%n_mode + do iatom = 1, uc%na + do ix = 1, 3 + if (norm2(qpoint%r) .lt. lo_sqtol) then + ! Get consistent Gamma + if (gpoint%omega(imode) .gt. lo_freqtol) then + f0 = 1.0_r8 / sqrt(gpoint%omega(imode)) + else + f0 = 0.0_r8 + end if + ialpha = (iatom -1) * 3 + ix + f1 = uc%invsqrtmass(iatom) + fh%ugv1(ialpha, imode) = gpoint%egv(ialpha, imode) * f0 * f1 + else + ! Just normal copy + if (ompoint%omega(imode) .gt. lo_freqtol) then + f0 = 1.0_r8/sqrt(ompoint%omega(imode)) + else + f0 = 0.0_r8 + end if + ialpha = (iatom - 1) * 3 + ix + f1 = uc%invsqrtmass(iatom) + fh%ugv1(ialpha, imode) = ompoint%egv(ialpha, imode) * f0 * f1 + end if + end do + end do + end do + + ! Pre-multiply eigenvectors with masses and frequencies + ctr = 0 + do iq = 1, qp%n_full_point + ctr = ctr + 1 + if (mod(ctr, mw%n) .ne. mw%r) cycle + do imode = 1, dr%n_mode + ! First get the thingy for q' + if (iq .eq. ind_gamma_full) then + ! Pick consistent Gamma + if (gpoint%omega(imode) .gt. lo_freqtol) then + f0 = 1.0_r8 / sqrt(gpoint%omega(imode)) + else + f0 = 0.0_r8 + end if + do iatom=1, uc%na + f1 = uc%invsqrtmass(iatom) + do ix = 1, 3 + ialpha = (iatom - 1)*3 + ix + fh%ugv2(ialpha, imode, iq) = gpoint%egv(ialpha, imode)*f0*f1 + end do + end do + else + ! Not at Gamma, don't have to care + if (dr%aq(iq)%omega(imode) .gt. lo_freqtol) then + f0 = 1.0_r8/sqrt(dr%aq(iq)%omega(imode)) + else + f0 = 0.0_r8 + end if + do iatom=1, uc%na + f1 = uc%invsqrtmass(iatom) + do ix=1, 3 + ialpha = (iatom - 1)*3+ix + fh%ugv2(ialpha, imode, iq) = dr%aq(iq)%egv(ialpha, imode)*f0*f1 + end do + end do + end if + ! Now we can get the things at q'' and q''' + do jq=1, qp%n_full_point + ! And finally for q''' + if (sr%omega4(imode, iq, jq) .gt. lo_freqtol) then + f0 = 1.0_r8 / sqrt(sr%omega4(imode, iq, jq)) + else + f0 = 0.0_r8 + end if + do iatom = 1, uc%na + f1 = uc%invsqrtmass(iatom) + do ix = 1, 3 + ialpha = (iatom - 1) * 3 + ix + fh%ugv4(ialpha, imode, iq, jq) = egv4(ialpha, imode, iq, jq) * f0 * f1 + end do + end do + end do + end do + end do + call mw%allreduce('sum', fh%ugv2) + call mw%allreduce('sum', fh%ugv4) + ! Space for the pre-transformed phi + allocate (fh%ptf_phi(dr%n_mode**4)) + fh%ptf_phi = 0.0_r8 + + ! Space for outer-producted eigenvectors ? + allocate (fh%evp1(dr%n_mode**2)) + allocate (fh%evp2(dr%n_mode**3)) + allocate (fh%evp3(dr%n_mode**4)) + fh%evp1 = 0.0_r8 + fh%evp2 = 0.0_r8 + fh%evp3 = 0.0_r8 + + if (verbosity .gt. 0) then + t1 = walltime() + write (*, *) "... pretransformed and made temporary space (", tochar(t1 - t0), "s)" + t0 = t1 + end if + end block pretransform + + call tmr%tock("4ph eigenvector scaling") + + ! Get the fourphonon matrix elements + fourphsc2: block + real(r8), dimension(3) :: qv2, qv3, qv4 + integer :: q2, q3, b1, b2, b3, b4, ctr, l + + ! Some space + t0 = walltime() + if (verbosity .gt. 0) call lo_progressbar_init() + + ! Space for matrix elements + allocate (sr%psi_4ph(dr%n_mode, dr%n_mode, dr%n_mode, dr%n_mode, qp%n_full_point, qp%n_full_point)) + sr%psi_4ph = 0.0_r8 + + ctr = 0 + l = 0 + do q2 = 1, qp%n_full_point + do q3 = 1, qp%n_full_point + ! fetch q-vectors + qv2 = qp%ap(q2)%r + qv3 = qp%ap(q3)%r + qv4 = qvec4(:, q2, q3) + ! do the half transform + call pretransform_phi4(fcf, qv2, qv3, qv4, fh%ptf_phi) + ! Get matrix elements for all modes + do b1=1, dr%n_mode + do b2=1, dr%n_mode + fh%evp1 = 0.0_r8 + call zgeru(dr%n_mode, dr%n_mode, (1.0_r8, 0.0_r8), fh%ugv2(:, b2, q2), 1, fh%ugv1(:, b1), 1, fh%evp1, dr%n_mode) + do b3=1, dr%n_mode + ! MPI division + l = l + 1 + if (mod(l, mw%n) .ne. mw%r) cycle + + fh%evp2 = 0.0_r8 + call zgeru(dr%n_mode, dr%n_mode*dr%n_mode, (1.0_r8, 0.0_r8), fh%ugv2(:, b3, q3), 1, fh%evp1, 1, fh%evp2, dr%n_mode) + do b4=1, dr%n_mode + fh%evp3 = 0.0_r8 + call zgeru(dr%n_mode, dr%n_mode*dr%n_mode*dr%n_mode, (1.0_r8, 0.0_r8), fh%ugv4(:, b4, q2, q3), 1, fh%evp2, 1, fh%evp3, dr%n_mode) + fh%evp3 = conjg(fh%evp3) + sr%psi_4ph(b1, b2, b3, b4, q2, q3) = dot_product(fh%evp3, fh%ptf_phi) + end do + end do + end do + end do + + ! Progressbar for four phonon + if (verbosity .gt. 0) then + ctr = ctr + 1 + if (lo_trueNtimes(ctr, 127, qp%n_full_point*qp%n_full_point)) then + call lo_progressbar('... fourphonon matrixelements', ctr, qp%n_full_point*qp%n_full_point, walltime() - t0) + end if + end if + end do + end do + ! sync across ranks + call mw%allreduce('sum', sr%psi_4ph) + + if (verbosity .gt. 0) then + t1 = walltime() + call lo_progressbar('... fourphonon matrixelements', qp%n_full_point*dr%n_mode, qp%n_full_point*dr%n_mode, t1 - t0) + t0 = t1 + end if + end block fourphsc2 + call tmr%tock("four-phonon matrix elements") + + ! Some cleanup + call mem%deallocate(egv4, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%deallocate(qvec4, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + call mem%tock(__FILE__, __LINE__, mw%comm) +end subroutine + +!> pre-transform to get half of the matrix elements for psi4 +subroutine pretransform_phi4(fcf, q2, q3, q4, ptf) + !> fourth order forceconstant + type(lo_forceconstant_fourthorder), intent(in) :: fcf + !> q-vectors + real(r8), dimension(3), intent(in) :: q2, q3, q4 + !> flattened, pretransformed matrix element + complex(r8), dimension(:), intent(out) :: ptf + + integer i, j, k, l, m + + complex(r8) :: expiqr + real(r8), dimension(3) :: rv2, rv3, rv4 + real(r8) :: iqr + integer :: a1, a2, a3, a4, ia, ib, ic, id, t, nb + + nb = fcf%na*3 + ptf = 0.0_r8 + do a1=1, fcf%na + do t = 1,fcf%atom(a1)%n + a2 = fcf%atom(a1)%quartet(t)%i2 + a3 = fcf%atom(a1)%quartet(t)%i3 + a4 = fcf%atom(a1)%quartet(t)%i4 + + rv2 = fcf%atom(a1)%quartet(t)%lv2 + rv3 = fcf%atom(a1)%quartet(t)%lv3 + rv4 = fcf%atom(a1)%quartet(t)%lv4 + + iqr = dot_product(q2, rv2) + dot_product(q3, rv3) + dot_product(q4, rv4) + iqr = -iqr*lo_twopi + expiqr = cmplx(cos(iqr), sin(iqr), r8) + do i=1, 3 + do j=1, 3 + do k=1, 3 + do l=1, 3 + ia = (a1 - 1)*3 + i + ib = (a2 - 1)*3 + j + ic = (a3 - 1)*3 + k + id = (a4 - 1)*3 + l + ! Now for the grand flattening scheme, to be consistent with the zgeru operations + m = (ia - 1)*nb*nb*nb + (ib - 1)*nb*nb + (ic - 1)*nb + id + ptf(m) = ptf(m) + fcf%atom(a1)%quartet(t)%m(i, j, k, l)*expiqr + end do + end do + end do + end do + end do + end do +end subroutine + end module