vgen_coulomb.F90 9.78 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
MODULE m_vgen_coulomb
8

9
  USE m_juDFT
10

11
CONTAINS
12 13 14

  SUBROUTINE vgen_coulomb( ispin, mpi, DIMENSION, oneD, input, field, vacuum, sym, stars, &
             cell, sphhar, atoms, den, yukawa_residual, vCoul, results )
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
    !     ***********************************************************
    !     FLAPW potential generator                           
    !     ***********************************************************
    !     Generates the Coulomb potential and optionally the density-potential integrals
    !     is takes a spin variable to indicate in which spin-channel the charge resides...
    !     
    !     ***********************************************************
    USE m_constants
    USE m_vmts
    USE m_intnv
    USE m_vvac
    USE m_vvacis
    USE m_vvacxy
    USE m_vintcz
    USE m_checkdopall
    USE m_types
    USE m_od_vvac
    USE m_od_vvacis
    USE m_convol
    USE m_psqpw
    USE m_cfft
    IMPLICIT NONE

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
    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
    logical,            intent(in)               :: yukawa_residual
    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
    integer                                      :: lh, n, nzst1
    INTEGER                                      :: imz, imzxy, ichsmrg, ivfft
    INTEGER                                      :: l, nat
    REAL                                         :: ani, g3, z, sig1dh, vz1dh
    COMPLEX, ALLOCATABLE                         :: alphm(:,:), psq(:)
    REAL,    ALLOCATABLE                         :: af1(:), bf1(:)
63 64 65 66 67
#ifdef CPP_MPI
    include 'mpif.h'
    integer:: ierr
#endif

68
    ALLOCATE ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3)  )
69 70 71 72 73 74 75 76


    vCoul%iter = den%iter

    !     ************** coulomb potential ***********************
    !     ----> create pesudo-charge density coefficients


77 78 79 80
    CALL timestart( "psqpw" )      
    CALL psqpw( mpi, atoms, sphhar, stars, vacuum, DIMENSION, cell, input, sym, oneD, &
         den%pw(:,ispin), den%mt(:,:,:,ispin), den%vacz(:,:,ispin), .FALSE., yukawa_residual, psq )
    CALL timestop( "psqpw" )
81 82 83 84 85
    IF (mpi%irank == 0) THEN
       !First the vacuum
       !     ------------------------------------------
       IF (oneD%odi%d1) THEN
          !-odim
86
          CALL timestart( "Vacuum" )
87 88

          !---> generates the m=0,gz=0 component of the vacuum potential
89
          CALL od_vvac( stars, vacuum, cell, psq, den%vacz(:,:,ispin), vCoul%vacz(:,:,ispin) )
90 91 92 93 94

          !---> generation of the vacuum warped potential components and
          !---> interstitial pw potential
          !---> vvacxy_5.F is a symmetrized potential generator

95 96 97 98 99
          CALL od_vvacis( oneD%odi%n2d, dimension, vacuum, oneD%odi%nq2, &
               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" )
100 101 102 103 104

          !+odim
       ELSEIF (input%film .AND. .NOT.oneD%odi%d1) THEN
          !     ----> potential in the  vacuum  region
          !       
105 106 107
          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) )
108

109 110
          CALL vvacxy( stars, vacuum, cell, sym, input, field, den%vacxy(:,:,:,ispin), vCoul%vacxy(:,:,:,ispin), alphm )
          CALL timestop( "Vacuum" )
111 112 113
       END IF
       !     ------------------------------------------
       !     ----> potential in the  interstitial  region
114
       CALL timestart( "interstitial" )
115 116 117 118 119 120 121 122 123 124 125 126 127
       WRITE (6,FMT=8010)
8010   FORMAT (/,5x,'coulomb potential in the interstitial region:')
       IF (input%film .AND. .NOT.oneD%odi%d1) THEN
          !           -----> create v(z) for each 2-d reciprocal vector
          ivfft =  3*stars%mx3 
          !         ivfft = 2*mx3 - 1
          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.GT.cell%amat(3,3)/2.) z = z - cell%amat(3,3)
128 129 130 131 132
                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 )
133 134 135 136 137 138 139 140
             ENDDO
             !                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
             !        --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
141
             CALL cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
             !            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.GE.-stars%mx3).AND.(k3.LE.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.NE.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
                ENDIF
             ENDDO
          ENDDO
166 167 168 169 170 171 172 173
        else if ( .not. input%film ) then
          if ( yukawa_residual ) then
            vCoul%pw(1:stars%ng3,ispin) = fpi_const * psq(1:stars%ng3) / ( stars%sk3(1:stars%ng3) ** 2 + input%preconditioning_param ** 2 )
            if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
          else
            vCoul%pw(1,ispin) = cmplx(0.0,0.0)
            vCoul%pw(2:stars%ng3,ispin) = fpi_const * psq(2:stars%ng3) / ( stars%sk3(2:stars%ng3) * stars%sk3(2:stars%ng3) )       
          end if
174 175 176 177 178 179 180 181
       END IF

       CALL timestop("interstitial")

    ENDIF ! mpi%irank == 0
    !     --------------------------------------------
    !     ---> potential in the muffin-tin spheres

182
    CALL timestart( "MT-spheres" )
183
#ifdef CPP_MPI
184
    CALL MPI_BCAST( den%mt, atoms%jmtd*(1+sphhar%nlhd) * atoms%ntype * dimension%jspd, MPI_DOUBLE_PRECISION, 0, mpi%mpi_comm, ierr )
185
#endif
186
    CALL vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vCoul%pw(:,ispin), den%mt(:,0:,:,ispin), &
187 188
         yukawa_residual, vCoul%mt(:,0:,:,ispin) )
    CALL timestop( "MT-spheres" )
189

190 191
    if( yukawa_residual ) return

192 193 194
    IF (mpi%irank == 0) THEN
       !     ---> check continuity of coulomb potential
       IF (input%vchk) THEN
195 196 197 198
          CALL timestart( "checking" )
          CALL checkDOPAll( input, DIMENSION, sphhar, stars, atoms, sym, vacuum, oneD, &
               cell, vCoul, ispin )
          CALL timestop( "checking" )
199 200 201 202 203
       END IF

       !Calculate potential-density integrals if results-variable is present!
       IF (PRESENT(results)) THEN
          IF (input%total) THEN
204
             CALL timestart( "den-pot integrals" )
205 206 207 208 209 210
             !     CALCULATE THE INTEGRAL OF n*Vcoulomb
             WRITE (6,FMT=8020)
             WRITE (16,FMT=8020)
8020         FORMAT (/,10x,'density-coulomb potential integrals',/)
             !       interstitial first
             !       convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
211
             CALL convol( stars, vCoul%pw_w(:,ispin), vCoul%pw(:,ispin), stars%ufft )
212
             results%te_vcoul = 0.0
213 214
             CALL int_nv( ispin, stars, vacuum, atoms, sphhar, cell, sym, input, oneD, &
                  vCoul, den, results%te_vcoul )
215 216 217 218 219

             WRITE (6,FMT=8030) results%te_vcoul
             WRITE (16,FMT=8030) results%te_vcoul
8030         FORMAT (/,10x,'total density-coulomb potential integral :', t40,f20.10)

220
             CALL timestop( "den-pot integrals" )
221 222 223 224 225
          END IF
       ENDIF
    ENDIF !irank==0

  END SUBROUTINE vgen_coulomb
226

227
END MODULE m_vgen_coulomb