Atomic potentials in Computem

11 minute read

Overview

This post contains my notes on how to get the correct unit of atomic potential from Computem, also some background theory of atomic potential and scattering factor, and some code snippets.

Last modified: 2022/12/27 by Chia-Hao Lee

Reference

Earl J. Kirkland, (2014), Advanced Computing in Electron Microscopy, 2nd edition
Earl J. Kirkland, (2020), Advanced Computing in Electron Microscopy, 3rd edition
- 3rd edition adds discussion of incostem, ABF, GPUs
Computem C++ source codes https://sourceforge.net/projects/computem/files/ (source.zip)
- Different versions might have different bugs, I typically use the 2018 versions but we should switch to the 2022 Aug. version soon.

How to get my projected object potential function ASAP ?

Create the crystal structure file (.dat) in the correct format (See Kirkland sec. 8.4.1)
Run atompot, use enough real and reciprocal space sampling
1. Large px number (small px size) in real space, and large super cell
Load the correct page of output.tif from Computem
1. In Python, use img = imageio.mimread('output.tif')[1]
2. In Matlab, use img = imread('output.tif', 2);
  - Computem output.tif has 3 pages (see sec. 8.3.1), 1st page is uint8, 2nd page is float32, 3rd page is just some parameters (max and min px values, px sizes). Note that ImageJ and tifffile in Python will only read the 1st page of the output, which is normalized to [0, 255] by Computem for display. The actual quantitative values are always stored in the 2nd page.
Multiply the img with a unit conversion factor 47.94 to get the projected object potential function in unit of volt - Ang.
1. projected_object_potential = img * 47.94 # Unit: volt-Ang
2. conversion factor = $\frac{h^2}{2\pi m_0 e}$= 0.4794*10^-18 volt-m^2 = 47.94 volt-Ang^2
  1. h: Planck’s constant, 6.6261*10^-34 J-s
  2. m0: electron rest mass, 9.1094*10^-31 kg
  3. e: Fundamental charge of electron, 1.6022*10^-19 C
3. The original unit of output.tif is 1/Ang, see the following section.

Notes:

The projected object potential function is kV independent. In order to get the object transmission function $t(x)$, one must multiply it by the interaction parameter $\sigma$ to convert it into phase shift. Therefore, $t(x) = exp[i \sigma v_z(x)]$, where $v_{z} (x)$ is the projected object potential.

Why the `atompot` output is in such (1/Ang) unit ?

atompot.cpp, line 56

The units of the output are such that when multiplied by Lambda*(M/M0) the result is the phase shift of the transmitted electron wave (i.e. a simple sum of the unadulterated Born approx. scattering amplitudes).

In actual Computem mulslice simulation, it takes the atompot output (in 1/Ang) and multiply by a scaling factor scale = wavelen * mm0 (Ang) to get the phase shift in radians. In other words, the kV dependence is included in the scaling factor.
```
  scale = wavlen * mm0;
  for( iy=0; iy<ny; iy++) {
      for( ix=0; ix<nx; ix++) {
          vz= scale * myFile(ix, iy);
          trans[ilayer].re(ix,iy) = (float) cos(vz);
          trans[ilayer].im(ix,iy) = (float) sin(vz);
      }
  }
```
- wavlen # Relativistic wavelength in angstrom
- mm0 = 1.0F + v0/511.0F; = $m/m0 = \gamma$ # dimensionless, $\gamma = 1.39$ for 200kV electron
  - v0 is in kV
- Despite it’s called vz after scaling in the source code, it’s actually the phase shift, or $\sigma v_z(x)$ following our previous notation.

Quantification test with a single Au atom

A single Au is simulated with 100 x 100 Ang supercell at [0.5, 0.5] using 4096 x 4096 pixels, px size = 2.44 pm
The atompot output has px range [-0.003886, 128.266], the px value @ 4 px away (~ 0.1 Ang) is 30.3904.
After scaling with 47.94 volt-Ang^2, the projected object potential of Au at 0.1 Ang radius is 30.3904 * 47.94 = 1457 volt-Ang, this value mathes the value (1.45 kV-Ang, see Kirkland Fig. 5.5) very well so we know the scaling is done correctly.

Overall STEM image/CBED simulation process

Input crystal structure (cell size, atom position, occupancy, Debye-Waller factors)
Create object slices if needed
In each slice, the projected object potential is calculated by the superposition of individual projected atomic potentials (Unit: V - Ang)
Calculate the object transmission function with $t(x) = exp[i \sigma v_z(x)]$
Calculate the incoming illumination wavefunction (focused probe, essentially lots of plane waves) $\psi_i(x)$
Multiply the illumination wavefunction by the objject transmission function to get the transmitted wavefunction $\psi_t(x)$ = $\psi_i(x) \times t(x)$ , which is the exit wave.
1. For focused probe, it has to be placed on each probe position w.r.t. the transmission function $\psi_i(x) \rightarrow \psi_i(x, x_p)$
Fourier transform the transmitted wavefunction (exit wave) to propagate from specimen plane to back focal plane, where diffraction patterns locate.

For CBED patterns, the square modulus of $\psi_t(k)$, that is, $

\psi_t(k)

^2$ will be your CBED pattern at the specific probe position.

For STEM images, integrate (sum) the square modulus of the wavefunction, $\int_{k_i}^{k_f}

\psi_t(k)

^2 dk$ at the k range specified by your detector will give you the pixel intensity at that probe position.

Repeat Step 6- 9 with a shifted probe for every probe positions

Notes on the theory of atomic potential and scattering factor

Atomic potential or charge distributions can be calculated from first principles using relativistic Hartree-Fock theory. Although we can tabulate all the atomic potentials, it’s not practical to use the tabulated numbers for simulation, so people use the parametrization approach instead.
The atom size is actually relatively constant with Z because the more electrons, the stronger Z will attracts it, keeping the size relatively constant to 1 Ang. Note that Fig 5.3 is only considering the electron charge density, while the atomic potential is mostly determined by the Z, or nuclei.
See Fig. 5.3
The atomic potential V(r) can be reasonably well approximated by 12 parameters (eqn. 5.9)
- a0: Bohr radius
- ai, bi, ci, di: Fitting parameters obtained by fitting the electron scattering factor with Hatree-Fock results using 3 Lorentzians and 3 Gaussians (Appendix C. 4)
  \[Va(x,y,z) = 2 \pi^2 a_0 e \sum_{i=1}^3 \frac{a_i}{r}exp(-2 \pi r \sqrt{b_i} + 2\pi^{5/2} a0 e \sum_{i=1}^3 c_i d_i^{-3/2} exp(-\pi^2 r^2 / d_i)\]
This form of $V_a(x,y,z)$ is analytically derived from the parameterized scattering factor
When integrating along Z direction, the projected atomic potential is: (eqn. 5.10)
- K0: modified Bessel function of zeroth order
- This is the form we actually used in the simulation process

\[v_z(x,y) = \int_{-\infin}^{+\infin}V_a(x,y,z)dz = 4 \pi^2 a_0 e \sum_{i=1}^3 {a_i}K_0(2 \pi r \sqrt{b_i}) + 2\pi^2 a_0 e \sum_{i=1}^3 \frac{c_i}{d_i} exp(-\pi^2 r^2 / d_i)\]

3D atomic potential (kV). See Fig. 5.4
Projected (integrated) atomic potential (V - Ang). See Fig. 5.5
Multiplying the projected potential with interaction parameter $\sigma$ you’ll get the phase shift. For example, the projected potential of Au at 0.1 Ang is 1.45 kV-Ang, and the $\sigma$ = 0.92/(kV-Ang), so the phase shift is 1.45*0.92~1.34 radians. A single Au is definitely not a weak phase object for 100 kV electron.
One might define rms (root-mean-square) radius as the atom size (See eqn. 5.11, 5.12)
The rms radisu of isolated single atom (See Fig. 5.7)
Practically, people use the 1st Born approximated atomic scattering factor to calculate the atomic potentials because:
1. It’s the Fourier transform of the atomic potential
2. It’s convenient to calculate projected potential $V_z(x, y)$ with Fourier slice/projection theorem
3. It’s independent to the acceleration voltage, unlike the complex Moliere version
4. It is well parametrized and the functional form can be analytically Fourier transformed into the atomic potential (eqn. 5.17)
  - ai, bi, ci, di: Fitting parameters obtained by fitting the electron scattering factor with Hatree-Fock results using 3 Lorentzians and 3 Gaussians (Appendix C. 4), a total of 12 parameters per atom.
  \[f_e(q) = \sum_{i=1}^3 \frac{a_i}{q^2+b_i} + \sum_{i=1}^3 c_i exp(-d_i q^2)\]
See Fig. 5.8

Under 1st Born approximation, the scattered wave amplitude (scattering factor) is the Fourier transform of the scatterer potential, which is atomic potential in this case.

The atomic potential is a real-valued function, therefore the 1st Born approximated scattering factor will also be a real-valued function

A common usage of scattering factor is to used it to constructe the complex structure factor $F_{hkl}$, where the kinematic diffraction intensity $I_{hkl} =

^2$

This formalism was developed for X-ray and their scattering factors are pretty well-described by kinematic scattering.
However, for electron diffraction, the scattered wave amplitude is simply not well-described by the 1st Born approximated scattering factor, because electrons are strongly scattered by even an single atom.

Sec. 5.2.4 Scattering Factor

The ﬁrst Born approximation is totally inadequate for directly calculating electron scattering in the electron microscope image.

The scattering factor should generally be a complex value for electron scatterings to preserve the number of electrons (optical theorem), however, in the first Born approximation, it’s been simply formulated as the Fourier transform of the atomic potential (real valued), so the first Born approximated scattering factor must also be a real valued function, this is essentially setting the phase of scattered wave to 0, and would fail to keep the number of particles. See Fig. 5.9

So one must use the atomic potential to calculate the object transmission function to better describe the electron scattering process. The electron scattering factor is introduced only to explain how we get the parametrized form of atomic potential in the simulation.

Proper calculation

1st Born approximated Scattering factor → Atomic potential → transmission function → transmitted wave function $\psi_t(x)$ → diffraction wave function $\psi_t(k)$→ diffraction pattern $

\psi_t(k)

^2$

Improper calculation (only works for kinematic scattering with weak phase object)

1st Born approximated Scattering factor → Structure factor $F$→ diffraction pattern $

^2$

Notes on Computem implementations of atomic potentials

Appatently Computem implements multiple approaches for the object potential function

scamp + featom for atompot, stemslic, mulslice and they’re all obsolete and will be discontinued soon (Stops updating since 2017)
vzatomLUT for autoslic, autostem
feMoliere + vzatom for incostem

Computem 2018 Oct. version

atompot: 2017 Aug. 5
- scamp + featom in Fourier space and then iFFT (This is probably why it’s been superceeded and obsolete)
- Warning: this program is obsolete and may be discontinued soon
- calculate projected atomic potentials (to use in multislice)
autoslic: 2018 Aug. 30
- vzatomLUT
- perform CTEM multislice with automatic slicing and FFTW and multithreaded using openMP
autostem: 2018 Oct. 4
- vzatomLUT
- Calculate STEM images using FFTW and multithreaded using openMP
incostem: 2016 Jul. 6
- feMoilere + vzatom
- calculate ADF-STEM images for very thin specimens using the incoherent image model
mulslice: 2018 May 4
- Read object potential from atompot
- Warning: this program is obsolete and may be discontinued soon
- perform traditional CTEM multislice calculation using FFTW with 4 thread(s)
stemslic: 2017 Aug. 10
- Read object potential from atompot
- Warning: this program is obsolete and may be discontinued soon
- perform STEM multislice using FFTW with 4 thread(s)

atompot approach: scamp + featom (obsolete but keep if for reference)

include slicelib.hpp
- slicelib.hpp and slicelib.cpp contains lots of relevant functions. hpp files contain declarations, while . cpp files contain implementation of what is declared in the . hpp file
- featom() : return scattering factor for a given atom
- ReadfeTable() : read fe scattering factor table
- sigma() : return the interaction parameter
- vatom() : return real space atomic potential (NOT projected)
- vzatom() : return real space projected atomic potential
- vzatomLUT() : same as vzatom() but with a look-up-table (faster)
- wavelength() : return electron wavelength in A for given keV

main workflow:

Load atom positions from .dat file
Calculate featom for each atom in k space
Call scamp to get complex scattering amplitude scampr, scampi. Note that scamp takes the real space position of the atom to generate the corresponding scampr and scampi
Accumulatively fill the complex diffraction wave function array with scamp and scampi for each atom
1. Think of each atom is accumulatively contributing to the final diffraction pattern

iFFT to get the real space object potential function

 /*  calculate scattering amplitude in upper half plane (FFTW)
    NOTE zero freg is in the bottom left corner and
      expands into all other corners - not in the center
      this is required for FFT, high freq is in the center

     - add extra complex conjugation to correct for the opposite
     sign convention used by FFTW

 */
             ncoeff = 0;
             for( iy=0; iy<=iymid; iy++) {
                 ky2 = ky[iy]*ky[iy] * ry2;
                 for( ix=0; ix<nx; ix++) {
                     k2 = kx[ix]*kx[ix] * rx2 + ky2;
                     if( k2 <= k2max) {
                         fe = scale * featom( iz, k2 );
                         scamp( kx[ix], ky[iy], &scampr, &scampi ) ;
                         cpix.re(ix,iy) += (float) (scampr * fe);  // real  === bad here ???
                         cpix.im(ix,iy) += (float) (-scampi * fe); // imag
                         ncoeff++;
                     }

                 }  /* end for(ix */
             } /* end for(iy */

             total2 = total2 + total1;
             //goto More;

         }  /* end top if( iz>= 1) */

     /*  end loop over types */
     } while( iz>= 1 );  // end do{  Readline....

     cout << "\n   for a grand total of " << total2 << " atoms" << endl;
     fp.close( );

     cpix.ifft();   /*    inverse fft */

     //  output results and find min and max to echo

scamp(): Complex specimen scattering amplitude, scampr, scampi = real and imag. parts of scattering amplitude returned by this routine

  /*--------------------- scamp() -----------------------*/
  /*
    Complex specimen scattering amplitude

    kx,ky = real scattering vectors
    scampr, scampi = real and imag. parts of scattering amplitude
             returned by this routine

    global variables used:
    x[natom], y[natom] = real array with atomic coord.
    occ[natom] = real array with occupations at each x,y
    natom      = integer number of atoms (all with same Z)
    nsx, nsy = integer number of cells to replicate in x,y

  */
  void scamp( float kx, float ky, double *scampr, double *scampi )
  {
      int i, ixc, iyc;
      double pi2, w1, w;
      double  scalex, scaley;

      double phasexr, phasexi, phaseyr, phaseyi;

      pi2 = 8.0 * atan( 1.0 );

      scalex = pi2*kx / ((double)nsx);
      scaley = pi2*ky / ((double)nsy);

      *scampr = *scampi = 0.0;

  /*  use a trick to sum nsx + nsy terms
      instead of nsx * nsy
      - cos()/sin() could be done recursive for speed
  */
      for( i=0; i<natom; i++) {

          phasexr = phasexi = 0.0;
          w1 = pi2 * kx * x[i];
          for( ixc=0; ixc<nsx; ixc++) {
              w = w1 + ixc*scalex;
              phasexr += cos(w);
              phasexi += sin(w);
          }

          phaseyr = phaseyi = 0.0;
          w1 = pi2 * ky * y[i];
          for( iyc=0; iyc<nsy; iyc++) {
              w = w1 + iyc*scaley;
              phaseyr += cos(w);
              phaseyi += sin(w);
          }

          *scampr += ( phasexr*phaseyr - phasexi*phaseyi ) * occ[i];
          *scampi += ( phasexr*phaseyi + phasexi*phaseyr ) * occ[i];

      } /* end for(i=0... */

      return;

   }  // end scamp()

Share on

Twitter Facebook LinkedIn

Chia-Hao Lee

Atomic potentials in Computem

Overview

How to get my projected object potential function ASAP ?

Why the `atompot` output is in such (1/Ang) unit ?

Quantification test with a single Au atom

Overall STEM image/CBED simulation process

Notes on the theory of atomic potential and scattering factor

Notes on Computem implementations of atomic potentials

Computem 2018 Oct. version

atompot approach: scamp + featom (obsolete but keep if for reference)

Share on

Leave a comment

You may also enjoy

Campus Cluster Usage Guide

7 Ways to Convert Velox emd Files

Chia-Hao Lee

Overview

How to get my projected object potential function ASAP ?

Why the atompot output is in such (1/Ang) unit ?

Quantification test with a single Au atom

Overall STEM image/CBED simulation process

Notes on the theory of atomic potential and scattering factor

Notes on Computem implementations of atomic potentials

Computem 2018 Oct. version

atompot approach: scamp + featom (obsolete but keep if for reference)

Share on

Leave a comment

You may also enjoy

Campus Cluster Usage Guide

7 Ways to Convert Velox emd Files

Why the `atompot` output is in such (1/Ang) unit ?