!! Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch !! !! This program is free software; you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation; either version 2, or (at your option) !! any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program; if not, write to the Free Software !! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA !! 02111-1307, USA. !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! MODULE POISSON ! ============== ! !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module poisson use mesh use cube_function implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! There are four public procedures in this module; all the rest is private. private public :: poisson_init, & poisson_solve, & poisson_sum, & poisson_fft !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! The next set of parameters fix the way in which the Hartree term is ! calculated: ! (i) performing directly the sum: set mode to HARTREE_SUM ! (ii) the FFT technique: set mode to HARTREE_FFT integer, parameter :: HARTREE_SUM = 1, & HARTREE_FFT = 2 integer, parameter :: mode = HARTREE_FFT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! The next parameters set the type of interaction that defines the system: ! (i) Coulombic: set interaction to COULOMB ! (ii) Yukawa: set interaction to YUKAWA ! In this latter case, gamma is the parameter that defines the Yukawa potential. integer, parameter :: COULOMB = 1, & YUKAWA = 2 integer, parameter :: interaction = COULOMB real(8), parameter :: gamma = 2._8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! These variables are only meant for internal use within the module. type(dcf) :: fft_cf real(8), pointer :: fft_coulb_FS(:,:) real(8), parameter :: pi = 3.141592653589793_8 contains !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE POISSON_SOLVE ! ======================== ! ! INPUT: ! rho [real(8), dimension(n,n)] : input density. ! -------- ! OUTPUT: ! v [real(8), dimension(n,n)] : potential that corresponds to input density rho. ! ! Notice that this routine just calls the real workers, poisson_sum, or poisson_fft, ! depending on how "mode" is set in the beginning of the module. !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine poisson_solve(rho, v) use mesh implicit none real(8), intent(in) :: rho(N, N) ! the density real(8), intent(out) :: v(N, N) ! the solution of the poisson equation select case(mode) case(HARTREE_SUM) call poisson_sum(rho, v) case(HARTREE_FFT) call poisson_fft(rho, v) end select end subroutine poisson_solve !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE POISSON_SUM ! ====================== ! ! INPUT: ! rho [real(8), dimension(n,n)] : input density. ! -------- ! OUTPUT: ! v [real(8), dimension(n,n)] : potential that corresponds to input density rho. ! ! This routine calculates the Hartree potential v created by the charge distribution ! rho, using the direct-sum method. !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine poisson_sum(rho, v) use mesh implicit none real(8), intent(in) :: rho(N, N) ! the density real(8), intent(out) :: v(N, N) ! the solution of the poisson equation !!!!!! MISSING CODE 9 !!!!!! END OF MISSING CODE end subroutine poisson_sum !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE POISSON_FFT ! ====================== ! ! INPUT: ! rho [real(8), dimension(n,n)] : input density. ! -------- ! OUTPUT: ! v [real(8), dimension(n,n)] : potential that corresponds to input density rho. ! ! This routine calculates the Hartree potential v created by the charge distribution ! rho, using the FFT based plane-wave method. ! ! You may find some strange calls to "dcf_alloc_RS" routines in some "cf" module. ! This is just because I essentially borrowed this from the octopus code, and ! this is the way it is done there. One does not need so much sophistication ! for this QD code, but it was faster to code it this way. ! ! In fact, this routine just does the Fourier transforms and calculates ! the product (second Eq. 27) in the notes. The "real" work of calculating the Fourier ! transform of the Coulomb or Yukawa interaction (first Eq. 27, or Eq. 29) ! is done by poisson_init. !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine poisson_fft(rho, pot) implicit none real(8), intent(in) :: rho(n, n) real(8), intent(out) :: pot(n, n) integer :: k, ix, iy call dcf_alloc_RS(fft_cf) ! allocate the cube in real space call dcf_alloc_FS(fft_cf) ! allocate the cube in Fourier space fft_cf%rs = 0.0_8 do ix = 1, n do iy = 1, n fft_cf%rs(fft_cf%n(1)/2+1-n/2+ix,fft_cf%n(2)/2+1-n/2+iy) = rho(ix, iy) enddo enddo call dcf_RS2FS(fft_cf) ! Fourier transform ! multiply by the FS of the Coulomb interaction fft_cf%FS(1:fft_cf%nx,1:fft_cf%n(2)) = fft_cf%FS(1:fft_cf%nx,1:fft_cf%n(2))*& fft_Coulb_FS(1:fft_cf%nx,1:fft_cf%n(2)) call dcf_FS2RS(fft_cf) ! Fourier transform back do ix = 1, n do iy = 1, n pot(ix, iy) = fft_cf%rs(fft_cf%n(1)/2+1-n/2+ix,fft_cf%n(2)/2+1-n/2+iy) enddo enddo call dcf_free_RS(fft_cf) ! memory is no longer needed call dcf_free_FS(fft_cf) end subroutine poisson_fft !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE POISSON_INIT ! ======================= ! ! No arguments. On output, the module variable fft_cf should be inialized, ! and ready to be used by poisson_fft. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine poisson_init() implicit none integer :: ix, iy, iz, ixx(2), db(2) real(8) :: temp(2), vec real(8) :: gpar,gperp,gx,gz,r_c real(8) :: DELTA_R = 1.0e-12_8 ! double the box to perform the fourier transforms db(1) = nint((1.0_8+sqrt(2.0_8))*n) db(2) = nint((1.0_8+sqrt(2.0_8))*n) call dcf_new(db, fft_cf) ! allocate cube function where we will perform call dcf_fft_init(fft_cf) ! the ffts r_c = sqrt(2.0_8)*n*delta ! store the fourier transform of the Coulomb interaction allocate(fft_Coulb_FS(fft_cf%nx, fft_cf%n(2))) fft_Coulb_FS = 0.0_8 ! Warning: this has to be set.... temp(1) = 2.0_8*pi/(db(1)*delta) temp(2) = 2.0_8*pi/(db(2)*delta) do iy = 1, db(2) ixx(2) = pad_feq(iy, db(2), .true.) do ix = 1, fft_cf%nx ixx(1) = pad_feq(ix, db(1), .true.) vec = sqrt(sum((temp(:)*ixx(:))**2)) select case(interaction) case(COULOMB) fft_coulb_fs(ix, iy) = r_c*besselint(vec*r_c) case(YUKAWA) fft_coulb_fs(ix, iy) = yukawa_fs(vec) end select end do end do end subroutine poisson_init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Next procedures are interal. They calculate the Fourier transform of the ! electronic interaction. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real(8) function yukawa_fs(x) result(y) real(8), intent(in) :: x y = 2.0_8*pi/(gamma*sqrt(1.0_8+x**2/gamma**2)) end function yukawa_fs real(8) function besselint(x) result(y) real(8), intent(in) :: x real(8), external :: bessel integer :: k real(8) :: z if(x < 0.2_8) then y = 2*pi - (pi/6.0_8)*x**2 return endif y = 0.0_8 k = 1 do z = bessel(k, x)/x y = y + z if(abs(z)<1.0e-9_8) exit k = k + 2 enddo y = 4*pi*y end function end module poisson