Skip to content
Permalink
Browse files

Doc and cleanup

  • Loading branch information
antoine-levitt committed Aug 7, 2019
1 parent 12cd2b9 commit b7613fb93150853e8aa3e79307f78dcabb4bb9d9
Showing with 59 additions and 53 deletions.
  1. +25 −19 README.md
  2. +13 −12 codegen/xc_fallback/README.md
  3. +1 −0 deps/build.jl
  4. +17 −19 src/core/energy_ewald.jl
  5. +3 −3 src/utils/guess_gaussian_sad.jl
@@ -5,39 +5,45 @@
[![Build Status on Linux](https://travis-ci.org/mfherbst/DFTK.jl.svg?branch=master)](https://travis-ci.org/mfherbst/DFTK.jl)
[![Coverage Status](https://coveralls.io/repos/mfherbst/DFTK.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/mfherbst/DFTK.jl?branch=master)

DFTK, short for the **density-functional toolkit** is a library of Julia routines
for experimentation with plane-wave-based density-functional theory (DFT).
The main aim at the moment is to provide a platform to facilitate numerical
analysis of algorithms and techniques related to DFT, such as self-consistent
field (SCF) schemes or pseudopotentials. For this we want to leverage as much
of the existing developments in plane-wave DFT and the related ecosystems
of Julia python or C codes as possible.
DFTK, short for the **density-functional toolkit** is a library of
Julia routines for experimentation with plane-wave-based
density-functional theory (DFT), as implemented in much larger
production codes such as Abinit, Quantum Espresso and VASP. The main
aim at the moment is to provide a platform to facilitate numerical
analysis of algorithms and techniques related to DFT. For this we want
to leverage as much of the existing developments in plane-wave DFT and
the related ecosystems of Julia python or C codes as possible.

The library is at a very early stage of development and the supported feature set
is thus very limited. Current features include:
is thus limited. Current features include:
- Lattice construction and problem setup based on [pymatgen](https://pymatgen.org/)
- Plane-wave discretisations building on top of
[FFTW](https://github.com/JuliaMath/FFTW.jl) and [fftw](http://fftw.org/).
[FFTW](https://github.com/JuliaMath/FFTW.jl).
- SCF routine based on [NLsolve](https://github.com/JuliaNLSolvers/NLsolve.jl)
and [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl).
- LDA and GGA functionals from [Libxc](https://github.com/unkcpz/Libxc.jl)
and [libxc](https://tddft.org/programs/libxc).
- LDA and GGA functionals from [Libxc](https://github.com/unkcpz/Libxc.jl).
- Band structure generation
- Support for both `Float64` (double precision) and `Float32` (single precision)
throughout the library (Only for selected DFT functionals at the moment).

**Note:** The code has not been properly verified against a standard DFT package
and might likely contain bugs. Use for production calculations
is not yet recommended.
and might contain bugs. Use for production calculations is not yet recommended.

## Installation
The package is not yet registered, so you need to install it from the github url:
The package is not yet registered, so you need to install it from the github url: use
```
(v1.1) pkg> add https://github.com/mfherbst/DFTK.jl
] add https://github.com/mfherbst/DFTK.jl.
```
from a Julia command line (version at least 1.1). You can then run the code in the `examples/` directory.

Some parts of the code require a working Python installation with the libraries `scipy` and `spglib`. The examples require `matplotlib` and `pymatgen`. Check out which version of python is used by the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package, and use the correponding package manager (usually `apt`, `pip` or `conda`) to install these libraries.

## Perspective
Despite the current focus on numerics, the intention is too keep the project
rather general, however, such that this platform could be useful for other
research in the context of plane-wave discretisations or DFT in the future.
Any contribution or discussion is welcome!
Despite the current focus on numerics, the intention is to keep the
project rather general, so that this platform could be useful for
other research in the future.

## Contact
Feel free to contact us (Michael Herbst and Antoine Levitt) directly,
to open issues and to submit pull requests. Any contribution or
discussion is welcome!
@@ -1,13 +1,14 @@
# Code generators for the functionals in src/core/xc_fallback
This directory provides routines that are used to generate the code in
`src/core/xc_fallback`. The intention is to be able to quickly add a
pure-Julia fallback implementation for simple functionals, such that
calculations in deviating precision or using non-standard Array types
can be performed. For the standard case DFTK still relies on libxc via Libxc.jl.

The generators in this folder are based upon the Maple code generators used
by [libxc](https://tddft.org/programs/libxc) in order to implement
its exchange-correlation functionals. The intention is to be able to quickly
add a pure-Julia fallback implementation for simple functionals, such that
calculations in deviating precision or using non-standard Array types can
be performed. For the standard case DFTK still relies on libxc via Libxc.jl.

Currently the generation procedure only generates rough code, which needs to
be manually edited in order to compile and to work. Thus a proper mass-generation
of functionals would require a more sophisticated solution, but the
implementation of selected important cases should be possible like so.
The generators in this folder are based upon the Maple code generators
used by [libxc](https://tddft.org/programs/libxc) in order to
implement its exchange-correlation functionals. Currently the
generation procedure only generates rough code, which needs to be
manually edited in order to compile and to work. Thus a proper
mass-generation of functionals would require a more sophisticated
solution, but the implementation of selected important cases should be
possible like so.
@@ -1,3 +1,4 @@
# This script is called automatically upon package installation.
using PyCall
if PyCall.conda
using Conda
@@ -3,11 +3,12 @@ using SpecialFunctions: erfc
"""
energy_ewald(lattice, [recip_lattice, ]charges, positions; η=nothing)
Compute electrostatic interaction energy per unit cell between point charges in a
uniform background of compensating charge to yield net neutrality. the `lattice`
and `recip_lattice` should contain the lattice and reciprocal lattice vectors as columns.
`charges` and `positions` are the point charges and their positions in fractional
coordinates.
Compute the electrostatic interaction energy per unit cell between point
charges in a uniform background of compensating charge to yield net
neutrality. the `lattice` and `recip_lattice` should contain the
lattice and reciprocal lattice vectors as columns. `charges` and
`positions` are the point charges and their positions (as an array of
arrays) in fractional coordinates.
"""
function energy_ewald(lattice, charges, positions; η=nothing)
T = eltype(lattice)
@@ -50,18 +51,16 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)

# Function to return the indices corresponding
# to a particular shell
# TODO switch to an O(N) implementation
function shell_indices(ish)
[[i,j,k] for i in -ish:ish for j in -ish:ish for k in -ish:ish
if maximum(abs.([i,j,k])) == ish]
end

# Loop over reciprocal-space shells
gsh = 0
gsh = 1 # Exclude G == 0
any_term_contributes = true
while any_term_contributes
# Notice that the first gsh this loop processes is 1
# in other words G == 0 is implicitly excluded.
gsh += 1
any_term_contributes = false

# Compute G vectors and moduli squared for this shell patch
@@ -82,6 +81,7 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
any_term_contributes = true
sum_recip += sum_strucfac * exp(-exponent) / Gsq
end
gsh += 1
end
# Amend sum_recip by proper scaling factors:
sum_recip = sum_recip * 4π / det(lattice)
@@ -93,24 +93,21 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
sum_real::T = -2η / sqrt(π) * sum(Z -> Z^2, charges)

# Loop over real-space shells
rsh = -1
rsh = 0 # Include R = 0
any_term_contributes = true
while any_term_contributes || rsh <= 0
# In this loop the first rsh, which is processed is rsh == 0
rsh += 1
while any_term_contributes || rsh <= 1
any_term_contributes = false

# Loop over R vectors for this shell patch
for R in shell_indices(rsh)
for (ti, Zi) in zip(positions, charges), (tj, Zj) in zip(positions, charges)
# Compute norm of the distance
dist = norm(lattice * (ti - tj - R))

# Avoid zero denominators
if dist <= eps(T)^(3/2)
# Avoid self-interaction
if rsh == 0 && ti == tj
continue
end

dist = norm(lattice * (ti - tj - R))

# erfc decays very quickly, so cut off at some point
if η * dist > max_erfc_arg
continue
@@ -120,6 +117,7 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing)
sum_real += Zi * Zj * erfc* dist) / dist
end # iat, jat
end # R
rsh += 1
end
(sum_recip + sum_real) / 2 # Amend by 1/2 (because of double counting)
(sum_recip + sum_real) / 2 # Divide by 1/2 (because of double counting)
end
@@ -20,21 +20,21 @@ function guess_gaussian_sad(basis, composition...)
end


## TODO give the formula in real space in the doc, and clarify that the atomic density is that of the valence electrons
@doc raw"""
Get the atomic decay length for an atom with `n_elec_core` core
and `n_elec_valence` valence electrons. The returned length parameter can be used
to generate an approximate atomic gaussian density in reciprocal space:
```math
\hat{ρ}(G) = Z \exp\left(-(2π \text{length} ⋅ G)^2)
\hat{ρ}(G) = Z \exp\left(-(2π \text{length} |G|)^2\right)
```
"""
function atom_decay_length(n_elec_core, n_elec_valence)
# Adapted from ABINIT/src/32_util/m_atomdata.F90,
# from which also the data has been taken.

# Round the number of valence electrons and return early if zero
n_elec_valence = round(Int, n_elec_valence)
if n_elec_valence == 0.0
if n_elec_valence == 0
return 0.0
end

0 comments on commit b7613fb

Please sign in to comment.
You can’t perform that action at this time.