The purpose of this package is to facilitate the computation of electrochemical (EC) interfaces with VASP. In particular, it allows computing grand potential variation curves.
It is inspired by this code written by Arthur Hagopian and Jean-Sebastien Filhol and published with 10.1021/acs.jctc.1c01237. The procedure is similar, with the following differences:
- The number of electrons, charge added and charge removed are stored in a file, called
ec_interface.yml
. - To get more precision, the code looks into
vaspout.h5
for some results (number of electrons, free and fermi energies), so it requires to compile VASP with-DVASP_HDF5
. - The procedure to locate the reference potential is slightly different (this should not affect the results).
- A way to plot the graph is not provided, but the
ec_results.csv
file will allow you to plot the results and make further analysis.
To use this package, you need to have a running VASP version, compiled with HDF5 support (-DVASP_HDF5
).
This code has been tested with VASP 6.3.2 and 6.4.1.
If you want to perform liquid phase calculations (either with the Poisson-Boltzmann or Homogeneous Background approaches), you also need to compile VASP with VASPsol (by following the instructions in their repository).
Finally, to install this package, you need a running Python 3 installation (Python >= 3.10 recommended), and
pip3 install git+https://github.com/pierre-24/ec-interface.git
Note: as this script install programs, you might need to add them (such as $HOME/.local/bin
, if you use --user
) to your $PATH
.
The following paragraph contains the standard way of getting a grand potential curve, with some tips and tricks along the way.
To perform an EC interface calculation, you need the following files in the same directory:
-
An
INCAR
file. In order to compute the information required for the EC calculation, it should at least contain the following instructions:LCHARG = .TRUE. LVHAR = .TRUE. LVTOT = .TRUE.
In practice, the following parameters are a good starting point:
!! convergence and precision both needs to be strict EDIFF = 1E-6 PREC = Acurate !! parameters required for EC calculations: LCHARG = .TRUE. ! Create CHGCAR LVHAR = .TRUE. ! Create LOCPOT LVTOT = .TRUE. ! Include VXC !! solvatation (optional) LSOL = .TRUE. ! Enable solvatation with VASPsol EB_K = 80.0 ! dielectric constant of the solvent TAU = 0.0 ! neglect cavitation energy NC_K = 1E-5 ! critical density parameter SIGMA_K = 0.6 ! Cavity shape parameter (default) LAMBDA_D_K = 1.5 ! for PBM, remove if HBM.
Refer to the VASP manual, 10.1021/acs.jctc.1c01237, and the VASPsol documentation for more details on these parameters and their values.
You might also want to increase the value of
NELM
(it may be more difficult to converge those calculations, especially with PCM),NBANDS
(all bands might get parrtially occupied as you add electrons), andLREAL
(as VASP is generally complaining). -
A
POSCAR
file (which contains a slab geometry) and its correspondingPOTCAR
(that you can generate withei-create-potcar
, see below). PAW potentials are strongly recommended.It is assumed that the C lattice vector matches the Z axis.
-
A
KPOINTS
file. -
A
ec_interface.yml
file (see below).
You can get info about a slab by runnning ei-check-slab
:
ei-check-slab POSCAR
Among others, the interslab distance (i.e., the vacuum between two repetition of the slab) should be adjusted (see 10.1021/acs.jctc.5b00170). See below for a tool to do so. Note that there should be enough vacuum to get an accurate value for the reference (vacuum) potential, which is used to compute the work function.
You can generate the POTCAR
file using ei-create-potcar
:
ei-create-potcar POSCAR -p /path/to/potpaw/
where /path/to/potpaw
is the directory containing all PAW potentials (H
, H_AE
, etc).
If your POSCAR
contains atoms for which you want to use an alternate potential, use -P
with a comma separated list of symbol=potential
, e.g., -P Li=Li_sv,C=C_h
.
You can get the number of electron that your system contains with ei-get-nzc
:
ei-get-nzc POSCAR POTCAR
Then, create a ec_interface.yml
.
You can start from the following:
ne_zc: 21 # this is the number of electrons in your system
ne_added: 0.05 # number of electron (charge) to add
ne_removed: 0.05 # number of electron to remove
step: 0.01 # step for adding/removing electrons
additional: [] # (optional) additional points, if required
The value of ne_zc
is the number of electrons that your system normally contains (zero charge).
Check, e.g., for NELECT
in a preliminary VASP output, or the output of ei-get-nzc
.
Adjust step
: it should be ~10⁻⁴ e Å⁻², to be multiplied by (twice) the slab surface, as recommended in 10.1021/acs.jctc.5b00170.
Then, pick a value for ne_added
and ne_removed
, which should be a multiple of step
.
Note that their values depend on the capacitance of your system... But they should not be large, since the corresponding change of potential should be within the acceptable range (see, e.g., 10.1088/1361-648X/ac0207 for more details).
Finally, you can create the inputs for the calculation using:
ei-make-directories
Note that by default, the POSCAR
, POTCAR
and KPOINTS
files are referred to by using symlinks, as they should be the same for all calculations.
The parameters are read from ec_interface.yml
. Additional options are:
-v
: to get details about the creation of the different directories,-f
: force the re-creation of directories if existing, and-c
: copy file inside created directories, instead of using symlinks to do so.
Run VASP in each and every directory that was created in the previous step.
If you are curious, you will notice that the difference between directories is the value of NELECT
, which sets the number of electrons in the calculation.
If every calculation went well, just run:
ei-extract-data
The parameters are read from ec_interface.yml
.
Use the -v
option to get details about extraction.
At the end of the procedure, a ec_results.h5
file should be created, containing the different data (number of electrons, free and fermi energies, vacuum and average potentials) in binary form.
Furthermore, in each directory, you'll find a charge_density_xy_avg.csv
and a local_potential_xy_avg.csv
file, which contains a XY-average of the CHGCAR
and LOCPOT
file, respectively: the first column contains the Z coordinates, the second correspond to the XY-averaged value times unit volume, while the third contains the plain XY-averaged value.
Note that you can also use ei-xy-average
to get the same information:
ei-xy-average CHGCAR > chg.csv
For CHGCAR
, to count the electrons in certain regions, you can also use:
ei-charge-intg CHGCAR
This program will select regions (in the z-direction) of low and high electron occupations and integrate over those.
Option -t
(threshold) allow to change the threshold for the detection of low/high regions.
Finally, run:
ei-compute-fee > results.csv
The parameters are read from ec_interface.yml
, while the data are read from ec_results.h5
.
Different options can be used to compute the FEE with different methods, see below.
At the end of the procedure, a results.csv
file (containing tab-separated data) should be created, which contains the following dataset, with:
- The number of electrons in the calculation (
NELECT
inINCAR
); - The free energy (equal to
E0
inOSZICAR
); - The vacuum potential (almost equal to the value of
FERMI_SHIFT
reported by VASPsol); - Fermi energy (equal to the value of
E-Fermi
inOUTCAR
); - The average potential of the unit cell.
And then another dataset, with:
-
The charge added or removed to the system;
-
The (absolute) work function, as
vacuum_potential - fermi_energy
; -
The relative potential versus reference,
work_function - ref
(ref=4.5V
by default, your can change this value with--ref
); -
The corresponding grand potential value. By default, it is computed as
free_energy - dn * fermi_energy
.Note: this might not be the correct free electrochemical energy you are looking for, and other methods are available: see this document for more information (TL;DR: use either
--pbm
or--hbm xxx
, wherexxx
is the fraction of active electrons).Note: the average potential should be about 0 at PZC. If it is not the case, you might want to use
--shift-avg
to set average to 0.
Please refer to 10.1039/c9cp06684e (and reference therein) for different information that you can extract from those data, such as the surface capacitances, the fukui functions, etc.
See this archive, which contains an example of calculation on a Li (100) slab using the PBM approach, inspired by 10.1021/acs.jctc.1c01237. It gives the following curve:
A capacitance of 0.0535 e/V is extracted from this curve using its second derivative. Due to the strong anharmonicity (the curve is clearly not symmetric around PZC), the actual value should be a little smaller.
To manipulate the geometries, it is generally recommended to use tools such as ASE or pymatgen to manipulate the geometries.
However, ec-interface
comes with a few handy (tough simple) tools to perform routine operations.
Use them with care!
- To adjust the interslab distance, you can use
ei-set-vacuum
, which creates a new geometry while enforcing vacuum size (i.e., the size of the last lattice vector). For example, to adjust the vacuum size to 25 Å:The new geometry is saved inmv POSCAR POSCAR_old ei-set-vacuum POSCAR_old -v 25.0 -o POSCAR
POSCAR
. Note that the slab is z-centered by the procedure. - To turn an XYZ geometry into POSCAR, you can use
ei-to-vasp-geometry
:It comes with a few options, such asei-to-vasp-geometry molecule.xyz --lattice=10,10,10 -o POSCAR
--lattice
to set the lattice vectors and--sort
to group atoms types (so that it is easier to create the POTCAR). - To merge two POSCARs, you can use
ei-merge-poscar
:It allows to merge two geometries. The lattice of the first geometry (hereei-to-vasp-geometry POSCAR_cell POSCAR_substrate --shift=5,5,7 -o POSCAR
POSCAR_cell
) is used in the final one. The--shift
option allows to reposition the second molecule in the first. Note thatSelective dynamics
information are kept.
Furthermore, since the idea is to compute the properties for different number of electrons, an insightful byproduct are the Fukui functions.
In particular, CHGCAR
files with:
ei-fukui EC_20.000/CHGCAR EC_20.050/CHGCAR -d 0.05 -o CHGCAR_fukui
where the -d
argument gives the value of CHGCAR
is the reference density ( -s
option is used, a symmetric difference formula is used, and the first file must then contain CHGCAR_fukui
file contains the Fukui function for
Contributions, either with issues or pull requests are welcomed.
If you want to contribute, this is the usual deal: start by forking, then clone your fork and use the following install procedures instead.
cd ec-interface
# definitely recommended in this case: use a virtualenv!
python -m venv virtualenv
source venv/bin/activate
# install also dev dependencies
make install
-
A good place to start is the list of issues. In fact, it is easier if you start by filling an issue, and if you want to work on it, says so there, so that everyone knows that the issue is handled.
-
Don't forget to work on a separate branch. Since this project follows the git flow, you should base your branch on
main
, not work in it directly:git checkout -b new_branch origin/main
-
Don't forget to regularly run the linting and tests:
make lint make test
Indeed, the code follows the PEP-8 style recommendations, checked by
flake8
. Having an extensive test suite is also a good idea to prevent regressions. -
Pull requests should be unitary, and include unit test(s) and documentation if needed. The test suite and lint must succeed for the merge request to be accepted.
My name is Pierre Beaujean, and I have a Ph.D. in quantum chemistry from the University of Namur (Belgium). I'm the main (and only) developer of this project, used in our lab. I use VASP in the frame of my post-doctoral research in order to study batteries and solid electrolyte interphrase, and I developed this project to ease my life.