Skip to content

Commit

Permalink
add README
Browse files Browse the repository at this point in the history
  • Loading branch information
@ArrogantGao
ArrogantGao committed Jan 26, 2024
1 parent b3d52c0 commit d98226e
Show file tree
Hide file tree
Showing 4 changed files with 98 additions and 21 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,12 @@ LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

[compat]
julia = "1"
julia = "1.7"
FFTW = "1.8"
FastTransforms = "0.15.16"
LoopVectorization = "0.12.166"
SpecialFunctions = "2.3.1"
LinearAlgebra = "1.10"
LinearAlgebra = "1.7"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Expand Down
87 changes: 86 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# ChebParticleMesh
# ChebParticleMesh.jl

<!-- [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://HPMolSim.github.io/ChebParticleMesh.jl/stable/) -->
<!-- [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://HPMolSim.github.io/ChebParticleMesh.jl/dev/) -->
Expand All @@ -7,3 +7,88 @@
[![Coverage](https://codecov.io/gh/HPMolSim/ChebParticleMesh.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/HPMolSim/ChebParticleMesh.jl)


`ChebParticleMesh.jl` is a package based on `Julia`, which provide a set of tools for the widely used **Particle-Mesh** methods.
Using the tools provided by this package, you can easily create a uniform grid with arbitrary dims and periodicity, and then interpolate the particles onto the grid with provided/self-defined window functions so that the problem can be transformed into reciprocal space via **Fast Fourier Transform** and then being solved there.

## Getting Started

Add this package in julia by typing `]` in Julia REPL and then type
```julia
pkg> add ChebParticleMesh
```
to install the package.

## Introduction

In this package, the process of Particle-Mesh is divided into SIX different steps, including

1. Pre-computating: Generating the grid, the window function and scaling factors
2. Interpolatation: Interpolating the point sources onto the grids
3. FFT: Using FFT to transform the density on grids into the reciprocal space
4. Scaling: Scale the fourier modes with the Green's function
5. IFFT: Using IFFT to transform the grid back to the real space
6. Gathering: Integrate the potential on the grid back to the particles

Specially, this package is called **Cheb**ParticleMesh.jl because we provided tools to approximate the window functions via piecewise Chebyshev poly, so that the calculation can be much faster when using some complicated kernels, such as the KB-kernel.

## Examples of Use

### Implentation of PPPM

In the following example, we will show how to use the `ChebParticleMesh.jl` to implentation of the long-range part of the well known [P3M method](https://en.wikipedia.org/wiki/P3M).

The target is to compute the double summation
$$
\sum_{\mathbf{m} \in \mathbb{Z}^3}\sum_{i = 1}^N \sum_{j = 1}^N q_i q_j \frac{\mathrm{erf}(\alpha |\mathbf{r}_{ij} + \mathbf{L_m}|)}{|\mathbf{r}_{ij}|}
$$
where $\mathrm{erf}(\cdot)$ for the error function and $\mathbf{L_m} = (m_x L_x, m_y L_y, m_z L_z)$ for the periodic images.

First we set up the system:
```Julia
L = 10.0, 10.0, 10.0 # size of the box
periodicity = (true, true, true) # fully periodic
extra_pad = (0, 0, 0) # no extra padding needed
N_real = (16, 16, 16) # 16 grids in each direction in real space
w = N_real ./ 4 # cutoff of interpolating
α = 0.5
```

Precomputation:
```julia
# Grid generation
gridinfo = GridInfo(N_real, w, periodicity, extra_pad, L)
gridbox = GridBox(gridinfo)

# using the Wkb window function, and generate its Chebyshev approximation
f_window = [x -> Wkb(x, (w[i] + 0.5) * gridinfo.h[i], 5.0 * w[i]) for i in 1:3]
cheb_coefs = tuple([ChebCoef(f_window[i], gridinfo.h[i], w[i], 10) for i in 1:3]...)

# using the Fourier transform of the window function and Green's function in Fourier space to generate the scaling factors
F_f_window = [x -> FWkb(x, (w[i] + 0.5) * gridinfo.h[i], 5.0 * w[i]) for i in 1:3]
func_scale = (kx, ky, kz) -> iszero(kx^2 + ky^2 + kz^2) ? zero(T) : (F_f_window[1](kx) * F_f_window[2](ky) * F_f_window[3](kz))^(-2) * exp(-(kx^2 + ky^2 + kz^2) / (4*α^2)) / (kx^2 + ky^2 + kz^2)
scalefactor = ScalingFactor(func_scale, gridinfo)
```

Generate some randomly allocated particles as input
```julia
n_atoms = 200
qs = 2.0 .* rand(n_atoms) .- 1.0
qs .-= sum(qs) / n_atoms
poses = [L .* (rand(), rand(), rand()) for i in 1:n_atoms]
```

Compute:
```julia
interpolate!(qs, poses, gridinfo, gridbox, cheb_coefs)
fft!(gridbox)
scale!(gridbox, scalefactor)
ifft!(gridbox)
E = gather(qs, poses, gridinfo, gridbox, cheb_coefs) / 8π
```

The full script is given in `example/PPPM.jl`, the result is compared with a mesh-free method, `Ewald3D` provided by `EwaldSummations.jl` and shows great convergency as the number of grid points increase.
![convergency](./example/PPPM.png)

## How to Contribute

If you find any bug or have any suggestion, please open an issue.
28 changes: 10 additions & 18 deletions example/PPPM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ qs = 2.0 .* rand(n_atoms) .- 1.0
qs .-= sum(qs) / n_atoms
L = 10.0, 10.0, 10.0
poses = [L .* (rand(), rand(), rand()) for i in 1:n_atoms]
α = 0.5
α = 0.1

function PPPM(qs::Vector{T}, poses::Vector{NTuple{3, T}}, N_real::NTuple{3, Int}, w::NTuple{3, Int}, L::NTuple{3, T}, α::T) where{T}
periodicity = (true, true, true)
Expand All @@ -15,10 +15,10 @@ function PPPM(qs::Vector{T}, poses::Vector{NTuple{3, T}}, N_real::NTuple{3, Int}
gridinfo = GridInfo(N_real, w, periodicity, extra_pad, L)
gridbox = GridBox(gridinfo)

f_window = [x -> Wkb(x, 6.5 * gridinfo.h[i], 30.0) for i in 1:3]
cheb_coefs = tuple([ChebCoef(f_window[i], gridinfo.h[i], w[i], 10) for i in 1:3]...)
f_window = [x -> Wkb(x, (w[i] + 0.5) * gridinfo.h[i], 8.0 * w[i]) for i in 1:3]
cheb_coefs = tuple([ChebCoef(f_window[i], gridinfo.h[i], w[i], 20) for i in 1:3]...)

F_f_window = [x -> FWkb(x, 6.5 * gridinfo.h[i], 30.0) for i in 1:3]
F_f_window = [x -> FWkb(x, (w[i] + 0.5) * gridinfo.h[i], 8.0 * w[i]) for i in 1:3]

func_scale = (kx, ky, kz) -> iszero(kx^2 + ky^2 + kz^2) ? zero(T) : (F_f_window[1](kx) * F_f_window[2](ky) * F_f_window[3](kz))^(-2) * exp(-(kx^2 + ky^2 + kz^2) / (4*α^2)) / (kx^2 + ky^2 + kz^2)

Expand All @@ -36,26 +36,18 @@ end

E_PPPM = []

for Nr in [8, 16, 32, 64, 128]
for Nr in [8, 16, 32, 64, 128, 256]
N_real = tuple([Nr for i in 1:3]...)
w = Int.(N_real ./ 4)
Ei = PPPM(qs, poses, N_real, w, L, α)
@show Ei
push!(E_PPPM, Ei)
end

E_Ewald = [];

Ewald_pos = Point{3, Float64}[Point{3, Float64}(pos) for pos in poses];
for s in 1.0:10.0
interaction = Ewald3DInteraction(n_atoms, s, α, L, ϵ = 1.0, ϵ_inf = Inf)
E_Ewald_i = Ewald3D_long_energy(interaction, Ewald_pos, qs)
push!(E_Ewald, E_Ewald_i)
end

@show E_PPPM
@show E_Ewald
s = 10.0
interaction = Ewald3DInteraction(n_atoms, s, α, L, ϵ = 1.0, ϵ_inf = Inf)
E_Ewald = Ewald3D_long_energy(interaction, Ewald_pos, qs)

plot(ylabel = "long range energy")
plot!(E_PPPM, label = "PPPM", marker = :circle)
plot!(E_Ewald, label = "Ewald", marker = :square)
plot([3, 4, 5, 6, 7, 8], log10.(abs.(E_PPPM .- E_Ewald)), xlabel = "log2(N real)", ylabel = "log10(error)", marker = :circle, label = "PPPM", dpi = 500)
savefig("./PPPM.png")
Binary file added example/PPPM.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

2 comments on commit d98226e

@ArrogantGao
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/99601

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.0 -m "<description of version>" d98226e75b4367d13efbf7f79c06e24b434341bc
git push origin v0.1.0

Please sign in to comment.