README.md

@@ -3,12 +3,7 @@ More on the theoretical background can be seen on [arXiv](https://arxiv.org/abs/
 
 # TODO
 ## Testing
-- Test on lat3_791
+- Run tests on different magnetic materials and compare it to Grogu Matlab --> ran on Jij_for_Marci_6p45ang, but I could not compare data
-- Test on varcell
-- Test on varcell_myrun
-- Test on varcell_myrun_2
-- Test scaling
-- Test rotations to x-y 30 and 60
 
 ## Developing
 - Magnetic entities cannot handle orbitals, only shells
@@ -17,7 +12,6 @@ More on the theoretical background can be seen on [arXiv](https://arxiv.org/abs/
 - Check if exchange field has scalar part
 - Add more tests!!
 - io stuff
-- logging
 
 # Building wheel
 See detailed documentation on [PYPI](https://packaging.python.org/en/latest/tutorials/packaging-projects/).

src/grogupy/__init__.py

@@ -23,7 +23,6 @@
 
 
 from grogupy.core import *
-from grogupy.globals import *
 from grogupy.io import *
 from grogupy.magnetism import *
 from grogupy.utilities import *

src/grogupy/globals.py

@@ -1,61 +0,0 @@
-# Copyright (c) [2024] []
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
-
-import numpy as np
-
-# Pauli matrices
-TAU_X = np.array([[0, 1], [1, 0]])
-TAU_Y = np.array([[0, -1j], [1j, 0]])
-TAU_Z = np.array([[1, 0], [0, -1]])
-TAU_0 = np.array([[1, 0], [0, 1]])
-
-
-# list of accepted input parameters
-ACCEPTED_INPUTS = [
-    "infile",
-    "outfile",
-    "scf_xcf_orientation",
-    "ref_xcf_orientations",
-    "kset",
-    "kdirs",
-    "ebot",
-    "eset",
-    "esetp",
-    "parallel_solver_for_Gk",
-    "padawan_mode",
-]
-
-DEFAULT_ARGUMENTS = dict(
-    infile=None,
-    outfile=None,
-    scf_xcf_orientation=np.array([0, 0, 1]),
-    ref_xcf_orientations=[
-        dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),
-        dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),
-        dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),
-    ],
-    kset=2,
-    kdirs="xyz",
-    ebot=None,
-    eset=42,
-    esetp=1000,
-    parallel_solver_for_Gk=False,
-    padawan_mode=True,
-)

src/grogupy/grogu.py

@@ -72,6 +72,10 @@ def parse_command_line():
         type=str,
         help="Output file name",
     )
+
+    #    parser.add_argument('--scf-orientation', dest = 'scf_xcf_orientation',    default = None, help = 'Output file name')
+    #    parser.add_argument('--ref-orientation', dest = 'ref_xcf_orientations',   default = None, help = 'Output file name')
+
     parser.add_argument(
         "--kset",
         dest="kset",
@@ -153,27 +157,21 @@ def main(simulation_parameters, magnetic_entities, pairs):
         except:
             print("numpy version unknown.")
 
-    # rename outfile
+        # rename outfile
-    if not simulation_parameters["outfile"].endswith(".pickle"):
+        if not simulation_parameters["outfile"].endswith(".pickle"):
-        simulation_parameters["outfile"] += ".pickle"
+            simulation_parameters["outfile"] += ".pickle"
-
+
-    # if ebot is not given put it 1 eV under the smallest energy
+        # if ebot is not given put it 0.1 eV under the smallest energy
-    if simulation_parameters["ebot"] is None:
+        if simulation_parameters["ebot"] is None:
-        try:
+            try:
-            eigfile = simulation_parameters["infile"][:-3] + "EIG"
+                eigfile = simulation_parameters["infile"][:-3] + "EIG"
-            simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 1
+                simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1
-            simulation_parameters["automatic_ebot"] = True
+                simulation_parameters["automatic_ebot"] = True
-        except:
+            except:
-            print("Could not determine ebot.")
+                print("Could not determine ebot.")
-            print("Parameter was not given and .EIG file was not found.")
+                print("Parameter was not given and .EIG file was not found.")
-    else:
+        else:
-        simulation_parameters["automatic_ebot"] = False
+            simulation_parameters["automatic_ebot"] = False
-
-    # add third direction for off-diagonal anisotropy
-    for orientation in simulation_parameters["ref_xcf_orientations"]:
-        o1 = orientation["vw"][0]
-        o2 = orientation["vw"][1]
-        orientation["vw"].append((o1 + o2) / np.sqrt(2))
 
     # read sile
     fdf = sisl.get_sile(simulation_parameters["infile"])
@@ -218,7 +216,7 @@ def main(simulation_parameters, magnetic_entities, pairs):
         ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2
 
     # identifying TRS and TRB parts of the Hamiltonian
-    TAUY = np.kron(np.eye(NO), TAU_Y)
+    TAUY = np.kron(np.eye(NO), tau_y)
     hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())])
     hTRS = (hh + hTR) / 2
     hTRB = (hh - hTR) / 2
@@ -296,7 +294,7 @@ def main(simulation_parameters, magnetic_entities, pairs):
         R = RotMa2b(simulation_parameters["scf_xcf_orientation"], orient["o"])
         rot_XCF = np.einsum("ij,jklm->iklm", R, XCF)
         rot_H_XCF = sum(
-            [np.kron(rot_XCF[i], tau) for i, tau in enumerate([TAU_X, TAU_Y, TAU_Z])]
+            [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]
         )
         rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]
 
@@ -468,8 +466,8 @@ def main(simulation_parameters, magnetic_entities, pairs):
 
         # calculate magnetic parameters
         for mag_ent in magnetic_entities:
-            K, consistency = calculate_anisotropy_tensor(mag_ent)
+            Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent)
-            mag_ent["K"] = K * sisl.unit_convert("eV", "meV")
+            mag_ent["K"] = np.array([Kxx, Kyy, Kzz]) * sisl.unit_convert("eV", "meV")
             mag_ent["K_consistency"] = consistency
 
         for pair in pairs:
@@ -512,17 +510,17 @@ def main(simulation_parameters, magnetic_entities, pairs):
 if __name__ == "__main__":
     # loading parameters
     # it is not clear how to give grogu.fdf path...
-    # command_line_arguments = parse_command_line()
+    command_line_arguments = parse_command_line()
-    # fdf_arguments, magnetic_entities, pairs = read_fdf(command_line_arguments["infile"])
+    fdf_arguments, magnetic_entities, pairs = read_fdf(command_line_arguments["infile"])
 
     # right now we do not use any of these input, but it shows
     # the order of priority when processing arguments
     default_arguments = False
     fdf_arguments = False
     command_line_arguments = False
-    # simulation_parameters = process_input_args(
+    simulation_parameters = process_input_args(
-    #    default_arguments, fdf_arguments, command_line_arguments, ACCEPTED_INPUTS
+        default_arguments, fdf_arguments, command_line_arguments, ACCEPTED_INPUTS
-    # )
+    )
 
     ####################################################################################################
     #                                   This is the input file for now                                 #
@@ -544,19 +542,21 @@ if __name__ == "__main__":
         dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])),
     ]
     simulation_parameters = dict()
-    simulation_parameters["infile"] = ""
+    simulation_parameters["infile"] = (
-    simulation_parameters["outfile"] = "./"
+        "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf"
+    )
+    simulation_parameters["outfile"] = "./Fe3GeTe2_notebook"
     simulation_parameters["scf_xcf_orientation"] = np.array([0, 0, 1])
     simulation_parameters["ref_xcf_orientations"] = [
         dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),
         dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),
         dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),
     ]
-    simulation_parameters["kset"] = 0
+    simulation_parameters["kset"] = 20
-    simulation_parameters["kdirs"] = ""
+    simulation_parameters["kdirs"] = "xy"
     simulation_parameters["ebot"] = None
-    simulation_parameters["eset"] = 0
+    simulation_parameters["eset"] = 600
-    simulation_parameters["esetp"] = 0
+    simulation_parameters["esetp"] = 10000
     simulation_parameters["parallel_solver_for_Gk"] = False
     simulation_parameters["padawan_mode"] = True
     ####################################################################################################

src/grogupy/io.py

@@ -24,9 +24,41 @@
 from pickle import dump, load
 
 import numpy as np
-from globals import ACCEPTED_INPUTS, DEFAULT_ARGUMENTS
 from sisl.io import fdfSileSiesta
 
+# list of accepted input parameters
+ACCEPTED_INPUTS = [
+    "infile",
+    "outfile",
+    "scf_xcf_orientation",
+    "ref_xcf_orientations",
+    "kset",
+    "kdirs",
+    "ebot",
+    "eset",
+    "esetp",
+    "parallel_solver_for_Gk",
+    "padawan_mode",
+]
+
+default_arguments = dict(
+    infile=None,
+    outfile=None,
+    scf_xcf_orientation=np.array([0, 0, 1]),
+    ref_xcf_orientations=[
+        dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),
+        dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),
+        dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),
+    ],
+    kset=2,
+    kdirs="xyz",
+    ebot=None,
+    eset=42,
+    esetp=1000,
+    parallel_solver_for_Gk=False,
+    padawan_mode=True,
+)
+
 
 def read_fdf(path):
     """It reads the simulation parameters, magnetic entities and pairs from the fdf.
@@ -365,8 +397,7 @@ def print_atoms_and_pairs(magnetic_entities, pairs):
             # coordinates and tag
             print(f"{tag}        {xyz[0]}        {xyz[1]}        {xyz[2]}")
             print("Consistency check: ", mag_ent["K_consistency"])
-            print("K:        # Kxx, Kxy, Kxz, Kyx, Kyy, Kyz, Kzx, Kzy, Kzz")
+            print("Anisotropy diag: ", mag_ent["K"])
-            print(mag_ent["K"])
         print("")
 
     print(

src/grogupy/magnetism.py

@@ -117,6 +117,8 @@ def parse_magnetic_entity(dh, atom=None, l=None, **kwargs):
 def calculate_anisotropy_tensor(mag_ent):
     """Calculates the renormalized anisotropy tensor from the energies.
 
+    It uses the grogu convention for output.
+
     Args:
         mag_ent : dict
             An element from the magnetic entities
@@ -124,34 +126,22 @@ def calculate_anisotropy_tensor(mag_ent):
     Returns:
         K : np.array_like
             elements of the anisotropy tensor
-        consistency_check : float
-            Absolute value of the difference from the consistency check
     """
 
     # get the energies
     energies = mag_ent["energies"]
 
-    K = np.zeros((3, 3))
-
     # calculate the diagonal tensor elements
-    K[0, 0] = energies[1, 1] - energies[1, 0]
+    Kxx = energies[1, 1] - energies[1, 0]
-    K[1, 1] = energies[0, 1] - energies[0, 0]
+    Kyy = energies[0, 1] - energies[0, 0]
-    K[2, 2] = 0
+    Kzz = 0
-
-    # calculate the off-diagonal tensor elements
-    K[0, 1] = (energies[2, 0] + energies[2, 1]) / 2 - energies[2, 2]
-    K[1, 0] = K[0, 1]
-    K[0, 2] = (energies[1, 0] + energies[1, 1]) / 2 - energies[1, 2]
-    K[2, 0] = K[0, 2]
-    K[1, 2] = (energies[0, 0] + energies[0, 1]) / 2 - energies[0, 2]
-    K[2, 1] = K[1, 2]
 
     # perform consistency check
-    calculated_diff = K[1, 1] - K[0, 0]
+    calculated_diff = Kyy - Kxx
     expected_diff = energies[2, 0] - energies[2, 1]
     consistency_check = abs(calculated_diff - expected_diff)
 
-    return K, consistency_check
+    return Kxx, Kyy, Kzz, consistency_check
 
 
 def calculate_exchange_tensor(pair):

src/grogupy/utilities.py

@@ -22,10 +22,15 @@
 """
 
 import numpy as np
-from globals import TAU_0, TAU_X, TAU_Y, TAU_Z
 from scipy.special import roots_legendre
 from sisl.io.siesta import eigSileSiesta
 
+# Pauli matrices
+tau_x = np.array([[0, 1], [1, 0]])
+tau_y = np.array([[0, -1j], [1j, 0]])
+tau_z = np.array([[1, 0], [0, -1]])
+tau_0 = np.array([[1, 0], [0, 1]])
+
 
 def commutator(a, b):
     """Shorthand for commutator.
@@ -184,7 +189,7 @@ def tau_u(u):
     # u is force to be of unit length
     u = u / np.linalg.norm(u)
 
-    return u[0] * TAU_X + u[1] * TAU_Y + u[2] * TAU_Z
+    return u[0] * tau_x + u[1] * tau_y + u[2] * tau_z
 
 
 #