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 1763 : 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 , ncsh
72 : real :: ani, g3, z , sigmaa(2)
73 : complex :: sig1dh, vz1dh, vmz1dh, vmz1dh_is, constantShift
74 : complex :: sigma_loc(2), sigma_loc2(2)
75 1763 : complex, allocatable :: alphm(:,:), psq(:)
76 1763 : 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 , l_gradientEfield
80 :
81 : #ifdef CPP_MPI
82 : integer:: ierr
83 : #endif
84 :
85 1763 : l_dfptvgen = PRESENT(stars2)
86 1763 : l_2ndord = PRESENT(iDir2)
87 1763 : l_corr = .FALSE. !ALL(ABS(den%vac)<1e-12)
88 1763 : vmz1dh_is = cmplx(0.0,0.0)
89 1763 : sigma_loc = CMPLX(0.0,0.0) !sigma_disc
90 1763 : sigma_loc2 = CMPLX(0.0,0.0) !MERGE(sigma_disc,cmplx(0.0,0.0),PRESENT(sigma_disc2))
91 1763 : l_gradientEfield = .FALSE.
92 1763 : if (present(iDtype)) then
93 2500290 : if (iDtype==0 .and. .not. ALL(ABS(den%pw)<1e-12)) l_gradientEfield = .TRUE.
94 : end if
95 :
96 1763 : vintcza = cmplx(0.0,0.0)
97 1763 : sig1dh = cmplx(0.0,0.0)
98 1763 : vz1dh = cmplx(0.0,0.0)
99 1763 : vmz1dh = cmplx(0.0,0.0)
100 1763 : vslope = cmplx(0.0,0.0)
101 1763 : constantShift = cmplx(0.0,0.0)
102 :
103 1763 : allocate ( alphm(stars%ng2,2), af1(3*stars%mx3), bf1(3*stars%mx3), psq(stars%ng3) )
104 1763 : sigma = (cell%amat(3,3) - 2*cell%z1)/8.0
105 1763 : vCoul%iter = den%iter
106 :
107 : ! PSEUDO-CHARGE DENSITY COEFFICIENTS
108 1763 : call timestart( "psqpw" )
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 1763 : & dfptdenimag, stars2, iDtype, iDir, dfptden0, iDir2 )
115 1763 : call timestop( "psqpw" )
116 :
117 : ! VACUUM POTENTIAL
118 1763 : if ( fmpi%irank == 0 ) then
119 1241 : if ( input%film ) then
120 : ! ----> potential in the vacuum region
121 475 : call timestart( "Vacuum" )
122 1900 : if ((.not.l_dfptvgen).or.norm2(stars%center)<1e-8) then
123 : ! If we do DPFT AND q/=0, there is no G_||+q_||=0 part! So all components are
124 : ! handled by the G_||/=0 parts in vvacis/vvacxy, that are told to explicitly
125 : ! start at star 1 instead of 2 for this!
126 259 : 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,l_dfptvgen )
127 : end if
128 475 : call vvacis( stars, vacuum, cell, psq, input, field, vCoul%vac(:vacuum%nmzxyd,:,:,ispin), l_dfptvgen, l_corr )
129 475 : call vvacxy( stars, vacuum, cell, sym, input, field, den%vac(:vacuum%nmzxyd,:,:,ispin), vCoul%vac(:vacuum%nmzxyd,:,:,ispin), alphm, l_dfptvgen )
130 475 : call timestop( "Vacuum" )
131 475 : if ( l_dfptvgen .AND. juphon%l_symVacLevel ) constantShift = (vCoul%vac(vacuum%nmzd,1,1,ispin) - vCoul%vac(vacuum%nmzd,1,2,ispin)) / 2
132 : end if
133 :
134 : ! INTERSTITIAL POTENTIAL
135 1241 : call timestart( "interstitial" )
136 1241 : write(oUnit, fmt=8010 )
137 : 8010 format (/,5x,'coulomb potential in the interstitial region:')
138 : ! in case of a film:
139 1241 : if ( input%film ) then
140 : ! create v(z) for each 2-d reciprocal vector
141 475 : ivfft = 3 * stars%mx3
142 475 : ani = 1.0 / real( ivfft )
143 66679 : do irec2 = 1, stars%ng2
144 : ! If we do DFPT, we want to fix the second vacuum at infinity to 0. This is WIP,
145 : ! as to how we want to do it eventually. Here, we calculate the necessary offset.
146 : !IF (l_dfptvgen.AND.irec2 == 1) vmz1dh_is = vintcz( stars, vacuum, cell, input, field, -cell%z1+1e-15, irec2, psq, &
147 : ! vCoul%vac(:,:,:,ispin), &
148 : ! rhobar, sig1dh, vz1dh, alphm, vslope, .TRUE., CMPLX(0.0,0.0) )
149 66204 : i = 0
150 4167723 : do i3 = 0, ivfft - 1
151 4101519 : i = i + 1
152 4101519 : z = cell%amat(3,3) * i3 * ani
153 4101519 : if (l_dfptvgen) gaussian = 1 - exp( -(z-cell%amat(3,3)/2.0)**2 / ( 2 * sigma**2 ) )
154 4101519 : if ( z > cell%amat(3,3) / 2. ) z = z - cell%amat(3,3)
155 : vintcza = vintcz( stars, vacuum, cell, input, field, z, irec2, psq, &
156 : vCoul%vac(:,:,:,ispin), &
157 4101519 : rhobar, sig1dh, vz1dh, alphm, vslope, sigma_disc, l_dfptvgen, l_corr, vmz1dh-vmz1dh_is )
158 4101519 : if (l_dfptvgen) THEN
159 3547926 : if (juphon%l_symVacLevel .AND. (irec2 == 1) ) vintcza = vintcza + constantShift
160 3547926 : vintcza = vintcza * gaussian
161 : end if
162 4101519 : af1(i) = real( vintcza )
163 4167723 : bf1(i) = aimag( vintcza )
164 : end do
165 : ! z = (i_sm-1)*ani
166 : ! IF (z > 0.5) z = z - 1.0
167 : ! af1(i_sm) = af1(i_sm) + z * delta
168 : ! bf1(i_sm) = bf1(i_sm) + z * deltb
169 : ! ENDDO
170 : ! ENDIF
171 :
172 :
173 : ! --> 1-d fourier transform and store the coefficients in vTot%pw( ,1)
174 66204 : call cfft( af1, bf1, ivfft, ivfft, ivfft, -1 )
175 : ! delta = ivfft * delta * 2 / fpi ! * amat(3,3)**2 * ani
176 66204 : i = 0
177 4168198 : do i3 = 0, ivfft - 1
178 4101519 : k3 = i3
179 4101519 : if ( k3 > floor( ivfft / 2. ) ) k3 = k3 - ivfft
180 4101519 : i = i + 1
181 4167723 : if ( ( k3 >= -stars%mx3 ) .and. ( k3 <= stars%mx3 ) ) then
182 2800550 : irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
183 : ! IF ( (irec2 == 1).AND.(i3 > 0) ) THEN ! smooth potential
184 : ! corr = 2.0*mod(abs(k3),2) - 1.0
185 : ! bf1(i) = bf1(i) + delta * corr / k3
186 : ! ENDIF
187 : ! ----> only stars within g_max sphere (shz oct.97)
188 2800550 : if ( irec3 /= 0 ) then
189 2477248 : xint = cmplx( af1(i), bf1(i) ) * ani
190 2477248 : nzst1 = stars%nstr(irec3) / stars%nstr2(irec2)
191 2477248 : vCoul%pw(irec3,1) = vCoul%pw(irec3,1) + xint / nzst1
192 : end if
193 : end if
194 : end do
195 : end do
196 475 : if (l_dfptvgen .AND. juphon%l_symVacLevel) vCoul%vac(:,1,:,ispin) = vCoul%vac(:,1,:,ispin) + constantShift
197 475 : sigma_disc = sigma_loc
198 475 : if (l_gradientEfield) then
199 18 : if (iDir == 3 ) then
200 6 : sigmaa(1) = ( field%efield%sigma + field%efield%sig_b(1) ) / cell%area
201 6 : sigmaa(2) = ( field%efield%sigma + field%efield%sig_b(2) ) / cell%area
202 :
203 6 : if ( (sigmaa(1) + sigmaa(2)) .lt. 1E-8 ) then
204 : ! Asymmetric setup in which an Efield contribution exists inside the film
205 12 : vCoul%pw(1,:) = vCoul%pw(1,:) - fpi_const * sigmaa(1) !/ cell%omtil
206 : end if
207 :
208 : ! vacuum contribution
209 : ! obtain mesh point (ncsh) of charge sheet for external electric field
210 6 : ncsh = field%efield%zsigma / vacuum%delz + 1.01
211 612 : do i = 1 , ncsh
212 : ! vacuum 1
213 1212 : vCoul%vac(i,1,1,:) = vCoul%vac(i,1,1,:) - fpi_const * sigmaa(1)
214 : ! vacuum 2
215 1218 : vCoul%vac(i,1,2,:) = vCoul%vac(i,1,2,:) + fpi_const * sigmaa(2)
216 : end do
217 : end if
218 : end if
219 : ! in case of a bulk system:
220 766 : else if ( .not. input%film ) then
221 766 : if ( vCoul%potdenType == POTDEN_TYPE_POTYUK ) then
222 : vCoul%pw(1:stars%ng3,ispin) = fpi_const * psq(1:stars%ng3) &
223 3840 : / ( stars%sk3(1:stars%ng3) ** 2 + input%preconditioning_param ** 2 )
224 : ! if( abs( real( psq(1) ) ) * cell%omtil < 0.01 ) vCoul%pw(1,ispin) = 0.0
225 : ! there is a better option now using qfix in mix
226 : else
227 763 : vCoul%pw(1,ispin) = cmplx(0.0,0.0)
228 763 : first_star = MERGE(2,1,stars%sk3(1)< 1E-9)
229 1020717 : vCoul%pw(first_star:stars%ng3,ispin) = fpi_const * psq(first_star:stars%ng3) / stars%sk3(first_star:stars%ng3) ** 2
230 : end if
231 : end if
232 1241 : call timestop("interstitial")
233 : end if ! fmpi%irank == 0
234 :
235 : ! MUFFIN-TIN POTENTIAL
236 1763 : call timestart( "MT-spheres" )
237 : #ifdef CPP_MPI
238 1763 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr) !should be totally useless, but needed anyway????
239 5289 : call MPI_BCAST( vcoul%pw, size(vcoul%pw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
240 1763 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr) !should be totally useless, but ...
241 : #endif
242 :
243 : call vmts( input, fmpi, stars, sphhar, atoms, sym, cell, juphon, dosf, vCoul%pw(:,ispin), &
244 : den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin), ispin, &
245 1763 : dfptdenimag, dfptvCoulimag, iDtype, iDir, iDir2 )
246 1763 : call timestop( "MT-spheres" )
247 :
248 1763 : if( vCoul%potdenType == POTDEN_TYPE_POTYUK ) return
249 :
250 1757 : if ( fmpi%irank == 0 ) then
251 1238 : CHECK_CONTINUITY: if ( input%vchk ) then ! TODO: We could use this for DFPT as well if we
252 5 : call timestart( "checking" ) ! passed an optional to checkDOPAll and modded
253 : call checkDOPAll( input, sphhar, stars, atoms, sym, vacuum, & ! slightly.
254 5 : cell, vCoul, ispin )
255 5 : call timestop( "checking" )
256 : end if CHECK_CONTINUITY
257 :
258 1238 : CALCULATE_DENSITY_POTENTIAL_INTEGRAL: if ( present( results ) ) then
259 349 : call timestart( "den-pot integrals" )
260 : ! CALCULATE THE INTEGRAL OF n*Vcoulomb
261 349 : write(oUnit, fmt=8020 )
262 : 8020 format (/,10x,'density-coulomb potential integrals',/)
263 : ! interstitial first
264 : ! convolute ufft and pot: F(G) = \sum_(G') U(G - G') V(G')
265 349 : call convol( stars, vCoul%pw_w(:,ispin), vCoul%pw(:,ispin))
266 349 : results%te_vcoul = 0.0
267 : call int_nv( ispin, stars, vacuum, atoms, sphhar, cell, sym, input, &
268 349 : vCoul, den, results%te_vcoul )
269 :
270 349 : write(oUnit, fmt=8030 ) results%te_vcoul
271 : 8030 format (/,10x,'total density-coulomb potential integral :', t40,f20.10)
272 :
273 349 : call timestop( "den-pot integrals" )
274 : end if CALCULATE_DENSITY_POTENTIAL_INTEGRAL
275 : end if !irank==0
276 5289 : end subroutine vgen_coulomb
277 : end module m_vgen_coulomb
|
