Parallel Poisson solver

This library provides a parallel Poission solver implemented in Vmol. A notable feature is that the non-zero Dirichlet boundary condition can be used (solved using Ricker's boundary erasure).

Prerequisites

This library requires FFTW and PFFT.

Usage

First, use MPI_Cart_create to generate the 3-dimensional mapping of MPI ranks to the 3-D cartesian coordinate.

   use MPI
   use poission_solver
   ...

   call MPI_Cart_create(MPI_COMM_WORLD, 3, &      ! 3 dimensions
         nprocxyz_corder, &                       ! Desired # of processors at each dimension, in a reversed order. Typically (PZ, PY, PX). PX*PY*PZ must be equal to total ranks.
         [.false., .false., .false.], &           ! Periodic = .false.
         .true., &                                ! Allow reordering (.true. for better performance, .false. may very slightly slow down but MPI_COMM_WORLD's node order matches the mappping)
         mpi_comm_cart, ierr)                     ! output communicator
	call MPI_Comm_rank(mpi_comm_cart, rank, ierr) ! When reorder = .true., this rank might not be equal to MPI_COMM_WORLD's.
    call MPI_Cart_coords(mpi_comm_cart, rank, 3, coords_corder, ierr)
    myposition = coords_corder(3:1:-1)            ! (px, py, pz) of my rank, 0 <= px < PX and such like.

Then, initialize the library by calling poisson_init.

    call poisson_init(ndim, mpi_comm_cart, &      ! ndim = total system size (NX, NY, NZ). (NX, NY, NZ) must be a multiple of (PX, PY, PZ).
         gridlen, &                               ! grid length (LX, LY, LZ)
         myposition, &                            ! my rank, only used to do sanity check
		 eigenvalue_type, &                       ! Use "eigenvalue_finite_difference_2" to use non-zero Dirichlet BC
         boundary_condition, &                    ! Use "poisson_bc_dirichlet_zero_staggered" to use non-zero Dirichlet BC
         verbose = .true., &                      ! Be loud and noisy
         measure = initialize_measure)            ! Use initial benchmark to choose faster routine (in PFFT and FFTW). Other choices are "initialize_estimate" (fast init, slow runtime) and "initialize_exhaustive" (slow init, fast runtime).

After the initialization is complete, the buffer phi_or_rho within poisson_solver library is allocated. This variable is a 3-dimensional matrix of real(8). Set this variable to the RHS of the Poisson equation, typically -4 * pi * rho (see also notes below). Then, if you want to use non-zero Dirichlet boundary condition, run the solver by

    call run_with_dirichlet_boundary(xmplane, xpplane, ymplane, ypplane, zmplane, zpplane) ! xmplane is 2-dimensional (y,z) plane at (xmin - 1/2, y, z), of the callee's rank.

Or, if you want the boundary condition with zeros (Dirichlet, Neumann etc.), just do call run().

After executing run_with_dirichlet_boundary or run, phi_or_rho is set to the electrostatic potential. If you want to recalculate with different charge densities, just set phi_or_rho and call run_with_dirichlet_boundary or run again.

If you want to release all the resources allocated by this library, call finalize.

See also an example in test_poisson_solver.f90.

Compile & link

The typical compilation commandline is:

mpif90 -O -c -I/path/to/FFTW/include -I/path/to/pfft/include -std=f2003 poission_solver.F90 -o poisson_solver.o

The typical linking order is:

mpif90 -o executable_file the_object_files_using_solver.o poisson_solver.o -L/path/to/pfft/lib -lpfft -L/path/to/fftw/lib -lfftw3_mpi -lfftw3

And please add as many optimization options you desire.

Notes

Size of phi_or_rho

The size of phi_or_rho might be larger than the size of (NX/PX, NY/PY, NZ/PZ). This is a feature of PFFT, required to achieve the better performance. Thus, for example,

    phi_or_rho(:, :, :) = -4 * pi * your_rho_variable(:, :, :)

may fail due to the size difference. You need to explicitly specify the size you want to save:

    CX = NX / PX
	CY = NY / PY
	CZ = NZ / PZ
    phi_or_rho(1:CX, 1:CY, 1:CZ) = -4 * pi * your_rho_variable(:, :, :)

Fortran 90 or 2003?

The file name is suffixed with .F90 but we use Fortran 2003. This is because Intel Fortran does not recognize .F03 nor .f03 as a suffix.

Asking to cite our works

If you use this library in your academic work, we humbly ask you to cite our work:

H. Takahashi, S. Sakuraba & A. Morita.
Large-Scale Parallel Implementation of Hartree-Fock Exchange Energy on the Real-Space Grids Using 3D-Parallel FFT.
Journal of Chemical Information and Modeling, in press. DOI: https://doi.org/10.1021/acs.jcim.9b01063 .

See also

• Kudos to PoisFFT, a library for solving the Poisson equation in parallel.
• BigDFT contains a Poisson solver library Psolver.

TODOs

• Move ugly module global variables into a single type.