Atomic-scale capacitance

In this tutorial you will perform an atomic-scale simulation of a parallel plate nano-capacitor, which means that the gap between the two parallel metal surfaces is in the nanometer range. This is currently not a realistic device, but the purpose is to show a methodology that can be used to compute the capacitance of other structures, including metal/semiconductor/metal structures.

../../_images/capacitor.png ../../_images/ppcap4.gif


Capacitor physics
Applying a bias \(V\) across the capacitor leads to a build-up of oppositely signed charge \(Q\) on the plates, end therefore an electric field \(E\) between them. The electric field strength can be written as \(E=\sigma/\varepsilon = V/d\), where \(\sigma = Q/A\) is the charge density and \(\varepsilon = k \varepsilon_0\) is the permittivity of the spacer layer between the plates (vacuum or dielectric). The capacitance can then be calculated from \(C = Q/V = Q / Ed = A \varepsilon/d\), where \(d\) is the distance between the plates.

The equations given above are valid for an ideal macroscopic parallel plate capacitor, but a nano-scale capacitor may in general be quite different. The main method used for computing the capacitance from first principles is to evaluate the electrostatic energy, and from that the capacitance. We will also show a slightly simpler approach based on the Mulliken charge populations of the atoms.

Build the parallel plate capacitor

Open the Builder builder_icon, and do the following:

  • Use Add ‣ From Database to add gold to the Stash.
  • Open the Builders ‣ Cleave plugin, and create a Au(111) electrode by entering the Miller indices as [111].
  • Repeat the structure 6 times in the C direction using Bulk Tools ‣ Repeat.
  • Select and delete the two atoms in the middle – they have atom indices 3 and 14.
  • Select the surface atoms (indices 7 and 8), and turn them into ghost atoms using the toolbar icon ghost_icon.
  • The distance between the two surfaces is still a bit too small. To change it, open Coordinate Tools ‣ Translate, and select all the atoms in the right-hand surface (including the nearest ghost atom) and move the atoms +2 Å in the Z direction.
  • Turn the bulk structure into a device configuration via Device Tools ‣ Device From Bulk and use the default settings.

The result will be the device as shown below.



The ghost atoms are used to improve the description of the electron density between the two surfaces. They are simply extra basis sets placed in the empty region. This technique, which is sometimes referred to as a vacuum basis, is also used for calculating the work function of a surface, as explained in a tutorial: Computing the work function of a metal surface using ghost atoms.


Zero bias

Send the device configuration to the Script Generator script_generator_icon by clicking the sendto_icon icon, and set up the script as follows:

  • Add a New Calculator calculator_icon and use the DFT-LDA engine with a SingleZetaPolarized basis set. Then set a 9x9x100 k-point grid.

  • Add the following Analysis objects analysis_icon: MullikenPopulation, ElectrostaticDifferencePotential, and ElectronDifferenceDensity.

  • Change the output filename to

  • Save the script as

Send the script to the Job Manager job_manager_icon and execute the calculation. It takes 15-20 minutes on a laptop.

Finite Bias

To compute the capacitance it is necessary to compute \(dq/dV\), where \(dq\) is the induced charge on the surfaces and \(V\) the bias voltage. Use the following script to perform the calculations for bias voltages of 0.25, 0.5, 0.75, and 1.0 V:

device_configuration = nlread('', DeviceConfiguration)[-1]
calculator = device_configuration.calculator()
for voltage in [0.25, 0.5, 0.75, 1.0]:
    # Set new calculator with modified electrode voltages on the configuration
    # use the self consistent state of  the old calculation as starting input.
          calculator(electrode_voltages=(0.5*voltage*Volt, -0.5*voltage*Volt)),
    nlsave('', device_configuration, object_id="SCF %s" % voltage)
    # -------------------------------------------------------------
    # Mulliken population
    # -------------------------------------------------------------
    mulliken_population = MullikenPopulation(device_configuration)
    nlsave('', mulliken_population, object_id="Mulliken %s" % voltage)
    # -------------------------------------------------------------
    # Electrostatic difference potential
    # -------------------------------------------------------------
    electrostatic_difference_potential = ElectrostaticDifferencePotential(device_configuration)
    nlsave('', electrostatic_difference_potential, object_id="EDP %s" % voltage)
    # -------------------------------------------------------------
    # Electron difference density
    # -------------------------------------------------------------
    electron_difference_density = ElectronDifferenceDensity(device_configuration)
    nlsave('', electron_difference_density, object_id="EDD %s" % voltage)

You can download the scipt here: Run it using the Job Manager or execute it from a terminal:

$ atkpython > bias_loop.out

It takes about an hour to finish this calculation on a laptop.


The previously converged calculation is used as an initial guess for the next calculation to speed up the SCF convergence. The object_id makes it easier to identify the data in post-SCF analysis.


We will introduce two different ways to calculate the capacitance:

  • a simple approach based on the induced charge;
  • a more general method based on evaluating the electrostatic energy.

Approach Based on the Induced Charge

The finite bias will build up a net positive/negative charge on the two gold surfaces, as shown in the plots below.


Now, the capacitance can in principle be computed as \(C=dq/dV\), where \(q\) is the total accumulated charge on one surface as a function of the bias. From the slope of the straight line you can therefore obtain the capacitance of the nano-capacitor:


The plots above, as well as the entire data analysis based on the Mulliken populations, can be achieved by running this script: It is a long script but a great opportunity to learn about custom data analysis in Python and ATK.


The script reports a capacitance of 5.77*10-22 Farad, as computed from fitting a straight line to the data points. To put this number into perspective, consider the capacitance of an idel parallel plate capacitor, written as \(C = \varepsilon A/d\), where \(\varepsilon=1\) for the vacuum permittivity and the area \(A\) of each plate can be calculated from the device structure. You can then calculate the parallel plate distance \(d=A/C\). This is also done in the script and it reports \(d\)=11.05 Å.

Comparing this value with the distance between the surface gold atoms (13.77 Å) the position of the so-called image plane is estimated to be 1.36 Å above the gold surface.

Electrostatic energy analysis

To compute the capacitance using a more general and accurate method, you will need to first calculate the electrostatic energy \(E\) as a function of the applied bias \(V\):

\[E(V)=\frac{1}{2}\int \delta v(\mathbf{r},V) \delta n(\mathbf{r}, V)d\mathbf{r}\]

where \(\delta v\) is the induced electrostatic potential and \(\delta n\) is the induced density.

The script extracts the potential and density at 1 V bias from the existing NetCDF file and produces the plots of the induced density and electrostatic potential. The script also reports two parameters:

  • the electric field in the middle of the vacuum gap \(E\)= 0.0905 V/Å;
  • the distance between the two image planes \(d\)=11.88 Å.


The position of the image planes is defined as the point where the electric field is zero, that is, where the electrostatic potential becomes flat.


The blue circles correspond to the atom positions. The main charge accumulation takes place in a very narrow region right above the surface, between the outer surface atoms and the ghost atoms.

From this, the capacitance can be extracted by fitting a parabola to the relation \(E(V) = CV^2/2\). The script calculates the electrostatic energy and fits the capacitance as shown below.



The script reports the capacitance \(C\)=5.45*10-22 Farad, in qualitative agreement with the value obtained from the simple approach using Mulliken charge analysis.

This script also reports \(d\)=11.69 Å, corresponding to an image plane at 1.04 Å above the gold surface, in very close accordance with the position of the peak of induced density.


What is the difference between the two methods?
In the Mulliken charge analysis, it is assumed that all electrons on the left surface have the same potential, and similarly for the right surface. When we convolute the electron density and the electrostatic potential we account for the fact that some part of the electron density is located in the region with a finite electric field. The electrostatic energy analysis is therefore more accurate for calculating the capacitance.

The electrostatic energy analysis also has the advantage that it does not require the researcher to separate the atoms as uniquely “left” and “right” – this is trivial for the two gold surfaces, but that may not always be the case.

Bias-dependent Capacitance

The analysis outlined above is based on the assumption that the capacitance is independent of the voltage. However, this is not always the case. A more advanced spline curve fitting allows you to investigate systems whose capacitance does depend on the bias voltage.

Provided that you have already run the scripts and in the previous section, you can now download and run the script to perform the fitting.


The data saved by the two scripts are used by the present script to accelerate the analysis.

The computed capacitance is shown below as a function of the bias. As expected, the capacitance changes in this case very little with the bias voltage.



The present example uses rather few data points to illustrate the method. In general, you should use more than five data points to get an accurate analysis.

Dielectric spacer material

The electric field in the vacuum gap is extremely high (over 900*106 V/m at 1 V bias) and would immediately lead to a breakdown discharge in the device. Replacing the vacuum gap between the capacitor plates with a dielectric material can partly solve this.

There are two different ways to do this using VNL-ATK:

  1. Insert a dielectric material, and explicitly form the two metal–semiconductor–metal interfaces.
  2. Use a cheaper and less rigorous approach, where you simply modify the dielectric constant of the region between the gold surfaces.

For demonstration purposes, we here go for the second option. There are in fact two different ways to modify the dielectric constant:

  • Implicit solvent method: Assign a dielectric constant different than 1 to the vacuum region by modifying settings in the Poisson solver.
  • Spatial region: Use the VNL Builder to place a dielectric spatial region between the surfaces.


In both cases it is necessary to use the Multigrid Poisson solver, but you can still have periodic boundary conditions along A and B.

Implicit solvent method

Follow the steps that led to calculation of the capacitance from electrostatic energy analysis, but change the Poisson solver to Multigrid. Do it like this:

  1. Copy the script to
  2. Send the new script to the Editor editor_icon and find the statement defining the Poisson solver. Change it to
# Poisson Solver Settings
device_poisson_solver = MultigridSolver(
  1. Change the output filename to throughout the script, and execute the script.
  2. Adapt the analysis scripts leading to the electrostatic energy analysis such that they use as input, and run them.

Repeating the simulations and analysis you will find the capacitance is increased by roughly a factor of 7. You can clearly see the expected discontinuity of the electric field indicated by the red circles below.


Dielectric spatial region

To try the same simulation using a dielectric spatial region, go to the Builder builder_icon and use the Miscellaneous ‣ Spatial Regions tool to place between the surfaces a box with a relative dielectric constant of 10:



The dielectric region is a rectangular box, while the unit cell is hexagonal for the present [111] fcc surface. To solve the mismatch, make sure the rectangular box is large enough in the AB plane to extend outside the unit cell. The parts outside the cell will be cut off during the calculation.

Then send the new geometry to the Script Generator script_generator_icon and repeat all the steps. Or simply open the structure in the Editor editor_icon and copy/paste the relevant parts of the scripts you already have, and change the names of the output files.

The results are illustrated below. The position of the dielectric region is indicated by the grey shading.



In both dieletric cases, the obtained capacitance is 2.7-3.8 zF, depending on the method used to compute it. The difference mainly arises from the fact that you for the spatial region approach need to make an explicit choice of how close the spatial region should be placed to the surface. Tthis is implicitly decided automatically with the implicit solvent approach, since the dielectric constant is set to 1 as close to each atom core as possible, also for ghost atoms.


  • The scripts used in the electrostatic energy analysis are self-contained and generic. They can be used to analyse other systems as well, after suitable modification of the file names.
  • The scripts used in the Mulliken charge analysis need to be modified before they can be used for other systems, to account for geometry changes.
  • These approaches are both valid only if the current through the device is very small. For example, the transmission is extremely small (it’s a capacitor – not a diode!). You can verify that the transmission is extremely small for this large vacuum gap, of the order 10-23.


What next?
Another interesting example of a nano-scale capacitor would be a once formed by Carbon Nanotube Junctions.