Calculating Molecular Magnetizabilities

Background

The magnetizability of an atom can formally defined as the second derivative of the electronic energy with respect to an applied magnetic field (other theoretical formulations exist, see Dimitrova et al.). It describes the diamagnetic response of the electrons in a molecule to an external magnetic field. This portion is always part of the magnetic response but it is not the full response (e.g., the paramagnetic response in systems with spins is not included and could be larger than the diamagnetic response, with opposite sign).

Tutorial

This tutorial will demonstrate how to calculate the Magnetizabilities for C₂H₂O

H₂C₂O (Geometry from Lutnaes et al.)

drawing

Step 1: Setup

The geometry of H₂C₂O was taken from Lutnaes et al. and converted to the geometry.in format for FHI-aims. If magnetizabilities are the only magnetic response observable being computed, then no further alterations to the geometry.in file are necessary.

# geometry.in
atom       0.000000    0.000000   -1.29234546954059    C
atom       0.000000    0.000000    .02048127473584    C
atom       0.000000    0.000000    1.18431817553677    O
atom       0.000000    0.93777557175219   -1.81423220511842    H
atom       0.000000   -0.93777557175219   -1.81423220511842    H

The control.in file can be set up as follows:

# control.in
xc                      pbe
relativistic            none
basis_threshold         0.0
magnetic_response       m

+ [basis sets for C, H, and O]

For this calculation the generalized gradient approximation PBE is sufficient. To appropriately compare to the non-relativistic reference, we set relativistic to none here. If relativity is important to your system, this can be enabled with relativistic atomic_zora scalar for a scalar relativistic treatment. The basis threshold must be set to 0.0 to allow the use of the basis sets tailored for NMR observables. Finally, the magnetic_response keyword must be set to m to compute only the magnetizabilities.

Specific to FHI-aims, the work Laasner et al. 2024 discusses in depth the impact of basis sets on the computation of various magnetic response parameters. For the purposes of this tutorial, we will compute chemical shieldings using 10 different basis set choices as can be found in the following table. All pathes are defined relative to the FHIaims base directory.

Basis Set	Path to species default files
FHI-aims "light"	./species_defaults/defaults_2020/light
FHI-aims "tight"	./species_defaults/defaults_2020/tight
NAO-VCC-2Z	./species_defaults/NAO-VCC-nZ/NAO-VCC-2Z
NAO-VCC-3Z	./species_defaults/NAO-VCC-nZ/NAO-VCC-3Z
NAO-VCC-4Z	./species_defaults/NAO-VCC-nZ/NAO-VCC-4Z
NAO-VCC-5Z	./species_defaults/NAO-VCC-nZ/NAO-VCC-5Z
NAO-J-2	./species_defaults/NAO-J-n/NAO-J-2
NAO-J-3	./species_defaults/NAO-J-n/NAO-J-3
NAO-J-4	./species_defaults/NAO-J-n/NAO-J-4
NAO-J-5	./species_defaults/NAO-J-n/NAO-J-5

Note that, in actual production calculations, it will not be necessary to probe that many basis sets.

In general, the FHI-aims "tight" species defaults are sufficient to produce accurate magnetizabilities at relatively low computational cost.

Step 2: Simulation and Locating Results

Now we can run the simulations! When a calculation has completed, open the aims.out file and search for:

  Magnetizability (10^-30 J/T^2)
  ------------------------------

Scrolling down you should be able to see the Magnetizability tensor broken down into the diamagnetic, paramagnetic, and total versions. For the purposes of this tutorial, only the total magnetizability tensor, and in particular the isotropic value immediately above this tensor is relevant.

Example:

    Diamagnetic:  -4.749157E+02 (isotropic)
        -4.285979E+02   0.000000E+00   0.000000E+00
         0.000000E+00  -4.845609E+02   0.000000E+00
         0.000000E+00   0.000000E+00  -5.115881E+02
    Paramagnetic:  5.153807E+01 (isotropic)
         3.636573E+01   1.097854E-12   6.289897E-13
         2.478904E-13   6.691885E+01  -3.826831E-13
         4.706894E-13  -5.455470E-13   5.132962E+01
    Total:        -4.233776E+02 (isotropic)
        -3.922322E+02   0.000000E+00   0.000000E+00
         0.000000E+00  -4.176421E+02   0.000000E+00
         0.000000E+00   0.000000E+00  -4.602585E+02

In this case (FHI-aims "light" species defaults), the magnetizability value is -423.3776*10^-30 J/T².

Step 3: Full Results and Comparison to Literature Values

The computed values for all 10 basis sets considered, as well as reference values from Lutnaes et al. can be found in the following table. The reference value given computed using "CCSD(T)" which is considered the gold standard for reference calculations. Note the short runtimes for all basis sets as compared to shieldings calculations (See shieldings tutorial). All calculations were computed locally on my MacBook pro with an M1 Max processor.

Basis Set	Runtime (s)	Magnetizability (\(10^{-30}\) \(J/T^2\))
Reference: CCSD(T)	unknown	−423.9
FHI-aims "light"	0.335 s	-423.3776
FHI-aims "tight"	2.175 s	-421.9135
NAO-VCC-2Z	1.770 s	-429.4391
NAO-VCC-3Z	3.253 s	-424.0457
NAO-VCC-4Z	6.078 s	-423.1789
NAO-VCC-5Z	13.118 s	-421.9167
NAO-J-2	9.495 s	-429.4176
NAO-J-3	14.635 s	-424.0307
NAO-J-4	27.163 s	-423.1682
NAO-J-5	60.879 s	-421.9080

From these results it is clear that the FHI-aims basis sets produce accurate results at much lower cost than other basis sets. For a further discussion of this see Laasner et al. 2024.