vgen_coulomb.F90 9.61 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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    ! PSEUDO-CHARGE DENSITY COEFFICIENTS
    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., 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
shift    
Matthias Redies committed
170
          vCoul%pw(1,ispin) = cmplx(-2.0,0.0) ! that's my line
Matthias Redies's avatar
Matthias Redies committed
171
          write (*,*) "shifting vtot by", vCoul%pw(1,ispin)
172
173
174
175
176
177
178
179
180
181
          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" )
182
#ifdef CPP_MPI
183
184
185
    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 ...
186
#endif
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
228
229
230
    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 )
          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')
          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
          write( 16, 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