LCOV - code coverage report
Current view: top level - vgen - vmts.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 74 75 98.7 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : module m_vmts
       2             : 
       3             : contains
       4             : 
       5         346 :   subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, potdenType, vr )
       6             : 
       7             :   !-------------------------------------------------------------------------
       8             :   ! This subroutine calculates the lattice harmonics expansion coefficients 
       9             :   ! of the Coulomb / Yukawa potential for all atom types.                     
      10             :   !
      11             :   ! In more detail:
      12             :   ! the radius-dependent coefficient of the potential is
      13             :   ! V_{lm}(r) = \int_0^R G_l(r,r') \rho_{lm}(r') r'^2 dr' + V_{lm}(R) P_l(r)
      14             :   ! where 
      15             :   ! [sphere boundary contribution - first part of code:]
      16             :   ! V_{lm}(R) is the spherical harmonics expansion of the interstitial 
      17             :   ! potential,
      18             :   ! P_l(r) is the derivative of the spherical Green's function on the sphere 
      19             :   ! boundary (times a constant),
      20             :   ! [sphere interior contribution - second part of code:]
      21             :   ! G_l(r,r') ) = green_factor u_1(r_<) u_2(r_>) is the spherical Green's 
      22             :   ! function with homogeneous solutions u_1 and u_2 and 
      23             :   ! r_< = min(r,r') and r_> = max(r,r')
      24             :   ! The integral is split in a part where r'=r_< and a part where r'=r_> and
      25             :   ! the integral from r to R is split in \int_0^R - \int_0^r. Resulting 
      26             :   ! terms depending solely on R (not on r) are summarised in the variable 
      27             :   ! termsR.
      28             :   !
      29             :   ! More information and equations can be found in
      30             :   ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118(2011) 
      31             :   !-------------------------------------------------------------------------
      32             : 
      33             : #include"cpp_double.h"
      34             :     use m_constants
      35             :     use m_types
      36             :     use m_intgr, only : intgr2
      37             :     use m_phasy1
      38             :     use m_sphbes
      39             :     use m_od_phasy
      40             :     use m_SphBessel
      41             :     implicit none
      42             : 
      43             :     type(t_input),  intent(in)        :: input
      44             :     type(t_mpi),    intent(in)        :: mpi
      45             :     type(t_stars),  intent(in)        :: stars
      46             :     type(t_sphhar), intent(in)        :: sphhar
      47             :     type(t_atoms),  intent(in)        :: atoms
      48             :     type(t_sym),    intent(in)        :: sym
      49             :     type(t_cell),   intent(in)        :: cell
      50             :     type(t_oneD),   intent(in)        :: oneD
      51             :     complex,        intent(in)        :: vpw(:)!(stars%ng3,input%jspins)
      52             :     real,           intent(in)        :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      53             :     integer,        intent(in)        :: potdenType
      54             :     real,           intent(out)       :: vr(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      55             : 
      56             :     complex                           :: cp, sm
      57             :     integer                           :: i, jm, k, l, lh, n, nd, lm, n1, nat, m, imax, lmax
      58         692 :     complex                           :: vtl(0:sphhar%nlhd, atoms%ntype)
      59         346 :     complex                           :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
      60             :     real                              :: green_factor, termsR
      61         692 :     real                              :: green_1    (1:atoms%jmtd), green_2    (1:atoms%jmtd)
      62         692 :     real                              :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
      63         692 :     real                              :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)
      64         346 :     real                              :: sbf(0:atoms%lmaxd)
      65         692 :     real, allocatable, dimension(:,:) :: il, kl
      66             :     
      67         346 :     !$ complex, allocatable :: vtl_loc(:,:)
      68             : #ifdef CPP_MPI
      69             :     include 'mpif.h'
      70             :     integer                       :: ierr(3)
      71         346 :     complex, allocatable          :: c_b(:)
      72             : 
      73             :     external MPI_REDUCE
      74             : #endif
      75             :     integer :: OMP_GET_NUM_THREADS, OMP_GET_THREAD_NUM
      76             : 
      77             : 
      78             : 
      79             :     ! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values
      80             :     ! of the interstitial Coulomb / Yukawa potential on the sphere boundary
      81             : 
      82        1034 :     vtl(:,:) = cmplx( 0.0, 0.0 )
      83             : 
      84             :     
      85             :     ! q=0 component
      86         346 :     if ( mpi%irank == 0 ) then
      87         173 :       vtl(0,1:atoms%ntype) = sfp_const * vpw(1)
      88             :     end if
      89             : 
      90             :     ! q/=0 components
      91             :     !$omp parallel default( none ) &
      92             :     !$omp& shared( mpi, stars, vpw, oneD, atoms, sym, cell, sphhar, vtl ) &
      93             :     !$omp& private( k, cp, pylm, nat, n, sbf, nd, lh, sm, jm, m, lm, l ) &
      94         692 :     !$omp& private( vtl_loc )
      95         346 :     !$ allocate(vtl_loc(0:sphhar%nlhd,atoms%ntype)) 
      96        1722 :     !$ vtl_loc(:,:) = cmplx(0.0,0.0)
      97         346 :     !$omp do
      98             :     do k = mpi%irank+2, stars%ng3, mpi%isize
      99      174129 :       cp = vpw(k) * stars%nstr(k)
     100      174129 :       if ( .not. oneD%odi%d1 ) then
     101      174129 :         call phasy1( atoms, stars, sym, cell, k, pylm )
     102             :       else
     103             :         call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, &
     104           0 :             atoms%neq, atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm )
     105             :       end if
     106      174129 :       nat = 1
     107      463910 :       do n = 1, atoms%ntype
     108      289781 :         call sphbes( atoms%lmax(n), stars%sk3(k) * atoms%rmt(n), sbf )
     109      289781 :         nd = atoms%ntypsy(nat)
     110    10937855 :         do lh = 0, sphhar%nlh(nd)
     111    10648074 :           l = sphhar%llh(lh,nd)
     112    10648074 :           sm = (0.0,0.0)
     113    32795186 :           do jm = 1, sphhar%nmem(lh,nd)
     114    22147112 :             m = sphhar%mlh(jm,lh,nd)
     115    22147112 :             lm = l * ( l + 1 ) + m + 1 
     116    32795186 :             sm = sm + conjg( sphhar%clnu(jm,lh,nd) ) * pylm(lm,n)
     117             :           end do
     118             :           !$ if (.false.) then
     119             :           vtl(lh,n) = vtl(lh,n) + cp * sbf(l) * sm
     120             :           !$ end if
     121    10937855 :           !$ vtl_loc(lh,n) = vtl_loc(lh,n) + cp * sbf(l) * sm
     122             :         end do
     123      463910 :         nat = nat + atoms%neq(n)
     124             :       end do
     125             :     end do
     126             :     !$omp end do
     127         692 :     !$omp critical
     128        1722 :     !$ vtl = vtl + vtl_loc
     129             :     !$omp end critical
     130         346 :     !$ deallocate(vtl_loc) 
     131             :     !$omp end parallel
     132             : #ifdef CPP_MPI
     133         346 :     n1 = ( sphhar%nlhd + 1 ) * atoms%ntype
     134         346 :     allocate( c_b(n1) )
     135         346 :     call MPI_REDUCE( vtl, c_b, n1, CPP_MPI_COMPLEX, MPI_SUM, 0, mpi%mpi_comm, ierr )
     136         346 :     if ( mpi%irank == 0 ) vtl = reshape( c_b, (/sphhar%nlhd+1,atoms%ntype/) )
     137         346 :     deallocate( c_b )
     138             : #endif
     139             : 
     140             : 
     141             : 
     142             :     ! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the 
     143             :     ! values of the sphere Coulomb/Yukawa potential on the sphere boundary
     144             : 
     145         346 :     if( mpi%irank == 0 ) then
     146         173 :     if ( potdenType == POTDEN_TYPE_POTYUK ) then
     147           3 :       allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
     148             :     end if
     149             : 
     150         173 :     nat = 1
     151         517 :     do n = 1, atoms%ntype
     152         344 :       nd = atoms%ntypsy(nat)
     153         344 :       imax = atoms%jri(n)
     154         344 :       lmax = sphhar%llh(sphhar%nlh(nd), nd)
     155         344 :       if ( potdenType == POTDEN_TYPE_POTYUK ) then
     156             :         !do concurrent (i = 1:imax)
     157       22725 :         do i = 1,imax
     158       11370 :           call ModSphBessel( il(0:,i), kl(0:,i), input%preconditioning_param * atoms%rmsh(i,n), lmax )
     159             :         end do
     160             :       end if
     161        5867 :       do lh = 0, sphhar%nlh(nd)
     162        5523 :         l = sphhar%llh(lh,nd)
     163        5523 :         if ( potdenType == POTDEN_TYPE_POTYUK ) then
     164         207 :           green_1(1:imax) = il(l,1:imax)
     165         207 :           green_2(1:imax) = kl(l,1:imax)
     166         207 :           green_factor    = fpi_const * input%preconditioning_param
     167             :         else
     168        5316 :           green_1(1:imax) = atoms%rmsh(1:imax,n) ** l
     169     3275192 :           green_2(1:imax) = 1.0 / ( green_1(1:imax) * atoms%rmsh(1:imax,n) )
     170        5316 :           green_factor    = fpi_const / ( 2 * l + 1 )
     171             :         end if
     172        5523 :         integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
     173     3432098 :         integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
     174        5523 :         call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
     175        5523 :         call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
     176        5523 :         termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
     177             :         vr(1:imax,lh,n) = green_factor * (   green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
     178        5867 :                                            + green_2(1:imax) *            integral_1(1:imax)   )
     179             :       end do
     180         517 :       nat = nat + atoms%neq(n)
     181             :     end do
     182         173 :     if ( potdenType == POTDEN_TYPE_POTYUK ) then
     183           3 :       deallocate( il, kl )
     184             :     end if
     185             : 
     186             :     
     187             : 
     188         173 :     if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     189         499 :       do n = 1, atoms%ntype
     190         499 :         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)
     191             :       end do
     192             :     end if
     193             :     end if
     194             : 
     195        1384 :   end subroutine vmts
     196             : 
     197             : end module m_vmts

Generated by: LCOV version 1.13