added globals, off-diagonal anysotropy and todo

README.md

@@ -3,7 +3,12 @@ More on the theoretical background can be seen on [arXiv](https://arxiv.org/abs/
 
 # TODO
 ## Testing
-- 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 lat3_791
+- 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
@@ -12,6 +17,7 @@ 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,6 +23,7 @@
 
 
 from grogupy.core import *
+from grogupy.globals import *
 from grogupy.io import *
 from grogupy.magnetism import *
 from grogupy.utilities import *

src/grogupy/globals.py

@@ -0,0 +1,61 @@
+# 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,10 +72,6 @@ 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",
@@ -161,11 +157,11 @@ def main(simulation_parameters, magnetic_entities, pairs):
     if not simulation_parameters["outfile"].endswith(".pickle"):
         simulation_parameters["outfile"] += ".pickle"
 
-        # if ebot is not given put it 0.1 eV under the smallest energy
+    # if ebot is not given put it 1 eV under the smallest energy
     if simulation_parameters["ebot"] is None:
         try:
             eigfile = simulation_parameters["infile"][:-3] + "EIG"
-                simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1
+            simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 1
             simulation_parameters["automatic_ebot"] = True
         except:
             print("Could not determine ebot.")
@@ -173,6 +169,12 @@ def main(simulation_parameters, magnetic_entities, pairs):
     else:
         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"])
 
@@ -216,7 +218,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
@@ -294,7 +296,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]
 
@@ -466,8 +468,8 @@ def main(simulation_parameters, magnetic_entities, pairs):
 
         # calculate magnetic parameters
         for mag_ent in magnetic_entities:
-            Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent)
+            K, consistency = calculate_anisotropy_tensor(mag_ent)
-            mag_ent["K"] = np.array([Kxx, Kyy, Kzz]) * sisl.unit_convert("eV", "meV")
+            mag_ent["K"] = K * sisl.unit_convert("eV", "meV")
             mag_ent["K_consistency"] = consistency
 
         for pair in pairs:
@@ -510,17 +512,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                                 #
@@ -552,7 +554,7 @@ if __name__ == "__main__":
         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"] = 20
+    simulation_parameters["kset"] = 10
     simulation_parameters["kdirs"] = "xy"
     simulation_parameters["ebot"] = None
     simulation_parameters["eset"] = 600

src/grogupy/io.py

@@ -24,41 +24,9 @@
 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.
@@ -397,7 +365,8 @@ 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("Anisotropy diag: ", mag_ent["K"])
+            print("K:        # Kxx, Kxy, Kxz, Kyx, Kyy, Kyz, Kzx, Kzy, Kzz")
+            print(mag_ent["K"])
         print("")
 
     print(

src/grogupy/magnetism.py

@@ -117,8 +117,6 @@ 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
@@ -126,22 +124,34 @@ 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
-    Kxx = energies[1, 1] - energies[1, 0]
+    K[0, 0] = energies[1, 1] - energies[1, 0]
-    Kyy = energies[0, 1] - energies[0, 0]
+    K[1, 1] = energies[0, 1] - energies[0, 0]
-    Kzz = 0
+    K[2, 2] = 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 = Kyy - Kxx
+    calculated_diff = K[1, 1] - K[0, 0]
     expected_diff = energies[2, 0] - energies[2, 1]
     consistency_check = abs(calculated_diff - expected_diff)
 
-    return Kxx, Kyy, Kzz, consistency_check
+    return K, consistency_check
 
 
 def calculate_exchange_tensor(pair):

src/grogupy/utilities.py

@@ -22,15 +22,10 @@
 """
 
 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.
@@ -189,7 +184,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
 
 
 #