Source code for asr.c2db.magnetic_anisotropy

"""Magnetic anisotropy."""
import numpy as np
from ase import Atoms
import asr
from asr.core import (command, ASRResult, prepare_result,
                      option, AtomsFile)

from asr.database.browser import (
    table, make_panel_description, describe_entry, href)
from math import pi

# We don't have mathjax I think, so we should probably use either html
# or unicode.  But z does not exist as unicode superscript, so we mostly
# use html for sub/superscripts.
def equation():
    i = '<sub>i</sub>'
    j = '<sub>j</sub>'
    z = '<sup>z</sup>'
    return (f'E{i} = '
            f'−1/2 J ∑{j} S{i} S{j} '
            f'− 1/2 B ∑{j} S{i}{z} S{j}{z} '
            f'− A S{i}{z} S{i}{z}')

# This panel description actually assumes that we also have results for the
# exchange recipe.

panel_description = make_panel_description(
Heisenberg parameters, magnetic anisotropy and local magnetic
moments. The Heisenberg parameters were calculated assuming that the
magnetic energy of atom i can be represented as


where J is the exchange coupling, B is anisotropic exchange, A is
single-ion anisotropy and the sums run over nearest neighbours. The
magnetic anisotropy was obtained from non-selfconsistent spin-orbit
calculations where the exchange-correlation magnetic field from a
scalar calculation was aligned with the x, y and z directions.

        href("""D. Torelli et al. High throughput computational screening for 2D
ferromagnetic materials: the critical role of anisotropy and local
correlations, 2D Mater. 6 045018 (2019)""",

def spin_axis(theta, phi):
    if theta == 0:
        return 'z'
    elif np.allclose(phi, 90):
        return 'y'
        return 'x'

def webpanel(result, context):
    magstate = context.magstate().result
    if magstate['magstate'] == 'NM':
        return []

    magtable = table(result, 'Property',
                     ['magstate', 'magmom',
                      'dE_zx', 'dE_zy'], kd=context.descriptions)

    panel = {'title':
                 f'Basic magnetic properties ({context.xcname})',
             'columns': [[magtable], []],
             'sort': 11}
    return [panel]

[docs]@prepare_result class Result(ASRResult): spin_axis: str E_x: float E_y: float E_z: float theta: float phi: float dE_zx: float dE_zy: float key_descriptions = { "spin_axis": "Magnetic easy axis", "E_x": "Soc. total energy, x-direction [eV/unit cell]", "E_y": "Soc. total energy, y-direction [eV/unit cell]", "E_z": "Soc. total energy, z-direction [eV/unit cell]", "theta": "Easy axis, polar coordinates, theta [radians]", "phi": "Easy axis, polar coordinates, phi [radians]", "dE_zx": "Magnetic anisotropy energy between x and z axis [meV/unit cell]", "dE_zy": "Magnetic anisotropy energy between y and z axis [meV/unit cell]" } formats = {"webpanel2": webpanel} def spin_index(self): axis = self['spin_axis'] if axis == 'z': index = 2 elif axis == 'y': index = 1 else: index = 0 return index def spin_angles(self): return self['theta'] * 180 / pi, self['phi'] * 180 / pi
[docs]@command('asr.c2db.magnetic_anisotropy') @option('-a', '--atoms', help='Atomic structure.', type=AtomsFile(), default='structure.json') @asr.calcopt def main(atoms: Atoms, calculator: dict = { 'name': 'gpaw', 'mode': {'name': 'pw', 'ecut': 800}, 'xc': 'PBE', 'kpts': {'density': 12.0, 'gamma': True}, 'occupations': {'name': 'fermi-dirac', 'width': 0.05}, 'convergence': {'bands': 'CBM+3.0'}, 'nbands': '200%', 'txt': 'gs.txt', 'charge': 0 }) -> Result: """Calculate the magnetic anisotropy. Uses the magnetic anisotropy to calculate the preferred spin orientation for magnetic (FM/AFM) systems. Returns ------- theta: Polar angle in radians phi: Azimuthal angle in radians """ from asr.c2db.magstate import main as magstate from import calculate as calculategs from gpaw.spinorbit import soc_eigenstates from gpaw.occupations import create_occ_calc magstateresults = magstate( atoms=atoms, calculator=calculator) magstate = magstateresults['magstate'] # Figure out if material is magnetic results = {} if magstate == 'NM': results['E_x'] = 0 results['E_y'] = 0 results['E_z'] = 0 results['dE_zx'] = 0 results['dE_zy'] = 0 results['theta'] = 0 results['phi'] = 0 results['spin_axis'] = 'z' return Result(data=results) calculateresult = calculategs(atoms=atoms, calculator=calculator) calc = calculateresult.calculation.load() width = 0.001 occcalc = create_occ_calc({'name': 'fermi-dirac', 'width': width}) Ex, Ey, Ez = (soc_eigenstates(calc, theta=theta, phi=phi, occcalc=occcalc).calculate_band_energy() for theta, phi in [(90, 0), (90, 90), (0, 0)]) dE_zx = Ez - Ex dE_zy = Ez - Ey DE = max(dE_zx, dE_zy) theta = 0 phi = 0 if DE > 0: theta = 90 if dE_zy > dE_zx: phi = 90 axis = spin_axis(theta, phi) results.update({'spin_axis': axis, 'theta': theta / 180 * pi, 'phi': phi / 180 * pi, 'E_x': Ex * 1e3, 'E_y': Ey * 1e3, 'E_z': Ez * 1e3, 'dE_zx': dE_zx * 1e3, 'dE_zy': dE_zy * 1e3}) return Result(data=results)
if __name__ == '__main__': main.cli()