VYukawaFilm.f90 17.2 KB
 Miriam Hinzen committed Oct 01, 2018 1 2 ``````module m_VYukawaFilm `````` Miriam Hinzen committed Oct 17, 2018 3 4 5 6 7 8 9 10 11 12 13 `````` ! Computation of the film-case Yukawa potential for the preconditioning of the charge density residual ! 1. part: The pseudo-charge density is used for the vacuum and interstitial potentials. ! 4. part: The MT potential is calculated at the end. ! These two parts do not change compared to the bulk case. ! 2. and 3. part: The vacuum and interstitial potentials are calculated by integration of the product of the charge density and ! the 1D Green function with respect to z (from minus infinity to infinity). The integrals for both potentials are split in 4: ! the seperating points are the boundaries of the slab and the point where the integration variable z' equals the point z where we ! want to evaluate the potential for. For the contribution from the vacuum charge density we do a numerical integration. For the ! contribution from the slab charge density we have analytical expressions of the integrals. `````` Miriam Hinzen committed Oct 01, 2018 14 15 16 17 `````` contains `````` Miriam Hinzen committed Oct 01, 2018 18 `````` `````` Miriam Hinzen committed Oct 01, 2018 19 20 21 22 23 24 `````` subroutine VYukawaFilm( stars, vacuum, cell, sym, input, mpi, atoms, sphhar, dimension, oneD, den, & VYukawa ) use m_constants use m_types use m_psqpw `````` Miriam Hinzen committed Oct 01, 2018 25 `````` use m_vmts `````` Miriam Hinzen committed Oct 01, 2018 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 `````` implicit none type(t_stars), intent(in) :: stars type(t_vacuum), intent(in) :: vacuum type(t_cell), intent(in) :: cell type(t_sym), intent(in) :: sym type(t_input), intent(in) :: input type(t_mpi), intent(in) :: mpi type(t_atoms), intent(in) :: atoms type(t_sphhar), intent(in) :: sphhar type(t_dimension), intent(in) :: dimension type(t_oneD), intent(in) :: oneD type(t_potden), intent(in) :: den type(t_potden), intent(inout) :: VYukawa complex :: psq(stars%ng3) `````` Miriam Hinzen committed Oct 01, 2018 43 `````` complex :: alphm(stars%ng2,2) `````` Miriam Hinzen committed Oct 01, 2018 44 45 46 47 48 49 50 51 52 `````` ! PSEUDO-CHARGE DENSITY call psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, den%pw(:,1), den%mt(:,:,:,1), den%vacz(:,:,1), .false., VYukawa%potdenType, psq ) ! VACUUM POTENTIAL `````` Miriam Hinzen committed Oct 01, 2018 53 54 55 `````` call VYukawaFilmVacuum( stars, vacuum, cell, sym, input, & psq, den%vacxy(:,:,:,1), den%vacz(:,:,1), & VYukawa%vacxy, VYukawa%vacz, alphm ) `````` Miriam Hinzen committed Oct 01, 2018 56 57 58 59 `````` ! INTERSTITIAL POTENTIAL `````` Miriam Hinzen committed Oct 01, 2018 60 61 62 `````` call VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, & psq, VYukawa%vacxy, VYukawa%vacz, alphm, & VYukawa%pw(:,1) ) `````` Miriam Hinzen committed Oct 01, 2018 63 64 65 66 `````` ! MUFFIN-TIN POTENTIAL `````` Miriam Hinzen committed Oct 01, 2018 67 `````` call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, VYukawa%mt(:,0:,:,1) ) `````` Miriam Hinzen committed Oct 01, 2018 68 69 70 71 72 `````` end subroutine VYukawaFilm `````` Miriam Hinzen committed Oct 01, 2018 73 `````` `````` Miriam Hinzen committed Oct 01, 2018 74 75 76 77 `````` subroutine VYukawaFilmVacuum( stars, vacuum, cell, sym, input, & psq, rhtxy, rht, & VVxy, VVz, alphm ) `````` Miriam Hinzen committed Oct 17, 2018 78 79 `````` ! 1. part: Compute the contribution from the interstitial charge density to the vacuum potential as a function of q_xy and z (analytic expression for integral) ! 2. part: Compute the contribution from the vacuum charge density to the vacuum potential as a function of q_xy and z by numerical integration `````` Miriam Hinzen committed Oct 01, 2018 80 `````` `````` Miriam Hinzen committed Oct 01, 2018 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 `````` use m_ExpSave use m_constants use m_types use m_intgr, only: intgz1Reverse use m_qsf implicit none type(t_stars), intent(in) :: stars type(t_vacuum), intent(in) :: vacuum type(t_cell), intent(in) :: cell type(t_sym), intent(in) :: sym type(t_input), intent(in) :: input complex, intent(in) :: psq(stars%ng3) complex, intent(in) :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2) real, intent(in) :: rht(vacuum%nmzd,2) complex, intent(out) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) ! this is the qxy /= 0 part of the vacuum potential real, intent(out) :: VVz(vacuum%nmzd,2) ! this is the qxy = 0 part of the vacuum potential complex, intent(out) :: alphm(stars%ng2,2) ! these are the integrals in upper and lower vacuum, now including the first star---integral for star ig2 is in alphm(ig2,ivac) complex :: sum_qz(2,stars%ng2) complex :: c_ph(-stars%mx3:stars%mx3,stars%ng2) complex :: signIqz real :: g_damped(stars%ng2), qz, sign, vcons(stars%ng2) real :: exp_m(vacuum%nmzd,stars%ng2), exp_p(vacuum%nmzd,stars%ng2) real :: expDhg(stars%ng2), expDg(stars%ng2) real :: z(vacuum%nmzd) `````` Miriam Hinzen committed Oct 17, 2018 108 `````` integer :: iz, irec2, irec3, ivac, iqz `````` Miriam Hinzen committed Oct 01, 2018 109 110 111 112 113 114 `````` complex :: fa(vacuum%nmzxyd,2:stars%ng2), fb(vacuum%nmzxyd,2:stars%ng2) complex :: alpha(vacuum%nmzxyd,2:stars%ng2,2), beta(vacuum%nmzxyd,2:stars%ng2,2), gamma(vacuum%nmzxyd,2:stars%ng2) real :: ga(vacuum%nmzd), gb(vacuum%nmzd) real :: delta(vacuum%nmzd,2), epsilon(vacuum%nmzd,2), zeta(vacuum%nmzd) `````` Miriam Hinzen committed Oct 17, 2018 115 116 `````` ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS `````` Miriam Hinzen committed Oct 01, 2018 117 `````` do iz = 1, vacuum%nmz `````` Miriam Hinzen committed Oct 17, 2018 118 `````` z(iz) = ( iz - 1 ) * vacuum%delz `````` Miriam Hinzen committed Oct 01, 2018 119 120 121 122 123 `````` end do do irec2 = 1, stars%ng2 g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 ) vcons(irec2) = tpi_const / g_damped(irec2) do iz = 1, vacuum%nmz `````` Miriam Hinzen committed Oct 17, 2018 124 `````` exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) ) `````` Miriam Hinzen committed Oct 01, 2018 125 126 127 `````` exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) ) end do expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) ) `````` Miriam Hinzen committed Oct 17, 2018 128 `````` expDg(irec2) = exp_save( -2 * cell%z1 * g_damped(irec2) ) `````` Miriam Hinzen committed Oct 01, 2018 129 130 131 132 133 134 135 136 137 138 139 `````` end do sum_qz = (0.,0.) VVxy = (0.,0.) VVz = 0. ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY do irec2 = 1, stars%ng2 do ivac = 1, vacuum%nvac sign = 3. - 2. * ivac `````` Miriam Hinzen committed Oct 17, 2018 140 141 `````` do iqz = -stars%mx3, stars%mx3 irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz) `````` Miriam Hinzen committed Oct 01, 2018 142 `````` ! use only stars within the g_max sphere -> stars outside the sphere have per definition index ig3n = 0 `````` Miriam Hinzen committed Oct 01, 2018 143 `````` if ( irec3 /= 0 ) then `````` Miriam Hinzen committed Oct 17, 2018 144 145 146 147 `````` c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz) qz = iqz * cell%bmat(3,3) signIqz = sign * ImagUnit * qz sum_qz(ivac,irec2) = sum_qz(ivac,irec2) + c_ph(iqz,irec2) * psq(irec3) * ( exp( signIqz * cell%z1 ) - exp( - signIqz * cell%z1 ) * expDg(irec2) ) / ( signIqz + g_damped(irec2) ) `````` Miriam Hinzen committed Oct 01, 2018 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 `````` endif enddo if( irec2 /= 1 ) then VVxy(1:vacuum%nmzxy,irec2,ivac) = vcons(irec2) * sum_qz(ivac,irec2) * exp_m(1:vacuum%nmzxy,irec2) else VVz(1:vacuum%nmz,ivac) = vcons(1) * sum_qz(ivac,1) * exp_m(1:vacuum%nmz,1) end if enddo enddo ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY ! case irec2 > 1: do irec2 = 2, stars%ng2 do ivac = 1, vacuum%nvac `````` Miriam Hinzen committed Oct 01, 2018 164 `````` ! integrands: `````` Miriam Hinzen committed Oct 17, 2018 165 166 `````` fa(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_m(1:vacuum%nmzxy,irec2) * expDhg(irec2) fb(1:vacuum%nmzxy,irec2) = rhtxy(1:vacuum%nmzxy,irec2-1,ivac) * exp_p(1:vacuum%nmzxy,irec2) * expDhg(irec2) `````` Miriam Hinzen committed Oct 01, 2018 167 168 169 170 171 172 173 `````` ! integrals: ! alpha(z,q_xy,ivac) = int_z^infty rho(z',q_xy,ivac) exp(-sqrt(q_xy**2+prec_param**2)*z') dz' ! beta (z,q_xy,ivac) = int_{D/2}^z rho(z',q_xy,ivac) exp(+sqrt(q_xy**2+prec_param**2)*z') dz' ! where for z < 0 the lower vacuum charge density (ivac=2) is defined by rho(q_xy,z,ivac=2) := rho(q_xy,-z,ivac=2) call intgz1Reverse( fa(:,irec2), vacuum%delz, vacuum%nmzxy, alpha(:,irec2,ivac), .true. ) call qsfComplex( vacuum%delz, fb(:,irec2), beta(:,irec2,ivac), vacuum%nmzxy, 1 ) ! alphm(q_xy,ivac) = alpha(D/2,q_xy,ivac) --- these integrals are also needed for the interstitial potential `````` Miriam Hinzen committed Oct 01, 2018 174 175 `````` alphm(irec2,ivac) = alpha(1,irec2,ivac) end do `````` Miriam Hinzen committed Oct 01, 2018 176 `````` if ( vacuum%nvac == 1 ) then `````` Miriam Hinzen committed Oct 01, 2018 177 178 179 180 181 182 183 184 185 186 `````` if ( sym%invs ) then alphm(irec2,2) = conjg( alphm(irec2,1) ) else alphm(irec2,2) = alphm(irec2,1) end if end if do ivac = 1, vacuum%nvac gamma(1:vacuum%nmzxy,irec2) = exp_m(1:vacuum%nmzxy,irec2) * ( alphm(irec2,mod(ivac,2)+1) + beta(1:vacuum%nmzxy,irec2,ivac) ) & + exp_p(1:vacuum%nmzxy,irec2) * alpha(1:vacuum%nmzxy,irec2,ivac) ! mod(ivac,2)+1 outputs the other ivac value where ( 2. * gamma(:,irec2) == gamma(:,irec2) ) gamma(:,irec2) = cmplx( 0., 0. ) `````` Miriam Hinzen committed Oct 17, 2018 187 `````` VVxy(1:vacuum%nmzxy,irec2,ivac) = VVxy(1:vacuum%nmzxy,irec2,ivac) + vcons(irec2) * gamma(1:vacuum%nmzxy,irec2) * expDhg(irec2) `````` Miriam Hinzen committed Oct 01, 2018 188 189 190 191 192 `````` end do end do ! case irec2 = 1: do ivac = 1, vacuum%nvac `````` Miriam Hinzen committed Oct 17, 2018 193 194 `````` ga(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_m(1:vacuum%nmz,1) * expDhg(1) gb(1:vacuum%nmz) = rht(1:vacuum%nmz,ivac) * exp_p(1:vacuum%nmz,1) * expDhg(1) `````` Miriam Hinzen committed Oct 01, 2018 195 196 197 `````` call intgz1Reverse( ga(:), vacuum%delz, vacuum%nmz, delta(:,ivac), .true. ) ! integrals call qsf( vacuum%delz, gb(:), epsilon(:,ivac), vacuum%nmz, 1 ) alphm(1,ivac) = delta(1,ivac) `````` Miriam Hinzen committed Oct 01, 2018 198 `````` end do `````` Miriam Hinzen committed Oct 17, 2018 199 `````` if ( vacuum%nvac == 1 ) alphm(1,2) = alphm(1,1) `````` Miriam Hinzen committed Oct 01, 2018 200 201 202 203 `````` do ivac = 1, vacuum%nvac zeta(1:vacuum%nmz) = exp_m(1:vacuum%nmz,1) * ( alphm(1,mod(ivac,2)+1) + epsilon(1:vacuum%nmz,ivac) ) & + exp_p(1:vacuum%nmz,1) * delta(1:vacuum%nmz,ivac) where ( 2. * zeta == zeta ) zeta = 0. `````` Miriam Hinzen committed Oct 17, 2018 204 `````` VVz(1:vacuum%nmz,ivac) = VVz(1:vacuum%nmz,ivac) + vcons(1) * zeta(1:vacuum%nmz) * expDhg(1) `````` Miriam Hinzen committed Oct 01, 2018 205 206 207 208 209 210 `````` end do end subroutine VYukawaFilmVacuum `````` Miriam Hinzen committed Oct 01, 2018 211 `````` `````` Miriam Hinzen committed Oct 01, 2018 212 213 214 215 `````` subroutine VYukawaFilmInterstitial( stars, vacuum, cell, sym, input, & psq, VVxy, VVz, alphm, & VIq ) `````` Miriam Hinzen committed Oct 17, 2018 216 217 218 219 220 221 222 223 224 `````` ! main parts: ! 1. part: Compute the contribution from the interstitial charge density to the interstitial potential as a function of q_xy and z (analytic expression for integral) ! 2. part: Add the contribution from the vacuum charge density to the interstitial potential, which had already been computed earlier for the vacuum potential ! 4. part: Compute the coefficients V^I(q_xy,q_z) from the function V^I(q_xy,z) by a 1D Fourier transform: ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z ) ! In order to be able to match the interstitial and vacuum potentials smoothly at the interface, the Fourier transform is done ! in a slightly larger region. -> 3. part ! 3. part: Interpolate the vacuum potential in a small region surrounding the slab `````` Miriam Hinzen committed Oct 01, 2018 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 `````` use m_ExpSave use m_constants use m_types use m_cfft implicit none type(t_stars), intent(in) :: stars type(t_vacuum), intent(in) :: vacuum type(t_cell), intent(in) :: cell type(t_sym), intent(in) :: sym type(t_input), intent(in) :: input complex, intent(in) :: psq(stars%ng3) complex, intent(in) :: VVxy(vacuum%nmzxyd,2:stars%ng2,2) real, intent(in) :: VVz(vacuum%nmzd,2) complex, intent(in) :: alphm(stars%ng2,2) complex, intent(out) :: VIq(stars%ng3) `````` Miriam Hinzen committed Oct 17, 2018 244 245 246 `````` real :: partitioning, rz, qz, q integer :: irec2, irec3, iz, jz, ivac, iqz, nfft, nzmax, nzmin, nzdh, nLower, nUpper, jvac complex, allocatable :: VIz(:,:), eta(:,:) `````` Miriam Hinzen committed Oct 01, 2018 247 248 `````` complex :: VIqz(-stars%mx3:stars%mx3,stars%ng2), c_ph(-stars%mx3:stars%mx3,stars%ng2) complex :: vcons1(stars%ng3) `````` Miriam Hinzen committed Oct 17, 2018 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 `````` real :: VIzReal(3*stars%mx3,stars%ng2), VIzImag(3*stars%mx3,stars%ng2) real, allocatable :: exp_m(:,:), exp_p(:,:) real, allocatable :: z(:) real :: g_damped(stars%ng2), vcons2(stars%ng2), expDhg(stars%ng2) ! DEFINITIONS / ALLOCATIONS / INITIALISATIONS ! grid points z_i: nfft = 3 * stars%mx3 ! number of grid points for Fourier transform partitioning = 1. / real( nfft ) nzmax = nfft / 2 ! index of maximal z below D~/2 nzmin = nzmax - nfft + 1 ! index of minimal z above -D~/2 nzdh = ceiling( cell%z1 / cell%amat(3,3) * nfft ) - 1 ! index of maximal z below D/2 allocate( z(nzmin:nzmax) ) ! definition of z_i: ! indexing: z_0 = 0; positive indices for z > 0; negative indices for z < 0 do iz = nzmin, nzmax z(iz) = cell%amat(3,3) * iz * partitioning `````` Miriam Hinzen committed Oct 01, 2018 267 `````` end do `````` Miriam Hinzen committed Oct 17, 2018 268 269 270 `````` ! other variables: allocate( VIz(nzmin:nzmax,stars%ng2), eta(nzmin:nzmax,stars%ng2) ) allocate( exp_m(nzmin:nzmax,stars%ng2), exp_p(nzmin:nzmax,stars%ng2) ) `````` Miriam Hinzen committed Oct 01, 2018 271 272 273 `````` do irec2 = 1, stars%ng2 g_damped(irec2) = sqrt( stars%sk2(irec2) ** 2 + input%preconditioning_param ** 2 ) vcons2(irec2) = -1. / ( 2. * g_damped(irec2) ) `````` Miriam Hinzen committed Oct 17, 2018 274 275 276 `````` do iz = nzmin, nzmax exp_m(iz,irec2) = exp_save( - g_damped(irec2) * z(iz) ) exp_p(iz,irec2) = exp_save( g_damped(irec2) * z(iz) ) `````` Miriam Hinzen committed Oct 01, 2018 277 `````` end do `````` Miriam Hinzen committed Oct 17, 2018 278 `````` expDhg(irec2) = exp_save( - cell%z1 * g_damped(irec2) ) `````` Miriam Hinzen committed Oct 01, 2018 279 280 281 282 283 284 285 286 287 `````` end do do irec3 = 1, stars%ng3 vcons1(irec3) = fpi_const * psq(irec3) / ( stars%sk3(irec3) ** 2 + input%preconditioning_param ** 2 ) end do VIz = (0.,0.) VIq = (0.,0.) ! CONTRIBUTION FROM THE INTERSTITIAL CHARGE DENSITY `````` Miriam Hinzen committed Oct 17, 2018 288 `````` `````` Miriam Hinzen committed Oct 01, 2018 289 `````` ! compute V^I(q_xy,z) as a function of q_xy and z `````` Miriam Hinzen committed Oct 17, 2018 290 291 292 293 294 295 296 297 298 299 300 301 302 `````` do irec2 = 1, stars%ng2 do iz = -nzdh, nzdh do iqz = -stars%mx3, stars%mx3 irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz) if ( irec3 /= 0 ) then ! use only stars within the g_max sphere c_ph(iqz,irec2) = stars%rgphs(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz) qz = iqz * cell%bmat(3,3) VIz(iz,irec2) = VIz(iz,irec2) & + vcons1(irec3) * c_ph(iqz,irec2) * & ( exp( ImagUnit * qz * z(iz) ) & + vcons2(irec2) * expDhg(irec2) * & ( ( g_damped(irec2) + ImagUnit * qz ) * exp_p(iz,irec2) * exp( ImagUnit * qz * cell%z1 ) & + ( g_damped(irec2) - ImagUnit * qz ) * exp_m(iz,irec2) * exp( - ImagUnit * qz * cell%z1 ) ) ) `````` Miriam Hinzen committed Oct 01, 2018 303 `````` end if `````` Miriam Hinzen committed Oct 17, 2018 304 `````` enddo `````` Miriam Hinzen committed Oct 01, 2018 305 `````` end do `````` Miriam Hinzen committed Oct 17, 2018 306 `````` ! irec2 loop continues `````` Miriam Hinzen committed Oct 01, 2018 307 308 `````` `````` Miriam Hinzen committed Oct 17, 2018 309 `````` ! CONTRIBUTION FROM THE VACUUM CHARGE DENSITY `````` Miriam Hinzen committed Oct 01, 2018 310 `````` `````` Miriam Hinzen committed Oct 17, 2018 311 312 `````` ! irec2 loop continues eta(-nzdh:nzdh,irec2) = exp_m(-nzdh:nzdh,irec2) * alphm(irec2,2) + exp_p(-nzdh:nzdh,irec2) * alphm(irec2,1) `````` Miriam Hinzen committed Oct 01, 2018 313 `````` where ( 2.0 * eta(:,irec2) == eta(:,irec2) ) eta(:,irec2) = cmplx( 0.0, 0.0 ) `````` Miriam Hinzen committed Oct 17, 2018 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 `````` VIz(-nzdh:nzdh,irec2) = VIz(-nzdh:nzdh,irec2) + tpi_const / g_damped(irec2) * eta(-nzdh:nzdh,irec2) ! irec2 loop continues ! INTERPOLATION IN VACUUM REGION ! use Lagrange polynomials of order 3 to interpolate the vacuum potential outside I ! q, q-1 and q-2 (scaled with +/-1 or +/-0.5) are the factors of the Lagrange basis polynomials ! irec2 loop continues do ivac = 1, 2 select case( ivac ) case( 1 ) nUpper = nzmax; nLower = nzdh + 1; jvac = 1 case( 2 ) nUpper = nzmin; nLower = -nzdh - 1; jvac = 2; if ( sym%invs .or. sym%zrfs ) jvac = 1 end select do iz = nLower, nUpper rz = ( abs( z(iz) ) - cell%z1 ) / vacuum%delz + 1.0 jz = rz ! index of maximal vacuum grid point below z_i q = rz - jz ! factor in Lagrange basis polynomials if ( irec2 == 1 ) then VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVz(jz, jvac) & - q * ( q - 2. ) * VVz(jz+1,jvac) & + 0.5 * q * ( q - 1. ) * VVz(jz+2,jvac) else if ( jz + 2 <= vacuum%nmzxy ) then VIz(iz,irec2) = 0.5 * ( q - 1. ) * ( q - 2. ) * VVxy(jz, irec2,jvac) & - q * ( q - 2. ) * VVxy(jz+1,irec2,jvac) & + 0.5 * q * ( q - 1. ) * VVxy(jz+2,irec2,jvac) if ( ( sym%invs .and. .not. sym%zrfs ) .and. ivac == 2 ) VIz(iz,irec2) = conjg( VIz(iz,irec2) ) end if end do end do end do ! irec2 `````` Miriam Hinzen committed Oct 01, 2018 348 `````` `````` Miriam Hinzen committed Oct 17, 2018 349 `````` ! 1D FOURIER TRANSFORM TO FIND THE COEFFICIENTS V^I(q_xy,q_z) `````` Miriam Hinzen committed Oct 01, 2018 350 `````` `````` Miriam Hinzen committed Oct 17, 2018 351 352 353 354 355 356 357 358 `````` ! change the indexing for the subroutine cfft, and split real and imaginary parts VIzReal(1:nzmax+1,:) = real( VIz(0:nzmax,:) ); VIzReal(nzmax+2:nfft,:) = real( VIz(nzmin:-1,:) ) VIzImag(1:nzmax+1,:) = aimag( VIz(0:nzmax,:) ); VIzImag(nzmax+2:nfft,:) = aimag( VIz(nzmin:-1,:) ) ! V^I(q_xy,z) = sum_{q_z} V^I(q_xy,q_z) * exp( ImagUnit * q_z * z ) do irec2 = 1, stars%ng2 call cfft( VIzReal(:,irec2), VIzImag(:,irec2), nfft, nfft, nfft, -1 ) ! irec2 loop continues `````` Miriam Hinzen committed Oct 01, 2018 359 `````` `````` Miriam Hinzen committed Oct 17, 2018 360 `````` ! reorder `````` Miriam Hinzen committed Oct 01, 2018 361 `````` VIqz(0,irec2) = cmplx( VIzReal(1,irec2), VIzImag(1,irec2) ) `````` Miriam Hinzen committed Oct 17, 2018 362 363 364 `````` do iqz = 1, stars%mx3 VIqz( iqz,irec2) = cmplx( VIzReal(iqz+1, irec2), VIzImag(iqz+1, irec2) ) VIqz(-iqz,irec2) = cmplx( VIzReal(nfft+1-iqz,irec2), VIzImag(nfft+1-iqz,irec2) ) `````` Miriam Hinzen committed Oct 01, 2018 365 `````` end do `````` Miriam Hinzen committed Oct 17, 2018 366 367 368 369 370 `````` ! add the computed components to V^I(q_xy,q_z): do iqz= -stars%mx3, stars%mx3 irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),iqz) if ( irec3 /= 0 ) VIq(irec3) = VIq(irec3) + VIqz(iqz,irec2) * partitioning / ( stars%nstr(irec3) / stars%nstr2(irec2) ) `````` Miriam Hinzen committed Oct 01, 2018 371 372 373 374 375 376 377 `````` end do end do end subroutine VYukawaFilmInterstitial `````` Miriam Hinzen committed Oct 01, 2018 378 ``end module m_VYukawaFilm``