3.44 ESP charges

This sections describes how to calculate partial charges by fitting to the electro static potential (ESP). These charges are widely used in the context of force fields ([222], [66], [280], [33], CHELP: [62], CHELPG: [43], RESP-charges: [24], CHELP-BOW/CHELMO: [279], ESP-charge from multi-pole-moments: [146]). FHI-aims implements a simple method for cluster calculation (molecules) as well as two methods for solids (periodic boundary conditions) [51], [61].
The starting point for these methods is the calculation of the electro static potential at a sufficiently high number of grid points outside the vdw radius of the atoms. To define a space region for the grid two parameters are necessary: a minimal radius and a maximal radius around the atoms. These radii are defined as multiples of the vdw-radius of the atoms, see figure 3.13 for details. The values for the vdw radii of most atoms in the periodic table have been taken from "http://de.wikipedia.org/wiki/Van-der-Waals-Radius" ([38], [266], [209]). For the generation of the points cubic grids are used. For cluster calculations points within a cube encapsulating the spheres with the maximal radius (multiple of the vdw radius) around all atoms are generated. For periodic boundary conditions the provided unitcell is used. The points within the superposition of the spheres with the minimal radius (minimal multiple of the vdw radius) are excluded. The atom-centered radial grids are also available, mainly for test purposes. Here, all points lie on N spheres with radii between the minimal and maximal multiple of the vdw radius. The spacing between the radii of the N spheres is equidistant or logarithmic.

Refer to caption — Figure 3.13: Definition of the volume used for the creation of grid points at which the potential is evaluated.

For the cluster case the function to fit to is a sum of Coulomb potentials with charges $q_{i}$ , the ESP-charges, at the atomic position $\mathbf{R}_{i}$ :

V_{ESP}(\mathbf{r})=\sum_{i=1}^{N_{at}}\frac{q_{i}}{|\mathbf{r}-\mathbf{R}_{i}|}

(3.203)

The $q_{i}$ are calculated by a least squares fit with the additional constraint of constant total charge $q_{tot}=\sum_{i=1}^{N_{at}}q_{i}$ . We use the method of Lagrange multipliers to minimize the function:

F=\sum_{k=1}^{N_{grid}}\left(V_{DFT}(\mathbf{r_{k}})-V_{ESP}(\mathbf{r_{k}})% \right)^{2}-\lambda\left(q_{tot}-\sum_{i=1}^{N_{at}}q_{i}\right)^{2}.

(3.204)

This can be translated into a system of linear equations:

\mathbf{\hat{A}}\mathbf{q}=\mathbf{B}.

(3.205)

with the $N_{at+1}$ x $N_{at+1}$ matrix $\mathbf{\hat{A}}$ :

	$\displaystyle A_{ij}=$	$\displaystyle\sum_{k=1}^{N_{grid}}\frac{1}{\|\mathbf{r_{k}}-\mathbf{R}_{i}\|}% \frac{1}{\|\mathbf{r_{k}}-\mathbf{R}_{j}\|}\quad\mathrm{with}\quad i,j\leq N_{at}$		(3.206)
	$\displaystyle A_{i=N_{at}+1,j}=$	$\displaystyle A_{i,j=N_{at}+1}=1;\qquad A_{i=N_{at}+1,j=N_{at}+1}=0$

and the $N_{at}+1$ vector $\mathbf{B}$ :

	$\displaystyle B_{i}=$	$\displaystyle\sum_{k=1}^{N_{grid}}\frac{V_{DFT}(\mathbf{r_{k}})}{\|\mathbf{r_{k% }}-\mathbf{R}_{i}\|}\quad\mathrm{with}\quad i\leq N_{at}$		(3.207)
	$\displaystyle B_{i=N_{at}+1}=q_{tot}$			(3.208)

$\mathbf{q}$ are the $N_{at}$ charges.
For solids (periodic boundary conditions) the situation is more complicated because all charges are repeated infinitely and the potential is only defined up-to an arbitrary offset. The methods to solve this problem are based on Ewald summation. Details on method 1 can be found here [51] and on method two here [61]. The function for the potential generated by the ESP charges centered at the atoms of the unit-cell now reads:

V_{ESP}(\mathbf{r})=\sum_{i=1,\mathbf{T}}^{N_{at}}q_{i}\frac{\mathrm{erfc}(% \alpha|\mathbf{r}-\mathbf{R_{i,\mathbf{T}}}|)}{|\mathbf{r}-\mathbf{R_{i,% \mathbf{T}}}|}+\frac{4\pi}{V_{cell}}\sum_{i=1,\mathbf{k}}^{N_{at}}q_{i}\mathrm% {cos}(\mathbf{k}(\mathbf{r}-\mathbf{R}_{i}))\frac{\mathrm{exp}^{-\frac{k^{2}}{% 4\alpha^{2}}}}{k^{2}}

(3.209)

with $\mathbf{T}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3}$ mapping the lattice positions. $\mathbf{a_{i}}$ the lattice vectors, $n_{i}\in\mathbb{Z}$ . $\mathbf{k}=m_{1}\mathbf{b}_{1}+m_{2}\mathbf{b}_{2}+m_{3}\mathbf{b}_{3}$ mapping the reciprocal space. $\mathbf{b_{i}}$ the reciprocal lattice vectors, $m_{i}\in\mathbb{Z}$ . $V_{cell}$ is the volume of the unit-cell. The parameter $\alpha$ is defined as $\alpha=\frac{\sqrt{\pi}}{R_{c}}$ with $R_{c}$ the cutoff radius of the Ewald summation. For the second method Wolf summation is implemented as well:

	$\displaystyle V_{ESP}(\mathbf{r})=$	$\displaystyle\sum_{i=1}^{N_{at}}q_{i}\left[\frac{\mathrm{erfc}(\sqrt{\alpha}\|% \mathbf{r}-\mathbf{R_{i}}\|)}{\|\mathbf{r}-\mathbf{R_{i}}\|}-\frac{\mathrm{erfc}(% \sqrt{\alpha}R_{c})}{R_{c}}+\left(\frac{\mathrm{erfc}(\sqrt{\alpha}R_{c})}{R_{% c}^{2}}+\frac{2\sqrt{\alpha}}{\sqrt{\pi}}\frac{\mathrm{exp}(\alpha R_{c}^{2})}% {R_{c}}\right)\right]$		(3.210)
		$\displaystyle\mathrm{x}(\|\mathbf{r}-\mathbf{R_{i}}\|-R_{c})$

The function to minimize for method 1 is:

	$\displaystyle F_{1}^{PBC}=$	$\displaystyle\sum_{k=1}^{N_{grid}}\left(V_{DFT}(\mathbf{r_{k}})-V_{ESP}(% \mathbf{r_{k}})+\frac{1}{N_{grid}}\sum_{j=1}^{N_{grid}}\left(V_{DFT}(\mathbf{r% _{j}})-V_{ESP}(\mathbf{r_{j}})\right)\right)^{2}$		(3.211)
		$\displaystyle-\lambda\left(q_{tot}-\sum_{i=1}^{N_{at}}q_{i}\right)+\sum_{i=1}^% {N_{at}}w_{i}\left(E_{i}^{0}+\chi_{i}q_{i}+\frac{1}{2}J_{i}^{00}q_{i}^{2}% \right).$

The $\chi_{i}$ and $J_{i}^{00}$ are electro-negativity and self-Coulomb interaction of the respective elements, which can be used to constrain the desired charges. $w_{i}$ are weighting factors for the constraints. This results in $\hat{A}$ and $\mathbf{B}$ for the linear equations system:

$\displaystyle A_{ij}=\sum_{k=1}^{N_{grid}}$	$\displaystyle\left(\frac{\partial V_{ESP}(\mathbf{r_{k}})}{\partial q_{i}}-% \frac{1}{N_{grid}}\sum_{j=1}^{N_{grid}}\frac{\partial V_{ESP}(\mathbf{r_{j}})}% {\partial q_{i}}\right)~{}\mathrm{x}$	(3.212)
	$\displaystyle\left(\frac{\partial V_{ESP}(\mathbf{r_{k}})}{\partial q_{j}}-% \frac{1}{N_{grid}}\sum_{m=1}^{N_{grid}}\frac{\partial V_{ESP}(\mathbf{r_{m}})}% {\partial q_{j}}\right)+\frac{w_{i}}{2}J_{i}^{00}\delta_{ij}\quad\mathrm{;}% \quad i,j\leq N_{at}$
$\displaystyle A_{i=N_{at}+1,j}=$	$\displaystyle A_{i,j=N_{at}+1}=1;\qquad A_{i=N_{at}+1,j=N_{at}+1}=0$
$\displaystyle B_{i}=\sum_{k=1}^{N_{grid}}$	$\displaystyle\left(V_{DFT}(\mathbf{r_{k}})-\frac{1}{N_{grid}}\sum_{j=1}^{N_{% grid}}V_{DFT}(\mathbf{r_{j}})\right)~{}\mathrm{x}$
	$\displaystyle\left(\frac{\partial V_{ESP}(\mathbf{r_{k}})}{\partial q_{i}}-% \frac{1}{N_{grid}}\sum_{m=1}^{N_{grid}}\frac{\partial V_{ESP}(\mathbf{r_{m}})}% {\partial q_{i}}\right)-w_{m}\frac{\chi_{m}}{2}\quad\mathrm{;}\quad i\leq N_{at}$
$\displaystyle B_{i=N_{at}+1}$	$\displaystyle=q_{tot}$

The function to minimize for method 2 is:

	$\displaystyle F_{2}^{PBC}=$	$\displaystyle\sum_{k=1}^{N_{grid}}\left(V_{DFT}(\mathbf{r_{k}})-\left(V_{ESP}(% \mathbf{r_{k}})+V_{DFT}^{offset}\right)\right)^{2}$		(3.213)
		$\displaystyle-\lambda\left(q_{tot}-\sum_{i=1}^{N_{at}}q_{i}\right)+\beta\sum_{% i=1}^{N_{at}}\left(q_{i}-q_{i0}\right)^{2}.$

The constraint charges $q_{i0}$ can be determined with other methods (e.g. Mulliken charge analysis [225]), $\beta$ is the weighing factor. This gives $\hat{A}$ and $\mathbf{B}$ :

	$\displaystyle A_{ij}=\sum_{k=1}^{N_{grid}}\left(\frac{\partial V_{ESP}(\mathbf% {r_{k}})}{\partial q_{i}}\frac{\partial V_{ESP}(\mathbf{r_{k}})}{\partial q_{j% }}\right)+\beta\delta_{ij}\quad\mathrm{;}\quad i,j\leq N_{at}$		(3.214)
	$\displaystyle A_{i=N_{at}+1,j}=A_{i,j=N_{at}+1}=\sum_{k=1}^{N_{grid}}\frac{% \partial V_{ESP}(\mathbf{r_{k}})}{\partial q_{j}}\quad j\leq N_{at}$
	$\displaystyle A_{i=N_{at}+2,j}=A_{i,j=N_{at}+2}=1$
	$\displaystyle A_{i=N_{at}+2,j=N_{at}+2}=A_{i=N_{at}+1,j=N_{at}+2}=A_{i=N_{at}+% 2,j=N_{at}+1}=0$
	$\displaystyle B_{i}=\sum_{k=1}^{N_{grid}}\left(V_{DFT}(\mathbf{r_{k}})\frac{% \partial V_{ESP}(\mathbf{r_{k}})}{\partial q_{i}}\right)-\beta q_{0i}\quad% \mathrm{;}\quad i\leq N_{at}$
	$\displaystyle B_{i=N_{at}+1}=\sum_{j=1}^{N_{grid}}V_{DFT}(\mathbf{r_{j}})$
	$\displaystyle B_{i=N_{at}+2}=q_{tot}$

Here the arbitrary offset of the potential $V_{offset}$ is an additional fitting parameter, the matrix $\hat{A}$ is of dimension $N_{at+2}$ x $N_{at+2}$ and $\mathbf{B}$ of dimension $N_{at+2}$ . For method 1 $V_{offset}$ can be calculated as:

V_{offset}=\sum_{k=1}^{N_{grid}}\left(V_{DFT}(\mathbf{r_{k}})-V_{ESP}(\mathbf{% r_{k}})\right)

(3.215)

from the fitted charges $q_{i}$ . As a measure for the quality of the fit the root-mean-square ( $R R M S$ ) is defined as:

RRMS=\left\{\frac{\sum_{k=1}^{N_{grid}}\left(\left(V_{ESP}(\mathbf{r_{k}})+V_{% offset}\right)-V_{DFT}(\mathbf{r_{k}})\right)^{2}}{\sum_{k=1}^{N_{grid}}\left(% V_{DFT}(\mathbf{r_{k}})\right)^{2}}\right\}^{2}

(3.216)

The current implementation is quite sensitive to the points chosen for the calculation of the electrostatic potential. Thorough studies regarding the parameters for the grid are advised! For periodic boundary conditions the ESP-charges calculated for transition densities are experimental, caution!! The ESP-charges from transition densities can be benchmarked against the dipole-moments calculated with compute_dipolematrix.
Example for a periodic system:

      output esp
      esp n_radius 10
      esp radius 1.0 2.
      esp pbc_method 1
      esp R_c 10
      output esp
      esp n_radius 10
      esp radius 1.0 2.
      esp pbc_method 2
      esp R_c 30
      output esp
      esp radius 1.0 2.
      esp pbc_method 2
      esp R_c 20
      esp equal_grid_n 10 10 10

The ESP-charges for the full potential will be calculated for a periodic system. At first with points within once the vdw radius and twice the vdw radius of the atoms, with 10 shells in between and a cutoff radius of 10 $\AA$ for the Ewald summation, method 1 (Ewald summation) is used. Secondly with points within once the vdw radius and twice the vdw radius of the atoms, with 10 shells in between and a cutoff radius of 30 $\AA$ , method 2 (Ewald summation) is used. Thirdly with points within once the vdw radius and twice the vdw radius of the atoms, with 10 $\times$ 10 $\times$ 10 initial points on a cubic grid, method 2 is used.

Warning: For periodic boundary conditions the atoms in the supplied geometry.in must be within the first (central) Wigner-Seitz-cell. For a two layer system it would normally be possible to have one layer within the (central) Wigner-Seitz-cell and the second one sticking out at one side. Periodic boundary conditions will take care. The atoms sticking out would be shifted back into the (central) Wigner-Seitz-cell. This might lead to different ESP-charges and Dipole matrix elements. They are only invariant for collective translations of all atoms. The program will stop if a non-valid geometry is detected.

Tags for general section of control.in:

output sub-tag: esp(control.in)

Usage: output esp
Purpose: Calculate the ESP-charges at the atomic position from the full density with default settings. The default settings are to use a radial grid with equidistant radii and without mapping back to the unit-cell in case of periodic boundary conditions. The radial multipliers for the exclusion radius are 3 (Minimum) and 8 (Maximum) times the vdw-radius for the cluster case and 1 (Minimum) and 2 (Maximum) times the vdw-radius for pbc. In both cases 5 radial-shells are used. For pbc a cutoff radius for the Ewald summation of 10 $\mathrm{\AA}$ is used. The maximal multiples of the lattice vectors ( $r_{m}$ ) and the reciprocal lattice vectors ( $k_{m}$ ) used in the summation are 7 (for both). No constraints are applied to the charges.
The keyword can be used multiple times to calculate esp charges with different settings. The optional sub-sub keywords (see below) have to be set after output esp for every calculation separately, otherwise defaults are used.

Optional keywords for output esp:

•

esp spin spin
Select spin spin state (1 or 2) for the calculation of the esp charges from the transition density $\rho_{ij}$ with states i -> j (default 1) at k-point kpoint (default 1).
•

esp state state_one state_two
Select states i=state_one -> j=state_two (i,j $\in$ [1,n_states]) for the calculation of the esp charges from the transition density $\rho_{ij}$ with spin spin (default 1) at k-point kpoint (default 1).
•

esp kpoint k-point
Select states the kpoint (k-point $\in$ [1,n_k_points]) for the calculation of the esp charges from the transition density $\rho_{ij}$ states i -> j (default 1) with spin spin (default 1).
•

esp radius min_radius max_radius
Select the minimal (min_radius) and maximal (max_radius) multiple of the vdw radius to select the volume where the potential is evaluated and the esp charges fitted. Defaults: 3 and 8 for clusters and 1 and 2 for pbc. max_radius should not be larger than smallest lattice vector for pbc and radial gird.
•

esp n_radius n_radius
Select the n_radius number of of radial shells that are used to create points in the selected volume. Default: 5
•

esp equal_grid_n n_grid_x n_grid_y n_grid_z
Select the number of points for the cubic grid in x - n_grid_x, y - n_grid_y and z - n_grid_z direction. Default: 10 10 10
•

esp rm r_m
PBC only. Maximum number r_m of multiples of the lattice vectors that are used in the real space sum of the Ewald summation. Default: 7
•

esp km k_m
PBC only. Maximum number k_m of multiples of the reciprocal lattice vectors that are used in the reciprocal space sum of the Ewald summation. Default: 7
•

esp R_c R_c
PBC only. Real space cutoff radius for the Ewald/Wolf-summation. Default: 20 $\mathrm{\AA}$ .
•

esp pbc_method method
PBC only. Switch between methods for periodic boundary conditions. Choices for method: 1 - method 1 with Ewald summation (3.44, 3.209); 2 - method 2 with Ewald summation (3.44, 3.209); 3 - method 3 with Wolf summation (3.44, 3.210). Default: 2
•

esp grid type
Switch between radial gird - 1, logarithmic radial grid - 2, cubic grid from lattice vectors - 3 and cubic grid from lattice vectors, but truncated in z-direction by maximal vdw-radius (exclude vacuum region for surface slabs) - 4. Default: 3.
•

esp output_cube type
Instead of fitting the ESP-charges to the potential, the potential is written to a cube-file (potential_esp_i.cube). A cubic grid is required. Switch between Hartree potential - 1, XC-potential - 2, X-potential - 3, C-potential - 4, Density - 5, and Hartree potential with the coordinates of voxels - 6.
•

esp output_fit
The ESP-charges are fitted to the potential and the potential calculated from the fitted ESP-charges is written to a cube-file (potential_esp_i.cube). A cubic grid is required.
•

esp use_dip_for_cube
The dipole correction is added to the cube output of the (fitted) potential. Only applied if a dipole correction was used during SCF.

Tag: esp_constraint(control.in)

Usage: esp_constraint method
Purpose: PBC only. Constrain the calculated ESP charges with periodic boundary conditions for method 1 or 2. The constraints will be used for all calculated ESP-charges. Due to technical reasons (array allocations) you have to specify the method you are using in control.in with this keyword. Choices for method are 1 and 2. The actual constraints have to be specified for each atom in geometry.in.

Tag: compute_esp_charges(control.in)

Usage: compute_esp_charges $E_{min}$ $E_{max}$ min_vdw_radius grid_type max_vdw_radius n_radius k-point
Purpose: Calculate and output the ESP-charges for the transition states within the energy window [ $E_{min}$ , $E_{max}$ ] at k-point k-point. Data will be written to file esp_element_k-point.dat. The grid type grid_type (1 - radial, 2 - logarithmic radial, 3 - cubic, 4 - cubic from radii) will be used. The potential will be calculated on the n_radius radial shells within the volume created by min_vdw_radius and max_vdw_radius times the vdw radius of the atoms or in a cube with n_radiusxn_radiusxn_radius oints. The dipole moment calculated from the fitted charges $q_{i}$ at atomic positions $R_{i}$ is also written to file:

\mathbf{d_{ij}}=\sum_{k=1}^{N_{at}}q_{k}\mathbf{R_{k}}

(3.217)

for transition states $i\rightarrow j$ .

Tags for geometry.in:

Tag: esp_constraint(geometry.in)

Usage: esp_constraint constraint_1 constraint_2 constraint_3
Purpose: PBC only. Define the constraints for the fit of the ESP-charges for each atom. Depending on the chosen method (esp_constraint method in control.in). Method 1 needs three parameters ( $\chi,J^{00},w$ ) 3.44 and method 2 needs two parameters ( $q_{0},\beta$ ) 3.44 as defined above.