Accuracy Testing

There are several ways on how to test accuracy of the calculations, which is nessary in order to know how good the results are.

Accuracy of Poisson equation Greens function

To test accuracy of Poisson equation and Helmholtz equation use test_accuracy function

HKQM.test_accuracy — Function

test_accuracy(a::Real, ne::Int, ng::Int, nt::Int; kwords) -> Dict

Test accuracy on Gaussian charge distribution self energy.

Arguments

a::Real : Simulation box size
ne::Int : Number of elements per dimension
ng::Int : Number of Gausspoints per dimension for r
nt::Int : Number of Gausspoints per dimension for t

Keywords

tmax=300 : Maximum t-value for integration
tmin=0 : Minimum t-value for integration
mode=:combination : Integration type - options :normal, :log, :loglocal, :local, :combination and :optimal
δ=0.25 : Localization parameter for local integration types
correction=true : Correction to tmax-> ∞ integration - true calculate correction, false do not calculate
tboundary=20 : Parameter for :combination mode. Switch to :loglocal for t>tboundary else use :log
α1=1 : 1st Gaussian is exp(-α1*r^2)
α2=1 : 2nd Gaussian is exp(-α2*r^2)
d=0 : Distance between two Gaussian centers

source

It can be called with

test_accuracy(5u"Å", 4, 24, 96)

Dict{String, Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(Eₕ,), 𝐋^2 𝐌 𝐓^-2, nothing}}} with 6 entries:
  "calculated integral"    => 24.7491 Eₕ
  "correction"             => 1.42426e-5 Eₕ
  "reference tail"         => 1.26221e-5 Eₕ
  "calculated"             => 24.7491 Eₕ
  "total reference energy" => 24.7394 Eₕ
  "reference integral"     => 24.7394 Eₕ

Which gives accuracy for 5 Å box with 4 elements per dimension and 24^3 Gauss-Lagrange points per element.

The default mode uses optimal_coulomb_tranformation to calculate Poisson equation

optimal_coulomb_tranformation

This is usually good enough.

Integration tuning

Work in progress

This is work in progress. More coming later...

To perform integration there is a so called t variable that has to be integrated from zero to infinity. The main contribution comes from small values.

Normal mode will use normal Gauss-Lagrange itegration
Logarithmic will distribute the points in logarithmic fashion
Local mean that average value around the points is calculated. This can be used with normal and logarithmic spacing

Fanally there is correction for very large values.

To test the integral accuracy calculate without correction and choose how many points are used for t-integrarion nt.

In example to use logarithmic spacing to integrate from t=20 to t=70 one can use

test_accuracy(5u"Å", 4, 24, 24; mode=:log, tmin=20, tmax=70, correction=false)

Dict{String, Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(Eₕ,), 𝐋^2 𝐌 𝐓^-2, nothing}}} with 6 entries:
  "calculated integral"    => 0.121273 Eₕ
  "correction"             => 0.00142426 Eₕ
  "reference tail"         => 0.00126212 Eₕ
  "calculated"             => 0.121273 Eₕ
  "total reference energy" => 24.7394 Eₕ
  "reference integral"     => 0.0141855 Eₕ

From the output we can see that integral is heavily overestimated, see "integration error" from output.

To test the same aree with loglocal mode gives

test_accuracy(5u"Å", 4, 24, 24; mode=:loglocal, tmin=20, tmax=70, correction=false)

Dict{String, Quantity{Float64, 𝐋^2 𝐌 𝐓^-2, Unitful.FreeUnits{(Eₕ,), 𝐋^2 𝐌 𝐓^-2, nothing}}} with 6 entries:
  "calculated integral"    => 0.0147545 Eₕ
  "correction"             => 0.00142426 Eₕ
  "reference tail"         => 0.00126212 Eₕ
  "calculated"             => 0.0147545 Eₕ
  "total reference energy" => 24.7394 Eₕ
  "reference integral"     => 0.0141855 Eₕ

and we can see that the accuracy was considerably improved.

To build integral tensor from parts you can use AbstractHelmholtzTensor types

HelmholtzTensorLinear
HelmholtzTensorLog
HelmholtzTensorLocalLinear
HelmholtzTensorLocalLog
HelmholtzTensorCombination

To overload the default tensor for calculations you need to create (=redefine) optimal_coulomb_tranformation function that returns HelmholtzTensor.