Properties at the nucleus¶
Observables near the atomic nucleus.¶
The purpose of this tutorial is to show how to compute several observables of interest in Mössbauer, NMR, and NQR spectroscopy, namely:
- the electric field gradient,
- the isomer shift,
- chemical shielding
This tutorial should take about 2 hours.
Note
Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.
In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:
export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/ # Do not change this line: path to tests dir
export ABI_PSPDIR=$ABI_TESTS/Pspdir/ # Do not change this line: path to pseudos dir
Examples in this tutorial use these shell variables: copy and paste
the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively,
source the set_abienv.sh script located in the ~abinit directory:
source ~abinit/set_abienv.sh
The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.
To execute the tutorials, create a working directory (Work*) and
copy there the input files of the lesson.
Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.
Electric field gradient¶
Various spectroscopies, including nuclear magnetic resonance and nuclear quadrupole resonance (NMR and NQR), as well as Mössbauer spectroscopy, show spectral features arising from the electric field gradient at the nuclear sites. Note that the electric field gradient (EFG) considered here arises from the distribution of charge within the solid, not due to any external electric fields.
The way that the EFG is observed in spectroscopic experiments is through its coupling to the nuclear electric quadrupole moment. The physics of this coupling is described in various texts, for example [Slichter1978]. Abinit computes the field gradient at each site, and then reports the gradient and its coupling based on input values of the nuclear quadrupole moments.
The electric field and its gradient at each nuclear site arises from the distribution of charge, both electronic and ionic, in the solid. The gradient especially is quite sensitive to the details of the distribution at short range, and so it is necessary to use the PAW formalism to compute the gradient accurately. For charge density \(n({\mathbf r})\), the potential \(V\) is given by $$ V({\mathbf r})=\int \frac{n({\mathbf r’})}{ |\mathbf{r}-\mathbf{r’}| } d{\mathbf r’} $$ and the electric field gradient is $$ V_{ij} = -\frac{\partial^2}{\partial x_i\partial x_j}V({\mathbf r}). $$ The gradient is computed at each nuclear site, for each source of charge arising from the PAW decomposition (see the tutorial PAW1 ). This is done in the code as follows [Profeta2003],[Zwanziger2008]:
- Valence space described by planewaves: expression for gradient is Fourier-transformed at each nuclear site.
- Ion cores: gradient is computed by an Ewald sum method
- On-site PAW contributions: moments of densities are integrated in real space around each atom, weighted by the gradient operator
The code reports each contribution separately if requested.
The electric field gradient computation is performed at the end of a ground-state calculation, and takes almost no additional time. The tutorial file is for stishovite, a polymorph of SiO\(_2\). In addition to typical ground state variables, only two additional variables are added:
nucefg 2
quadmom 0.0 -0.02558
# Computation of the electric field gradient at the atomic sites in stishovite # Stishovite is a polymorph of SiO2 and has the rutile structure. # Computation of EFG requires the use of PAW, Norm-conserving pseudopotentials are # not appropriate for this calculation. # here are the two variables needed to trigger computation of the field gradient nucefg 2 # NUClear properties EFG, value of 2 means to also print its components quadmom 0.0 -0.02558 # quadrupole moments of the nuclei (Si-29, O-17 here) # so that couplings as measured in NMR or NQR are also # reported, in MHz frequency units #Definition of the unit cell # stishovite has a primitive tetragonal cell, so two sides are equal, # the third different, and the angles all 90 degrees acell 8.0298189198E+00 8.0298189198E+00 5.1240093322E+00 Bohr rprim 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 #Definition of the atom types and pseudopotentials ntypat 2 # two types of atoms znucl 14 8 # atomic number of atoms, will be cross checked against pseudopotential files # here Silicon (14) and Oxygen (8). Pseudopotentials must be listed in this order. # this order defines type 1 as silicon, type 2 as oxygen pp_dirpath "$ABI_PSPDIR" # pseudopotential directory pseudos "Si.GGA-PBE-rpaw-1.55.abinit, O.GGA-PBE-rpaw-1.45.abinit" #Definition of the atoms natom 6 # 6 atoms in the unit cell typat 1 1 2 2 2 2 # atoms listed as silicons, then oxygens # next are the 6 atomic positions, in order of typat (silicon, then oxygen) # xred gives them in units of the cell translation vectors, as is commonly # done in crystallography xred 3.1301524738E-20 1.2463357829E-20 -8.3403774579E-17 5.0000000000E-01 5.0000000000E-01 5.0000000000E-01 3.0660149474E-01 3.0660149474E-01 -8.3403774579E-17 8.0660149474E-01 1.9339850526E-01 5.0000000000E-01 1.9339850526E-01 8.0660149474E-01 5.0000000000E-01 6.9339850526E-01 6.9339850526E-01 -8.3403774579E-17 #Numerical parameters of the calculation : planewave basis set and k point grid ecut 15 # plane wave cutoff energy, always need to check convergence of ecut in real runs pawecutdg 20 # cutoff energy for PAW fine grid ngkpt 8 8 6 # Parameters of Monkhorst-Pack kpoint grid. Convergence should be checked for # this parameter as well. #Parameters for the SCF procedure nstep 10 tolvrs 1.0D-18 # suppress printing of density, wavefunctions, etc for this tutorial prtden 0 prtwf 0 prteig 0 ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tnuc_1.abo, tolnlines= 12, tolabs= 1.2e-8 , tolrel= 3.0e-3 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = J. Zwanziger #%% keywords = PAW #%% description = #%% topics = EFG, SmartSymm #%%<END TEST_INFO>
The first variable instructs Abinit to compute and print the electric field gradient, and the second gives the quadrupole moments of the nuclei, in barns, one for each type of atom. A standard source for quadrupole moments is [Pyykko2008]. Here we are considering silicon and oxygen, and in particular Si-29, which has zero quadrupole moment, and O-17, the only stable isotope of oxygen with a non-zero quadrupole moment.
After running the file tnuc_1.abi through Abinit, you can find the following near the end of the output file:
Electric Field Gradient Calculation
atom : 1 typat : 1
Nuclear quad. mom. (barns) : 1.0000 Cq (MHz) : 0.0000 eta : 0.0000
efg eigval (au) : -0.152323 ; (V/m^2) : -1.48017693E+21
- eigvec : -0.000000 -0.000000 -1.000000
efg eigval (au) : -0.054274 ; (V/m^2) : -5.27401886E+20
- eigvec : 0.707107 -0.707107 0.000000
efg eigval (au) : 0.206597 ; (V/m^2) : 2.00757882E+21
- eigvec : 0.707107 0.707107 -0.000000
total efg : 0.076161 0.130436 -0.000000
total efg : 0.130436 0.076161 -0.000000
total efg : -0.000000 -0.000000 -0.152323
efg_el : 0.095557 0.004024 -0.000000
efg_el : 0.004024 0.095557 -0.000000
efg_el : -0.000000 -0.000000 -0.191114
efg_ion : -0.099183 0.005966 0.000000
efg_ion : 0.005966 -0.099183 0.000000
efg_ion : 0.000000 0.000000 0.198365
efg_paw : 0.079787 0.120445 0.000000
efg_paw : 0.120445 0.079787 0.000000
efg_paw : 0.000000 0.000000 -0.159574
This fragment gives the gradient at the first atom, which was silicon. Note that the gradient is not zero, but the coupling is—that’s because the quadrupole moment of Si-29 is zero, so although there’s a gradient there’s nothing in the nucleus for it to couple to.
Atom 3 is an oxygen atom, and its entry in the output is:
atom : 3 typat : 2
Nuclear quad. mom. (barns) : -0.0256 Cq (MHz) : 6.6150 eta : 0.1403
efg eigval (au) : -1.100599 ; (V/m^2) : -1.06949233E+22
- eigvec : 0.707107 -0.707107 0.000000
efg eigval (au) : 0.473085 ; (V/m^2) : 4.59714112E+21
- eigvec : 0.000000 0.000000 -1.000000
efg eigval (au) : 0.627514 ; (V/m^2) : 6.09778216E+21
- eigvec : 0.707107 0.707107 0.000000
total efg : -0.236543 0.864057 0.000000
total efg : 0.864057 -0.236543 0.000000
total efg : 0.000000 0.000000 0.473085
efg_el : -0.036290 -0.075078 0.000000
efg_el : -0.075078 -0.036290 0.000000
efg_el : 0.000000 0.000000 0.072579
efg_ion : -0.016807 0.291185 -0.000000
efg_ion : 0.291185 -0.016807 -0.000000
efg_ion : -0.000000 -0.000000 0.033615
efg_paw : -0.183446 0.647950 0.000000
efg_paw : 0.647950 -0.183446 0.000000
efg_paw : 0.000000 0.000000 0.366891
Now we see the electric field gradient coupling, in frequency units, along with the asymmetry of the coupling tensor, and, finally, the three contributions to the total. Note that the valence part, efg_el, is small, while the ionic part and the on-site PAW part are larger. In fact, the PAW part is largest; this is why these calculations give very poor results with norm-conserving pseudopotentials, and need the full accuracy of PAW to capture the behavior near the nucleus. Experimentally, the nuclear quadrupole coupling for O-17 in stishovite is reported as \(6.5\pm 0.1\) MHz, with asymmetry \(0.125\pm 0.05\) [Xianyuxue1994]. It is not uncommon for PAW-based EFG calculations to give coupling values a few percent too large; often this can be improved by using PAW datasets with smaller PAW radii, at the expense of more expensive calculations [Zwanziger2016].
Fermi contact interaction¶
The Fermi contact interaction arises from overlap of the electronic wavefunctions with the atomic nucleus, and is an observable for example in Mössbauer spectroscopy [Greenwood1971]. In Mössbauer spectra, the isomer shift \(\delta\) is expressed in (SI) velocity units as $$ \delta = \frac{2\pi}{3}\frac{c}{E_\gamma}\frac{Z e^2}{ 4\pi\epsilon_0} ( |\Psi (R)_A|^2 - |\Psi (R)_S|^2 )\Delta\langle r^2\rangle $$ where \(\Psi(R)\) is the electronic wavefunction at nuclear site \(R\), for the absorber (A) and source (S) respectively; \(c\) is the speed of light, \(E_\gamma\) is the nuclear transition energy, and \(Z\) the atomic number; and \(\Delta\langle r^2\rangle\) the change in the nuclear size squared. All these quantities are assumed known in the Mössbauer spectrum of interest, except \(|\Psi(R)|^2\), the Fermi contact term.
Abinit computes the Fermi contact term in the PAW formalism by using as observable \(\delta(R)\), that is, the Dirac delta function at the nuclear site [Zwanziger2009]. Like the electric field gradient computation, the Fermi contact calculation is performed at the end of a ground- state calculation, and takes almost no time. There is a tutorial file for SnO\(_2\), which, like stishovite studied above, has the rutile structure. In addition to typical ground state variables, only one additional variable is needed:
nucfc 1
# Computation of the Fermi contact term at the atomic sites in SnO2 # SnO2 has the rutile structure nucfc 1 # requests Fermi contact at each nuclear site #Definition of the unit cell acell 3*8.9535222618 # each primitive vector in rprim is scaled # by its value in acell, to give internally # cell vectors dimensioned as length in atomic units # SnO2 has a primitive tetragonal cell, so two sides are equal, # the third different, and the angles all 90 degrees rprim 1.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 1.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.672541156606163 #Definition of the atom types and pseudopotentials ntypat 2 # two types of atoms znucl 50 8 # atomic number of atoms, will be cross checked against pseudopotential files # here Tin (50) and Oxygen (8). Pseudopotentials must be listed in this order. # this order defines type 1 as tin, type 2 as oxygen pp_dirpath "$ABI_PSPDIR" # pseudopotential directory pseudos "Psdj_paw_pbe_std/Sn-sp.xml, Psdj_paw_pbe_std/O.xml" #Definition of the atoms natom 6 # 6 atoms in the unit cell typat 1 1 2 2 2 2 # atoms listed as silicons, then oxygens # next are the 6 atomic positions, in order of typat (silicon, then oxygen) # xred gives them in units of the cell translation vectors, as is commonly # done in crystallography xred 0.000000000000000 0.000000000000000 0.000000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.307100000000000 0.307100000000000 0.000000000000000 0.807100000000000 0.192900000000000 0.500000000000000 0.192900000000000 0.807100000000000 0.500000000000000 0.692900000000000 0.692900000000000 0.000000000000000 #Numerical parameters of the calculation : planewave basis set and k point grid ecut 15 # plane wave cutoff energy, always need to check convergence of ecut in real runs pawecutdg 20 # cutoff energy for PAW fine grid ngkpt 8 8 6 # Parameters of Monkhorst-Pack kpoint grid. Convergence should be checked for # this parameter as well. #Parameters for the SCF procedure nstep 10 tolvrs 1.0D-18 # suppress printing of density, wavefunctions, etc for this tutorial prtden 0 prtwf 0 prteig 0 ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tnuc_2.abo, tolnlines= 12, tolabs= 4.1e-09, tolrel= 3.000e-03 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = J. Zwanziger #%% keywords = PAW #%% description = #%% topics = EFG, SmartSymm #%%<END TEST_INFO>
After running this file, inspect the output and look for the phrase “Fermi-contact Term Calculation”. There you’ll find the FC output for each atom; in this case, the Sn atoms, typat 1, yield a contact term of 71.6428 (density in atomic units, \(a^{-3}_0\)).
To interpret Mössbauer spectra you need really both a source and an absorber; in the tutorial we provide also a file for \(\alpha\)-Sn (grey tin, which is non-metallic).
# Computation of the Fermi contact term at the atomic sites in \alpha-Sn # Nonmetallic form of tin nucfc 1 # requests Fermi contact at each nuclear site #Definition of the unit cell acell 3*12.262810608 # each primitive vector in rprim is scaled # by its value in acell, to give internally # cell vectors dimensioned as length in atomic units # a-Sn has a cubic conventional cell, so the primitive cell # has 60 deg angles rprim 0.000000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.000000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.000000000000000 #Definition of the atom types and pseudopotentials ntypat 1 # one types of atom znucl 50 # atomic number of atoms, will be cross checked against pseudopotential files # here Tin (50) is the atom type pp_dirpath "$ABI_PSPDIR" # pseudopotential directory pseudos "Psdj_paw_pbe_std/Sn-sp.xml" #Definition of the atoms natom 2 # 2 atoms in the primitive unit cell typat 1 1 # atoms listed as Sn Sn # next are the atomic positions # xred gives them in units of the cell translation vectors, as is commonly # done in crystallography xred 0.000000000000000 0.000000000000000 0.000000000000000 0.250000000000000 0.250000000000000 0.250000000000000 #Numerical parameters of the calculation : planewave basis set and k point grid ecut 15 # plane wave cutoff energy, always need to check convergence of ecut in real runs pawecutdg 20 # cutoff energy for PAW fine grid ngkpt 4 4 4 # Parameters of Monkhorst-Pack kpoint grid. Convergence should be checked for # this parameter as well. nshiftk 4 # the default M-P shift of 1/2, 1/2, 1/2 breaks the symmetry of the # reciprocal lattice but these shifts preserve it. See the abinit # documentation for the shiftk input variable shiftk 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 #Parameters for the SCF procedure nstep 10 tolvrs 1.0D-18 # suppress printing of density, wavefunctions, etc for this tutorial prtden 0 prtwf 0 prteig 0 ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tnuc_3.abo, tolnlines= 3, tolabs= 6.000e-10, tolrel= 3.000e-10 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = J. Zwanziger #%% keywords = PAW #%% description = #%% topics = EFG, SmartSymm #%%<END TEST_INFO>
If you run this file, you should find a contact term of 102.0748.
To check your results, you can use experimental data for the isomer shift \(\delta\) for known compounds to compute \(\Delta\langle r^2\rangle\) in the above equation (see [Zwanziger2009]). Using our results above together with standard tin Mössbauer parameters of \(E_\gamma = 23.875\) keV and an experimental shift of 2.2 mm/sec for \(\alpha\)-Sn relative to SnO\(_2\), we find \(\Delta\langle r^2\rangle = 5.67\times 10^{-3}\mathrm{fm}^2\), in decent agreement with other calculations of 6–7\(\times 10^{-3}\mathrm{fm}^2\) [Svane1987], [Svane1997].
Chemical shielding¶
The chemical shielding is routinely measured in NMR spectroscopy, and arises from the orbital currents induced by an external magnetic field. From an energetic point of view, it can be understood as the joint response between an external magnetic field, and the nuclear magnetic dipole moment at an atomic site, that is, $$ \sigma_{ij} = \frac{\partial^2 E}{\partial m_i\partial B_j} $$ for shielding \(\sigma\) and magnetic dipole \(m\). The total energy could then be viewed as $$ E = E^0 - m\cdot B + m_i\sigma_{ij}B_j, $$ where the first term is the unperturbed energy, the second is the Zeeman interaction, and the third involves the shielding. The effect can thus be viewed either as a bare dipole interacting with a shielded magnetic field (the traditional view), or conversely [Thonhauser2009] [Ceresoli2010], a shielded dipole interacting with the bare magnetic field. In Abinit we adopt the second view point, and so compute the perturbed energy \(\partial E/\partial B\) in the presence of a nuclear dipole moment. The necessary expressions are found in [Zwanziger2023], see also [Gonze2011a] [Ceresoli2006], and are quite complex; in highly abbreviated form we evaluate $$ \partial E/\partial B_\alpha = -M_{\alpha} = \frac{i}{2}\epsilon_{\alpha\beta\gamma}\sum_n^{\mathrm{occ}} \int d\mathbf{k} \langle \partial_{k_\beta} u_{n\mathbf{k}}|(H_{\mathbf{k}}+E_{\mathbf{k}})|\partial_{k_\gamma} u_{n\mathbf{k}}\rangle $$ for the induced magnetic moment \(M\) in direction \(\alpha\). Notice the implied summation over the antisymmetric unit tensor \(\epsilon_{\alpha\beta\gamma}\); the structure is hence that of a cross product, which yields the circulation induced by the Hamiltonian. In the full PAW treatment implemented in Abinit there are a number of other related terms as well. Note in this expression the appearance of the derivatives of the wavefunctions; these are obtained from the DDK perturbation in Abinit, as was examined in tutorial DFPT1 ).
Carrying out the calculation in Abinit is a two-step process: first, the ground state wavefunctions must be calculated in the presence of a nuclear magnetic dipole moment, on the atom of interest in the direction of interest; then these wavefunctions are used to compute the DDK wavefunctions under the same conditions, followed by assembly of the terms making up \(\sigma_{ij}\). Here is a sample input file:
# example of orbmag calculation of chemical shielding. system is # Ne atom at center of an empty box. # # for an atom, use a single k pt at Gamma kptopt 0 nkpt 1 kpt 3*0.0 # for a solid, one would use kptopt 3 and ngkpt 8 8 8 (say) nband 4 occopt 1 # dipole of strength 1, in direction x, on atom 1 nucdipmom 1.0 0.0 0.0 # no symmetries--in general a single dipole in the unit cell # destroys TR symmetry and most if not all spatial symmetries nsym 1 symrel 1 0 0 0 1 0 0 0 1 symmorphi 0 istwfk *1 ndtset 2 # two datasets, first is ground state tolvrs1 1.0D-18 nstep1 25 # added to force stop on automatic tests, in general should be larger # second is DDK perturbation, in all 3 directions, # with orbmag 2 (or 1) to trigger orbmag calculation at # the end of the DDK rfddk2 1 rfdir2 1 1 1 iscf2 -3 getwfk2 -1 orbmag2 2 tolwfr2 1.0D-20 nstep2 10 # lamb shielding for core electrons of this PAW set lambsig 341.0D-6 # ecut 30 gives better convergence ecut 15 pawecutdg 30 # some flags that orbmag and nucdipmom require optforces 0 optstress 0 pawcpxocc 2 usexcnhat 0 paral_atom 0 paral_kgb 0 # Structural parameters # cubic box 20 bohr/side acell 3*20.0 natom 1 ntypat 1 typat 1 znucl 10 # Ne atom at center of box xred 3*1/2 pp_dirpath="$ABI_PSPDIR/Psdj_paw_pbe_std/" pseudos="Ne.xml" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tnuc_4.abo, tolnlines= 65, tolabs= 1.1e-5, tolrel= 1.01e0 #%% [paral_info] #%% max_nprocs = 1 #%% [extra_info] #%% authors = J. Zwanziger #%% keywords = PAW #%% description = #%% topics = NMR, DFPT #%%<END TEST_INFO>
This file is for an isolated neon atom, that is, a neon atom at the center of a large box. There are two chained calculations: first, the ground state, and then, the DDK perturbation (rfelfd = 2 or equivalently rfddk = 1). In both cases, a nuclear magnetic dipole moment is applied to the neon atom with the input variable nucdipmom. The orbital magnetic moment calculation is triggered at the end of the DDK calculation with
orbmag2 2
Both calculations must be done with the nuclear dipole moment on the atom of interest, in the direction of interest. Here, there is only a single atom, and due to spherical symmetry, all directions are equivalent, so it is sufficient to use
nucdipmom 1.0 0.0 0.0
Notice that we use a dipole moment of size 1; this is purely for convenience as we will see below. Running this file should take 3 minutes or so. After completion you will find in the output file:
Orbital magnetic moment computed with DFPT derivative wavefunctions
Orbital magnetic moment, Cartesian directions :
-5.54877537E-04 6.30837909E-15 4.28342968E-15
Chern vector, Cartesian directions :
-2.52162344E-09 1.10115145E-17 6.96089655E-18
The first result is the induced magnetic moment, which is a vector–note that it points in the original dipole direction. To get the shielding, divide the dipole strength (here, 1), and multiply by (-1) to get the normal NMR sign convention–the result is 554.9~ppm. This result is underconverged–the fully converged result, with ecut=30, is 552.1 ppm. For comparison, the all-electron, wavefunction-based calculation using post-Hartree-Fock CCSD(T) of [Vaara2003relativistic] yields 551.9 ppm, so clearly the present method is capable of excellent accuracy. The second line is the Chern vector [Ceresoli2006], that is, the integral of the Berry curvature, which outside of exotic topological insulators, should be zero. We include it as a convergence check. In this fast test calculation we obtained \(-2.5\times 10^{-9}\); in fully converged calculations one would like numbers about 3 orders of magnitude smaller.
A final detail to note: the induced moment due to the filled core orbitals, is due to what is called the Lamb shielding. Usual PAW datasets do not include this value, but you can add it to your input file with the variable lambsig. To compute it you must obtain the core wavefunctions, which is a normal option in Atompaw (see tutorial PAW2). Then, you can integrate the core wavefunctions to generate the Lamb shielding, through
where \(\alpha\) is the fine structure constant.
Moving forward, should you want to compute the shielding tensor for an atomic site in a solid, the procedure is identical, except that you will have to run the sequence three times, once for each direction of the applied dipole. Then you will have obtained three component vectors \(M_\alpha\), \(M_\beta\), and \(M_\gamma\), which you can assemble into a \(3\times 3\) matrix, and diagonalize to find the principal values and directions. The result is the shielding tensor at the site of interest, and orientated with respect to the Cartesian directions.