Advanced charge density and bond order analysis¶

NBO analysis of EFG, NMR chemical shifts, NMR spin-spin coupling¶

For certain molecular properties it is possible to perform detailed analyses in terms of Natural Bond Orbitals (NBOs) and Natural Localized Molecular Orbitals (NLMOs). These features generally require a sequence of ADF and/or property code runs. An initial non-relativistic or scalar relativistic ADF run, followed by the generation of NBO and NLMO data, is required, and the resulting data files need to be present in subsequent property calculations, along with a keyword indicating that the NBO analysis is requested in the property module.

We have noted in the past some slight loss of numerical accuracy of the results after going through the various orbital transformations in the NBO - NLMO sequence. It is important that the user verifies in each case that the total contributions from the analysis are in agreement with the total calculated property, within the numerical integration accuracy limits. In order to assist the user with this, the analysis program always print the total analysis contributions, including small non-printed values.

Moreover, there appears to be a problem with the analysis of the Fock matrix in the NBO program in conjunction with ADF calculations. Therefore please do NOT use the Fock matrix second order perturbation theory analysis in NBO at this time. We will remove this disclaimer once the issue has been fixed. Applications of the NBO-NLMO property analysis codes have so far given no indication that the Fock matrix issue interferes with the analysis.

Important note: If properties are analyzed from within spin-orbit relativistic computations, the NBO/NLMO analysis is performed in terms of scalar (spin-free) relativistic orbitals, as detailed in the technical references. The results from these analyses are exact in the sense that they fully reproduce the final spin-orbit property result, and they allow to dissect the property in terms of more intuitive one component real scalar relativistic localized orbitals. Typically, the property analysis in a spin-orbit calculation involves contributions from unoccupied scalar NLMOs, whereas there are no such contributions if a non-relativistic or scalar relativistic property is analyzed.

Available properties for NBO analysis: EFG, NMR chemical shifts and NMR spin-spin coupling.

NBO analysis of EFG

EFGs: non-relativistic and scalar ZORA, in ADF/AOResponse. Requires initial ADF run with

AOresponse
 donothing
End

in order to generate orbitals that re equivalent to those generated in the subsequent ADF run where the EFG is calculated. Alternatively, simply calculate the EFG twice, once before the NBO generation step, and once afterward.

The next step (see below) is to create the NBOs and the required data files for the analysis. Afterward, in the second ADF run, use

Aoresponse
 efg NUC nbo
end

efg NUC nbo

Here NUC is the number of the nucleus at which the EFG is to be computed (ADF internal atom ordering). Example: efg 1 nbo.

The threshold for printing the EFG-NBO contributions to output can be adjusted via the Tresh options. The default is 0.05, which means that only orbitals with absolute value contribution larger than 5% of the total EFG are printed. To increase the number of contributions printed, specify a smaller threshold. For example:

Aoresponse
   efg 1 nbo Tresh 0.01
end

In addition to the optional NBO analysis, the EFG program in AOResponse prints a Mulliken type analysis of the EFG principal components, and an analysis in terms of canonical MOs.

WARNING: The ordering of the principal components is lowest to highest including the sign. That is, we have \(V11 \leq V22 \leq V33\). This does not conform to the usual convention of \(|V11| \leq |V22| \leq |V33|\). Please make sure you select the right component for your analysis.

Example job: $AMSHOME/examples/adf/AlCl3_efgnbo. For an explanation of the output and a general usage tutorial, see [7]. Further references and recommended citations, see [8].

NBO analysis of NMR Chemical shift

An implementation is currently available for spin-orbit ZORA computations. If scalar ZORA calculations are to be analyzed, provide the input keyword

FAKESO

in the NMR input (outside of the ‘nmr’ keyword). If this feature is requested one should restrict the calculation to a single shielding tensor per NMR run. It would be good practice to check the results against regular NMR calculations where the analysis feature is not requested. No ZORA scaling is applied in the analysis results. The data should be equivalent to a regular computation in the NMR input with

NMR
 u1k best
 calc all
END

Depending on whether scalar or spin-orbit calculations are to be analyzed, the sequence of calculations is different:

scalar:

1. ADF, scalar ZORA, symmetry NOSYM, closed shell

2. generate NBOs and required data files for analysis

3. NMR with FAKESO and analysis keywords, use adf.rkf (previously known as TAPE21), TAPE10 from step 1.

spin-orbit:

1. ADF, scalar ZORA, closed shell

2. generate NBOs and required data files for analysis

3. delete adf.rkf, TAPE10, TAPE15

4. ADF, spin-orbit ZORA, symmetry NOSYM, closed shell

5. NMR with analysis keywords, using adf.rkf, TAPE10 from step 4

In the NMR run, in addition to the NMR keyword, provide the following

analysis
 print 0.01
 canonical
 nbo
 components
end

The optional canonical keyword can be used independently from the NBO analysis features. It enables an analysis of the shielding in terms of the canonical MOs. The components keyword is optional and enables an analysis not only of the isotropic shielding but also of each principal component of the tensor. The print keyword selects printout of contributions relative to the total diamagnetic, paramagnetic. In the example, only contributions greater than 1% are printed. Set to zero to print ALL contributions.

Example job: $AMSHOME/examples/adf/CH4_nmrnbo. References [40] [41] [42].

NBO analysis of NMR spin-spin coupling (J-coupling)

Non-relativistic, scalar ZORA, spin-orbit ZORA

The sequence of jobs is similar to those in the NMR section.

scalar or non-relativistic:

1. ADF, scalar ZORA or nonrel.

2. generate NBOs and required data files for analysis

3. CPL with analysis keyword, use adf.rkf (previously known as TAPE21), TAPE10 from step 1.

spin-orbit:

1. ADF, scalar ZORA

2. generate NBOs and required data files for analysis

3. delete adf.rkf, TAPE10, TAPE15

4. ADF, spin-orbit ZORA

5. CPL with analysis keyword, using adf.rkf, TAPE10 from step 4

In the CPL run provide the following ‘contributions’ keyword to enable the analysis

nmrcoupling
   ... other options
   contributions 1E19 nbo
end

The numerical value selects a print threshold in SI units of T**2/J for the analysis. Increase the value to obtain less detail in the analysis. By default, ‘contributions’ triggers an analysis of the J-coupling in terms of canonical MOs. The nbo keyword enables in addition the NBO-NLMO analysis.

Please note that due to the history of how the program was developed the output from the scalar/nrel. analysis and from the spin-orbit calculations differs somewhat. The qualitative content is the same.

In scalar ZORA or non-relativistic CPL calculations without the SD term an orbital based analysis is only performed for the Fermi-contact mechanism. If you also need an analysis for the PSO and SD mechanisms but do not want to run a spin-orbit calculation with ADF please use the SD or NOSD keywords which will cause the spin-orbit branch of the CPL code to be used. In ZORA spin-orbit calculations the FC, SD, PSO, and cross terms are analyzed together by default. You can selectively switch them on or off in order to get individual mechanism analyses. The DSO mechanism is often negligible. An analysis tool for this mechanism has therefore not yet been developed.

Example job: $AMSHOME/examples/adf/CPL_CH3OH_NBO. References NMR spin-spin couplings with NBO analysis [42] [43] [44] [45]:

Generation of NBOs

How to generate the NBOs, NLMOs, and the data files needed for these calculations (step 2 below is step 2 in the examples above):

1. run ADF with scalar ZORA or non-relativistic options, and keep adf.rkf and TAPE15.

AMS_JOBNAME=Scalar $AMSBIN/ams <<eor
  ...
eor

# run adfnbo in WRITE  mode to create the gennbo input file FILE47
# and one of the required property analysis files, adfnbo.kf

$AMSBIN/adfnbo << eor
 ADFfile Scalar.results/adf.rkf
 TAPE15file Scalar.results/TAPE15
 write
 spherical
eor

rm -f adfnbo.37 adfnbo.39 adfnbo.49 adfnbo.48
$AMSBIN/gennbo6 << FILE47

# run adfnbo in COPY mode to create the second property analysis
# file, adfnbo2.kf

$AMSBIN/adfnbo <<  eor
 ADFfile Scalar.results/adf.rkf
 TAPE15file Scalar.results/TAPE15
 spherical
 copy
eor

# run adfnbo in READ mode: prepare locorb on TAPE21

$AMSBIN/adfnbo <<  eor
 ADFfile scalar.results/adf.rkf
 TAPE15file scalar.results/TAPE15
 spherical
 read
eor

rm -f adfnbo.37 adfnbo.39 adfnbo.49 adfnbo.48

# keep the adf.rkf after this sequence in order to
# be able to plot the NBOs and NLMOs with amsview

# clean up, keep adfnbo*.kf for any NBO property analyses.

Local, atomic, and non-local properties¶

The QTAIM input block replaces the Bader keyword. The former Bader Reactivity option now corresponds to the new ConceptualDFT block (see below).

The QTAIM functionality encompasses two different but related features: a topological analysis of the electron density gradient field [10], and a calculation of condensed atomic properties by integrating various local functions over atomic basins [46] [47]. The topological analysis is performed at all levels of the analysis (normal, extended or full) and it produces a list of the molecule’s electron density critical points and bond paths that can be visualized by the AMSview GUI module. Besides, the following properties at the critical points are calculated: electron density, its Laplacian and Hessian matrix, ellipticity, Jenkins’ metallicity index [12], reduced density gradient, DFT energy density variation rates [13], and Abramov’s local energy components [14]. At the Extended analysis level the following condensed atomic properties are calculated: charges, density Laplacian, dipole and quadrupole moments, Shannon and Renyi entropies, spatial distribution of the molecular orbitals over atomic basins. If the Energy keyword is specified, the atomic energies [11] are also calculated. At the Full level, the atom-atom matrix elements of the localization (LI) and delocalization (DI) indices [48] [49] are computed. In ADF 2023, we also added the (optional) calculation of the Source Function introduced by Gatti and Bader [54] [55]. This function is evaluated for each bond or ring critical point of the molecule, allowing a decomposition of the values of the electron density at these points as a sum of atomic contributions. Since this calculation requires the integration of the laplacian of the electron density, we strongly recommend a numerical quality set (at least) to ‘good’ (or higher).

QTAIM

Type:

Block

Description:

This block is used to request a topological analysis of the gradient field of the electron density, also known as the Bader’s analysis. If this block is specified without any sub-key, only local properties are calculated.

AnalysisLevel

Type:

Multiple Choice

Default value:

Normal

Options:

[Normal, Extended, Full]

Description:

Set the level of the QTAIM analysis: Normal - topology analysis and properties at the density critical points, Extended - same as Normal plus condensed atomic descriptors, Full - same as Extended plus non-local descriptors.

Source

Type:

Bool

Default value:

No

GUI name:

Source Function

Description:

Calculate the Source Function at BCPs and RCPs.

AtomsToDo

Type:

Integer List

GUI name:

Include atoms

Description:

List of atoms for which condensed descriptors are to be calculated. By default all atoms are included.

Spacing

Type:

Float

Default value:

0.5

Unit:

Bohr

Description:

Specifies spacing of the initial Cartesian grid when searching for critical points. It may be useful to specify a smaller value than the default if some critical points are missed. This will result in a more accurate but slower calculation.

Note that the accuracy of condensed atomic properties can be estimated by integrating the Laplacian of the electron density over a given basin, which, ideally, must vanish. The accuracy of the method can be improved by using a larger integration grid (NumericalQuality). Usually, the default grid is sufficient to achieve the accuracy of 10^-3 a.u. (differences of milliHartree in the energies). Unfortunately, the convergence of the electron density Laplacian integral with the grid size is not monotonous. Therefore this type of Bader atomic property calculation should be considered for applications where computational efficiency is critical and moderate accuracy is sufficient.

ADF2AIM¶

The ADF utility adf2aim (original name rdt21) developed by Xavi Lopez, Engelber Sans and Carles Bo converts an ADF adf.rkf (formerly known as TAPE21) to WFN format (for Bader analysis)

The program rdt21 is now called adf2aim and is part of the ADF package.

The WFN file is an input file for the third party program Xaim (see http://www.quimica.urv.es/XAIM for details), which is a graphical user interface to programs that can perform the Bader analysis. Usage of adf2aim can be found in the Examples Document.

Aromaticity index with QTAIM¶

The cyclic delocalization of mobile electrons in two or three dimensions is probably one of the key aspects that characterize aromatic compounds. However, aromaticity cannot be measured directly by any physical or chemical experiment because it is not a well-defined magnitude. One can define indicators of aromaticity, some of which are implemented in ADF. The matrix of localization/delocalization indices (LI-DI) and the following from it Iring (ring index) and MCI (multi center index) aromaticity indices (see Ref. [16]) can be computed for rings specified by the user. The rings are specified using the AROMATICITY input block, one line per ring.

Aromaticity
   RingAtoms integer_list
End

Aromaticity

Type:

Block

Description:

Calculate aromaticity indicators, i.e. the matrix of localization/delocalization indices (LI-DI), Iring (ring index) and MCI (multi center index) aromaticity indices.

RingAtoms

Type:

Integer List

Recurring:

True

Description:

List of the atom numbers that form a ring. For each ring a separate line.

For example, for a benzene molecule where atoms 1-6 are carbon atoms the input would look like this:

AROMATICITY
  ringatoms 1 2 3 4 5 6
END

The program does not check whether the atoms in the AROMATICITY input block are connected so it is responsibility of the user to specify them correctly. Moreover, for the Iring index to be calculated correctly the atoms must be listed in the order they are connected in the ring.

Calculation of the aromaticity indices invokes the QTAIM analysis automatically and uses its results so all accuracy considerations for the QTAIM key above are also valid here. This feature requires that the calculation is performed without symmetry (Symmetry NOSYM).

Results of the calculation are printed to the output under the “Aromaticity Index” header. The first table presents the total number of electrons per atom and their breakdown in “localized” and “delocalized”. Note: the total number of electrons per atom \(N_{total} = N_{loc} + N_{deloc} / 2\) since delocalized electrons are counted twice. The second table presents the Iring and MCI ring indices for each ring specified in the AROMATICITY block. These results (and the original LiDi matrix from which they are computed) are stored in the Properties section in the t21 file.

Global, atomic, and non-local descriptors¶

The ConceptualDFT block replaces and extends the set of features behind the Bader Reactivity keyword found in the ADF versions before ADF2019. The frontier molecular orbital (FMO) approximation is used for all calculations. Orbital degeneracy is taken into account following Martinez [17]. All these descriptors are described in Ref. [18]. Which exactly descriptors are calculated depends on the AnalysisLevel key:

• Normal: global conceptual DFT descriptors, i.e. electronic chemical potential, electronegativity (its opposite), hardness and softness, hyperhardness [19], electrophilicity index [20], dissociation energies (nucleofuge and electrofuge, see Ref. [21]), electrodonating and electroaccepting powers [22] and the net electrophilicity [23]. Also, new global descriptors called global dual descriptors \(\Delta f+\) and \(\Delta f-\) are calculated [24], which correspond to the spatial integral of the positive (electrophilic) or negative (nucleophilic) regions of the original dual descriptor, respectively. These descriptors allow for comparing total reactivity of molecules, in contrast with Fukui functions or original dual descriptor that integrate over the whole space to 1 or 0 for all molecules.

• Extended: both global and atomic properties, including condensed (QTAIM) descriptors in both canonical and grand canonical ensembles. More precisely, the following properties are computed: condensed Fukui functions [26], dual descriptor [25] [27] and some of their composite functions [50] [51] [52], as well as condensed local electrophilicity and nucleophilicity [28].

• Full: all properties are calculated, including non-local ones, such as the condensed linear response function [29] and the softness kernel [30].

Please note that atomic electronegativities [31] are calculated separately with an optional Electronegativity keyword. In contrast with other properties, it requires all-electron basis sets and triggers the TotalEnergy keyword.

ConceptualDFT

Type:

Block

Description:

Conceptual DFT Properties

AnalysisLevel

Type:

Multiple Choice

Default value:

Normal

Options:

[Normal, Extended, Full]

Description:

Set the level of the ConceptualDFT analysis: Normal - global descriptors only, Extended - both global and condensed (QTAIM) local descriptors, Full - all descriptors including non local ones.

Electronegativity

Type:

Bool

Default value:

No

GUI name:

Atomic electronegativities

Description:

Calculate atomic electronegativities. Requires an all-electron calculation (no frozen core), triggers the TotalEnergy and increases the [AnalysisLevel] to at least Extended.

AtomsToDo

Type:

Integer List

GUI name:

Include atoms

Description:

Define a subset of atoms for which properties are calculated. If the [Domains] block is present then this list specifies which atoms are used to define the domains bounding box.

For more information about Conceptual DFT, we suggest the reading of these excellent reviews: Refs. [33] and [32].

Domains of the dual descriptor¶

In a recent paper [34], Tognetti et al. proposed a real-space partitioning scheme of the dual descriptor (DD) into chemically meaningful reactivity domains. Basically, a predominantly electrophilic (nucleophilic) DD domain constitutes a region of space where this function remains positive (negative). They can be visualized using the GUI and quantitative information can be obtained using the additional Domains key block.

ConceptualDFT

Type:

Block

Description:

Conceptual DFT Properties

Domains

Type:

Block

Description:

Calculate integrated properties for the domains (same sign) of the dual descriptor.

Border

Type:

Float

Default value:

7.0

Unit:

Bohr

Description:

Set the extent of the Cartesian grid. Extent is the distance between a face of the grid’s bounding box and the most outlying atom in the corresponding direction. If the [AtomsToDo] key is present, the bounding box is created around the specified atoms.

Display

Type:

Float

Default value:

0.005

Description:

Domains for which the integrated DD value is smaller (in magnitude) than the specified value are omitted from the printed output.

Enabled

Type:

Bool

Default value:

No

GUI name:

Properties of reactivity domains

Description:

Calculate properties of reactivity domains.

Ensemble

Type:

Multiple Choice

Default value:

Canonical

Options:

[Canonical, GrandCanonical]

Description:

Statistical ensemble for DD domains. Canonical: DD values are calculated using the statistical canonical ensemble. GrandCanonical: DD values are calculated using the statistical grand canonical ensemble. The grand canonical DD corresponds to (S^2 f(2) - (gamma/eta^3) f^0), where f(2) is the canonical DD, gamma and eta - the hyper-hardness and hardness of the chemical system, respectively, and f^0 is the mean Fukui function. This statistical ensemble is a natural choice when comparing two chemical systems with a different number of electrons.

Radius

Type:

Float

Default value:

0.0

Description:

This option adds a sphere around each nucleus, excluding all points inside it. This can help to separate domains around an atom or to exclude core electrons. Be careful when using this option. In particular, the radius of the sphere should exceed two or three times the [Spacing] value to be effective. By default, no spheres are added.

Spacing

Type:

Float

Default value:

0.1

Unit:

Bohr

Description:

Specifies spacing (distance between neighboring points) of the rectangular Cartesian grid used when searching for DD domains. It may be useful to specify a smaller value (or increase the size of the grid, see [Border] key) if a substantial part of the electronic density is accounted for.

Threshold

Type:

Float

Default value:

0.001

Description:

Arbitrary value of dual descriptor used to separate DD domains (values below this threshold are ignored).