Calculate magnetic parameters
grogupy is specifically written in a modular way, to fully take advantage of python classes. In this simple example we show how to use grogupy to simulate the magnetic interaction between two atoms from a density functional theory (DFT) calculation and create an output file for magnopy. In this example we will use a non-collinear Siesta calculation of the CrI3 system with spin-orbit coupling.
First import grogupy, which will import the most important classes, functions and variables in its namespace.
[1]:
import grogupy
/Users/danielpozsar/Documents/Documents_MacBook Air/studies/elte/phd/grogu/.venv/lib/python3.13/site-packages/grogupy/_tqdm.py:41: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from tqdm.autonotebook import tqdm
Then we start to create the basic classes based on the desired parameters. First create the Kspace class which contains the parameters fot the Brillouin zone integration. Because in this case we only take a single layer of CrI3 and there is a large vacuum in the perpendicular direction it is sufficient to only integrate in the plane of the material. Furthermore for a fast initial calculation a 10x10 grid of k-points should be enough, but this is not fully converged.
[2]:
CrI3_kspace = grogupy.Kspace(kset=[10, 10, 1])
CrI3_kspace
[2]:
<grogupy.Kspace kset=[10 10 1], NK=100>
Then we would like to set the parameters for the complex integral for the Green’s function calculation. This can be done through the Contour class. Because CrI3 is an insulator 50 sample points should be enough. The energy minimum should be set below the smallest eigenvalue from the Siesta Hamiltonian.
[3]:
CrI3_contour = grogupy.Contour(emin=-30, eset=50, esetp=100)
CrI3_contour
[3]:
<grogupy.Contour emin=-35, emax=0, eset=50, esetp=100>
Then we can create the Hamiltonian class, which extracts and stores the Hamiltonian and geometrical data from the Siesta files using the sisl library. Furthermore we must provide the exchange field orientation in the DFT calculation, which is usually the perpendicular direction from the 2D material.
[4]:
CrI3_hamiltonian = grogupy.Hamiltonian(
infile="./../../../benchmarks/CrI3/CrI3.fdf",
scf_xcf_orientation=[0, 0, 1],
)
CrI3_hamiltonian
Setting up Hamiltonian: 100%|██████████| 147/147 [00:00<00:00, 618.10it/s]
[4]:
<grogupy.Hamiltonian scf_xcf_orientation=[0. 0. 1.], orientation=[0. 0. 1.], NO=216>
After these initial steps we can create the Builder class and initialize with the above information. We have to provide the directions of the rotated exchange field and the two perpendicular directions are calculated automatically.
[5]:
orientations = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
CrI3 = grogupy.Builder(ref_xcf_orientations=orientations)
CrI3.add_kspace(CrI3_kspace)
CrI3.add_contour(CrI3_contour)
CrI3.add_hamiltonian(CrI3_hamiltonian)
CrI3
[5]:
<grogupy.Builder npairs=0, numk=100, kset=[10 10 1], eset=50>
Then we have to define the magnetic entities, which can be a list of orbitals from the Siesta Hamiltonian. We use sisl functions to extract this information based on a more human like definition. For example take two shells from two Cr atoms.
[6]:
magnetic_entities = [dict(atom=0, l=2), dict(atom=1, l=2)]
With this information we can create two MagneticEntity instances.
[7]:
CrI3.add_magnetic_entities(magnetic_entities)
Add magnetic entities: 100%|██████████| 2/2 [00:00<00:00, 95.47it/s]
This would be enough to calculate the anisotropy on both sites, but it is not enough to build a spin model, for this we need the exchange tensor between all the pairs. Generally this can be done in a similar way. First a list of dictionaries must be created, where each dictionary contains ai, aj and Ruc, where the first two is the index from the magnetic_entities and the third is the supercell shift of the second magnetic entity.
[8]:
pairs = [
dict(ai=0, aj=1, Ruc=[0, 0, 0]),
dict(ai=0, aj=1, Ruc=[1, 0, 0]),
]
CrI3.add_pairs(pairs)
CrI3
Add pairs: 100%|██████████| 2/2 [00:00<00:00, 4857.33it/s]
[8]:
<grogupy.Builder npairs=2, numk=100, kset=[10 10 1], eset=50>
Now every information is contained in the CrI3 instance to run the simulation. There are multiple parameters that can be changed to tune the runtime and the precision of the simulation, but for now let us use the default parameters.
[9]:
CrI3.solve()
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 646.84it/s]
Rotating Exchange field: 100%|██████████| 3/3 [00:00<00:00, 78.04it/s]
Rotation 1: 100%|██████████| 100/100 [00:08<00:00, 12.04it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 599.39it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 7936.24it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 1261.91it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 6311.97it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 1386.36it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 10754.63it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 871.12it/s]
Rotating Exchange field: 100%|██████████| 3/3 [00:00<00:00, 92.28it/s]
Rotation 2: 100%|██████████| 100/100 [00:07<00:00, 12.63it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 505.36it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 4396.55it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 995.65it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 6875.91it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 1169.20it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 357.51it/s]
Rotation 3: 100%|██████████| 100/100 [00:07<00:00, 13.30it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 815.44it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 6927.01it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 892.97it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 5194.18it/s]
Extracting exchange field: 100%|██████████| 101/101 [00:00<00:00, 1119.77it/s]
Setup perturbations for rotated hamiltonian: 100%|██████████| 2/2 [00:00<00:00, 7430.12it/s]
Now the instance parameters are filled with the calculated exchane and on-site anisotropy parameters. We can print the instance CrI3 to get the information of the simulation and use the grogupy.save_magnopy() function to return the results in magnopy’s input format.
[10]:
print(CrI3)
--------------------------------------------------------------------------------
Metadata
grogupy version: 0.3.1
Architecture: CPU
SLURM job ID: Could not be determined.
Input file: /Users/danielpozsar/Documents/Documents_MacBook Air/studies/elte/phd/grogu/benchmarks/CrI3/CrI3.fdf
--------------------------------------------------------------------------------
Hamiltonian
Spin model: generalised-grogu
Spin mode: SPIN-ORBIT
Number of orbitals: 216
Unique ID: 8252670334522874585 3429370498709139588
--------------------------------------------------------------------------------
Solver
Number of threads in the parallel cluster: 1
Parallelization is over: No parallelization
Solver used for Greens function calculation: parallel
Maximum number of Greens function samples per batch: 50
Low memory mode: False
--------------------------------------------------------------------------------
Cell [Ang]
7.00000000e+00 0.00000000e+00 0.00000000e+00
-3.50000000e+00 6.06217782e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 3.00000000e+01
--------------------------------------------------------------------------------
Exchange field rotations
DFT axis: [0. 0. 1.]
Quantization axis and perpendicular directions:
1.00000000e+00 0.00000000e+00 0.00000000e+00 -->
0.00000000e+00 0.00000000e+00 -1.00000000e+00
0.00000000e+00 1.00000000e+00 0.00000000e+00
0.00000000e+00 7.07106781e-01 -7.07106781e-01
0.00000000e+00 1.00000000e+00 0.00000000e+00 -->
1.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 0.00000000e+00 -1.00000000e+00
7.07106781e-01 0.00000000e+00 -7.07106781e-01
0.00000000e+00 0.00000000e+00 1.00000000e+00 -->
1.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 1.00000000e+00 0.00000000e+00
7.07106781e-01 7.07106781e-01 0.00000000e+00
--------------------------------------------------------------------------------
Kspace
Total number of k points: 100
K points in each directions: [10 10 1]
--------------------------------------------------------------------------------
Contour
Eset: 50
Esetp: 100
Ebot: -35
Etop: 0
--------------------------------------------------------------------------------
[11]:
print(grogupy.save_magnopy(CrI3, comments=False))
================================================================================
GROGU INFORMATION
================================================================================
Hamiltonian convention
Double counting true
Normalized spins true
Intra-atomic factor +1
Exchange factor +0.5
================================================================================
Cell (Ang)
7.00000000 0.00000000 0.00000000
-3.50000000 6.06217782 0.00000000
0.00000000 0.00000000 30.00000000
================================================================================
Magnetic sites
Number of sites 2
Name x (Ang) y (Ang) z (Ang) s sx sy sz
0Cr(l:2) -0.00005220 0.00003006 15.00000000 1.74787087 0.00002581 0.02469071 0.99969514
1Cr(l:2) 3.50005256 2.02069554 14.99999997 1.74787735 0.00002145 -0.02470149 0.99969487
================================================================================
Intra-atomic anisotropy tensor (meV)
--------------------------------------------------------------------------------
0Cr(l:2)
Matrix
0.42677015 -0.04533856 -0.05362904
-0.04533856 0.51904088 -0.06215431
-0.05362904 -0.06215431 0.00000000
--------------------------------------------------------------------------------
1Cr(l:2)
Matrix
0.41333470 0.04549295 -0.04741143
0.04549295 0.51906709 0.06198310
-0.04741143 0.06198310 0.00000000
--------------------------------------------------------------------------------
================================================================================
Exchange tensor (meV)
Number of pairs 2
--------------------------------------------------------------------------------
Name1 Name2 i j k d (Ang)
--------------------------------------------------------------------------------
0Cr(l:2) 1Cr(l:2) 0 0 0 4.041512368971766
Matrix
-1.82204027 -0.28722064 -0.31246363
-0.72951690 -1.26385375 0.29209800
-0.21129252 0.51775736 -1.66394270
--------------------------------------------------------------------------------
0Cr(l:2) 1Cr(l:2) 1 0 0 10.69276805834447
Matrix
-0.02851866 -0.00233102 0.00177980
0.00205790 -0.02548159 0.00216193
-0.00146106 -0.00290484 -0.02175032
--------------------------------------------------------------------------------
================================================================================
More information can be found in the Tutorials.