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 682 : 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, constantShift
74 : complex :: mat2ord(5,3,3), sigma_loc(2), sigma_loc2(2)
75 682 : complex, allocatable :: alphm(:,:), psq(:)
76 682 : real, allocatable :: af1(:), bf1(:)
77 : real :: gaussian, sigma ! smoothing function in case of DFPT + Film
78 : LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
79 : LOGICAL :: l_2ndord, l_corr
80 :
81 : #ifdef CPP_MPI
82 : integer:: ierr
83 : #endif
84 :
85 682 : l_dfptvgen = PRESENT(stars2)
86 682 : l_2ndord = PRESENT(iDir2)
87 682 : l_corr = .FALSE. !ALL(ABS(den%vac)<1e-12)
88 682 : vmz1dh_is = cmplx(0.0,0.0)
89 682 : sigma_loc = CMPLX(0.0,0.0) !sigma_disc
90 682 : sigma_loc2 = CMPLX(0.0,0.0) !MERGE(sigma_disc,cmplx(0.0,0.0),PRESENT(sigma_disc2))
91 :
92 682 : vintcza = cmplx(0.0,0.0)
93 682 : sig1dh = cmplx(0.0,0.0)
94 682 : vz1dh = cmplx(0.0,0.0)
95 682 : vmz1dh = cmplx(0.0,0.0)
96 682 : vslope = cmplx(0.0,0.0)
97 682 : constantShift = cmplx(0.0,0.0)
98 :
99 682 : allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3) )
100 682 : sigma = (cell%amat(3,3) - 2*cell%z1)/8.0
101 682 : vCoul%iter = den%iter
102 :
103 : ! PSEUDO-CHARGE DENSITY COEFFICIENTS
104 682 : call timestart( "psqpw" )
105 682 : if (.not.l_dfptvgen) then
106 : call psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, &
107 682 : & juphon, den, ispin, .false., vCoul%potdenType, psq, sigma_loc )
108 0 : else if (.not.l_2ndord) then
109 : ! If we do DFPT, the MT density perturbation has an imaginary part that needs to be explicitly carried
110 : ! as another variable dfptdenimag%mt and results in the same component for the Coulomb potential later on.
111 : ! Also, the ionic qlm behave differently.
112 : call psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, &
113 : & juphon, den, ispin, .false., vCoul%potdenType, psq, sigma_loc,&
114 0 : & dfptdenimag%mt(:,:,:,ispin), stars2, iDtype, iDir, dfptden0%mt(:,:,:,ispin), dfptden0%pw(:,ispin) )
115 : else
116 0 : call make_mat_2nd(mat2ord)
117 : call psqpw( fmpi, atoms, sphhar, stars, vacuum, cell, input, sym, &
118 : & juphon, den, ispin, .false., vCoul%potdenType, psq, sigma_loc,&
119 0 : & dfptdenimag%mt(:,:,:,ispin), stars2, iDtype, iDir, dfptden0%mt(:,:,:,ispin), dfptden0%pw(:,ispin), iDir2, mat2ord )
120 : end if
121 682 : call timestop( "psqpw" )
122 :
123 : ! VACUUM POTENTIAL
124 682 : if ( fmpi%irank == 0 ) then
125 341 : if ( input%film ) then
126 : ! ----> potential in the vacuum region
127 37 : call timestart( "Vacuum" )
128 148 : if ((.not.l_dfptvgen).or.norm2(stars%center)<1e-8) then
129 : ! If we do DPFT AND q/=0, there is no G_||+q_||=0 part! So all components are
130 : ! handled by the G_||/=0 parts in vvacis/vvacxy, that are told to explicitly
131 : ! start at star 1 instead of 2 for this!
132 37 : if (.not.l_2ndord) then
133 : !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 )
134 37 : 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 )
135 : else
136 : !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 )
137 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 )
138 : end if
139 : end if
140 37 : call vvacis( stars, vacuum, cell, psq, input, field, vCoul%vac(:vacuum%nmzxyd,:,:,ispin), l_dfptvgen, l_corr )
141 37 : call vvacxy( stars, vacuum, cell, sym, input, field, den%vac(:vacuum%nmzxyd,:,:,ispin), vCoul%vac(:vacuum%nmzxyd,:,:,ispin), alphm, l_dfptvgen )
142 37 : call timestop( "Vacuum" )
143 37 : if ( l_dfptvgen .AND. juphon%l_symVacLevel ) constantShift = (vCoul%vac(vacuum%nmzd,1,1,ispin) - vCoul%vac(vacuum%nmzd,1,2,ispin)) / 2
144 : end if
145 :
146 : ! INTERSTITIAL POTENTIAL
147 341 : call timestart( "interstitial" )
148 341 : write(oUnit, fmt=8010 )
149 : 8010 format (/,5x,'coulomb potential in the interstitial region:')
150 : ! in case of a film:
151 341 : if ( input%film ) then
152 : ! create v(z) for each 2-d reciprocal vector
153 37 : ivfft = 3 * stars%mx3
154 37 : ani = 1.0 / real( ivfft )
155 13487 : do irec2 = 1, stars%ng2
156 : ! If we do DFPT, we want to fix the second vacuum at infinity to 0. This is WIP,
157 : ! as to how we want to do it eventually. Here, we calculate the necessary offset.
158 : !IF (l_dfptvgen.AND.irec2 == 1) vmz1dh_is = vintcz( stars, vacuum, cell, input, field, -cell%z1+1e-15, irec2, psq, &
159 : ! vCoul%vac(:,:,:,ispin), &
160 : ! rhobar, sig1dh, vz1dh, alphm, vslope, .TRUE., CMPLX(0.0,0.0) )
161 13450 : i = 0
162 532435 : do i3 = 0, ivfft - 1
163 518985 : i = i + 1
164 518985 : z = cell%amat(3,3) * i3 * ani
165 518985 : if (l_dfptvgen) gaussian = 1 - exp( -(z-cell%amat(3,3)/2.0)**2 / ( 2 * sigma**2 ) )
166 518985 : if ( z > cell%amat(3,3) / 2. ) z = z - cell%amat(3,3)
167 518985 : if (.not.l_2ndord) then
168 : vintcza = vintcz( stars, vacuum, cell, input, field, z, irec2, psq, &
169 : vCoul%vac(:,:,:,ispin), &
170 518985 : rhobar, sig1dh, vz1dh, alphm, vslope, sigma_disc, l_dfptvgen, l_corr, vmz1dh-vmz1dh_is )
171 : else
172 : vintcza = vintcz( stars, vacuum, cell, input, field, z, irec2, psq, &
173 : vCoul%vac(:,:,:,ispin), &
174 0 : rhobar, sig1dh, vz1dh, alphm, vslope, sigma_disc, l_dfptvgen, l_corr, vmz1dh-vmz1dh_is, sigma_disc2 )
175 : end if
176 518985 : if (l_dfptvgen) THEN
177 0 : if (juphon%l_symVacLevel .AND. (irec2 == 1) ) vintcza = vintcza + constantShift
178 0 : vintcza = vintcza * gaussian
179 : end if
180 518985 : af1(i) = real( vintcza )
181 532435 : bf1(i) = aimag( vintcza )
182 : end do
183 : ! z = (i_sm-1)*ani
184 : ! IF (z > 0.5) z = z - 1.0
185 : ! af1(i_sm) = af1(i_sm) + z * delta
186 : ! bf1(i_sm) = bf1(i_sm) + z * deltb
187 : ! ENDDO
188 : ! ENDIF
189 :
190 :
191 : ! --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
192 13450 : call cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
193 : ! delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
194 13450 : i = 0
195 532472 : do i3 = 0, ivfft - 1
196 518985 : k3 = i3
197 518985 : if ( k3 > floor( ivfft / 2. ) ) k3 = k3 - ivfft
198 518985 : i = i + 1
199 532435 : if ( ( k3 >= -stars%mx3 ) .and. ( k3 <= stars%mx3 ) ) then
200 359440 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
201 : ! IF ( (irec2 == 1).AND.(i3 > 0) ) THEN ! smooth potential
202 : ! corr = 2.0*mod(abs(k3),2) - 1.0
203 : ! bf1(i) = bf1(i) + delta * corr / k3
204 : ! ENDIF
205 : ! ----> only stars within g_max sphere (shz oct.97)
206 359440 : if ( irec3 /= 0 ) then
207 226326 : xint = cmplx( af1(i), bf1(i) ) * ani
208 226326 : nzst1 = stars%nstr(irec3) / stars%nstr2(irec2)
209 226326 : vCoul%pw(irec3,1) = vCoul%pw(irec3,1) + xint / nzst1
210 : end if
211 : end if
212 : end do
213 : end do
214 37 : if (l_dfptvgen .AND. juphon%l_symVacLevel) vCoul%vac(:,1,:,ispin) = vCoul%vac(:,1,:,ispin) + constantShift
215 37 : sigma_disc = sigma_loc
216 : ! in case of a bulk system:
217 304 : else if ( .not. input%film ) then
218 304 : if ( vCoul%potdenType == POTDEN_TYPE_POTYUK ) then
219 : vCoul%pw(1:stars%ng3,ispin) = fpi_const * psq(1:stars%ng3) &
220 3840 : / ( stars%sk3(1:stars%ng3) ** 2 + input%preconditioning_param ** 2 )
221 : ! if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
222 : ! there is a better option now using qfix in mix
223 : else
224 301 : vCoul%pw(1,ispin) = cmplx(0.0,0.0)
225 301 : first_star = MERGE(2,1,stars%sk3(1)< 1E-9)
226 608773 : vCoul%pw(first_star:stars%ng3,ispin) = fpi_const * psq(first_star:stars%ng3) / stars%sk3(first_star:stars%ng3) ** 2
227 : end if
228 : end if
229 341 : call timestop("interstitial")
230 : end if ! fmpi%irank == 0
231 :
232 : ! MUFFIN-TIN POTENTIAL
233 682 : call timestart( "MT-spheres" )
234 : #ifdef CPP_MPI
235 682 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr) !should be totally useless, but needed anyway????
236 2046 : call MPI_BCAST( vcoul%pw, size(vcoul%pw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
237 682 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr) !should be totally useless, but ...
238 : #endif
239 :
240 682 : IF (.NOT.l_dfptvgen) THEN
241 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, juphon, dosf, vCoul%pw(:,ispin), &
242 682 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin) )
243 0 : ELSE IF (.NOT.l_2ndord) THEN
244 : ! For DFPT there is a) an imaginary part to the potential and b) a different treatment
245 : ! for the ionic 1/r (now 1/r^2) contribution.
246 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, juphon, dosf, vCoul%pw(:,ispin), &
247 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin), &
248 0 : dfptdenimag%mt(:,0:,:,ispin), dfptvCoulimag%mt(:,0:,:,ispin), iDtype, iDir )
249 : ELSE
250 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, juphon, dosf, vCoul%pw(:,ispin), &
251 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin), &
252 0 : dfptdenimag%mt(:,0:,:,ispin), dfptvCoulimag%mt(:,0:,:,ispin), iDtype, iDir, iDir2, mat2ord )
253 : END IF
254 682 : call timestop( "MT-spheres" )
255 :
256 682 : if( vCoul%potdenType == POTDEN_TYPE_POTYUK ) return
257 :
258 676 : if ( fmpi%irank == 0 ) then
259 338 : CHECK_CONTINUITY: if ( input%vchk ) then ! TODO: We could use this for DFPT as well if we
260 4 : call timestart( "checking" ) ! passed an optional to checkDOPAll and modded
261 : call checkDOPAll( input, sphhar, stars, atoms, sym, vacuum, & ! slightly.
262 4 : cell, vCoul, ispin )
263 4 : call timestop( "checking" )
264 : end if CHECK_CONTINUITY
265 :
266 338 : CALCULATE_DENSITY_POTENTIAL_INTEGRAL: if ( present( results ) ) then
267 337 : call timestart( "den-pot integrals" )
268 : ! CALCULATE THE INTEGRAL OF n*Vcoulomb
269 337 : write(oUnit, fmt=8020 )
270 : 8020 format (/,10x,'density-coulomb potential integrals',/)
271 : ! interstitial first
272 : ! convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
273 337 : call convol( stars, vCoul%pw_w(:,ispin), vCoul%pw(:,ispin))
274 337 : results%te_vcoul = 0.0
275 : call int_nv( ispin, stars, vacuum, atoms, sphhar, cell, sym, input, &
276 337 : vCoul, den, results%te_vcoul )
277 :
278 337 : write(oUnit, fmt=8030 ) results%te_vcoul
279 : 8030 format (/,10x,'total density-coulomb potential integral :', t40,f20.10)
280 :
281 337 : call timestop( "den-pot integrals" )
282 : end if CALCULATE_DENSITY_POTENTIAL_INTEGRAL
283 : end if !irank==0
284 2046 : end subroutine vgen_coulomb
285 :
286 0 : subroutine make_mat_2nd(mat2ord)
287 : use m_constants
288 : complex, intent(out) :: mat2ord(5,3,3)
289 0 : mat2ord = cmplx(0.0,0.0)
290 :
291 0 : mat2ord(1,1,1) = sqrt(3.0/2.0)
292 0 : mat2ord(1,1,2) = sqrt(3.0/2.0)*ImagUnit
293 0 : mat2ord(1,2,1) = sqrt(3.0/2.0)*ImagUnit
294 0 : mat2ord(1,2,2) = -sqrt(3.0/2.0)
295 :
296 0 : mat2ord(2,1,3) = sqrt(3.0/2.0)
297 0 : mat2ord(2,2,3) = sqrt(3.0/2.0)*ImagUnit
298 0 : mat2ord(2,3,1) = sqrt(3.0/2.0)
299 0 : mat2ord(2,3,2) = sqrt(3.0/2.0)*ImagUnit
300 :
301 0 : mat2ord(3,1,1) = -1.0! - 0.5 ! TODO: What the hell is this value???
302 0 : mat2ord(3,2,2) = -1.0! - 0.5 ! TODO: What the hell is this value???
303 0 : mat2ord(3,3,3) = 2.0! - 0.5 ! TODO: What the hell is this value???
304 :
305 0 : mat2ord(4,1,3) = -sqrt(3.0/2.0)
306 0 : mat2ord(4,2,3) = sqrt(3.0/2.0)*ImagUnit
307 0 : mat2ord(4,3,1) = -sqrt(3.0/2.0)
308 0 : mat2ord(4,3,2) = sqrt(3.0/2.0)*ImagUnit
309 :
310 0 : mat2ord(5,1,1) = sqrt(3.0/2.0)
311 0 : mat2ord(5,1,2) = -sqrt(3.0/2.0)*ImagUnit
312 0 : mat2ord(5,2,1) = -sqrt(3.0/2.0)*ImagUnit
313 0 : mat2ord(5,2,2) = -sqrt(3.0/2.0)
314 :
315 0 : mat2ord = sqrt(4.0*pi_const/5.0) * mat2ord
316 0 : end subroutine
317 :
318 : end module m_vgen_coulomb
|
