vgen_coulomb.F90 9.51 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
             cell, sphhar, atoms, dosf, den, vCoul, results )
15 16 17 18 19 20 21 22 23 24 25
    !----------------------------------------------------------------------------
    ! 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
    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 
53
    LOGICAL,            INTENT(IN)               :: dosf
54 55 56 57 58 59
    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
60
    integer                                      :: lh, n, nzst1
61 62 63 64 65
    integer                                      :: imz, imzxy, ichsmrg, ivfft
    integer                                      :: l, nat
    real                                         :: ani, g3, z, sig1dh, vz1dh
    complex, allocatable                         :: alphm(:,:), psq(:)
    real,    allocatable                         :: af1(:), bf1(:)
66 67 68 69
#ifdef CPP_MPI
    include 'mpif.h'
    integer:: ierr
#endif
70 71 72 73
 
    
    allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3)  )
    vCoul%iter = den%iter
74 75 76



77 78
    ! PSEUDO-CHARGE DENSITY COEFFICIENTS
    call timestart( "psqpw" )      
79
    call psqpw( mpi, atoms, sphhar, stars, vacuum,  cell, input, sym, oneD, &
80 81 82 83 84 85 86 87 88 89 90 91 92 93
         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
94
        call od_vvacis( oneD%odi%n2d, input, vacuum, oneD%odi%nq2, &
95 96 97 98
             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" )
99 100
        !---> generation of the vacuum warped potential components and       
      elseif ( input%film .and. .not. oneD%odi%d1 ) then
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 137
        !     ----> 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
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 167
          !        --> 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 )
168 169
          ! 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
170
        else
Matthias Redies's avatar
Matthias Redies committed
171
          vCoul%pw(1,ispin) = cmplx(0.0,0.0) 
172 173 174 175 176 177 178 179
          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" )
180
#ifdef CPP_MPI
181 182 183
    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 ...
184
#endif
185
    call vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, dosf, vCoul%pw(:,ispin), &
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
               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