Contents

• 1   Introduction
  • 1.1   Algorithm
• 2   Quick start guide
  • 2.1   Example
• 3   Installation
• 4   Validation
  • 4.1   Test set
• 5   Known issues and limitations
• 6   Release history highlights
• 7   Copyright and attributions
  • 7.1   Included third-party components
  • 7.2   License

1   Introduction

SLABCC calculates an a posteriori energy correction for charged slab models under 3D periodic boundary conditions (PBC) based on the method proposed in:

Hannu-Pekka Komsa and Alfredo Pasquarello, Finite-Size Supercell Correction for Charged Defects at Surfaces and Interfaces, Physical Review Letters 110, 095505 (2013) DOI: 10.1103/PhysRevLett.110.095505 (Supplements)

This method estimates the error in the total energy of the charged models under 3D PBC due to the excess charge in the real system using Gaussian models. The model charge is assumed to be embedded in a medium with a dielectric-tensor profile ε(k) depending only on a single Cartesian space axis (k) that is orthogonal to the slab. The energy correction is calculated as:

E_corr = E_isolated - E_periodic - qΔV

where E_corr is the total energy correction for the model; E_periodic is the energy of the model charge, calculated by solving the periodic Poisson equation. E_isolated is the energy of the model charge embedded in the dielectric medium and can be determined by extrapolation. q is the total extra charge, and ΔV is the difference between the potential of the Gaussian model charge system and the DFT calculations.

The code can also calculate the charge correction for the 2D models under PBC. The isolated energies for the 2D models are calculated by extrapolation based on the method proposed in:

Hannu-Pekka Komsa, Natalia Berseneva, Arkady V. Krasheninnikov, and Risto M. Nieminen, Charged Point Defects in the Flatland: Accurate Formation Energy Calculations in Two-Dimensional Materials, Physical Review X 4, 031044 (2014) DOI: 10.1103/PhysRevX.4.031044 (Erratum)

And by the cylindrical Bessel expansion of the Poisson equation as proposed in:

Ravishankar Sundararaman, and Yuan Ping, First-principles electrostatic potentials for reliable alignment at interfaces and defects, The Journal of Chemical Physics 146, 104109 (2017) DOI: 10.1063/1.4978238

SLABCC have been initially developed for the Bremen Center for Computational Materials Science (BCCMS)

1.1   Algorithm

• slabcc reads the VASP output files (CHGCAR and LOCPOT), and calculates the extra charge distribution and the potential due to the extra charge. All the input files should correspond to the same geometry.
• The extra charge is modeled by the sum of Gaussian charge distributions as follows:

\rho(r) = \sum_{i}\frac{q_i}{\sigma_{i}^{3}(2\pi)^{3/2}} \exp \left ({- \frac{r_{i}^{2}}{2\sigma_{i}^{2}} } \right )

or trivariate Gaussian charge distributions as:

\rho(r) = \sum_{i}\frac{q_i}{\sigma_{i,x}\sigma_{i,y}\sigma_{i,z}(2\pi)^{3/2}} \exp \left ({- \frac{r_{i,x}^{2}}{2\sigma_{i,x}^{2}} - \frac{r_{i,y}^{2}}{2\sigma_{i,y}^{2}}- \frac{r_{i,z}^{2}}{2\sigma_{i,z}^{2}} } \right )

which gives a charge distribution normalized to q_i with a standard deviation of σ_i (charge_sigma) for each Gaussian distribution i centered at r_i (charge_position).

• The generated Gaussian model charge will be embedded in a dielectric medium with a profile of the form:

\epsilon (k) =  \frac{\epsilon_2-\epsilon_1}{2} \text{erf}\left(\frac{k-k_0 }{\beta}\right)+\frac{\epsilon_2+\epsilon_1}{2}

where k₀ is the interface position in Cartesian k-direction, ε₁ and ε₂ are dielectric tensors on either side of the interface (diel_in & diel_out) and β (diel_taper) defines the smoothness of transition assuming anisotropic dielectric tensor as:

\epsilon =
 \left | \begin{matrix}
   \epsilon_{11} & 0 & 0 \\
   0 & \epsilon_{22} & 0 \\
   0 & 0&  \epsilon_{33}
\end{matrix}  \right |

• The potential due to the charge distribution ρ under 3D PBC embedded in the dielectric medium ε(k) is calculated by solving the Poisson equation in Fourier space:

\epsilon(k) \nabla^2 V(r)+\frac{\partial}{\partial k} \epsilon(k)\frac{\partial}{\partial k}V(r) = -\rho(r)

• A non-linear optimization routine minimizes the difference between our calculated V(r) for the model charge and the V resulted from the VASP calculation by changing the position of the model Gaussian charge, its width, and the position of the slab interfaces.
• The E_periodic is calculated as:

E = \frac{1}{2} \int V(r) \rho(r) \, dr

• E_isolated is calculated the same way as E_periodic but with extrapolation of the fixed model charge embedded in an infinitely large dielectric medium. For the bulk and slab models, the extrapolation is done linearly. For the monolayer models (2D systems), the following equation is used for the extrapolation [10.1103/PhysRevX.8.039902]:

E = c_0 + c_1 x + c_2 x^2 + d e^{-c_3 x}

where c_i are the fitting parameters and

d =  \frac{c_1 - \frac{\partial E_M}{\partial x}}{c_3}

guarantees the correct energy gradient at x(=1/α)→0. E_M being the Madelung energy.

• ΔV is calculated at the position least affected by the model charge.

More information about the algorithms and the implementation details can be found `here`__.

2   Quick start guide

To calculate the charge correction slabcc needs the following files:

• Input parameters file (default: slabcc.in)
• CHGCAR of the neutral system from the VASP calculation (default: CHGCAR.N)
• CHGCAR of the charged system from the VASP calculation (default: CHGCAR.C)
• LOCPOT of the neutral system from the VASP calculation (default: LOCPOT.N)
• LOCPOT of the charged system from the VASP calculation (default: LOCPOT.C)

Input parameters file for a slab should minimally include (all in relative scale [0 1]):

• charge_position: position of the localized charge
• diel_in: dielectric tensor of the slab
• normal_direction: direction normal to the surface
• interfaces: position of the surfaces of the slab (in the normal direction)

2.1   Example

The following examples list the `input parameters`_ to be defined in slabcc.in file, assuming the VASP outputs (LOCPOT and CHGCAR files) are in the same directory.

1. Minimum input: The program models the extra charge with a Gaussian charge distribution localized around the position (charge_position= 0.24 0.56 0.65) in a slab model with a normal direction of (normal_direction = y) and surfaces at (interfaces = 0.25 0.75). The dielectric tensor inside the slab is assumed to be isotropic (diel_in = 4.8):

   charge_position = 0.24  0.56  0.65
   diel_in = 4.8
   normal_direction = y
   interfaces = 0.25 0.75
   

The program will use the default values for the other parameters to:

• Load the CHGCAR of charged and neutralized systems.
• Load the LOCPOT of charged and neutralized systems.
• Calculate the total extra charge from the difference between the charged and neutralized CHGCARs.
• Optimize the charge_position, interfaces and charge_sigma.
• Calculate the total energy correction for the charged system.
• Write all the input parameters used for calculation, the optimized parameters, and the results to the output file.

2. Correction with multiple localized Gaussian charges: If a single charge cannot represent your localized charge properly, you can use multiple Gaussian charges in your model. You have to define the positions of each Gaussian charge, as shown in the example below:

   LOCPOT_charged = CHARGED_LOCPOT
   LOCPOT_neutral = UNCHARGED_LOCPOT
   CHGCAR_charged = CHARGED_CHGCAR
   CHGCAR_neutral = UNCHARGED_CHGCAR
   charge_position = 0.24  0.56  0.65; 0.20  0.50  0.65
   diel_in = 4.8
   normal_direction = a
   interfaces = 0.25 0.75
   

3. Correction for the uniform dielectric medium, e.g., bulk models: You must have the same dielectric tensor inside and outside:

   LOCPOT_charged = CHARGED_LOCPOT
   LOCPOT_neutral = UNCHARGED_LOCPOT
   CHGCAR_charged = CHARGED_CHGCAR
   CHGCAR_neutral = UNCHARGED_CHGCAR
   charge_position = 0.24  0.56  0.65
   diel_in = 4.8
   diel_out = 4.8
   

4. Correction for the monolayers, i.e., 2D models (without extrapolation): To use the Bessel expansion of the Poisson equation for calculating the isolated energy of the 2D models, the in-plane dielectric constants must be equal and the model must be surrounded by a vacuum. Use the extrapolation method (extrapolate=yes) for more general cases:

   LOCPOT_charged = CHARGED_LOCPOT
   LOCPOT_neutral = UNCHARGED_LOCPOT
   CHGCAR_charged = CHARGED_CHGCAR
   CHGCAR_neutral = UNCHARGED_CHGCAR
   2D_model = yes
   charge_position = 0.5 0.4 0.56
   interfaces = 0.66 0.46
   normal_direction = z
   diel_in = 6.28 6.28 1.83
   diel_out = 1
   

5. Correction for the monolayers, i.e., 2D models (with extrapolation): To calculate the isolated energy by fitting the extrapolation results with the non-linear formula, extrapolation to relatively large cell sizes (1/α < 0.2) is necessary. To avoid large discretization errors, the size of the extrapolation grid will be automatically increased:

   LOCPOT_charged = CHARGED_LOCPOT
   LOCPOT_neutral = UNCHARGED_LOCPOT
   CHGCAR_charged = CHARGED_CHGCAR
   CHGCAR_neutral = UNCHARGED_CHGCAR
   2D_model = yes
   extrapolate = yes
   charge_position = 0.5 0.4 0.56
   interfaces = 0.66 0.46
   normal_direction = z
   diel_in = 6.28 6.28 1.83
   

3   Installation

1. Prerequisites:

1. Compiler: You need a C++ compiler with C++14 standard support (e.g. g++ 5.0 or later)
2. BLAS/OpenBLAS/MKL: You can use BLAS+LAPACK for the matrix operations inside the slabcc but it is highly recommended to use one of the high performance replacements, e.g., the OpenBLAS/MKL instead. If you don't have OpenBLAS installed on your system, follow the guide on the OpenBLAS website. Please refer to the Armadillo documentation for linking to other BLAS replacements.
3. FFTW: If you don't have FFTW installed on your system, follow the guide on the FFTW website. Alternatively, you can use the FFTW interface of the MKL.

2. Configuration: Set compilation parameters through environment variables.

1. $CC: C compiler (default: gcc)
2. $CXX: C++ compiler (default: g++)
3. $FFTW_HOME: path to FFTW library home
4. $FFTW_LIB: FFTW library flag (default: -lfftw3)
5. $BLAS_HOME: path to BLAS library home
6. $BLAS_LIB: BLAS library flags (default: -lblas -llapack -lpthread)
7. $EXTRA_FLAGS: extra compiler flags for CC and CXX
8. $LD_EXTRA_FLAGS: extra linker flags

3. Compilation: Run the command make in the bin/ to compile the slabcc.
4. Cleanup: You can run make clean to remove the compiled objects. make distclean additionally removes all the compiled objects of the bundled external libraries.

4   Validation

We are trying to keep the slabcc compatible with as many compilers as possible by using only the standard features of the C++ language. But it is not possible to guarantee this due to the dependency on third-party components. The current version of the slabcc has been build and validated on:

• Ubuntu Linux 16.04

• with GNU C/C++ compilers (5), OpenBLAS, FFTW

• Ubuntu Linux 22.04

• with GNU C/C++ compilers (9,11), OpenBLAS, FFTW
• with GNU C/C++ compilers (11), MKL (2023)
• with Intel oneAPI DPC++/C++ Compiler (2023), MKL (2023)
• with LLVM Clang (14), OpenBLAS, FFTW

• AlmaLinux 8.7

• with GNU C/C++ compilers (8), BLAS, FFTW

• openSUSE Leap 15.4

• with GNU C/C++ compilers (10), BLAS, FFTW

4.1   Test set

You can download a complete test set including input files, input parameters, and expected output here! You can also run the regression tests and verify their results with make test. You'll need numdiff for these tests.

5   Known issues and limitations

• Only orthogonal cells are supported.
• Maximum line length of the input file (slabcc.in) is 4000 bytes.
• Bessel expansion of the Poisson equation cannot be used for the calculation of isolated energies for the 2D models with anisotropic in-plane screening, trivariate Gaussian model change, or the models that are not surrounded by the vacuum (diel_out > 1). The extrapolation method must be used in these cases.

6   Release history highlights

• 2019-06-13: version 0.8 - OO redesign
• 2019-05-14: version 0.7 - Added discretization error mitigation
• 2019-04-04: version 0.6 - Added trivariate Gaussian model charge and selective charge optimization support
• 2019-03-18: version 0.5 - Added 2D model support
• 2018-10-10: version 0.4 - Added spdlog and several user interface and performance improvements
• 2018-07-29: version 0.3 - First public release

7   Copyright and attributions

Copyright (c) 2018-2023, University of Bremen, M. Farzalipour Tabriz

The source code and all the documentation are available under the 2-Clause BSD License. For more information, see license.

This code uses several open-source components, each of which is located under a separate sub-directory of src/. The copyrights of these libraries belong to their respective owners. Any modification made to those codes is also published under the same license. We acknowledge and are grateful to these developers and maintainers for their valuable contributions to this software and, more importantly, to the free software society.

The attributions are also present in the binary file and can be accessed by using --copyright flag.

7.1   Included third-party components

• Armadillo C++ Linear Algebra Library licensed under the Apache License 2.0

• Copyright 2008-2018, Conrad Sanderson
• Copyright 2008-2016, National ICT Australia (NICTA)
• Copyright 2017-2018, Arroyo Consortium
• Copyright 2017-2018, Data61, CSIRO
• This product includes software developed by Conrad Sanderson
• This product includes software developed at National ICT Australia (NICTA)
• This product includes software developed at Arroyo Consortium
• This product includes software developed at Data61, CSIRO

• inih (INI Not Invented Here) licensed under the 3-clause BSD license

• © 2009, Ben Hoyt, et al.

• clara licensed under the Boost Software License 1.0

• © 2014, Phil Nash, Martin Hořeňovský, et al.

• spline (Cubic Spline Interpolation) licensed under the Beer-Ware License 42

• © 2011, Devin Lane

• NLOPT licensed under the GNU LGPL

• © 2007-2014, Massachusetts Institute of Technology
• © 2007-2014, Steven G. Johnson et al.

• spdlog licensed under the MIT License

• © 2016, Gabi Melman, et al.

• Boost.Predef licensed under the Boost Software License 1.0

• © 2005-2018 Rene Rivera
• © 2015 Charly Chevalier
• © 2015 Joel Falcou, et al.

7.2   License

Copyright (c) 2018-2023, University of Bremen, M. Farzalipour Tabriz

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.