HeisenbergExchange

class HeisenbergExchange(configuration, kpoints=None, atom_indices=None, exchange_interaction_range=None, number_of_contour_points=None, integral_lower_bound=None, precalculate_eigenstates=None)

Class for calculating the inter-site Heisenberg exchange coupling matrix for a BulkConfiguration.

Parameters:

• configuration (BulkConfiguration | MoleculeConfiguration) – The BulkConfiguration for which to calculate the inter-site coupling matrix.

• kpoints (sequence (size 3) of int | MonkhorstPackGrid | KpointDensity) – The k-points with which to calculate the Heisenberg exchange constants.
  Default: MonkhorstPackGrid(1, 1, 1)

• atom_indices (All | list of ints.) – The atom indices for which to calculate the Heisenberg exchange constants.
  Default: All i.e. all atoms in the configuration.

• exchange_interaction_range (list of ints) – The range of the included interactions. Interactions are calculated between the primitive cell of the input configuration and translated cells given by the exchange_interaction_range. The exchange_interaction_range determines the number of translations of the configuration in the A, B, and C-directions given as a list of three positive integers, e.g. [3, 3, 3]. For a MoleculeConfigurtion [1, 1, 1] will always be used.
  Default: [5, 5, 5] for BulkConfiguration and [1, 1, 1] for a MoleculeConfiguration

• number_of_contour_points (int >= 2) – The number of points on the circle contour.
  Default: 30

• integral_lower_bound (PhysicalQuantity of type energy) – The distance between the lowest Fermi-level and the lowest energy circle contour point.
  Default: An energy determined on the chosen pseudopotentials or 1.5 Hartree for the semi-empirical calculators.

• precalculate_eigenstates (bool) – Boolean determining if the eigenstates and energies at all the k-points should be pre-calculated and stored for the real-space Green’s functions. Setting this option to True will reduce the computation time, but requires more memory.
  Default: False

DMIVectors()

Returns:

The Dzyaloshinskii-Moria interaction (DMI) vectors in the configuration. The DMI is only calculated if spin-orbit coupling is included. The shape of the array is [n_translations, n_atoms, n_atoms, 3], where n_translations is the number of translations, given by the product of the repetitions,

Return type:

PhysicalQuantity of type energy.

atomIndices()

Returns:

The atom indices.

Return type:

list of int.

atomSpins()

Returns:

The atomic spin polarizations.

Return type:

list of floats.

contour()

Returns:

The contour method to be used.

Return type:

OzakiContour | SemiCircleContour

couplingMatrix()

Returns:

The Heisenberg exchange coupling matrix corresponding the model \(H_{ij} = J_{ij} \hat{e}_i \cdot \hat{e}_j\). The shape of the array is [n_translations, n_atoms, n_atoms], where n_translations is the number of translations, given by the product of the repetitions, and n_atoms is the number of atoms in the configuration.

Return type:

PhysicalQuantity of type energy.

curieTemperature()

Calculates the Curie temperature based on the mean field approximation, \(T_c k_B = 2 J_0/3\), where \(J_0 = \sum_i J_{0i}\).

Returns:

The calculated Curie temperature.

Return type:

PhysicalQuantity of type temperature.

curieTemperatureRPA(q_grid=None, r_max=None)

Calculates the Curie temperature based on the random phase approximation (RPA), \(\frac{1}{T_c k_B} = \frac{6 \mu_B}{M} \sum_q \frac{1}{E(q)}\), where where \(E(q)\) is the magnon energy, and \(M\) is the magnetization.

Parameters:

• q_grid (MonkhorstPackGrid) – The q-points used for the RPA calculation.
  Default: MonkhorstPackGrid(8,8,8)

• r_max (PhysicalQuantity of type length.) – Maximum radius to include in the sum.
  Default: Include all calculated interactions.

Returns:

The calculated Curie temperature.

Return type:

PhysicalQuantity of type temperature.

distances()

Returns:

The distance between atoms corresponding to the coupling matrices. The shape of the array is [n_translations, n_atoms, n_atoms], where n_translations is the number of translations, given by the product of the repetitions, and n_atoms is the number of atoms in the configuration.

Return type:

PhysicalQuantity of type length.

exchangeInteractionRange()

Returns:

The interaction range of the configuration in the A, B, and C directions.

Return type:

list of int.

exchangeStiffness(r_max=None, damping_factor=None, magnetic_tolerance=None, atom_indices=None)

Calculate the exchange stiffness constant as \(A_{ex}=\frac{2 \mu_B}{3M} \sum_i J_{0i} \cdot R_{0i}^2 \cdot e^{-d\cdot R_{0i}}\), where \(M\) is the magnetization and \(d\) is a damping factor. If a damping factor is used, the actual exchange stiffness should be extrapolated to zero damping factor, as described in DOI:10.1103/PhysRevB.64.174402.

Parameters:

• r_max (PhysicalQuantity of type length.) – Maximum radius to include in the sum.
  Default: Include all calculated radii.

• damping_factor (PhysicalQuantity of type inverse length.) – Damping factor to be added to the sum.
  Default: 0*Ang**-1.

• magnetic_tolerance (float.) – Only include contribution from atoms with a magnetic moment larger than the magnetic tolerance.
  Default: 0.1

• atom_indices (list of int) – Include contributions to the exchange stiffness from atoms in atom_indices and their interaction with all the others.
  Default: All atom_indices defined upon construction

Returns:

The calculated exchange stiffness.

Return type:

PhysicalQuantity of type energy * distance^2.

extrapolatedExchangeStiffness(damping_factor_range=None, r_max=None, polynomial_order=None, magnetic_tolerance=None, atom_indices=None)

Calculate the extrapolated exchange stiffness constant as \(\lim_{\eta\rightarrow 0} \frac{2 \mu_B}{(3M} \sum_i J_{0i} R_{0i}^2 e^{-\eta R_{0i}}\), where \(M\) is the magnetization and \(\eta\) is a damping factor. If a damping factor is used, the actual exchange stiffness should be extrapolated to zero damping factor, as described in DOI:10.1103/PhysRevB.64.174402.

Parameters:

• damping_factor_range – List of damping factors used to extrapolate to zero damping_factor. The list must contain five or more entities.
  Default: A range of damping factors ranging from 0.1/d_0 to 1/d_0, where d_0 is the smallest nearest neighbor distance.

• r_max (PhysicalQuantity of type length.) – Maximum radius to include in the sum.
  Default: Include all calculated radii.

• polynomial_order (int.) – Order of the polynomial used to extrapolate to zero damping.
  Default: 5

• magnetic_tolerance (float.) – Only include contribution from atoms with a magnetic moment larger than the magnetic tolerance.
  Default: 0.1

• atom_indices (list of int) – Include contributions to the exchange stiffness from atoms in atom_indices and their interaction with all the others.
  Default: All atom_indices defined uppon construction

Returns:

The calculated exchange stiffness extrapolated to zero damping factor, a list of exchange stiffness values for each damping factor, the list of damping factors used, and the polynomial fit.

Return type:

PhysicalQuantity of type energy * distance^2, PhysicalQuantity of type energy * distance^2, PhysicalQuantity of type inverse length, numpy.polyfit

fourierTransformCouplings(qpoint, r_max=None, damping_factor=None)

Calculate the Fourier transform of the J_ij matrix.

Parameters:

• qpoint (list of three floats.) – Fractional q-point

• r_max (PhysicalQuantity of type inverse length.) – Maximum radius to include in the sum.
  Default: Include all calculated interactions.

• damping_factor – Exponential damping factor. A factor of exp(-damping_factor*R) is multiplied to the fourier transform.
  Default: Zero.

Returns:

The fourier transformed matrix.

Return type:

PhysicalQuantity of type energy.

kpoints()

Returns:

The k-points with which to calculate the Heisenberg exchange constants.

Return type:

class:~.MonkhorstPackGrid | RegularKpointGrid

magnonBandstructure(qpoints=None, points_per_segment=None, route=None, r_max=None, magnetic_tolerance=None)

Parameters:

• q_path (list of list of floats) – List of fractional q-points. This option is mutually exclusive to q_route and points_per_segment. The shape of the list is (number of q-points, 3), e.g. [[0.0, 0.0, 0.0], [0.25, 0.0, 0.0], [0.5, 0.0, 0.0]]

• route (list of str) – The route to take through the Brillouin zone as a list of symmetry points of the unit cell, e.g. ['G', 'X', 'G']. This option is mutually exclusive to q_path.

• points_per_segment (int) – The number of points per segment of the route. This option is mutually exclusive to q_path.
  Default: 20

• r_max (PhysicalQuantity of type length.) – Maximum interaction radius to include in the calculation.
  Default: Include all calculated interactions.

• magnetic_tolerance (float.) – Only include contribution from atoms with a magnetic moment larger than the magnetic tolerance.
  Default: 0.1

Returns:

A bandstructure object which can be saved and visualized.

Return type:

PhononBandstructure

metatext()

Returns:

The metatext of the object or None if no metatext is present.

Return type:

str | None

nlprint(stream=None)

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:

stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()

scaledCouplingMatrix()

Returns the coupling matrix elements scaled with the atomic spin polarizations. This corresponds to the matrix elements from the model \(H_{ij} = J_{ij} \hat{S}_i \cdot \hat{S}_j\).

Returns:

The coupling matrix elements. The shape of the array is [n_translations, n_atoms, n_atoms], where n_translations is the number of translations, given by the product of the repetitions, and n_atoms is the number of atoms in the configuration.

Return type:

PhysicalQuantity of type energy

setMetatext(metatext)

Set a given metatext string on the object.

Parameters:

metatext (str | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

spinComponentMatrices()

Returns:

Returns spin-component A-matrices as defined in e.g. Computer Physics Communications 264 (2021) 107938, eq. (14).

The shape of the array is [4, 4, n_translations, n_atoms, n_atoms], where the order of the first two indices is 0, x, y, z components, and where n_translations is the number of translations, given by the product of the repetitions, and n_atoms is the number of atoms in the configuration.

The A-matrices are not calculate for MoleculeConfigurations.

Return type:

PhysicalQuantity of type energy.

translations()

Returns:

The list of translations to be used.

Return type:

list of lists of int.

uniqueCouplingMatrixElementsAndDistances(atom_index=None, tolerance=None)

Parameters:

• atom_index (int) – The index of the atom for which to get the distances and coupling matrix elements.

• tolerance (PhysicalQuantity of type energy) – Energy tolerance for when two matrix elements are considered equal.
  Default: 0.01 meV

Returns:

Returns lists of distances and coupling matrix elements with unique elements as well as the multiplicity of the unique distances.

Return type:

PhysicalQuantity of type length, PhysicalQuantity of type energy, numpy.array of ints.

uniqueString()

Return a unique string representing the state of the object.

Usage Examples

Calculate Heisenberg exchange coupling for Fe (BCC):

# -*- coding: utf-8 -*-
# -------------------------------------------------------------
# Bulk Configuration
# -------------------------------------------------------------

# Set up lattice
lattice = BodyCenteredCubic(2.8665*Angstrom)

# Define elements
elements = [Iron]

# Define coordinates
fractional_coordinates = [[ 0.,  0.,  0.]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = SGGA.PBE

k_point_sampling = KpointDensity(
    density_a=7.0*Angstrom,
    )
numerical_accuracy_parameters = NumericalAccuracyParameters(
    density_mesh_cutoff=120.0*Hartree,
    k_point_sampling=k_point_sampling,
    )

calculator = LCAOCalculator(
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave('Fe-BCC.hdf5', bulk_configuration)

# -------------------------------------------------------------
# Heisenberg Exchange
# -------------------------------------------------------------
kpoints = KpointDensity(
    density_a=12.0*Angstrom,
    )

heisenberg_exchange = HeisenbergExchange(
    configuration=bulk_configuration,
    kpoints=kpoints,
    atom_indices=All,
    exchange_interaction_range=(9, 9, 9),
    number_of_contour_points=30,
    integral_lower_bound=Automatic,
    precalculate_eigenstates=False,
    )
nlsave('Fe-BCC.hdf5', heisenberg_exchange)

heisenberg_exchange.py

Notes

The Heisenberg exchange Hamiltonian is given by [1], [2]

(17)\[H = -\sum_{i\neq j}J_{ij}\hat{e}_i\hat{e}_j\]

where \(\hat{e}_i\) is the normalized local spin vector on atom \(i\) and \(J_{ij}\) is the Heisenberg exchange coupling constant.

We calculate the exchange coupling matrix elements \(J_{i\mathbf{0}, j\mathbf{R}}\) between a atom \(i\) in the central unit cell labeled \(\mathbf{0}\) and atom \(j\) in a different unit cell displaced from the central one with a lattice vector \(\mathbf{R}\). We use a real-space Green’s function implementation [3], [4] of the original Liechtenstein-Katsnelson-Antropov-Gubanov (LKAG) formula [1]:

(18)\[J_{0i,\mathbf{R}j} = -\frac{1}{4\pi}\int_{-\infty}^{E_F} d\epsilon {\rm Im Tr} \left[\Delta_i G_{i,j}(\epsilon, \mathbf{R}) \Delta_j G_{i,j}(\epsilon, -\mathbf{R}) \right]\]

where \(i\) and \(j\) label atomic indices within a unit cell, \(\mathbf{R}\) is a lattice vector, and where \(\Delta_i = H_i^\uparrow - H_i^\downarrow\) is the on-site difference between the up and down part of the Hamiltonian matrix. The real-space Green’s functions are calculated as

(19)\[G(\epsilon, \mathbf{R}) = \sum_{\mathbf{k}_i} w_{\mathbf{k}_i} G(\mathbf{k}_i, \epsilon) e^{-i\mathbf{k}_i\cdot \mathbf{R}}\]

where \(w_{\mathbf{k}_i}\) is the k-point weight and with the reciprocal space Green’s function given by:

\[G(\mathbf{k}, \epsilon) = (\epsilon \mathbf{S}(\mathbf{k}) - \mathbf{H}(\mathbf{k}))^{-1}\]

The energy integral in (18) is calculated using complex contour integration.

Curie temperature

The Curie temperature, \(T_C\), of a ferromagnet is the temperature at which the average magnetic moment becomes zero due to random fluctuations of the local magnetic moment directions. \(T_C\) can be estimated from the Heisenberg model in different ways. Within the mean field approximation (MFA) one have [2]

(20)\[T_C^{MFA} = \frac{2}{3k_B}\sum_{j\mathbf{R}\neq i\mathbf{0}}J_{i\mathbf{0}j\mathbf{R}}\]

The MFA Curie temperature is calculated as

curie_temperature_mfa = heisenberg_exchange.curieTemperature()

The MFA formula typically overestimates the Curie temperature when compared to experiments. A different approach is based on the random phase approximation (RPA) [2]:

(21)\[\frac{1}{k_B T_C^{RPA}} = \frac{6\mu_B}{M}\frac{1}{N_q}\sum_{\mathbf{q}}\frac{1}{E(\mathbf{q})}\]

where \(E(\mathbf{q})\) is the spin wave energy at q-point \(\mathbf{q}\) and \(M\) is the magnetic moment per atom and \(\mu_B\) is the Bohr magneton. The RPA Curie temperature is calculated as

q_grid = MonkhorstPackGrid(10,10,10)
curie_temperature_rpa = heisenberg_exchange.curieTemperatureRPA(q_grid=q_grid)

Since an average over the Brillouin zone is performed, the RPA calculation needs to be converged with respect to the number of q-points.

Exchange stiffness

The exchange stiffness describes the energy of a long wavelength spin wave excitation. In particular, if the energy as function of wave number \(q\) is written as

\[E(q) = E_0 + A_{ex} q^2 + \dots\]

the exchange stiffness constant, \(A_{ex}\) is seen to be given by the curvature of the \(E(q)\) vs. \(q\) curve. \(A_{ex}\) can be calculated from the Heisenberg exchange parameters as [2]

(22)\[A_{ex} = \frac{2\mu_B}{3M} \sum_{\mathbf{R}}\sum_{ij}J_{i\mathbf{0}, j\mathbf{R}} r_{i\mathbf{0}, j\mathbf{R}}^2\]

where \(r_{i\mathbf{0}, j\mathbf{R}}\) is the distance between atom \(i\) in cell \(\mathbf{0}\) and atom \(j\) in cell \(\mathbf{R}\). This can be calculated as

exchange_stiffness = heisenberg_exchange.exchangeStiffness()

In practice, (22) can be difficult to converge with respect to the sum over lattice vectors \(\mathbf{R}\). This is particularly pronounced in the case of iron where the exchange parameters are very long ranged. In this case it is possible to dampen the oscillations and extrapolate the result to zero damping as:

(23)\[\begin{split}A_{ex} &= \lim_{\eta\rightarrow 0} A_{ex}(\eta) \\ A_{ex}(\eta) &= \frac{2\mu_B}{3M} \sum_{\mathbf{R}}\sum_{ij}J_{i\mathbf{0}, j\mathbf{R}} r_{i\mathbf{0}, j\mathbf{R}}^2 e^{-\eta r_{i\mathbf{0}, j\mathbf{R}}}\end{split}\]

The extrapolation to \(\eta \rightarrow 0\) can be calculated as

extrapolated_exchange_stiffness, exchange_stiffness_vec, damping_factor_range, p_fit \
    = heisenberg_exchange.extrapolatedExchangeStiffness()

where the final extrapolated result is given by extrapolated_exchange_stiffness. The other returned quantities are a vector of calculated exchange stiffness values obtained for the damping factors given by damping_factor_range. Finally the polynomial fit is also returned.

Spin-scaling

The Heisenberg exchange Hamiltonian is sometimes written as

\[H = -\sum_{i\neq j}\tilde{J}_{ij}\hat{S}_i\hat{S}_j\]

where \(\hat{S}_i\) is the local spin vector on atom \(i\) and \(J_{ij}\) is the Heisenberg exchange coupling constant. Following [1], [5], in QuantumATK we calculate the coupling constants for the Hamiltonian

\[H = -\sum_{i \neq j}J_{ij}\hat{e}_i\hat{e}_j\]

where the spin operators are substituted by unit vectors \(\hat{e}_i\) pointing in the direction of the i’th site magnetization.

Both the coupling constants \(J_{ij}\) and \(\tilde{J}_{ij}\) are available in a HeisenbergExchange object and are simply relates as

\[\tilde{J}_{ij} = \frac{J_{ij}}{S_i\cdot S_j}\]

where \(S_i\) is calculated from the MullikenPopulation on each atom as \(S_i = (M_i^{up} - M_i^{down})/2\) where \(M_i^{up}\) is the Mulliken population for up-electrons on atom \(i\).

The coupling constants \(J_{ij}\) are obtained as:

coupling_matrix = heisenberg_exchange.couplingMatrix()

and the coupling constants \(\tilde{J}_{ij}\) are obtained as:

scaled_coupling_matrix = heisenberg_exchange.scaledCouplingMatrix()

[1] (1,2,3)

A.I. Liechtenstein, M.I. Katsnelson, V.P. Antropov, and V.A. Gubanov. Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys. Journal of Magnetism and Magnetic Materials, 67(1):65 – 74, 1987. URL: http://www.sciencedirect.com/science/article/pii/0304885387907219, doi:https://doi.org/10.1016/0304-8853(87)90721-9.

[2] (1,2,3,4)

M. Pajda, J. Kudrnovský, I. Turek, V. Drchal, and P. Bruno. Ab initio calculations of exchange interactions, spin-wave stiffness constants, and curie temperatures of fe, co, and ni. Phys. Rev. B, 64:174402, Oct 2001. URL: https://link.aps.org/doi/10.1103/PhysRevB.64.174402, doi:10.1103/PhysRevB.64.174402.

[3]

Asako Terasawa, Munehisa Matsumoto, Taisuke Ozaki, and Yoshihiro Gohda. Efficient algorithm based on liechtenstein method for computing exchange coupling constants using localized basis set. Journal of the Physical Society of Japan, 88(11):114706, 2019. URL: https://doi.org/10.7566/JPSJ.88.114706, doi:10.7566/JPSJ.88.114706.

[4]

Xu He, Nicole Helbig, Matthieu J. Verstraete, and Eric Bousquet. Tb2j: a python package for computing magnetic interaction parameters. Computer Physics Communications, 264:107938, 2021. URL: https://link.aps.org/doi/10.1103/PhysRevB.64.174402, doi:10.1103/PhysRevB.64.174402.

[5]

Myung Joon Han, Taisuke Ozaki, and Jaejun Yu. Electronic structure, magnetic interactions, and the role of ligands in $\mathrm Mn_n(n=4,12)$ single-molecule magnets. Phys. Rev. B, 70:184421, 2004. URL: https://link.aps.org/doi/10.1103/PhysRevB.70.184421, doi:10.1103/PhysRevB.70.184421.