[Physics] Implementing simple atom model using density functional theory (DFT)

atomic-physicscomputational physicsdensity-functional-theoryquantum mechanicssoftware

I am trying to write computer code which will find the energy and density function for an atom with $Z$ protons and $N$ electrons. I am working in 1D for simplicity and would like to make the overall code as simple as possible (I am using Hartree units aswell). I am using also using the Hartree approximation. This is my overall algorithm:

Finds hydrogen like wave-functions (i.e. no electron-election interaction)
Works out density function via
$n(z) = \sum_i \left | \phi_i(z) \right |^2$
Works out Hartree potential
$ V_H(r) = \int_0^\infty \frac{n(\vec{s})}{\left | \vec{r}-\vec{s} \right |} d\vec{s} $
Uses Kohn-Sham equation to find new individual wavefuncitons
$\left (-\tfrac{1}{2}\bigtriangledown^2 + \frac{Z}{\left |r \right |} + V_H(\vec{r}) + V_{EX}(\vec{r}) \right )\phi_i = \varepsilon_i \phi_i$
Go back to step 1. and repeat.

I am having a number of difficulty and would be very grateful for some help with any of my questions.

Firstly to compute the Hartree potential I need to turn it into a 1D problem, so we assume even spherical distribution of charge:

$ V_H(r) = \int_0^\infty 4\pi s^2 \frac{ n(s)}{\left |r-s\right |} ds $

However this blows up when $r = s$. I have been looking into Gauss's law but had no luck.

Additionally, I am unsure how to compute the exchange-correlation term. I am trying to find a function which I can easily implement.

Any help would be great.

Best Answer

The Hartree potential is typically calculated not through the integral you're giving, but by solving Poisson's equation. Look here or there for details about this calculation applied to DFT, and here for what I found to be a good lower-level introduction to Poisson solvers in general. Good luck for your quest, it's really a nontrivial thing to do!

As an aside: what basis set are you going to employ?

Related Solutions

[Physics] Density Functional Theory DFT-D, DFT-D2 and DFT-D3

Empirical dispersion accounts for the fact that Density Functional Theory systematically underestimates Van der Waals interactions. DFT-D and DFT-D2 are very simple. DFT-D3 is a bit more complex.

DFT-D and DFT-D2 both consist of a pairwise -1/r^6 attraction, like a Lennard-Jones attraction, scaled by a parameter C6 and a damping function which prevents the energy from going to -∞ at r=0. The distance parameter r0 is just the sum of the covalent radii of the interacting pair of atoms. The energy parameter C6 is the geometric mean of the two element-wise C6 parameters. Importantly, neither C6 or r0 depend on anything but the pair of elements involved. That makes DFT-D and DFT-D2 very easy to implement.

The following code implements the DFT-D2 energy as used in the wB97X-D functional ("Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections", Chai & Head-Gordon, 2008). DFT-D and variants use other choices of damping function but are otherwise the same. Parameters C6 and covalent radii are available in "Semiempirical GGA-type density functional constructed with a long-range dispersion correction", Grimme, 2006. Currently both of these papers are available as PDFs on Google Scholar. You can get the parameters c6_coefficients and covalent_radii there.

def dftd2_wB97XD(r, ai, aj):
    r6 = r**6
    c6 = sqrt(c6_coefficients[ai] * c6_coefficients[aj])
    r0 = covalent_radii[ai] + covalent_radii[aj]
    fdamp = 1.0 / (1.0 + 6.0 * (r / r0)**-12.0)
    return -fdamp * c6 / r6

Be careful about the units of the parameters - J/mol nm^6 for C6, and Å for the covalent radii. In a real use case, you probably want C6 in kcal/mol Å^6 or Hartree Å^6.

DFT-D3 is conceptually only a little more complex: -C6/r^6 * fdamp6 + C8/r^8 * fdamp8. However, the parameters C6 and C8 depend on the coordination number of each atom, estimated using the distances to other nearby atoms. This makes DFT-D3 more accurate than the earlier versions. A good description can be found at https://www.vasp.at/wiki/index.php/DFT-D3.

There is also an even newer variant, DFT-D4, which continues to add accuracy and complexity by first estimating the atomic partial charges, and then using them to estimate the dispersion parameters. In this way, each atom's dispersion parameters depend on the entire surrounding environment. The paper for the D4 correction is here: https://chemrxiv.org/articles/A_Generally_Applicable_Atomic-Charge_Dependent_London_Dispersion_Correction_Scheme/7430216

Related Question