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

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : module m_vgen_coulomb
       8             : 
       9             :   use m_juDFT
      10             : 
      11             : contains
      12             : 
      13         346 :   subroutine vgen_coulomb( ispin, mpi, dimension, oneD, input, field, vacuum, sym, stars, &
      14             :              cell, sphhar, atoms, den, vCoul, results )
      15             :     !----------------------------------------------------------------------------
      16             :     ! FLAPW potential generator                           
      17             :     !----------------------------------------------------------------------------
      18             :     ! Generates the Coulomb or Yukawa potential and optionally the 
      19             :     ! density-potential integrals
      20             :     ! vCoul%potdenType = POTDEN_TYPE_POTYUK -> Yukawa case
      21             :     ! It takes a spin variable to indicate in which spin-channel the charge 
      22             :     ! resides.
      23             :     !----------------------------------------------------------------------------
      24             : 
      25             :     use m_constants
      26             :     use m_types
      27             :     use m_vmts
      28             :     use m_intnv
      29             :     use m_vvac
      30             :     use m_vvacis
      31             :     use m_vvacxy
      32             :     use m_vintcz
      33             :     use m_checkdopall
      34             :     use m_od_vvac
      35             :     use m_od_vvacis
      36             :     use m_convol
      37             :     use m_psqpw
      38             :     use m_cfft
      39             :     implicit none
      40             : 
      41             :     integer,            intent(in)               :: ispin
      42             :     type(t_mpi),        intent(in)               :: mpi
      43             :     type(t_dimension),  intent(in)               :: dimension
      44             :     type(t_oneD),       intent(in)               :: oneD
      45             :     type(t_input),      intent(in)               :: input
      46             :     type(t_field),      intent(inout)            :: field
      47             :     type(t_vacuum),     intent(in)               :: vacuum
      48             :     type(t_sym),        intent(in)               :: sym
      49             :     type(t_stars),      intent(in)               :: stars
      50             :     type(t_cell),       intent(in)               :: cell
      51             :     type(t_sphhar),     intent(in)               :: sphhar
      52             :     type(t_atoms),      intent(in)               :: atoms 
      53             :     type(t_potden),     intent(in)               :: den
      54             :     type(t_potden),     intent(inout)            :: vCoul
      55             :     type(t_results),    intent(inout), optional  :: results
      56             : 
      57             :     complex                                      :: vintcza, xint, rhobar
      58             :     integer                                      :: i, i3, irec2, irec3, ivac, j, js, k, k3
      59             :     integer                                      :: lh, n, nzst1
      60             :     integer                                      :: imz, imzxy, ichsmrg, ivfft
      61             :     integer                                      :: l, nat
      62             :     real                                         :: ani, g3, z, sig1dh, vz1dh
      63         692 :     complex, allocatable                         :: alphm(:,:), psq(:)
      64         692 :     real,    allocatable                         :: af1(:), bf1(:)
      65             : #ifdef CPP_MPI
      66             :     include 'mpif.h'
      67             :     integer:: ierr
      68             : #endif
      69             :  
      70             :     
      71         346 :     allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3)  )
      72         346 :     vCoul%iter = den%iter
      73             : 
      74             : 
      75             : 
      76             :     ! PSEUDO-CHARGE DENSITY COEFFICIENTS
      77         346 :     call timestart( "psqpw" )      
      78             :     call psqpw( mpi, atoms, sphhar, stars, vacuum,  cell, input, sym, oneD, &
      79         346 :          den%pw(:,ispin), den%mt(:,:,:,ispin), den%vacz(:,:,ispin), .false., vCoul%potdenType, psq )
      80         346 :     call timestop( "psqpw" )
      81             : 
      82             : 
      83             : 
      84             :     ! VACUUM POTENTIAL
      85         346 :     if ( mpi%irank == 0 ) then
      86         173 :       if ( oneD%odi%d1 ) then
      87           0 :         call timestart( "Vacuum" )
      88             :         !---> generates the m=0,gz=0 component of the vacuum potential
      89           0 :         call od_vvac( stars, vacuum, cell, psq, den%vacz(:,:,ispin), vCoul%vacz(:,:,ispin) )
      90             :         !---> generation of the vacuum warped potential components and
      91             :         !---> interstitial pw potential
      92             :         !---> vvacxy_5.F is a symmetrized potential generator
      93             :         call od_vvacis( oneD%odi%n2d, input, vacuum, oneD%odi%nq2, &
      94             :              oneD%odi%kv, cell, oneD%odi%M, stars, oneD%odi%nst2, &
      95             :              oneD, den%vacz(:,:,ispin), den%vacxy(:,:,:,ispin), psq, &
      96           0 :              vCoul%vacz(:,:,ispin), sym, vCoul%vacxy(:,:,:,ispin), vCoul%pw(:,ispin) )
      97           0 :         call timestop( "Vacuum" )
      98             :         !---> generation of the vacuum warped potential components and       
      99         173 :       elseif ( input%film .and. .not. oneD%odi%d1 ) then
     100             :         !     ----> potential in the  vacuum  region
     101           5 :         call timestart( "Vacuum" ) 
     102           5 :         call vvac( vacuum, stars, cell, sym, input, field, psq, den%vacz(:,:,ispin), vCoul%vacz(:,:,ispin), rhobar, sig1dh, vz1dh )
     103           5 :         call vvacis( stars, vacuum, sym, cell, psq, input, field, vCoul%vacxy(:,:,:,ispin) )
     104           5 :         call vvacxy( stars, vacuum, cell, sym, input, field, den%vacxy(:,:,:,ispin), vCoul%vacxy(:,:,:,ispin), alphm )
     105           5 :         call timestop( "Vacuum" )
     106             :       end if
     107             : 
     108             : 
     109             : 
     110             :       ! INTERSTITIAL POTENTIAL
     111         173 :       call timestart( "interstitial" )
     112         173 :       write( 6, fmt=8010 )
     113             : 8010  format (/,5x,'coulomb potential in the interstitial region:')
     114             :       ! in case of a film:
     115         173 :       if ( input%film .and. .not.oneD%odi%d1 ) then
     116             :         ! create v(z) for each 2-d reciprocal vector
     117           5 :         ivfft = 3 * stars%mx3 
     118           5 :         ani = 1.0 / real( ivfft )
     119        1130 :         do irec2 = 1, stars%ng2
     120        1125 :           i = 0
     121       72075 :           do i3 = 0, ivfft - 1
     122       70950 :             i = i + 1
     123       70950 :             z = cell%amat(3,3) * i3 * ani
     124       70950 :             if ( z > cell%amat(3,3) / 2. ) z = z - cell%amat(3,3)
     125             :             vintcza = vintcz( stars, vacuum, cell, sym, input, field, z, irec2, psq, &
     126             :                               vCoul%vacxy(:,:,:,ispin), vCoul%vacz(:,:,ispin), &
     127       70950 :                               rhobar, sig1dh, vz1dh, alphm )
     128       70950 :             af1(i) = real( vintcza )
     129       72075 :             bf1(i) = aimag( vintcza )
     130             :           end do
     131             :           !                z = (i_sm-1)*ani
     132             :           !                IF (z > 0.5) z = z - 1.0
     133             :           !                af1(i_sm) = af1(i_sm) + z * delta
     134             :           !                bf1(i_sm) = bf1(i_sm) + z * deltb
     135             :           !              ENDDO
     136             :           !            ENDIF
     137             : 
     138             : 
     139             :           !        --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
     140        1125 :           call cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
     141             :           !            delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
     142        1125 :           i = 0
     143        1130 :           do i3 = 0, ivfft - 1
     144       70950 :             k3 = i3
     145       70950 :             if ( k3 > floor( ivfft / 2. ) ) k3 = k3 - ivfft
     146       70950 :             i = i + 1
     147       72075 :             if ( ( k3 >= -stars%mx3 ) .and. ( k3 <= stars%mx3 ) ) then
     148       48425 :               irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
     149             :               !                 IF ( (irec2 == 1).AND.(i3 > 0) ) THEN                 ! smooth potential
     150             :               !                   corr = 2.0*mod(abs(k3),2) - 1.0
     151             :               !                   bf1(i) = bf1(i) + delta * corr / k3
     152             :               !                 ENDIF
     153             :               !       ----> only stars within g_max sphere (shz oct.97)
     154       48425 :               if ( irec3 /= 0 ) then
     155       31463 :                 xint = cmplx( af1(i), bf1(i) ) * ani
     156       31463 :                 nzst1 = stars%nstr(irec3) / stars%nstr2(irec2)
     157       31463 :                 vCoul%pw(irec3,1) = vCoul%pw(irec3,1) + xint / nzst1
     158             :               end if
     159             :             end if
     160             :           end do
     161             :         end do
     162             :       ! in case of a bulk system:
     163         168 :       else if ( .not. input%film ) then
     164         168 :         if ( vCoul%potdenType == POTDEN_TYPE_POTYUK ) then
     165             :           vCoul%pw(1:stars%ng3,ispin) = fpi_const * psq(1:stars%ng3) &
     166           3 :             / ( stars%sk3(1:stars%ng3) ** 2 + input%preconditioning_param ** 2 )
     167             :           ! if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
     168             :           ! there is a better option now using qfix in mix
     169             :         else
     170         165 :           vCoul%pw(1,ispin) = cmplx(0.0,0.0) 
     171         165 :           vCoul%pw(2:stars%ng3,ispin) = fpi_const * psq(2:stars%ng3) / stars%sk3(2:stars%ng3) ** 2
     172             :         end if
     173             :       end if
     174         173 :     call timestop("interstitial")
     175             :     end if ! mpi%irank == 0
     176             : 
     177             : 
     178             : 
     179             :     ! MUFFIN-TIN POTENTIAL
     180         346 :     call timestart( "MT-spheres" )
     181             : #ifdef CPP_MPI
     182         346 :     CALL MPI_BARRIER(mpi%mpi_comm,ierr) !should be totally useless, but needed anyway????
     183         346 :     call MPI_BCAST( vcoul%pw, size(vcoul%pw), MPI_DOUBLE_COMPLEX, 0, mpi%mpi_comm, ierr )
     184         346 :     CALL MPI_BARRIER(mpi%mpi_comm,ierr) !should be totally useless, but ...
     185             : #endif
     186             :     call vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vCoul%pw(:,ispin), &
     187         346 :                den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin) )
     188         346 :     call timestop( "MT-spheres" )
     189             : 
     190             : 
     191             : 
     192         376 :     if( vCoul%potdenType == POTDEN_TYPE_POTYUK ) return
     193             : 
     194             : 
     195             : 
     196         340 :     if ( mpi%irank == 0 ) then
     197         170 :       CHECK_CONTINUITY: if ( input%vchk ) then
     198           5 :         call timestart( "checking" )
     199             :         call checkDOPAll( input, dimension, sphhar, stars, atoms, sym, vacuum, oneD, &
     200           5 :                           cell, vCoul, ispin )
     201           5 :         call timestop( "checking" )
     202             :       end if CHECK_CONTINUITY
     203             : 
     204         170 :       CALCULATE_DENSITY_POTENTIAL_INTEGRAL: if ( present( results ) ) then
     205         170 :         if ( input%total ) then
     206         170 :           call timestart( "den-pot integrals" )
     207             :           !     CALCULATE THE INTEGRAL OF n*Vcoulomb
     208         170 :           write( 6, fmt=8020 )
     209             : 8020      format (/,10x,'density-coulomb potential integrals',/)
     210             :           !       interstitial first
     211             :           !       convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
     212         170 :           call convol( stars, vCoul%pw_w(:,ispin), vCoul%pw(:,ispin), stars%ufft )
     213         170 :           results%te_vcoul = 0.0
     214             :           call int_nv( ispin, stars, vacuum, atoms, sphhar, cell, sym, input, oneD, &
     215         170 :                        vCoul, den, results%te_vcoul )
     216             : 
     217         170 :           write( 6, fmt=8030 ) results%te_vcoul
     218             : 8030      format (/,10x,'total density-coulomb potential integral :', t40,f20.10)
     219             : 
     220         170 :           call timestop( "den-pot integrals" )
     221             :         end if
     222             :       end if CALCULATE_DENSITY_POTENTIAL_INTEGRAL
     223             :     end if !irank==0
     224             : 
     225         382 :   end subroutine vgen_coulomb
     226             : 
     227             : end module m_vgen_coulomb

Generated by: LCOV version 1.13