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 : #ifdef CPP_MPI
11 : use mpi
12 : #endif
13 : contains
14 :
15 696 : subroutine vgen_coulomb( ispin, fmpi, input, field, vacuum, sym, juphon, stars, &
16 : cell, sphhar, atoms, dosf, den, vCoul, sigma_disc, results, dfptdenimag, dfptvCoulimag, dfptden0, stars2, iDtype, iDir, iDir2, sigma_disc2 )
17 : !----------------------------------------------------------------------------
18 : ! FLAPW potential generator
19 : !----------------------------------------------------------------------------
20 : ! Generates the Coulomb or Yukawa potential and optionally the
21 : ! density-potential integrals
22 : ! vCoul%potdenType = POTDEN_TYPE_POTYUK -> Yukawa case
23 : ! It takes a spin variable to indicate in which spin-channel the charge
24 : ! resides.
25 : !----------------------------------------------------------------------------
26 :
27 : use m_constants
28 : use m_types
29 : use m_vmts
30 : use m_intnv
31 : use m_vvac
32 : use m_vvacis
33 : use m_vvacxy
34 : use m_vintcz
35 : use m_checkdopall
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) :: fmpi
43 :
44 :
45 : type(t_input), intent(in) :: input
46 : type(t_field), intent(in) :: field
47 : type(t_vacuum), intent(in) :: vacuum
48 : type(t_sym), intent(in) :: sym
49 : type(t_juphon), intent(in) :: juphon
50 : type(t_stars), intent(in) :: stars
51 : type(t_cell), intent(in) :: cell
52 : type(t_sphhar), intent(in) :: sphhar
53 : type(t_atoms), intent(in) :: atoms
54 : LOGICAL, INTENT(IN) :: dosf
55 : type(t_potden), intent(in) :: den
56 : type(t_potden), intent(inout) :: vCoul
57 : COMPLEX, INTENT(INOUT) :: sigma_disc(2)
58 : type(t_results), intent(inout), optional :: results
59 :
60 : TYPE(t_potden), OPTIONAL, INTENT(IN) :: dfptdenimag, dfptden0
61 : TYPE(t_potden), OPTIONAL, INTENT(INOUT) :: dfptvCoulimag
62 : TYPE(t_stars), OPTIONAL, INTENT(IN) :: stars2
63 : INTEGER, OPTIONAL, INTENT(IN) :: iDtype, iDir ! DFPT: Type and direction of displaced atom
64 : INTEGER, OPTIONAL, INTENT(IN) :: iDir2 ! DFPT: 2nd direction for 2nd order VC
65 : COMPLEX, OPTIONAL, INTENT(IN) :: sigma_disc2(2)
66 :
67 : complex :: vintcza, xint, rhobar,vslope
68 : integer :: i, i3, irec2, irec3, ivac, j, js, k, k3
69 : integer :: lh, n, nzst1, first_star
70 : integer :: imz, imzxy, ichsmrg, ivfft
71 : integer :: l, nat
72 : real :: ani, g3, z
73 : complex :: sig1dh, vz1dh, vmz1dh, vmz1dh_is
74 : complex :: mat2ord(5,3,3), sigma_loc(2), sigma_loc2(2)
75 696 : complex, allocatable :: alphm(:,:), psq(:)
76 696 : real, allocatable :: af1(:), bf1(:)
77 : LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
78 : LOGICAL :: l_2ndord, l_corr
79 :
80 : #ifdef CPP_MPI
81 : integer:: ierr
82 : #endif
83 :
84 696 : l_dfptvgen = PRESENT(stars2)
85 696 : l_2ndord = PRESENT(iDir2)
86 4224 : l_corr = ALL(ABS(den%vac)<1e-12)
87 696 : vmz1dh_is = cmplx(0.0,0.0)
88 696 : sigma_loc = sigma_disc
89 2088 : sigma_loc2 = MERGE(sigma_disc,cmplx(0.0,0.0),PRESENT(sigma_disc2))
90 :
91 696 : vintcza = cmplx(0.0,0.0)
92 696 : sig1dh = cmplx(0.0,0.0)
93 696 : vz1dh = cmplx(0.0,0.0)
94 696 : vmz1dh = cmplx(0.0,0.0)
95 696 : vslope = cmplx(0.0,0.0)
96 :
97 696 : allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3) )
98 696 : vCoul%iter = den%iter
99 :
100 : ! PSEUDO-CHARGE DENSITY COEFFICIENTS
101 696 : call timestart( "psqpw" )
102 696 : if (.not.l_dfptvgen) then
103 : call psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, &
104 696 : & juphon, den, ispin, .false., vCoul%potdenType, psq, sigma_loc )
105 0 : else if (.not.l_2ndord) then
106 : ! If we do DFPT, the MT density perturbation has an imaginary part that needs to be explicitly carried
107 : ! as another variable dfptdenimag%mt and results in the same component for the Coulomb potential later on.
108 : ! Also, the ionic qlm behave differently.
109 : call psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, &
110 : & juphon, den, ispin, .false., vCoul%potdenType, psq, sigma_loc,&
111 0 : & dfptdenimag%mt(:,:,:,ispin), stars2, iDtype, iDir, dfptden0%mt(:,:,:,ispin), dfptden0%pw(:,ispin) )
112 : else
113 0 : call make_mat_2nd(mat2ord)
114 : call psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, &
115 : & juphon, den, ispin, .false., vCoul%potdenType, psq, sigma_loc,&
116 0 : & dfptdenimag%mt(:,:,:,ispin), stars2, iDtype, iDir, dfptden0%mt(:,:,:,ispin), dfptden0%pw(:,ispin), iDir2, mat2ord )
117 : end if
118 696 : call timestop( "psqpw" )
119 :
120 : ! VACUUM POTENTIAL
121 696 : if ( fmpi%irank == 0 ) then
122 348 : if ( input%film ) then
123 : ! ----> potential in the vacuum region
124 44 : call timestart( "Vacuum" )
125 176 : if ((.not.l_dfptvgen).or.norm2(stars%center)<1e-8) then
126 : ! If we do DPFT AND q/=0, there is no G_||+q_||=0 part! So all components are
127 : ! handled by the G_||/=0 parts in vvacis/vvacxy, that are told to explicitly
128 : ! start at star 1 instead of 2 for this!
129 44 : if (.not.l_2ndord) then
130 : !call vvac( vacuum, stars, cell, input, field, psq, den%vac(:,1,:,ispin), vCoul%vac(:,1,:,ispin), rhobar, sig1dh, vz1dh,vslope,.FALSE..AND.l_dfptvgen,vmz1dh,sigma_disc )
131 44 : call vvac( vacuum, stars, cell, input, field, psq, den%vac(:,1,:,ispin), vCoul%vac(:,1,:,ispin), rhobar, sig1dh, vz1dh,vslope,l_corr,vmz1dh,sigma_disc )
132 : else
133 : !call vvac( vacuum, stars, cell, input, field, psq, den%vac(:,1,:,ispin), vCoul%vac(:,1,:,ispin), rhobar, sig1dh, vz1dh,vslope,.TRUE.,vmz1dh,sigma_disc,sigma_disc2 )
134 0 : call vvac( vacuum, stars, cell, input, field, psq, den%vac(:,1,:,ispin), vCoul%vac(:,1,:,ispin), rhobar, sig1dh, vz1dh,vslope,l_corr,vmz1dh,sigma_disc,sigma_disc2 )
135 : end if
136 : end if
137 44 : call vvacis( stars, vacuum, cell, psq, input, field, vCoul%vac(:vacuum%nmzxyd,:,:,ispin), l_dfptvgen, l_corr )
138 44 : call vvacxy( stars, vacuum, cell, sym, input, field, den%vac(:vacuum%nmzxyd,:,:,ispin), vCoul%vac(:vacuum%nmzxyd,:,:,ispin), alphm, l_dfptvgen )
139 44 : call timestop( "Vacuum" )
140 : end if
141 :
142 : ! INTERSTITIAL POTENTIAL
143 348 : call timestart( "interstitial" )
144 348 : write(oUnit, fmt=8010 )
145 : 8010 format (/,5x,'coulomb potential in the interstitial region:')
146 : ! in case of a film:
147 348 : if ( input%film ) then
148 : ! create v(z) for each 2-d reciprocal vector
149 44 : ivfft = 3 * stars%mx3
150 44 : ani = 1.0 / real( ivfft )
151 17421 : do irec2 = 1, stars%ng2
152 : ! If we do DFPT, we want to fix the second vacuum at infinity to 0. This is WIP,
153 : ! as to how we want to do it eventually. Here, we calculate the necessary offset.
154 : !IF (l_dfptvgen.AND.irec2 == 1) vmz1dh_is = vintcz( stars, vacuum, cell, input, field, -cell%z1+1e-15, irec2, psq, &
155 : ! vCoul%vac(:,:,:,ispin), &
156 : ! rhobar, sig1dh, vz1dh, alphm, vslope, .TRUE., CMPLX(0.0,0.0) )
157 17377 : i = 0
158 665953 : do i3 = 0, ivfft - 1
159 648576 : i = i + 1
160 648576 : z = cell%amat(3,3) * i3 * ani
161 648576 : if ( z > cell%amat(3,3) / 2. ) z = z - cell%amat(3,3)
162 648576 : if (.not.l_2ndord) then
163 : vintcza = vintcz( stars, vacuum, cell, input, field, z, irec2, psq, &
164 : vCoul%vac(:,:,:,ispin), &
165 648576 : rhobar, sig1dh, vz1dh, alphm, vslope, sigma_disc, l_dfptvgen, l_corr, vmz1dh-vmz1dh_is )
166 : else
167 : vintcza = vintcz( stars, vacuum, cell, input, field, z, irec2, psq, &
168 : vCoul%vac(:,:,:,ispin), &
169 0 : rhobar, sig1dh, vz1dh, alphm, vslope, sigma_disc, l_dfptvgen, l_corr, vmz1dh-vmz1dh_is, sigma_disc2 )
170 : end if
171 648576 : af1(i) = real( vintcza )
172 665953 : bf1(i) = aimag( vintcza )
173 : end do
174 : ! z = (i_sm-1)*ani
175 : ! IF (z > 0.5) z = z - 1.0
176 : ! af1(i_sm) = af1(i_sm) + z * delta
177 : ! bf1(i_sm) = bf1(i_sm) + z * deltb
178 : ! ENDDO
179 : ! ENDIF
180 :
181 :
182 : ! --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
183 17377 : call cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
184 : ! delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
185 17377 : i = 0
186 665997 : do i3 = 0, ivfft - 1
187 648576 : k3 = i3
188 648576 : if ( k3 > floor( ivfft / 2. ) ) k3 = k3 - ivfft
189 648576 : i = i + 1
190 665953 : if ( ( k3 >= -stars%mx3 ) .and. ( k3 <= stars%mx3 ) ) then
191 449761 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
192 : ! IF ( (irec2 == 1).AND.(i3 > 0) ) THEN ! smooth potential
193 : ! corr = 2.0*mod(abs(k3),2) - 1.0
194 : ! bf1(i) = bf1(i) + delta * corr / k3
195 : ! ENDIF
196 : ! ----> only stars within g_max sphere (shz oct.97)
197 449761 : if ( irec3 /= 0 ) then
198 282753 : xint = cmplx( af1(i), bf1(i) ) * ani
199 282753 : nzst1 = stars%nstr(irec3) / stars%nstr2(irec2)
200 282753 : vCoul%pw(irec3,1) = vCoul%pw(irec3,1) + xint / nzst1
201 : end if
202 : end if
203 : end do
204 : end do
205 44 : sigma_disc = sigma_loc
206 : ! in case of a bulk system:
207 304 : else if ( .not. input%film ) then
208 304 : if ( vCoul%potdenType == POTDEN_TYPE_POTYUK ) then
209 : vCoul%pw(1:stars%ng3,ispin) = fpi_const * psq(1:stars%ng3) &
210 3840 : / ( stars%sk3(1:stars%ng3) ** 2 + input%preconditioning_param ** 2 )
211 : ! if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
212 : ! there is a better option now using qfix in mix
213 : else
214 301 : vCoul%pw(1,ispin) = cmplx(0.0,0.0)
215 301 : first_star = MERGE(2,1,stars%sk3(1)< 1E-9)
216 608773 : vCoul%pw(first_star:stars%ng3,ispin) = fpi_const * psq(first_star:stars%ng3) / stars%sk3(first_star:stars%ng3) ** 2
217 : end if
218 : end if
219 348 : call timestop("interstitial")
220 : end if ! fmpi%irank == 0
221 :
222 : ! MUFFIN-TIN POTENTIAL
223 696 : call timestart( "MT-spheres" )
224 : #ifdef CPP_MPI
225 696 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr) !should be totally useless, but needed anyway????
226 2088 : call MPI_BCAST( vcoul%pw, size(vcoul%pw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
227 696 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr) !should be totally useless, but ...
228 : #endif
229 :
230 696 : IF (.NOT.l_dfptvgen) THEN
231 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, dosf, vCoul%pw(:,ispin), &
232 696 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin) )
233 0 : ELSE IF (.NOT.l_2ndord) THEN
234 : ! For DFPT there is a) an imaginary part to the potential and b) a different treatment
235 : ! for the ionic 1/r (now 1/r^2) contribution.
236 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, dosf, vCoul%pw(:,ispin), &
237 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin), &
238 0 : dfptdenimag%mt(:,0:,:,ispin), dfptvCoulimag%mt(:,0:,:,ispin), iDtype, iDir )
239 : ELSE
240 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, dosf, vCoul%pw(:,ispin), &
241 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin), &
242 0 : dfptdenimag%mt(:,0:,:,ispin), dfptvCoulimag%mt(:,0:,:,ispin), iDtype, iDir, iDir2, mat2ord )
243 : END IF
244 696 : call timestop( "MT-spheres" )
245 :
246 696 : if( vCoul%potdenType == POTDEN_TYPE_POTYUK ) return
247 :
248 690 : if ( fmpi%irank == 0 ) then
249 345 : CHECK_CONTINUITY: if ( input%vchk ) then ! TODO: We could use this for DFPT as well if we
250 4 : call timestart( "checking" ) ! passed an optional to checkDOPAll and modded
251 : call checkDOPAll( input, sphhar, stars, atoms, sym, vacuum, & ! slightly.
252 4 : cell, vCoul, ispin )
253 4 : call timestop( "checking" )
254 : end if CHECK_CONTINUITY
255 :
256 345 : CALCULATE_DENSITY_POTENTIAL_INTEGRAL: if ( present( results ) ) then
257 344 : call timestart( "den-pot integrals" )
258 : ! CALCULATE THE INTEGRAL OF n*Vcoulomb
259 344 : write(oUnit, fmt=8020 )
260 : 8020 format (/,10x,'density-coulomb potential integrals',/)
261 : ! interstitial first
262 : ! convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
263 344 : call convol( stars, vCoul%pw_w(:,ispin), vCoul%pw(:,ispin))
264 344 : results%te_vcoul = 0.0
265 : call int_nv( ispin, stars, vacuum, atoms, sphhar, cell, sym, input, &
266 344 : vCoul, den, results%te_vcoul )
267 :
268 344 : write(oUnit, fmt=8030 ) results%te_vcoul
269 : 8030 format (/,10x,'total density-coulomb potential integral :', t40,f20.10)
270 :
271 344 : call timestop( "den-pot integrals" )
272 : end if CALCULATE_DENSITY_POTENTIAL_INTEGRAL
273 : end if !irank==0
274 2088 : end subroutine vgen_coulomb
275 :
276 0 : subroutine make_mat_2nd(mat2ord)
277 : use m_constants
278 : complex, intent(out) :: mat2ord(5,3,3)
279 0 : mat2ord = cmplx(0.0,0.0)
280 :
281 0 : mat2ord(1,1,1) = sqrt(3.0/2.0)
282 0 : mat2ord(1,1,2) = sqrt(3.0/2.0)*ImagUnit
283 0 : mat2ord(1,2,1) = sqrt(3.0/2.0)*ImagUnit
284 0 : mat2ord(1,2,2) = -sqrt(3.0/2.0)
285 :
286 0 : mat2ord(2,1,3) = sqrt(3.0/2.0)
287 0 : mat2ord(2,2,3) = sqrt(3.0/2.0)*ImagUnit
288 0 : mat2ord(2,3,1) = sqrt(3.0/2.0)
289 0 : mat2ord(2,3,2) = sqrt(3.0/2.0)*ImagUnit
290 :
291 0 : mat2ord(3,1,1) = -1.0! - 0.5 ! TODO: What the hell is this value???
292 0 : mat2ord(3,2,2) = -1.0! - 0.5 ! TODO: What the hell is this value???
293 0 : mat2ord(3,3,3) = 2.0! - 0.5 ! TODO: What the hell is this value???
294 :
295 0 : mat2ord(4,1,3) = -sqrt(3.0/2.0)
296 0 : mat2ord(4,2,3) = sqrt(3.0/2.0)*ImagUnit
297 0 : mat2ord(4,3,1) = -sqrt(3.0/2.0)
298 0 : mat2ord(4,3,2) = sqrt(3.0/2.0)*ImagUnit
299 :
300 0 : mat2ord(5,1,1) = sqrt(3.0/2.0)
301 0 : mat2ord(5,1,2) = -sqrt(3.0/2.0)*ImagUnit
302 0 : mat2ord(5,2,1) = -sqrt(3.0/2.0)*ImagUnit
303 0 : mat2ord(5,2,2) = -sqrt(3.0/2.0)
304 :
305 0 : mat2ord = sqrt(4.0*pi_const/5.0) * mat2ord
306 0 : end subroutine
307 :
308 : end module m_vgen_coulomb
|
