module m_vmts

contains

  subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, yukawa_residual, vr )

  !-------------------------------------------------------------------------
  ! This subroutine calculates the lattice harmonics expansion coefficients 
  ! of the Coulomb / Yukawa potential for all atom types.                     
  !
  ! In more detail:
  ! the radius-dependent coefficient of the potential is
  ! V_{lm}(r) = \int_0^R G_l(r,r') \rho_{lm}(r') r'^2 dr' + V_{lm}(R) P_l(r)
  ! where 
  ! [sphere boundary contribution - first part of code:]
  ! V_{lm}(R) is the spherical harmonics expansion of the interstitial 
  ! potential,
  ! P_l(r) is the derivative of the spherical Green's function on the sphere 
  ! boundary (times a constant),
  ! [sphere interior contribution - second part of code:]
  ! G_l(r,r') ) = green_factor u_1(r_<) u_2(r_>) is the spherical Green's 
  ! function with homogeneous solutions u_1 and u_2 and 
  ! r_< = min(r,r') and r_> = max(r,r')
  ! The integral is split in a part where r'=r_< and a part where r'=r_> and
  ! the integral from r to R is split in \int_0^R - \int_0^r. Resulting 
  ! terms depending solely on R (not on r) are summarised in the variable 
  ! termsR. 
  !-------------------------------------------------------------------------

#include"cpp_double.h"
    use m_constants
    use m_types
    use m_intgr, only : intgr2
    use m_phasy1
    use m_sphbes
    use m_od_phasy
    use m_SphBessel
    implicit none

    type(t_input),  intent(in)        :: input
    type(t_mpi),    intent(in)        :: mpi
    type(t_stars),  intent(in)        :: stars
    type(t_sphhar), intent(in)        :: sphhar
    type(t_atoms),  intent(in)        :: atoms
    type(t_sym),    intent(in)        :: sym
    type(t_cell),   intent(in)        :: cell
    type(t_oneD),   intent(in)        :: oneD
    complex,        intent(in)        :: vpw(:)!(stars%ng3,input%jspins)
    real,           intent(in)        :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
    logical,        intent(in)        :: yukawa_residual
    real,           intent(out)       :: vr(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)

    complex                           :: cp, sm
    integer                           :: i, jm, k, l, lh, n, nd, lm, n1, nat, m, imax, lmax
    complex                           :: vtl(0:sphhar%nlhd, atoms%ntype)
    complex                           :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
    real                              :: green_factor, termsR
    real                              :: green_1    (1:atoms%jmtd), green_2    (1:atoms%jmtd)
    real                              :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
    real                              :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)
    real                              :: sbf(0:atoms%lmaxd)
    real, allocatable, dimension(:,:) :: il, kl
    
    !$ COMPLEX vtl_loc(0:sphhar%nlhd,atoms%ntype)
#ifdef CPP_MPI
    include 'mpif.h'
    integer                       :: ierr(3)
    complex, allocatable          :: c_b(:)

    external MPI_REDUCE
#endif
    integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM


    !-------------------------------------------------------------------------
    ! sphere boundary contribution to the coefficients calculated from the values
    ! of the interstitial Coulomb / Yukawa potential on the sphere boundary

    vtl(:,:) = cmplx( 0.0, 0.0 )

#ifdef CPP_MPI
    call MPI_BCAST( vpw, size(vpw), CPP_MPI_COMPLEX, 0, mpi, ierr )
#endif

    ! q=0 component
    if ( mpi%irank == 0 ) then
      vtl(0,1:atoms%ntype) = sfp_const * vpw(1)
    end if

    ! q/=0 components
    !$omp parallel default( none ) &
    !$omp& shared( mpi, stars, vpw, oneD, atoms, sym, cell, sphhar, vtl ) &
    !$omp& private( k, cp, pylm, nat, n, sbf, nd, lh, sm, jm, m, lm, l ) &
    !$omp& private( vtl_loc ) 
    !$ vtl_loc(:,:) = cmplx(0.d0,0.d0)
    !$omp do
    do k = mpi%irank+2, stars%ng3, mpi%isize
      cp = vpw(k) * stars%nstr(k)
      if ( .not. oneD%odi%d1 ) then
        call phasy1( atoms, stars, sym, cell, k, pylm )
      else
        call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, &
            atoms%neq, atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm )
      end if
       
      nat = 1
      do n = 1, atoms%ntype
        call sphbes( atoms%lmax(n), stars%sk3(k) * atoms%rmt(n), sbf )
        nd = atoms%ntypsy(nat)
        do lh = 0, sphhar%nlh(nd)
          l = sphhar%llh(lh,nd)
          sm = (0.0,0.0)
          do jm = 1, sphhar%nmem(lh,nd)
            m = sphhar%mlh(jm,lh,nd)
            lm = l * ( l + 1 ) + m + 1 
            sm = sm + conjg( sphhar%clnu(jm,lh,nd) ) * pylm(lm,n)
          end do
        !$ if (.false.) then
        vtl(lh,n) = vtl(lh,n) + cp * sbf(l) * sm
        !$ end if
        !$ vtl_loc(lh,n) = vtl_loc(lh,n) + cp * sbf(l) * sm
        end do
        nat = nat + atoms%neq(n)
      end do
    end do
    !$omp end do
    !$omp critical
    !$ vtl = vtl + vtl_loc
    !$omp end critical
    !$omp end parallel
#ifdef CPP_MPI
    n1 = ( sphhar%nlhd + 1 ) * atoms%ntype
    allocate( c_b(n1) )
    call MPI_REDUCE( vtl, c_b, n1, CPP_MPI_COMPLEX, MPI_SUM, 0, mpi%mpi_comm, ierr )
    if ( mpi%irank == 0 ) vtl = reshape( c_b, (/sphhar%nlhd+1,atoms%ntype/) )
    deallocate (c_b)
#endif

    !-------------------------------------------------------------------------
    ! sphere interior contribution to the coefficients calculated from the 
    ! values of the sphere Coulomb/Yukawa potential on the sphere boundary
    
    if ( yukawa_residual ) then
      allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
    end if

    nat = 1
    do n = 1, atoms%ntype
      nd = atoms%ntypsy(nat)
      imax = atoms%jri(n)
      lmax = sphhar%llh(sphhar%nlh(nd), nd)
      if ( yukawa_residual ) then
        do concurrent (i = 1:imax)
          call ModSphBessel( il(0:,i), kl(0:,i), input%preconditioning_param * atoms%rmsh(i,n), lmax )
        end do
      end if
      do lh = 0, sphhar%nlh(nd)
        l = sphhar%llh(lh,nd)
        if ( yukawa_residual ) then
          green_1(1:imax) = il(l,1:imax)
          green_2(1:imax) = kl(l,1:imax)
          green_factor    = fpi_const * input%preconditioning_param
        else
          green_1(1:imax) = atoms%rmsh(1:imax,n) ** l
          green_2(1:imax) = 1.0 / ( green_1(1:imax) * atoms%rmsh(1:imax,n) )
          green_factor    = fpi_const / ( 2 * l + 1 )
        end if
        integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
        integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
        call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
        call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
        termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
        vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
                                           + green_2(1:imax) *            integral_1(1:imax)   )
      end do
      nat = nat + atoms%neq(n)
    end do
    if ( yukawa_residual ) then
      deallocate( il, kl )
    end if

    if ( .not. yukawa_residual ) then
      do n = 1, atoms%ntype
        vr(1:atoms%jri(n),0,n) = vr(1:atoms%jri(n),0,n) - sfp_const * ( 1.0 / atoms%rmsh(1:atoms%jri(n),n) - 1.0 / atoms%rmt(n) ) * atoms%zatom(n)
      end do
    end if

  end subroutine vmts

end module m_vmts