Developers Manual:Ion-Ion interaction

From OctopusWiki
Jump to: navigation, search

The following description works, but it is slow due to the interpolation required and less modular. So now this term is implemented using an standard Ewald sum; see here for the formalism.

Energy

The ion-ion interaction term is tricky to treat, as it is an infinite sum with a finite result. The standard way to handle it is through an Ewald sum; however in Octopus this is done in a slightly different way (formally equivalent).

The energy term (per volume unit) is:


E_{ii}=\frac12\sum_i^{\infty}\sum_j^{in}\frac{z_iz_j}{r_{ij}}

in is a sum over atoms inside the cell, out only over atoms outside and \infty is over all the atoms.

The first separation that is


E_{ii}=\frac12\sum_i^{in}\sum_j^{in}\frac{z_iz_j}{r_{ij}} + \frac12\sum_i^{out}\sum_j^{in}\frac{z_iz_j}{r_{ij}}

the first term is the energy of a finite system, the second can be seen as the interaction of the particles in the cell with the potential generated by the particles outside, we will call it E_{io}.

Now, as in the case of the Ewald sum we will separate the potential in two parts:


\frac{z_j}{r}=z_j\frac{1-\mbox{erf}(\alpha r)}{r} + z_j\frac{\mbox{erf}(\alpha r)}{r}=v^{sr}_j(r)+v^{lr}_j(r)\ ,

with \alpha chosen to match the separation of the pseudo-potential. Now


E_{io}=\frac12\sum_i^{out}\sum_j^{in}z_jv^{sr}_i(r_j) + \frac12\sum_j^{in}z_j\sum_i^{out}v^{lr}_i(r_j)

As v^{sr}_j(r) is short-ranged (i.e. exists r_c such that v^{sr}_j(r)=0 for r>r_c) the first term is a finite sum and can be calculated directly by including periodic copies of the atoms in the cell up to an r_c radius.

To obtain the second term we first see that


V(r)=\sum^{\infty}_iv^{lr}_i(r)

can be obtained by solving the Poisson equation


\nabla^2V(r)=-4\pi\sum_i^{\infty}n^{lr}(r)=-4\pi n(r)

and n(r) can be calculated directly because n^{lr}, the charge distribution that generates v^{lr}, is short ranged (it is a gaussian). Actually this is the same way that the long range part of the pseudopotential is calculated, so we already have to calculate V.

From here is direct that


V^{out}(r)=\sum_i^{out}v^{lr}_i(r) = V(r) - \sum_i^{in}v^{lr}_i(r)\ .

So the second term of E_{io} can be calculated by evaluating V^{out}\! in the atomic positions.


E^{lr}=\frac12\sum_j^{in}z_jV^{out}(r_j)

In the case of Octopus this means we have to interpolate the function, as the atoms are not always over a grid point. Currently this is quite slow, if we want to do dynamics with periodic systems this has to be improved.

Forces

The forces are the gradients of the previous terms with respect to the positions of the atoms in the cell.

The in-in term is as usual for finite systems:



\vec{F}_i=\sum_j^{in}\frac{z_iz_j}{r_{ij}^3}(\vec{r}_i-\vec{r}_j)\,,

the short range part of the in-out term is similar:



\vec{F}^{sr}_i=\frac12\sum_j^{out}z_i\frac{d v^{sr}(r_j)}{dr}\frac{(\vec{r}_i-\vec{r}_j)}{r_{ij}}\ ,

while the long range part is simply:


\vec{F}^{lr}_i=\frac12z_i\vec{\nabla}V^{out}(r_i)
.