poisson_test.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio,
!! G. Bertsch, M. Oliveira
!!
!! 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., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module poisson_test_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use io_function_oct_m
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use space_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                      &
    poisson_test
 
contains
 
  !-----------------------------------------------------------------
  subroutine poisson_test(this, space, mesh, latt, namespace, repetitions)
    type(poisson_t),         intent(in) :: this
    class(space_t),          intent(in) :: space
    class(mesh_t),           intent(in) :: mesh
    type(lattice_vectors_t), intent(in) :: latt
    type(namespace_t),       intent(in) :: namespace
    integer,                 intent(in) :: repetitions
 
    real(real64), allocatable :: rho(:), vh(:), vh_exact(:), xx(:, :), xx_per(:)
    real(real64) :: alpha, beta, rr, delta, ralpha, hartree_nrg_num, &
      hartree_nrg_analyt, lcl_hartree_nrg
    real(real64) :: total_charge
    integer :: ip, ierr, iunit, nn, n_gaussians, itime, icell
    real(real64) :: threshold
    type(lattice_iterator_t) :: latt_iter
 
    push_sub(poisson_test)
 
    if (mesh%box%dim == 1) then
      call messages_not_implemented('Poisson test for 1D case', namespace=namespace)
    end if
 
    !%Variable PoissonTestPeriodicThreshold
    !%Type float
    !%Default 1e-5
    !%Section Hamiltonian::Poisson
    !%Description
    !% This threshold determines the accuracy of the periodic copies of
    !% the Gaussian charge distribution that are taken into account when
    !% computing the analytical solution for periodic systems.
    !% Be aware that the default leads to good results for systems
    !% that are periodic in 1D - for 3D it is very costly because of the
    !% large number of copies needed.
    !%End
    call parse_variable(namespace, 'PoissonTestPeriodicThreshold', 1e-5_real64, threshold)
 
    ! Use two gaussians with different sign
    n_gaussians = 2
 
    safe_allocate(     rho(1:mesh%np))
    safe_allocate(      vh(1:mesh%np))
    safe_allocate(vh_exact(1:mesh%np))
    safe_allocate(xx(1:space%dim, 1:n_gaussians))
    safe_allocate(xx_per(1:space%dim))
 
    rho = m_zero
    vh = m_zero
    vh_exact = m_zero
 
    alpha = m_four*mesh%spacing(1)
    write(message(1),'(a,f14.6)')  "Info: The alpha value is ", alpha
    write(message(2),'(a)')        "      Higher values of alpha lead to more physical densities and more reliable results."
    call messages_info(2, namespace=namespace)
    beta = m_one / (alpha**space%dim * sqrt(m_pi)**space%dim)
 
    write(message(1), '(a)') 'Building the Gaussian distribution of charge...'
    call messages_info(1, namespace=namespace)
 
    ! Set the centers of the Gaussians by hand
    xx(1, 1) = m_one
    xx(2, 1) = -m_half
    if (space%dim == 3) xx(3, 1) = m_two
    xx(1, 2) = -m_two
    xx(2, 2) = m_zero
    if (space%dim == 3) xx(3, 2) = -m_one
    xx = xx * alpha
 
    ! Density as sum of Gaussians
    rho = m_zero
    do nn = 1, n_gaussians
      do ip = 1, mesh%np
        call mesh_r(mesh, ip, rr, origin = xx(:, nn))
        rho(ip) = rho(ip) + (-1)**nn * beta*exp(-(rr/alpha)**2)
      end do
    end do
 
    total_charge = dmf_integrate(mesh, rho)
 
    write(message(1), '(a,e23.16)') 'Total charge of the Gaussian distribution', total_charge
    call messages_info(1, namespace=namespace)
 
    write(message(1), '(a)') 'Computing exact potential.'
    call messages_info(1, namespace=namespace)
 
    ! This builds analytically its potential
    vh_exact = m_zero
    latt_iter = lattice_iterator_t(latt, m_one/threshold)
    do nn = 1, n_gaussians
      ! sum over all periodic copies for each Gaussian
      write(message(1), '(a,i2,a,i9,a)') 'Computing Gaussian ', nn, ' for ', latt_iter%n_cells, ' periodic copies.'
      call messages_info(1, namespace=namespace)
 
      do icell = 1, latt_iter%n_cells
        xx_per = xx(:, nn) + latt_iter%get(icell)
        !$omp parallel do private(rr, ralpha)
        do ip = 1, mesh%np
          call mesh_r(mesh, ip, rr, origin=xx_per)
          select case (space%dim)
          case (3)
            if (rr > r_small) then
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * loct_erf(rr/alpha)/rr
            else
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * (m_two/sqrt(m_pi))/alpha
            end if
          case (2)
            ralpha = rr**2/(m_two*alpha**2)
            if (ralpha < 100.0_real64) then
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * beta * (m_pi)**(m_three*m_half) * alpha * exp(-rr**2/(m_two*alpha**2)) * &
                loct_bessel_in(0, rr**2/(m_two*alpha**2))
            else
              vh_exact(ip) = vh_exact(ip) + (-1)**nn * beta * (m_pi)**(m_three*m_half) * alpha * &
                (m_one/sqrt(m_two*m_pi*ralpha))
            end if
          end select
        end do
      end do
    end do
 
    ! This calculates the numerical potential
    do itime = 1, repetitions
      call dpoisson_solve(this, namespace, vh, rho)
    end do
 
    ! Output results
    iunit = io_open("hartree_results", namespace, action='write')
    delta = dmf_nrm2(mesh, vh-vh_exact)
    write(iunit, '(a,e23.16)') 'Hartree test (abs.) = ', delta
    delta = delta/dmf_nrm2(mesh, vh_exact)
    write(iunit, '(a,e23.16)') 'Hartree test (rel.) = ', delta
 
    ! Calculate the numerical Hartree energy (serially)
    lcl_hartree_nrg = m_zero
    do ip = 1, mesh%np
      lcl_hartree_nrg = lcl_hartree_nrg + rho(ip) * vh(ip)
    end do
    lcl_hartree_nrg = lcl_hartree_nrg * mesh%spacing(1) * mesh%spacing(2) * mesh%spacing(3)/m_two
#ifdef HAVE_MPI
    call mpi_reduce(lcl_hartree_nrg, hartree_nrg_num, 1, &
      mpi_double_precision, mpi_sum, 0, mpi_comm_world, mpi_err)
    if (mpi_err /= 0) then
      write(message(1),'(a)') "MPI error in MPI_Reduce; subroutine poisson_test of file poisson.F90"
      call messages_warning(1, namespace=namespace)
    end if
#else
    hartree_nrg_num = lcl_hartree_nrg
#endif
 
    ! Calculate the analytical Hartree energy (serially, discrete - not exactly exact)
    lcl_hartree_nrg = m_zero
    do ip = 1, mesh%np
      lcl_hartree_nrg = lcl_hartree_nrg + rho(ip) * vh_exact(ip)
    end do
    lcl_hartree_nrg = lcl_hartree_nrg * mesh%spacing(1) * mesh%spacing(2) * mesh%spacing(3) / m_two
#ifdef HAVE_MPI
    call mpi_reduce(lcl_hartree_nrg, hartree_nrg_analyt, 1, &
      mpi_double_precision, mpi_sum, 0, mpi_comm_world, mpi_err)
    if (mpi_err /= 0) then
      write(message(1),'(a)') "MPI error in MPI_Reduce; subroutine poisson_test of file poisson.F90"
      call messages_warning(1, namespace=namespace)
    end if
#else
    hartree_nrg_analyt = lcl_hartree_nrg
#endif
 
    write(iunit, '(a)')
 
    if (mpi_world%rank == 0) then
      write(iunit,'(a,e23.16)') 'Hartree Energy (numerical)  =', hartree_nrg_num, 'Hartree Energy (analytical) =',hartree_nrg_analyt
    end if
 
    call io_close(iunit)
 
    call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_rho", namespace, space, &
      mesh, rho, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_exact", namespace, space, &
      mesh, vh_exact, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_numerical", namespace, space, &
      mesh, vh, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_rho", namespace, space, &
      mesh, rho, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_exact", namespace, space, &
      mesh, vh_exact, unit_one, ierr)
    call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_numerical", namespace, space, &
      mesh, vh, unit_one, ierr)
    ! not dimensionless, but no need for unit conversion for a test routine
 
    safe_deallocate_a(rho)
    safe_deallocate_a(vh)
    safe_deallocate_a(vh_exact)
    safe_deallocate_a(xx)
 
    pop_sub(poisson_test)
  end subroutine poisson_test
 
end module poisson_test_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: