SLABCC 1.0.0

Contents

1 Introduction
- 1.1 Algorithm
2 Quick start guide
- 2.1 Example
3 Installation
4 Validation
- 4.1 Test set
5 Command-line parameters
6 Input parameters
7 Results and the generated files
- 7.1 Output files
8 FAQ
9 Known issues and limitations
10 Release history highlights
11 Copyright and attributions
- 11.1 Included third-party components
- 11.2 License

1 Introduction

SLABCC calculates 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). The manual is accessible online.

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:

ρ(r) = ∑_i(q_i)/(σ³_i(2π)^3 ⁄ 2)exp⎛⎝ − (r²_i)/(2σ²_i)⎞⎠

or trivariate Gaussian charge distributions as:

ρ(r) = ∑_i(q_i)/(σ_i, xσ_i, yσ_i, z(2π)^3 ⁄ 2)exp⎛⎝ − (r²_i, x)/(2σ²_i, x) − (r²_i, y)/(2σ²_i, y) − (r²_i, z)/(2σ²_i, z)⎞⎠

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:

ϵ(k) = (ϵ₂ − ϵ₁)/(2) erf⎛⎝(k − k₀)/(β)⎞⎠ + (ϵ₂ + ϵ₁)/(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:

ϵ = ||||| ϵ₁₁ 0 0 0 ϵ₂₂ 0 0 0 ϵ₃₃ |||||

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:

ϵ(k)∇²V(r) + (∂)/(∂k)ϵ(k)(∂)/(∂k)V(r) = − ρ(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 = (1)/(2)∫V(r)ρ(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₀ + c₁x + c₂x² + de^− c₃x

where c_i are the fitting parameters and

d = (c₁ − (∂E_M)/(∂x))/(c₃)

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.

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. Afterwards, slabcc will:

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.

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
```

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

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
```
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

Prerequisites:

Source: Download the latest stable release and extract the files. You can also clone the git repository and use the latest commit on the master branch to get the latest changes.

Compiler: You need a C++ compiler with C++14 standard support (e.g. g++ 5.0 or later)

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.

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.

Configuration: Set compilation parameters through environment variables.

$CC: C compiler (default: gcc)

$CXX: C++ compiler (default: g++)

$FFTW_HOME: path to FFTW library home

$FFTW_LIB_FLAG: FFTW library flag (default: -lfftw3)

$BLAS_HOME: path to BLAS library home

$BLAS_LIB_FLAG: BLAS library flags (default: -lblas -llapack -lpthread)

$EXTRA_FLAGS: extra compiler flags for CC and CXX

$LD_EXTRA_FLAGS: extra linker flags

Compilation: Run the command make in the bin/ to compile the slabcc.
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. However, it is not possible to guarantee this due to the dependency on third-party components. The current version of the slabcc has been built and validated inside containers based on:

Ubuntu Linux 16.04 (dockerfile)

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

Ubuntu Linux 22.04 (dockerfiles 1 2)

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 (dockerfile)

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

openSUSE Leap 15.4 (dockerfile)

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

Older versions of the code were also being tested on MS Windows 10 (latest toolchains: Intel C/C++ compilers 19.0.4 and Microsoft C/C++ compilers 19.20.27508 linked to MKL 19.0.4 and FFTW 3.3.5), support for which is currently dropped.

4.1 Test set

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

5 Command-line parameters

You can run slabcc without any additional options. Alternatively, you can use the following options to modify its behavior:

`-h, --help`	display usage information (this list)
`-i, --input <input_file>`	slabcc input file name
`-o, --output <output_file>`	slabcc output file name
`-l, --log <log_file>`	slabcc log file name
`-d, --diff`	calculate charge and potential differences only
`-m, --man`	show quick start guide
`-v, --version`	show slabcc version
`-c, --copyright`	show copyright information

6 Input parameters

slabcc reads all its parameters from the input file (by default, slabcc.in). You can change the input file's name using the command-line parameters. The input file is processed as follows:

Lines starting with # will be treated as comments. Inline comments are also allowed.
Double quotation marks will be removed from the strings.
A warning will be issued for any unidentified parameter.
A warning will be issued for the use of the deprecated parameters.
All the coordinates must be in fractional form [0-1].
Boolean (True/False) parameters can also be declared as 0/1, on/off, yes/no, .true./.false.
Parameter names can be written in small or capital letters.
For vectors and matrices, columns are separated by a " "(space), while rows are separated by a ";" (semicolon).
Lines starting with a space " " will be treated as the continuation of the last parameter's value.
Subsequent definitions for any parameter will be concatenated with the existing definition.

Parameter	Description and the options / `example`	Default value
`2d_model`	Calculate the charge correction for a 2D model	false
`charge_fraction`	Fraction of the total extra charge in each localized Gaussian model charge (in the case of multiple Gaussian charges) `charge_fraction = 0.4 0.6`	The extra charge will be equally divided among all positions
`charge_position`	Center of the model Gaussian charges `charge_position = 0.2 0.5 0.3` `charge_position = 0.2 0.2 0.2; 0.3 0.4 0.3`
`charge_rotation`	Rotation angles around each axis for the trivariate Gaussian charges in arc degrees (-90, 90) `charge_rotation = 0 45 0` `charge_rotation = 45 0 0; 0 -10 0`	0
`charge_sigma`	width of each localized Gaussian charge. It can be 1 or, in the case of trivariate models, 3 parameters per localized Gaussian charge. For trivariate Gaussian models, defining a single parameter per charge, sets the sigma values to be equal in all directions. for a single Gaussian charge `charge_sigma = 1` for multiple Gaussian charges `charge_sigma = 1 1.5` for two trivariate Gaussian charges `charge_sigma = 1 2 3; 1.5 2.5 3.5;`	1 (for each charge in each direction)
`charge_trivariate`	Use trivariate Gaussian model along the main axis	false
`CHGCAR_charged`	Charge density file (CHGCAR) of the charged system `CHGCAR_charged = CHGCAR1`	CHGCAR.C
`CHGCAR_neutral`	Charge density file (CHGCAR) of the neutral system `CHGCAR_neutral = CHGCAR2`	CHGCAR.N
`diel_in`	Diagonal elements of the static dielectric tensor inside the slab. If only a single value is given, all of them will be assumed to be equal. `diel_in = 3` `diel_in = 3 4 5`	1
`diel_out`	Diagonal elements of the static dielectric tensor outside of the slab	1
`diel_taper`	The steepness of the transition between diel_in and diel_out (β in the dielectric profile formula)	1
`extrapolate`	Calculate the isolated energy using the extrapolation method	opposite of the `2d_model` parameter
`extrapolate_grid_x`	Extrapolation grid size multiplier. The number of grid points in each direction will be multiplied by this value. extrapolate_grid_x > 1 will use a larger grid in the extrapolations, which will increase its accuracy but will requires more memory and computational power. extrapolate_grid_x = 1 will use the same grid size as the VASP input files in the extrapolation. extrapolate_grid_x < 1 will use a smaller grid for the extrapolations, which increases the speed and decreases the memory usage, but the energies for the higher orders of the extrapolation may not be accurate! `extrapolate_grid_x = 1.8`	1
`extrapolate_steps_number`	Number of extrapolation steps in the calculation of E_isolated	10: for 2D models 4: for the rest
`extrapolate_steps_size`	Size of extrapolation steps with respect to the initial supercell size	1: for 2D models 0.5: for the rest
`interfaces`	Interfaces of the slab in normal direction `interfaces = 0.11 0.40`	0.25 0.75
`LOCPOT_charged`	Local potential file (LOCPOT) of the charged system `LOCPOT_charged = LOCPOT1`	LOCPOT.C
`LOCPOT_neutral`	Local potential file (LOCPOT) of the neutral system `LOCPOT_neutral = LOCPOT2`	LOCPOT.N
`normal_direction`	Normal direction of the slab: one of x/y/z or a/b/c corresponding to the 1st, 2nd and 3rd vectors in the input file's cell vectors `normal_direction = b`	z
`optimize`	Optimizer master switch. Upon deactivation, it takes precedence over all the other optimization options: `optimize_charge_fraction`, `optimize_charge_position`, `optimize_charge_rotation`, `optimize_charge_sigma`, and `optimize_interfaces` true: evaluate each of the optimization switches individually false: deactivate all optimization switches	true
`optimize_algorithm`	Optimization algorithm in the NLOPT library BOBYQA : Bound Optimization BY Quadratic Approximation [1] COBYLA: Constrained Optimization BY Linear Approximation [2] SBPLX: S.G. Johnson's implementation of the Subplex (subspace-searching simplex) algorithm [3] `optimize_algorithm = SBPLX`	BOBYQA
`optimize_charge_fraction`	true: find the optimal values for the model's charge_fraction parameter to construct the best model charge which mimics the potential obtained from the VASP calculation false: do not change the charge_fraction parameter	true
`optimize_charge_position`	true: find the optimal values for the model's charge_position parameter to construct the best model charge which mimics the potential obtained from the VASP calculation false: do not change the charge_position parameter	true
`optimize_charge_rotation`	true: find the optimal values for the model's charge_rotation parameter to construct the best model charge which mimics the potential obtained from the VASP calculation. This can only be used for the trivariate Gaussian models. false: do not change the charge_rotation parameter	false
`optimize_charge_sigma`	true: find the optimal values for the model's charge_sigma parameter to construct the best model charge which mimics the potential obtained from the VASP calculation false: do not change the charge_sigma parameter	true
`optimize_grid_x`	Optimization grid size multiplier. The number of the grid points in each direction will be multiplied by this value. optimize_grid_x > 1 will use a larger grid in the optimization which will increase its accuracy but will requires more memory and the computational power. [usually this is not necessary] optimize_grid_x = 1 will use the same grid as the VASP input files in the optimization optimize_grid_x < 1 will use a smaller grid for the optimization which increases the speed and decreases the memory usage but the parameters obtained using very small grid sizes may be inaccurate!	0.8
`optimize_interfaces`	true: find the optimal values for the model's interfaces to construct the best model which mimics the potential obtained from the VASP calculation false: do not change the position of interfaces in the model charge	true
`optimize_maxsteps`	Maximum number of optimization steps `optimize_maxsteps = 2000`
`optimize_maxtime`	Maximum time for optimization in minutes `optimize_maxtime = 1440`
`optimize_tolerance`	Relative optimization tolerance (convergence criteria) for root mean square error of the model potential	0.01
`slab_center`	Center of the slab. During the calculations, everything will be shifted to keep this point at the center. This point must be inside of the slab. `slab_center = 0.2 0.7 0.5`	0.5 0.5 0.5
`verbosity`	Verbosity of the program [4] 0: No extra info. Only write the output file. Logging is disabled. 1: Display calculated energy correction terms. Write the planar averaged potential and charge for the Gaussian model charge and the extra-charge of QM calculations in the direction normal to the slab surface. 2: Write extra-charge density, extra-charge potential and dielectric profiles. Display debug info including the compilation machine info and a few important enviroment variables. 3: Write the planar averaged files in all directions. 4: Display the time passed since the start of slabcc (in seconds) and a description of each calculation step (trace mode)	1

[1]	M.J.D. Powell, The BOBYQA algorithm for bound constrained optimization without derivatives, Department of Applied Mathematics and Theoretical Physics, Cambridge England, technical report NA2009/06 (2009).

[2]	M.J.D. Powell, Direct search algorithms for optimization calculations, Acta Numerica, Vol. 7(1998) pp. 287-336

[3]	T.H. Rowan, Functional Stability Analysis of Numerical Algorithms, Ph.D. thesis, Department of Computer Sciences, University of Texas at Austin, 1990.

[4]	Each verbosity level includes all the outputs from the lower verbosity options. Check the files table for complete list of the output files.

7 Results and the generated files

slabcc writes its calculated energy correction values to the standard output as well as the output file. All reported energy values are in eV.

Depending on the verbosity level of your choice, you may get additional reports from each part of the calculation in the standard output and/or extra output files.

7.1 Output files

The parsed input variables or their default values and the calculation results will be written to the output file (by default: slabcc.out) You can change this file’s name using the command-line parameters. A typical output file is shown below:

# Parameters read from the file or their default values:
2d_model = no
charge_fraction = 1
charge_position = 0.5 0.5 0.37;
charge_rotation = 0 0 0;
charge_sigma = 1;
charge_trivariate = no
CHGCAR_charged = ../03-V_Cl_pos/CHGCAR
CHGCAR_neutral = ../02-V_Cl/CHGCAR
diel_in = 2.45
diel_out = 1
diel_taper = 1
extrapolate = yes
extrapolate_grid_x = 1
extrapolate_steps_number = 4
extrapolate_steps_size = 0.5
interfaces = 0 0.375
LOCPOT_charged = ../03-V_Cl_pos/LOCPOT
LOCPOT_neutral = ../02-V_Cl/LOCPOT
normal_direction = z
optimize_algorithm = COBYLA
optimize_charge_fraction = yes
optimize_charge_position = yes
optimize_charge_rotation = no
optimize_charge_sigma = yes
optimize_grid_x = 0.8
optimize_interfaces = yes
optimize_maxsteps = 0
optimize_maxtime = 0
optimize_tolerance = 0.01
slab_center = 0.5 0.5 0.25
verbosity = 5

[Optimized_model_parameters]
interfaces_optimized = 0.942000748357 0.455672787711
charge_sigma_optimized = 1.4132676877
charge_position_optimized = 0.501460639345 0.50145532106 0.385476689493;

[Results]
dV = -0.00291385176718
E_periodic of the model charge = 2.0404453156f
E_isolated of the model charge = 2.59716677886
Energy correction for the model charge (E_iso-E_per-q*dV) = 0.559635314929

Planar average files are written as the double column in plain text format. The first column represents the coordinates along the axis (in Angstrom) and the second column is the planar average value. The files are named as: "slabcc_{1}{2}{XXX}.dat" where:

{1}: N: Neutral system, C: Charged system, D: Difference
{2}: X/Y/Z: Corresponds to the 1st, 2nd, and the 3rd axis in the input files
{XXX}: CHG: CHGCAR, POT: LOCPOT

All the possible output files and the minimum value of the verbosity parameter for activation of each are listed in the table below:

file name	Description	verbosity
slabcc_CXCHG.dat	Planar average of charged CHGCAR file in X direction	3
slabcc_CXPOT.dat	Planar average of charged LOCPOT file in X direction	3
slabcc_CYCHG.dat	Planar average of charged CHGCAR file in Y direction	3
slabcc_CYPOT.dat	Planar average of charged LOCPOT file in Y direction	3
slabcc_CZCHG.dat	Planar average of charged CHGCAR file in Z direction	3
slabcc_CZPOT.dat	Planar average of charged LOCPOT file in Z direction	3
slabcc_D.CHGCAR	Difference in the neutral and charged CHGCAR files	2
slabcc_D.LOCPOT	Difference in the neutral and charged LOCPOT files	2
slabcc_DIEL.dat	Generated dielectric profile (ε₁₁ ε₂₂ ε₃₃) along the normal axis to the surface	3
slabcc_DXCHG.dat	Planar average of extra charge (neutral and charged difference) CHGCAR file in X direction	3
slabcc_DXPOT.dat	Planar average of extra charge (neutral and charged difference) LOCPOT file in X direction	3
slabcc_DYCHG.dat	Planar average of extra charge (neutral and charged difference) CHGCAR file in Y direction	3
slabcc_DYPOT.dat	Planar average of extra charge (neutral and charged difference) LOCPOT file in Y direction	3
slabcc_DZCHG.dat	Planar average of extra charge (neutral and charged difference) CHGCAR file in Z direction	3
slabcc_DZPOT.dat	Planar average of extra charge (neutral and charged difference) LOCPOT file in Z direction	3
slabcc_M.CHGCAR	CHGCAR of the Gaussian model	2
slabcc_M.LOCPOT	LOCPOT of the Gaussian model	2
slabcc_MXCHG.dat	Planar average of model charge in X direction	3
slabcc_MXPOT.dat	Planar average of model potential in X direction	3
slabcc_MYCHG.dat	Planar average of model charge in Y direction	3
slabcc_MYPOT.dat	Planar average of model potential in Y direction	3
slabcc_MZCHG.dat	Planar average of model charge in Z direction	3
slabcc_MZPOT.dat	Planar average of model potential in Z direction	3
slabcc_NXCHG.dat	Planar average of neutral CHGCAR file in X direction	3
slabcc_NXPOT.dat	Planar average of neutral LOCPOT file in X direction	3
slabcc_NYCHG.dat	Planar average of neutral CHGCAR file in Y direction	3
slabcc_NYPOT.dat	Planar average of neutral LOCPOT file in Y direction	3
slabcc_NZCHG.dat	Planar average of neutral CHGCAR file in Z direction	3
slabcc_NZPOT.dat	Planar average of neutral LOCPOT file in Z direction	3

Note: The planar averaged potential and charge files corresponding to the normal direction will be written in verbosity = 1

8 FAQ

How do I obtain the CHGCAR and LOCPOT files from VASP calculations? You can add the following tags to your INCAR file to get the LOCPOT and CHGCAR files:
```
LVTOT = .TRUE.
LVHAR = .TRUE.
LCHARG = .TRUE.
```

After obtaining the files for your charged system, do the calculation again without relaxing (changing) the geometry to get the necessary files for the neutralized system.

Do I need to perform spin-polarized calculations in VASP? Although the slabcc only reads the sum of both spins, for a proper description of the charge distribution in your system, you may need to perform a spin-polarized calculation.
How can I speed up the model parameter optimization process? You can try using a different optimization algorithm or improving the initial guess for the model parameters to speed up the optimization. As a last resort, you can also use a smaller computation grid for the optimization (optimize_grid_x < 1), or increase the optimization convergence criteria (optimize_tolerance) to speed up the process, but the accuracy of the obtained results may not be adequate. Ideally, you should use the optimzied parameters on a coarser grid as the initial guess for optimization on a finer (initial) grid.
Why do I need to provide an initial guess for the parameters that will be optimized? The optimization algorithms used in slabcc are local error minimization algorithms. Their success and performance highly depend on the initial guess for the provided parameters.
How should I decide on the initial guess for the parameters that will be optimized? As a rule of thumb, start with a single Gaussian charge as your model. Set its position to the expected position of the charge localization. Use the location of the surface atoms as the interface position. You can use the "-d" switch in the command line (./slabcc -d) to just generate the CHGCAR and the LOCPOT files for the extra charge and their planar averages without shifting the input files to the slab_center. These files will guide you on how to provide the initial guess for the input parameters.
Can I turn off the optimization for the input parameters? Yes. But optimization ensures the model charge mimics the original localized charge at large distances as closely as possible. If you turn off the optimization, you must be aware of the possible side-effects and definitely check your results.

Can I run the slabcc on a computational cluster? Yes. BUT… Although slabcc hugely benefits from the multicore architecture of the computation nodes using OpenMP, it has not yet been parallelized using MPI. Therefore, it won’t use more than one machine at a time.
Is the slabcc free? Can I use its source code in my own software? slabcc is released under the 2-Clause BSD license which permits this software to be modified, redistributed, and/or used for commercial purposes provided that the source retains the original copyright owner's name (University of Bremen, M. Farzalipour Tabriz) and full text of the license (LICENSE.txt)
How accurate are the slabcc results? The accuracy of the final results depends on various factors, including the accuracy and grid size of the input files and the provided input parameters. The optimization algorithm used for parameter estimation is a non-linear local optimizer, which means that the result will highly depend on its initial conditions. Models with different numbers of Gaussian charges have different accuracy and may be compared with caution. In the case of models with multiple charges, the results must be vigorously checked. You must always do your own testing before using the results. There are a few known issues and limitations to the slabcc code and its algorithm. Also keep in mind that this is free software, and as the license explicitly mentions, there is absolutely no warranty for its fitness for any particular purpose.

How can I check the slabcc results? slabcc can calculate the planar averaged potential and charge files for the extra charge in the input files and the model Gaussian charge. You should compare the model charge distribution and potential, especially in the direction normal to the surface, and compare them to the original VASP results. For example, if z is the normal direction in your slab model (normal_direction = z), then you should compare slabcc_MZCHG.dat and slabcc_MZPOT.dat, with slabcc_DZCHG.dat and slabcc_DZPOT.dat, respectively. Check the files table for a complete list of the output files.

Another method to test the effectiveness of the charge correction is to increase the thickness of the vacuum in your slab model and check the charge-corrected total energies. If the charge correction is done properly, the energy values must be independent of the (adequately large) vacuum thickness.

How should I cite the slabcc? Please cite the slabcc as: (You can download the citation in the RIS format from here!)

Meisam Farzalipour Tabriz, Bálint Aradi, Thomas Frauenheim, Peter Deák, SLABCC: Total energy correction code for charged periodic slab models, Computer Physics Communications, Vol. 240C (2019), pp. 101-105, DOI: 10.1016/j.cpc.2019.02.018

How can I extract the files in slabcc_data.tar.xz? You can use Tar + XZ Utils as:

tar -xvf slabcc_data.tar.xz

Alternatively, you can use WinRAR or 7zip.

Something is not working! What should I do?

If you need help with compiling the code or running it on a cluster, please contact your system administrator.

If you have found a bug in the code, please report it here.

9 Known issues and limitations

Only orthogonal cells are supported.
The 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.

10 Release history highlights

2023-06-05: version 1.0 - Improved build, error mitigation, and correctness checks
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

11 Copyright and attributions

The source code and all 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.

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

11.1 Included third-party components

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

Copyright 2008-2023 Conrad Sanderson

Copyright 2008-2016 National ICT Australia (NICTA)

Copyright 2017-2023 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 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.

11.2 License

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

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
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.