Forces

The Probe-Particle Model (the AFM part - ppafm) calculates all the forces that affect the probe particle (PP):

1. Forces caused by the interaction with (metallic) tip

   • The lateral forces - x and y are holding the (PP) in its equilibrium position from its lateral sides. These forces are harmonic (spring) forces. For CO they are around 0.24 N/m. Note: For the Double Probe Particle model different forces are used, see its wiki.
   • The radial force - R is holding the PP in the appropriate distance bellow the last metallic-tip atom. !!! Unlike in the original Fortran code and the paper this force is also harmonic (spring) !!! - This allows for better computational stability, than the original Lennard-Jones (L-J) model. L-J could in extreme cases lead to the loss of PP, which is physically right, but the exact constant is unknown and is extremely tip-dependent. Anyway, if the PP (standing for the flexible tip-apex) is lost in theory or experiment it causes loss of high resolution and the results are generally not usable. From our experience the usage of less precise potential is not causing any additional problems with the usage of spring constant around 20.0 or 30.0 N/m.

2. Forces caused by the interaction with the sample

   • The interaction with the sample is explained in detail below in the section Force-field models

During the relaxation, all the forces are taking into account and once a new equilibrium is found the code calculates the forces acting on the metallic forces. These tip forces can be either plotted or used for computation of the frequency shift.

Force-field models

The force field governing the interaction between the PP and the sample atoms in ppafm is divided into three components:

$$E(\vec{r}) = E_\mathrm{Pauli}(\vec{r}) + E_\mathrm{vdW}(\vec{r}) + E_\mathrm{ES}(\vec{r}).$$

The three components are the Pauli repulsion $E_\mathrm{Pauli}$, the van der Waals (vdW) attraction $E_\mathrm{vdW}$, and the electrostatic interaction (ES) $E_\mathrm{ES}$. For each of these, ppafm provides more than one model to choose from:

• Pauli:
  • Lennard-Jones
  • Full-density based model
• vdW:
  • Lennard-Jones
  • DFT-D3
• ES:
  • Point-charge
  • Hartree

While in principle it is possible to combine these freely, there are a few combinations that make the most sense. In order of increasing physical accuracy:

• Lennard-Jones without electrostatics
• Lennard-Jones + point-charge electrostatics
• Lennard-Jones + Hartree electrostatics
• Full-density based model + DFT-D3 vdW + Hartree electrostatics

The different force-field models are explained in more detail in the following, with examples of how to use them in the ppafm CLI. For GPU-accelerated calculations, see the tutorial in the Python API documentation.

Lennard-Jones

The Lennard-Jones (LJ) force field is a commonly used pair-wise force field that models the Pauli repulsion and van der Waals attraction interactions between atoms. It has the form

$$E_\mathrm{LJ}(\vec{r}_{ab}) = \varepsilon_{ab}\left[ \left( \frac{\sigma_{ab}}{r_{ab}} \right)^{12} - 2\left( \frac{\sigma_{ab}}{r_{ab}} \right)^{6} \right]$$

where $\vec{r}_{ab}$ is a vector between a pair of atoms $a$ and $b$, and $\sigma_{ab}$ and $\varepsilon_{ab}$ are parameters that define the length scale and strength of interaction, respectively, for the given pair of elements. The $r^{-12}$-term models the Pauli repulsion and the $r^{-6}$ models the vdW attraction. The parameters are tabulated for every element of the periodic table, and the pair parameters are determined by the mixing rules:

$$\sigma_{ab} = \frac{\sigma_a + \sigma_{b}}{2}$$ $$\varepsilon_{ab} = \sqrt{\varepsilon_a \varepsilon_b}$$

In ppafm the parameters values from the OPLS force field are used. You can copy and adjust this file in the local directory of calculations. The model will read it automatically.

The force is pre-calculated as a sum of forces originating from all the sample atoms, and stored in FFLJ_?.xsf or FFLJ.npz files.

In the ppafm CLI, the Lennard-Jones force field can be computed using the ppafm-generate_ljff command. Give the input geometry using the --input/-i option. For example, in the Graphene example we run

ppafm-generate-ljff -i Gr6x6N3hole.xyz

This outputs the LJ force field in three files FFLJ_{x,y,z}.xsf for each of the spatial directions. Optionally one can choose to also output energy using the --energy/-E option. Also note the --ffModel option which when set to vdW only outputs the vdW contribution without the Pauli term.

Electrostatics

ppafm supports two different ways of computing the electrostatics term, using either point charges or the Hartree potential. The point-charge model is the simpler one, where the electrostatic force is obtained from the Coulomb force between the PP and the sample atoms, which all have a point charge associated with them. The force is

$$\vec{F}_\mathrm{ES}(\vec{r}_{ab}) = \frac{1}{4 \pi \varepsilon_0} \frac{q_a q_b}{r_{ab}^3} \vec{r}_{ab},$$

where $\varepsilon_0$ is the vacuum permittivity, and $q_a$ and $q_b$ are the charges. The point charges are typically obtained from for example Mulliken or Hirshfeld charge analysis.

The point-charge electrostatic force field can be obtained in the ppafm CLI by using the ppafm-generate-elff-point-charges command. For example, in the Graphene example we use

ppafm-generate-elff-point-charges -i Gr6x6N3hole.xyz --tip s

This outputs the electrostatic force field in three files FFel_{x,y,z}.xsf for each of the spatial directions. The point charges are contained in the input xyz file. Inspection of the Gr6x6N3hole.xyz file shows how the point charges are included as an additional fifth column in the xyz file. The Extended XYZ format of ASE is also supported.

The point-charge approximation works reasonably well in many places, but for a more physically accurate model we have to consider the full continuous charge density of the sample and the tip. Using the Hartree potential of the sample $V_\mathrm{S}$ and the charge distribution of the PP $\rho_\mathrm{PP}$, the electrostatic interaction energy can be calculated as

$$E_{ES}(\vec{r}) = \int_{\vec{r}' \in \mathbb{R}^3} \rho_\mathrm{PP}(\vec{r}' + \vec{r}) V_\mathrm{S}(\vec{r}') d \vec{r}'$$

The electrostatic force field from the Hartree potential can be generated in the ppafm CLI by using the ppafm-generate-elff command. For example, in the PTCDA_HARTREE_dz2 example we run

ppafm-generate-elff -i LOCPOT.xsf --tip dz2 -f npy

using the Hartree potential in LOCPOT.xsf generated with VASP.

Usage of multipoles is available for both of methods as well. Generally speaking it is important to consider the charge and the multipole character of the PP. The character of the charge distribution depends on what the functionalization of the tip is. Typically for a CO molecule we use the quadrupole charge distribution of the dz2 orbital, but for example for a Xe-tip the simpler monopole s-orbital is a better model. In the point-charge approximation the quadrupole is created by placing three point charges in a $-q$, $2q$, $-q$ pattern equidistantly on the z-axis. The width of the charge distribution or the distance of the point charges can be controlled using the sigma parameter. The multipole moment is calculated by Q*sigma for dipole or Q*sigma*sigma for quadrupole.

Full-density based model

The Pauli repulsion modelled by the $r^{-12}$-term in the LJ-potential is not physically well justified. In order the arrive at a better approximation of the Pauli repulsion, Ellner et al. introduced the full-density based model (FDBM). It models the repulsion as an overlap integral between the sample ($\rho_\mathrm{S}$) and PP ($\rho_\mathrm{PP}$) electron densities

$$E_\mathrm{Pauli}(\vec{r}) = A \int_{\vec{r}' \in \mathbb{R}^3} \left[ \rho_\mathrm{PP}(\vec{r}' + \vec{r}) \rho_\mathrm{S}(\vec{r}') \right]^{\beta} d\vec{r}',$$

where $A$ and $\beta$ are fitting parameters. In the original paper the fitting parameters were chosen individually for each molecule by comparison to DFT calculations, ranging from 12..22 for $A$ and 1.07..1.12 for $\beta$. However, reasonable parameter values can be found also by visual inspection, noting that the observed contrast is significantly more sensitive to the $\beta$ parameter, and that adjusting both parameters in the same direction along with adjustments in the scanning height usually results in similar-looking contrast.

The Pauli force field from the FDBM can be calculated in the ppafm CLI by using the ppafm-conv-rho command. For example, in the pyridineDensOverlap example we run

ppafm-conv-rho -s sample/CHGCAR.xsf -t tip/CHGCAR.xsf -B 1.0 -E

where we provide the sample electron density using the -s option and the tip electron density using the -t option. Both densities, sample/CHGCAR.xsf and tip/CHGCAR.xsf, must be defined for the same unit cell geometry. The value of the $\beta$-exponent is set with the --Bpauli/-B option. The output force field is saved in the three files FFpauli_{x,y,z}.xsf for each of the spatial directions. Also note that we have to subsequently call the ppafm-relaxed-scan with the --noLJ option to use the Pauli from the separate files as well as --Apauli to set the value of the $A$-parameter.

DFT-D3 vdW

In their implementation of the FDBM, Ellner et al. chose to use the DFT-D3 instead of Lennard-Jones vdW interaction to accompany the Pauli interaction. Following this choice, the DFT-D3 energy/force calculation is also implemented in ppafm, specifically in the Becke-Johnson damping form:

$$E_\mathrm{vdW}(\vec{r}_\mathrm{PP}) = - \sum_{a=1}^{N_\mathrm{at}} s_6 \frac{C_6^{a\mathrm{PP}}}{r_{a\mathrm{PP}}^6 + f(R_{a\mathrm{PP}}^0)^6} + s_8 \frac{C_8^{a\mathrm{PP}}}{r_{a\mathrm{PP}}^8 + f(R_{a\mathrm{PP}}^0)^8},$$

where the sum is over all of the atoms in the sample, $r_{a\mathrm{PP}}$ is the distance between the PP and atom $a$,

$$f(R_{a\mathrm{PP}}^0) = a_1 R_{a\mathrm{PP}}^0 + a_2$$

is a damping function,

$$R_{a\mathrm{PP}}^{0}= \sqrt{\frac{C_8^{a\mathrm{PP}}}{C_6^{a\mathrm{PP}}}}$$

is a cutoff radius, and the $C_6^{a\mathrm{PP}}$ and $C_8^{a\mathrm{PP}}$ are coefficients for the energy and length scale. The parameters $s_6$, $s_8$, $a_1$, and $a_2$ are additional scaling parameters for different DFT functionals. ppafm takes the tabulated values of all of these parameters from the original Fortran implementation of DFT-D3.

The $C_{6/8}^{a\mathrm{PP}}$-coefficients are different from the LJ-coefficients in that they depend on the environment of the atoms in order to account for the changing polarizabilities of the atoms depending on their bonding configuration. In principle the PP would be included in this environment. However, due to chemically inert nature of the functionalized AFM tip, in the ppafm the PP is left out of this calculation. This makes the calculation more efficient, since the coefficients can be calculated just once for each atom in the sample, independent of the position of the PP.

The DFT-D3 vdW force field can be calculated in the ppafm CLI using the ppafm-generate-dftd3 command. For example, in the pyridineDensOverlap example we run

ppafm-generate-dftd3 -i sample/LOCPOT.xsf --df_name PBE

Here we choose to use the $s_6$, $s_8$, $a_1$, and $a_2$ parameters corresponding to the PBE DFT-functional, since that is the functional used for generating the densities and potentials for the FDBM calculation. The parameters can also be specified manually using the --df_params option. The output force field is saved in the three files FFvdW_{x,y,z}.xsf for each of the spatial directions. Also note that we have to subsequently call the ppafm-relaxed-scan with the --noLJ option to use the vdW from the separate files.

Comparison of DFT vs (original) L-J vs FDBM-D3 implementation:

Total forces including van der Waals for 6 different model systems, 4 molecules (benzene, bromomethane, carbon monoxide, and water) and 2 atoms (He, Xe). The tip is always a CO molecule, with the oxygen atom pointing towards the sample. On the molecules, the tip is placed over the C (in benzene), Br (in bromomethane), and O (in CO and H2O) atoms, respectively. The benzene and H2O molecules are flat-lying, CH3Br points with the Br atom towards the tip, CO is in an upright position wit the O atom pointing towards the tip:

Van der Waals components only for comparison. This demonstrates that the D3 model of dispersion forces as implemented in PPAFM agrees reasonably well with the dispersion forces implemented in VASP:

Forces without the van der Waals component. This shows the discrepancy between FDBM of PPAFM and VASP while also demonstrating that the problem is not with the dispersion (vdW) forces but somewhere else: