-
Notifications
You must be signed in to change notification settings - Fork 7
Open
Description
Is your feature request related to a problem? Please describe.
No, just adding an example
Describe the solution you'd like
Probably something more general and cleaner than what is provided below.
Describe alternatives you've considered
NA
Additional context
Thanks for this awesome library! I'd like to use it for magnetic structure refinement.
I made a quick and dirty example reproducing the analysis for LaMnO3 appearing in this tutorial.
The code is quite ugly, I barely tested it and it only works for representations of dimension 1, but it should be easy to extend it. Maybe it turns out to be useful to someone out there or maybe it can be improved to become a more general magnetic example. If not, feel free to close this issue.
from ase.io import read
from spgrep import get_spacegroup_irreps
from spgrep.representation import get_character, check_spacegroup_representation, project_to_irrep
import numpy as np
# This implements Appendix B in
# https://neutrons2.ornl.gov/conf/2014/magstr/docs/Tutorial_Magnetic_Structures.pdf
# This can be obtained from https://legacy.materialsproject.org/materials/mp-1180700/cif?type=conventional_standard&download=true
atoms = read('LaMnO3_mp-1180700_conventional_standard.cif')
qpoint = [0., 0., 0.]
# We need to use all positions, otherwise a different spacegroup is found.
irreps, rotations, translations, mapping_little_group = get_spacegroup_irreps(
atoms.cell.array, atoms.get_scaled_positions(), atoms.get_atomic_numbers(), qpoint
)
print("Found ", len(irreps), " irreps")
little_rotations = rotations[mapping_little_group]
little_translations = translations[mapping_little_group]
little_order = len(little_rotations)
# Compute rotations of Axial vectors
axial_rep = np.zeros_like(little_rotations)
for i, l in enumerate(little_rotations):
axial_rep[i] = l * np.linalg.det(l)
# Compute permutations and consider only magnetic atoms
# In this case, Mn atoms are 4,5,6,7
# code below is taken from phonon example
positions = atoms.get_scaled_positions()[4:8]
nmagat = len(positions)
# Operation-`i` moves atom-`kappa` to `permutations[i, kappa]`
permutations = np.zeros((little_order, nmagat), dtype=int)
for i, (Ri, vi) in enumerate(zip(little_rotations, little_translations)):
for kappa, position in enumerate(positions):
new_pos = np.remainder(Ri @ position + vi, 1)
for kappa2, position2 in enumerate(positions):
if np.allclose(position2, new_pos):
permutations[i, kappa] = kappa2
break
shifts = np.zeros((little_order, nmagat, 3))
for i, (Ri, vi) in enumerate(zip(little_rotations, little_translations)):
perm_i = permutations[i]
shifts[i] = positions @ Ri.T + vi[None, :] - positions[perm_i]
perm_rep = np.zeros((little_order, 4, 4), dtype=np.complex128)
for i, Ri in enumerate(little_rotations):
for kappa in range(nmagat):
kappa2 = permutations[i, kappa]
perm_rep[i, kappa2, kappa] = np.exp(
-2j * np.pi * np.dot(Ri.T @ qpoint, shifts[i, kappa])
)
# Compute Magnetic operators
from spgrep.representation import get_direct_product
magop = get_direct_product(perm_rep, axial_rep)
# Compute characters of magnetic operators and its decomposition into irreps
magchi = get_character(magop)
for i, irrep in enumerate(irreps):
print (f'Irrep {i} contained', magchi @ irrep.reshape(8,-1) / little_order)
# Compute basis for a single component and a single irrep
r = np.zeros(3*nmagat, dtype=complex)
irreducible_representation_num = 0
a = np.zeros(3*nmagat)
a[0] = 1 # x component first atom
for g in range(little_order):
chi = get_character(irreps[irreducible_representation_num])[g]
r += np.conj(chi) * ( magop[g] @ a )
r /= little_order
for i in range(nmagat):
print(r[3*i:3*(i+1)])
# Same as above, but for all irreps and all components
for i, irrep in enumerate(irreps):
r = np.zeros([3*nmagat,3], dtype=complex)
a = np.zeros([3*nmagat, 3])
a[0,0] = 1
a[1,1] = 1
a[2,2] = 1
for g in range(little_order):
chi = get_character(irrep)[g]
r += np.conj(chi) * ( magop[g] @ a )
r /= little_order
print(f'Irreducible representation {i}\n Basis:')
for j in range(nmagat):
print(f'Atom {j}')
print(r[3*j:3*(j+1)])
Edit: added conjugate, does not change the result but it's formally correct.
Reactions are currently unavailable
Metadata
Metadata
Assignees
Labels
No labels