vgen_coulomb.F90 9.45 KB
Newer Older
1 2 3 4 5
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
6

7 8 9 10 11 12
module m_vgen_coulomb

  use m_juDFT

contains

13
  subroutine vgen_coulomb( ispin, mpi, dimension, oneD, input, field, vacuum, sym, stars, &
14 15 16 17 18 19 20 21 22 23 24 25
             cell, sphhar, atoms, den, vCoul, results )
    !----------------------------------------------------------------------------
    ! FLAPW potential generator                           
    !----------------------------------------------------------------------------
    ! Generates the Coulomb or Yukawa potential and optionally the 
    ! density-potential integrals
    ! vCoul%potdenType = POTDEN_TYPE_POTYUK -> Yukawa case
    ! It takes a spin variable to indicate in which spin-channel the charge 
    ! resides.
    !----------------------------------------------------------------------------

    use m_constants
26
    use m_types
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
    use m_vmts
    use m_intnv
    use m_vvac
    use m_vvacis
    use m_vvacxy
    use m_vintcz
    use m_checkdopall
    use m_od_vvac
    use m_od_vvacis
    use m_convol
    use m_psqpw
    use m_cfft
    implicit none

    integer,            intent(in)               :: ispin
    type(t_mpi),        intent(in)               :: mpi
    type(t_dimension),  intent(in)               :: dimension
    type(t_oneD),       intent(in)               :: oneD
    type(t_input),      intent(in)               :: input
    type(t_field),      intent(inout)            :: field
    type(t_vacuum),     intent(in)               :: vacuum
    type(t_sym),        intent(in)               :: sym
    type(t_stars),      intent(in)               :: stars
    type(t_cell),       intent(in)               :: cell
    type(t_sphhar),     intent(in)               :: sphhar
    type(t_atoms),      intent(in)               :: atoms 
    type(t_potden),     intent(in)               :: den
    type(t_potden),     intent(inout)            :: vCoul
    type(t_results),    intent(inout), optional  :: results

    complex                                      :: vintcza, xint, rhobar
    integer                                      :: i, i3, irec2, irec3, ivac, j, js, k, k3
59
    integer                                      :: lh, n, nzst1
60 61 62 63 64
    integer                                      :: imz, imzxy, ichsmrg, ivfft
    integer                                      :: l, nat
    real                                         :: ani, g3, z, sig1dh, vz1dh
    complex, allocatable                         :: alphm(:,:), psq(:)
    real,    allocatable                         :: af1(:), bf1(:)
65 66 67 68
#ifdef CPP_MPI
    include 'mpif.h'
    integer:: ierr
#endif
69 70 71 72
 
    
    allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3)  )
    vCoul%iter = den%iter
73 74 75



76 77
    ! PSEUDO-CHARGE DENSITY COEFFICIENTS
    call timestart( "psqpw" )      
78
    call psqpw( mpi, atoms, sphhar, stars, vacuum,  cell, input, sym, oneD, &
79 80 81 82 83 84 85 86 87 88 89 90 91 92
         den%pw(:,ispin), den%mt(:,:,:,ispin), den%vacz(:,:,ispin), .false., vCoul%potdenType, psq )
    call timestop( "psqpw" )



    ! VACUUM POTENTIAL
    if ( mpi%irank == 0 ) then
      if ( oneD%odi%d1 ) then
        call timestart( "Vacuum" )
        !---> generates the m=0,gz=0 component of the vacuum potential
        call od_vvac( stars, vacuum, cell, psq, den%vacz(:,:,ispin), vCoul%vacz(:,:,ispin) )
        !---> generation of the vacuum warped potential components and
        !---> interstitial pw potential
        !---> vvacxy_5.F is a symmetrized potential generator
93
        call od_vvacis( oneD%odi%n2d, input, vacuum, oneD%odi%nq2, &
94 95 96 97
             oneD%odi%kv, cell, oneD%odi%M, stars, oneD%odi%nst2, &
             oneD, den%vacz(:,:,ispin), den%vacxy(:,:,:,ispin), psq, &
             vCoul%vacz(:,:,ispin), sym, vCoul%vacxy(:,:,:,ispin), vCoul%pw(:,ispin) )
        call timestop( "Vacuum" )
98 99
        !---> generation of the vacuum warped potential components and       
      elseif ( input%film .and. .not. oneD%odi%d1 ) then
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
        !     ----> potential in the  vacuum  region
        call timestart( "Vacuum" ) 
        call vvac( vacuum, stars, cell, sym, input, field, psq, den%vacz(:,:,ispin), vCoul%vacz(:,:,ispin), rhobar, sig1dh, vz1dh )
        call vvacis( stars, vacuum, sym, cell, psq, input, field, vCoul%vacxy(:,:,:,ispin) )
        call vvacxy( stars, vacuum, cell, sym, input, field, den%vacxy(:,:,:,ispin), vCoul%vacxy(:,:,:,ispin), alphm )
        call timestop( "Vacuum" )
      end if



      ! INTERSTITIAL POTENTIAL
      call timestart( "interstitial" )
      write( 6, fmt=8010 )
8010  format (/,5x,'coulomb potential in the interstitial region:')
      ! in case of a film:
      if ( input%film .and. .not.oneD%odi%d1 ) then
        ! create v(z) for each 2-d reciprocal vector
        ivfft = 3 * stars%mx3 
        ani = 1.0 / real( ivfft )
        do irec2 = 1, stars%ng2
          i = 0
          do i3 = 0, ivfft - 1
            i = i + 1
            z = cell%amat(3,3) * i3 * ani
            if ( z > cell%amat(3,3) / 2. ) z = z - cell%amat(3,3)
            vintcza = vintcz( stars, vacuum, cell, sym, input, field, z, irec2, psq, &
                              vCoul%vacxy(:,:,:,ispin), vCoul%vacz(:,:,ispin), &
                              rhobar, sig1dh, vz1dh, alphm )
            af1(i) = real( vintcza )
            bf1(i) = aimag( vintcza )
          end do
          !                z = (i_sm-1)*ani
          !                IF (z > 0.5) z = z - 1.0
          !                af1(i_sm) = af1(i_sm) + z * delta
          !                bf1(i_sm) = bf1(i_sm) + z * deltb
          !              ENDDO
          !            ENDIF
137 138


139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166
          !        --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
          call cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
          !            delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
          i = 0
          do i3 = 0, ivfft - 1
            k3 = i3
            if ( k3 > floor( ivfft / 2. ) ) k3 = k3 - ivfft
            i = i + 1
            if ( ( k3 >= -stars%mx3 ) .and. ( k3 <= stars%mx3 ) ) then
              irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
              !                 IF ( (irec2 == 1).AND.(i3 > 0) ) THEN                 ! smooth potential
              !                   corr = 2.0*mod(abs(k3),2) - 1.0
              !                   bf1(i) = bf1(i) + delta * corr / k3
              !                 ENDIF
              !       ----> only stars within g_max sphere (shz oct.97)
              if ( irec3 /= 0 ) then
                xint = cmplx( af1(i), bf1(i) ) * ani
                nzst1 = stars%nstr(irec3) / stars%nstr2(irec2)
                vCoul%pw(irec3,1) = vCoul%pw(irec3,1) + xint / nzst1
              end if
            end if
          end do
        end do
      ! in case of a bulk system:
      else if ( .not. input%film ) then
        if ( vCoul%potdenType == POTDEN_TYPE_POTYUK ) then
          vCoul%pw(1:stars%ng3,ispin) = fpi_const * psq(1:stars%ng3) &
            / ( stars%sk3(1:stars%ng3) ** 2 + input%preconditioning_param ** 2 )
167 168
          ! if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
          ! there is a better option now using qfix in mix
169
        else
Matthias Redies's avatar
Matthias Redies committed
170
          vCoul%pw(1,ispin) = cmplx(0.0,0.0) 
171 172 173 174 175 176 177 178 179 180
          vCoul%pw(2:stars%ng3,ispin) = fpi_const * psq(2:stars%ng3) / stars%sk3(2:stars%ng3) ** 2
        end if
      end if
    call timestop("interstitial")
    end if ! mpi%irank == 0



    ! MUFFIN-TIN POTENTIAL
    call timestart( "MT-spheres" )
181
#ifdef CPP_MPI
182 183 184
    CALL MPI_BARRIER(mpi%mpi_comm,ierr) !should be totally useless, but needed anyway????
    call MPI_BCAST( vcoul%pw, size(vcoul%pw), MPI_DOUBLE_COMPLEX, 0, mpi%mpi_comm, ierr )
    CALL MPI_BARRIER(mpi%mpi_comm,ierr) !should be totally useless, but ...
185
#endif
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
    call vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vCoul%pw(:,ispin), &
               den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin) )
    call timestop( "MT-spheres" )



    if( vCoul%potdenType == POTDEN_TYPE_POTYUK ) return



    if ( mpi%irank == 0 ) then
      CHECK_CONTINUITY: if ( input%vchk ) then
        call timestart( "checking" )
        call checkDOPAll( input, dimension, sphhar, stars, atoms, sym, vacuum, oneD, &
                          cell, vCoul, ispin )
        call timestop( "checking" )
      end if CHECK_CONTINUITY

      CALCULATE_DENSITY_POTENTIAL_INTEGRAL: if ( present( results ) ) then
        if ( input%total ) then
          call timestart( "den-pot integrals" )
          !     CALCULATE THE INTEGRAL OF n*Vcoulomb
          write( 6, fmt=8020 )
8020      format (/,10x,'density-coulomb potential integrals',/)
          !       interstitial first
          !       convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
          call convol( stars, vCoul%pw_w(:,ispin), vCoul%pw(:,ispin), stars%ufft )
          results%te_vcoul = 0.0
          call int_nv( ispin, stars, vacuum, atoms, sphhar, cell, sym, input, oneD, &
                       vCoul, den, results%te_vcoul )

          write( 6, fmt=8030 ) results%te_vcoul
8030      format (/,10x,'total density-coulomb potential integral :', t40,f20.10)

          call timestop( "den-pot integrals" )
        end if
      end if CALCULATE_DENSITY_POTENTIAL_INTEGRAL
    end if !irank==0

  end subroutine vgen_coulomb

end module m_vgen_coulomb