Update

pyproject.toml

@@ -0,0 +1,44 @@
+[build-system]
+requires = [
+    "setuptools>=42",
+    "wheel",
+    "pybind11"
+]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "undi"
+description = "Muon-Nuclear Dipolar Interaction"
+version = "1.1"
+authors = [
+    { name = "Pietro Bonfa", email = "bonfus@gmail.com" }
+]
+readme = "README.rst"
+license = { file = "LICENSE" }
+classifiers = [
+    "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
+    "Programming Language :: Python :: 3 :: Only",
+    "Operating System :: OS Independent",
+    "Intended Audience :: Science/Research",
+    "Topic :: Scientific/Engineering :: Physics"
+]
+keywords = ["muSR spin nuclear"]
+requires-python = ">=3.12"
+dependencies = [
+    "numpy>=2.0",
+    "qutip",
+    "pybind11"
+]
+
+[tool.setuptools.packages.find]
+include = ["undi*"]
+
+[tool.setuptools.package-data]
+"undi.kubo_toyabe" = ["*.txt"]
+
+[project.optional-dependencies]
+fast_cpp = ["pybind11"]
+progress_bar = ["tqdm"]
+
+[project.urls]
+Source = "http://github.com/bonfus/undi"

setup.py

@@ -1,4 +1,4 @@
-from setuptools import setup
+from setuptools import setup, find_packages
 import os
 from glob import glob
 
@@ -62,11 +62,12 @@
       author='Pietro Bonfa',
       author_email='bonfus@gmail.com',
       license='GPLv3',
-      packages=['undi', 'undi.fast'],
+      packages=['undi', 'undi.fast', 'undi.kubo_toyabe'],
+      package_data={'undi.kubo_toyabe': ['*.txt']},
       include_package_data=True,
       install_requires = [
                           'numpy',
-                          'qutip'
+                          'qutip',
                           ],
       extras_require = {
         'Fast C++ implementation':  ["pybind11"],

undi/__main__.py

@@ -54,7 +54,7 @@ def main():
     - structure (str): Structure to be parsed.
     - --max-hdim (int, optional): Dimension of the Hilbert space. Increase to have better results. Default is 1000.
     - --atom-as-muon (str, optional): Atom to be considered as muon. Default is "H".
-    - --Bmod (float, optional): Magnetic field modulus. Default is 0.0.
+    - --B_mod (float, optional): Magnetic field modulus. Default is 0.0.
     - --convergence-check (bool, optional): Flag to check for convergence. Default is False.
     - --algorithm (str, optional): Algorithm to be used for the analysis. Default is "fast".
     - --sample_size_average (int, optional): Number of random samples to for powder average. Default is 1000.
@@ -67,15 +67,15 @@ def main():
     parser.add_argument('structure', type=str, help='Structure to be parsed')
     parser.add_argument('--max-hdim', type=int, default=1e3, help='Dimension of the Hilbert space. Increase to have better results. Default is 1000.')
     parser.add_argument('--atom-as-muon', type=str, default="H", help='Atom to be considered as muon. Default is "H".')
-    parser.add_argument('--Bmod', type=float, default=0., help='Magnetic field modulus. Default is 0.0.')
+    parser.add_argument('--B_mod', type=float, default=0., help='Magnetic field modulus. Default is 0.0.')
     parser.add_argument('--convergence-check', type=bool, default=False, help='Flag to check for convergence. Default is False.')
     parser.add_argument('--algorithm', type=str, default="fast", help='Algorithm to be used for the analysis. Default is "fast".')
     parser.add_argument('--sample_size_average', type=int, default=1000, help='Number of random samples to for powder average. Default is 1000.')
     parser.add_argument('--dump', type=bool, default=True, help='Flag to dump the results. Default is True.')
 
     args = parser.parse_args()
     structure = args.structure
-    Bmod = args.Bmod
+    B_mod = args.B_mod
     convergence_check = args.convergence_check
     atom_as_muon = args.atom_as_muon
     max_hdim = args.max_hdim
@@ -91,7 +91,7 @@ def main():
 
     results = execute_undi_analysis(
         structure=atms,
-        Bmod=Bmod,
+        B_mod=B_mod,
         atom_as_muon=atom_as_muon,
         max_hdim=max_hdim,
         convergence_check=convergence_check,

undi/isotopes.py

@@ -478,10 +478,10 @@ class Element:
     #    ],
     #)
     info = DotArray(np.loadtxt(easyspin_isotope_info, comments='%',
-                                dtype=[('Z', 'int32'), ('A', 'int32'), ('Stable', 'U1'),
-                                ('Symbol', 'U2'), ('Element', 'U10'), ('Spin', 'float'),
-                                ('G_factor', 'float'), ('Abundance', 'float'),
-                                ('Quadrupole', 'float') ]))
+                            dtype=[('Z', np.int32), ('A', np.int32), ('Stable', 'U1'),
+                            ('Symbol', 'U2'), ('Element', 'U10'), ('Spin', np.float64),
+                            ('G_factor', np.float64), ('Abundance', np.float64),
+                            ('Quadrupole', np.float64) ]))
 
     def __init__(self, symbol):
         self.isotopes = []

undi/kubo_toyabe/KT.py

@@ -0,0 +1,87 @@
+import numpy as np
+import pandas as pd
+from ase import neighborlist
+from importlib_resources import files
+
+from undi import kubo_toyabe as isotopedata
+
+file = files(isotopedata) / "isotopedata.txt"
+
+info = pd.read_table(
+    file,
+    comment="%",
+    sep="\s+",
+    names=[
+        "Z",
+        "A",
+        "Stable",
+        "Symbol",
+        "Element",
+        "Spin",
+        "G_factor",
+        "Abundance",
+        "Quadrupole",
+    ],
+)
+
+munhbar = 7.622593285e6 * 2 * np.pi  # mu_N/hbar, SI
+# (2/3)(μ_0/4pi)^2 (planck2pi 2pi × 135.5 MHz/T )^2 = 5.374 021 39 × 10^(−65) kg²·m^(6)·A^(−2)·s^(−4)
+factor = 5.37402139e-5  # angstrom instead of m
+
+
+def get_isotopes(Z):
+    return info[info.Z == Z][["Abundance", "Spin", "G_factor"]].to_numpy()
+
+
+def compute_second_moments(atms, cutoff_distances={}):
+    """
+    Compute second moments taking care of isotope averages
+    """
+    tot_H = np.count_nonzero(atms.get_atomic_numbers() == 1)
+
+    species_avg = {}
+    for e in np.unique(atms.get_atomic_numbers()):
+        if e == 1:
+            continue
+
+        species_avg[e] = 0.0
+        for a in get_isotopes(e):
+            species_avg[e] += (a[0] / 100) * a[1] * (a[1] + 1) * (munhbar * a[2]) ** 2
+
+    # compute second moments
+    specie_contribs = {}
+    for e in np.unique(atms.get_atomic_numbers()):
+        if e == 1:
+            continue
+        sum = 0.5 * np.sum(
+            neighborlist.neighbor_list(
+                "d", atms, cutoff={(1, e): cutoff_distances.get(e, 40)}
+            )
+            ** -6
+        )
+        specie_contribs[e] = species_avg[e] * sum * factor / tot_H
+
+    return specie_contribs
+
+
+def kubo_toyabe(tlist, Gmu_S2):
+    """Calculates the Kubo-Toyabe polarization for the nuclear arrangement
+    provided in input.
+
+    Parameters
+    ----------
+    tlist : numpy.array
+        List of times at which the muon polarization is observed.
+
+    Returns
+    -------
+    numpy.array
+        Kubo-Toyabe function, for a powder in zero field.
+    """
+    # this is gamma_mu times sigma^2
+    return 0.333333333333 + 0.6666666666 * (1 - Gmu_S2 * np.power(tlist, 2)) * np.exp(
+        -0.5 * Gmu_S2 * np.power(tlist, 2)
+    )
+
+
+#### end for KT