# Dielectric constant (via finite-field) tutorial

## Requirements

### Software components

- nanotools
- RESCU+

### Pseudopotentials

I will need the following pseudopotential file:

- C_LDA_TM_DZP.mat

Let's copy it from the pseudo database to the current working directory and `export RESCUPLUS_PSEUDO=$PWD`.

## Briefing

In this tutorial, I show how to calculate the dielectric constant of diamond looking at the screening of a finite external potential.

## Code

An xyz file named `diamond.fxyz`

```shell
8
# diamond
C   0.0000   0.0000   0.0000
C   0.0000   0.5000   0.5000
C   0.5000   0.0000   0.5000
C   0.5000   0.5000   0.0000
C   0.2500   0.2500   0.2500
C   0.2500   0.7500   0.7500
C   0.7500   0.2500   0.7500
C   0.7500   0.7500   0.2500
```

and the following Python script

```python
from nanotools import Atoms, Cell, System, TotalEnergy
from nanotools.dielectric import DielectricConstant
a = 3.562
cell = Cell(avec=[[a,0.,0.],[0.,a,0.],[0.,0.,a]], resolution=0.15)
atoms = Atoms(fractional_positions="diamond.fxyz")
sys = System(cell=cell, atoms=atoms)
sys.supercell([4,1,1])
ecalc = TotalEnergy(sys)
ecalc.solver.set_mpi_command("mpiexec -n 16")
calc = DielectricConstant(ecalc, axis=0, amplitude=0.1)
calc.solve()
die = calc.get_dielectric_constant()
print(f"epsilon = {die}")
fig = calc.plot_external_potential()
```

## Explanations

In this section, I show how to obtain the dielectric constant of diamond.
Here is a high level view of the calculation workflow:

1. Create a [`TotalEnergy` calculator](#total-energy-calculator) and specify calculation parameters like resolution, k-sampling, tolerance, etc.
2. Create a [`DielectricConstant` calculator](#dielectric-constant-calculator).
3. Solve the Kohn-Sham equation for a system in vacuum and a system in a finite external potential and compute the dielectric constant.

I will introduce a sawtooth potential which preserve periodic boundary conditions.
The electrons will react by screening the external potential.
The difference between the vacuum system and biased system will be a sawtooth potential scaled by the dielectric constant.
In other words, I will compute

$$\epsilon = {\partial V^{saw}(\mathbf{r}) \over \partial\mathbf{r}} / {\partial \Delta V^{eff}(\mathbf{r}) \over \partial\mathbf{r}}$$

where the derivatives are taken at the center of the linear sections of the sawtooth potential.

### Total energy calculator

I first define a diamond crystal system with a lattice constant of 3.562 Ang.
I create a $4\times 1\times 1$ supercell because I want the sawtooth potential to be linear over a relatively extended region.
It is good practice to test various supercell lengths.
I then create a total energy calculator (`ecalc`) from this `System` object.
For more detail on how to setup a total energy calculation, refer to [total energy tutorial](../etot/tutorial_etot.md).

### Dielectric constant calculator

I then initialize a DC calculator from the total energy calculator.
I also specify the parameters `axis` and `amplitude` which control along which axis the sawtooth potential varies and what is the amplitude of the sawtooth potential.
The amplitude should not be too large, otherwise the linear response approximation will break down, and not too small, otherwise numerical inaccuracies will creep in.
The left and right fit should give the same answer, as they do here.
Asymmetry may arise if the sawtooth potential does not break on atomic plane.
The vacuum and perturbed systems are solved invoking `solve`.
The calculator will save results to `nano_die_out.json` as the solving progresses.
The results can then be loaded with the class method `read`.

The value of the dielectric constant can be calculated calling

```python
die = calc.get_dielectric_constant()
print(f"epsilon = {die}")
```

yielding `epsilon = 5.578101320105986`.

The linear response approximation is verified plotting $V^{saw}(\mathbf{r})$ and $\Delta V^{eff}(\mathbf{r})$.
Calling `plot_external_potential` you should see

![Potential screening in diamond](../figs/die-0-C.png)