Locating Critical Points of Densities Given as Grids
Introduction
In critic2, the location of critical points (CPs) of the electron density and other fields is carried out by the AUTO keyword. AUTO works by applying Newton’s method repeatedly starting at a collection of starting points in the system, known as seeds. Newton requires knowledge of the first and second derivatives of the electron density. If the density is given analytically and the gradient and the Hessian can be calculated exactly, Newton’s method typically finds all critical points in the system and typically it only fails either because more seeds are required or because of the appearance of spurious critical points in the vacuum region. These two cases can be easily solved by appropriate use of the SEED and DISCARD subkeywords of AUTO.
If the electron density is given on a grid, analytical derivatives are not available because the density is known only at a discrete set of points, the nodes of a threedimensional uniform grid. The simplest way of trying to locate its critical points in this case is to employ threedimensional interpolation to extend the electron density to all points in the system, and to enable the calculation of derivatives. There are two problems with this naive approach:

Close to the nuclei, where the density is very steep, the interpolation methods create spurious oscillations that interfere with the CP location algorithm. The oscillations are more severe in the derivatives of the electron density, and worse the higher the derivative.

Density grids for periodic systems (solids) are often nonorthogonal. If the grid nodes used for the interpolation depend on the grid geometry, then the result of the interpolation depends on the particular basis chosen for representing the unit cell, which introduces spuriously privileged directions in the density. This is known as lattice bias.
In a 2022 article (OterodelaRoza, J. Chem. Phys. 156, 224116 (2022)) an interpolation method for the allelectron density was proposed that addresses these two problems. This method is based on polyharmonic spline interpolation combined with a smoothing function, and can be accessed by using the SMOOTHRHO option to the LOAD keyword. The number of grid nodes on which the interpolation is based (NENV) as well as the number of interpolants used in the linear combination (FDMAX) can be controlled by the user but in general touching these values is not required. Spurious critical points can still appear if the density is noisy or if the grid is not sufficiently fine.
Also, note that the SMOOTHRHO keyword is not the default interpolation method in critic2, and needs to be activated explicitly. This is because it is only suitable for interpolating allelectron densities. Other fields would require other smoothing functions (as described in the article). If you need to find critical points of other fields, please let me know.
Using VASP
The allelectron density in VASP can be obtained using the LAECHG
tag in the INCAR
file. VASP writes two files, AECCAR0
and
AECCAR2
, containing the core density and the reconstructed valence
electron density on the grid. The allelectron density is obtained as
the sum of the two fields. For example, for the covellite (CuS) system
in the package:
# Read the crystal structure either of the AECCARs
crystal covellite.AECCAR0
# Calculate the reconstructed allelectron density as the sum of core
# and valence densities
load covellite.AECCAR0
load covellite.AECCAR2
load as "$1+$2" id rhoae smoothrho
reference rhoae
# Find the critical points and write them to a vmd file
auto
cpreport covellite.vmd cell border graph
The CPREPORT keyword is used to write the set of critical points to a vmd file for visualization, together with the graph.
Using abinit
The allelectron density can be written from a PAW calculation in
abinit using the prtden 3
keyword, but this requires the separate
calculation of the core density contributions in a
core_density_atom_type1.dat
file. Instead, you can use prtden 2
and write the reconstructed valence electron density, and then employ
critic2’s core density tables to add the core contribution.
Assuming you went with prtden 2
, a _PAWDEN
file was generated that
contains the valence allelectron density. In a first step, we augment
this density with the core contribution and write the result to a new
grid file (rhoae.cube
):
## Read the crystal structure from the PAWDEN file
crystal quartz_o_PAWDEN
## Load the reconstructed valence density and augment it with the core
## contribution. The pseudopotential charge (the number of valence
## electrons) is 4 for Si and 6 for O.
load quartz_o_PAWDEN zpsp si 4 o 6 core
## Write the core+valence density to a cube file (rhoae.cube). We
## cannot interpolate from the above because it is not a pure grid
## field.
cube grid field 1 file rhoae.cube
The cube file written in this way can then be used to find the CPs:
## Read the structure from the cube file.
crystal rhoae.cube
## Load the allelectron density and activate the smoothrho
## interpolation.
load rhoae.cube smoothrho
## Find the critical points and write them to a vmd file
## for visualization
auto
cpreport quartz.vmd cell border graph
Using Quantum ESPRESSO
In Quantum ESPRESSO, the reconstructed allelectron density can be
generated using a PAW calculation and the plot_num=21
option in
pp.x
. This option writes a cube file (rhoae.cube
) containing the
density. Finding the critical points from this cube file is
straightforward. For instance, for the benzene crystal included in the
example package:
## Read the crystal structure from the cube file
crystal rhoae.cube
## Load the density with smoothrho interpolation
load rhoae.cube smoothrho
## Find the critical points and write them to a vmd file
auto
cpreport benzene.vmd border molmotif graph
Using molecular cube files (Gaussian,…)
The SMOOTHRHO interpolation can be applied to any allelectron density, periodic or nonperiodic. In particular, it can be used to reconstruct the density and find CPs in molecules whose electron density is given on a grid, for instance, as a cube file. The input is similar to the other examples, with two exceptions. First, the MOLECULE keyword is used to read the geometry. Second, a density cutoff is set so that critical points with very low density are discarded. This is used to prevent the appearance of spurious critical points in the vacuum region, where the density is flat and essentially zero. A value of \(10^{5}\) atomic units is a good choice. For the phenol dimer example included in the package:
## Read the molecular geometry from the cube file.
molecule phenol_phenol.cube
## Read the electron density from the cube file and activate the
## smoothrho interpolation.
load phenol_phenol.cube smoothrho
## Locate the critical points and discard the CPs with density ($1)
## less than 1e5 atomic units.
auto discard "$1 < 1e5"
## Write the results to a vmd file for visualization.
cpreport phenol_phenol.vmd graph
Visualization
The pictures above can be generated using vmd and the output of the CPREPORT keyword. The detailed description of how this keyword works and the visualization options is given in the manual. In short, to view the file, do:
vmd e bleh.vmd
and you can remove the atoms by writing:
label delete Atoms all
in the Tk console.
Example files package
Files: example_08_01.tar.xz.