Towards full Python interface to the core C code of CLEED.

1. New cmdline script called `cleedpy-leed` and provide initial implementation.
2. Small modifications to the config module.
3. Provide preprocess module to prepare the inputs.
4. Move `mk_cg_coef` and `mk_ylm_coef` functions from test_cleed.c to leed.c.
5. Provide a working example of yaml input for the code.

cleedpy/cleed/src/leed.c

@@ -36,9 +36,15 @@ real ** leed(
     struct eng_str eng;
 
     // Printing stuff
+    printf("TEST SASHA\n");
     inp_showbop(bulk, over, phs_shifts);
+    printf("END TEST SASHA\n");
     return 0;
 
+
+    mk_cg_coef (2*v_par->l_max); // Setting up Clebsh Gordan coefficients as global variables.
+    mk_ylm_coef(2*v_par->l_max); // Setting up spherical harmonics coefficients as global variables.
+
     eng.ini = energy_list[0];
     eng.stp = energy_list[1] - energy_list[0];
     eng.fin = energy_list[energy_list_size-1];

cleedpy/cleed/test_cleed.c

@@ -237,9 +237,6 @@ FILE *res_stream;
  Prepare some often used parameters.
 *********************************************************************/
 
-  mk_cg_coef (2*v_par->l_max);
-  mk_ylm_coef(2*v_par->l_max);
-
   // Construct energy list
   energy_list_size = (eng->fin - eng->ini)/eng->stp + 1;
   energy_list = (real *) malloc(energy_list_size * sizeof(real));

cleedpy/cli/generate_schema.py

@@ -3,12 +3,12 @@
 
 import typer
 
-from ..config import NonGeometricalParameters
+from ..config import InputParameters
 
 
 def generate_schema(path: Path):
     """Generate schema for the LEED calculations"""
-    schema = NonGeometricalParameters.model_json_schema()
+    schema = InputParameters.model_json_schema()
 
     with open(path, "w") as f:
         f.write(json.dumps(schema, indent=4))

cleedpy/cli/leed.py

@@ -0,0 +1,33 @@
+from pathlib import Path
+
+import typer
+
+from ..config import load_parameters
+from ..physics.constants import BOHR_TO_ANGSTROM
+from ..preprocess import extract_bulk_parameters, transform_cells
+
+
+def leed(inptut_file: Path):
+    """Leed CLI."""
+    config = load_parameters(inptut_file)
+    transformation_matrix = transform_cells(config)
+
+    bulk_parameters = extract_bulk_parameters(config, transformation_matrix)
+    print(
+        f"Optical potential: {bulk_parameters.vr} eV (Vr), {bulk_parameters.vi} eV (Vi)"
+    )
+    print(f"Temperature: {bulk_parameters.temp} K")
+    print(
+        f"""Bulk 2-dim. unit cell: {bulk_parameters.a[0]}
+    {bulk_parameters.a[1]*BOHR_TO_ANGSTROM} {bulk_parameters.a[3]*BOHR_TO_ANGSTROM}
+    {bulk_parameters.a[2]*BOHR_TO_ANGSTROM} {bulk_parameters.a[4]*BOHR_TO_ANGSTROM}"""
+    )
+
+
+def cli():
+    """Leed CLI."""
+    typer.run(leed)
+
+
+if __name__ == "__main__":
+    cli()

cleedpy/cli/validate_parameters.py

@@ -3,16 +3,16 @@
 import typer
 import yaml
 
-from ..config import NonGeometricalParameters
+from ..config import InputParameters
 
 
 def validate(path: Path):
     """Validate the input parameters"""
     with open(path) as f:
         data = yaml.load(f, Loader=yaml.FullLoader)
         print(data)
-        NonGeometricalParameters.model_validate(data)
-        obj = NonGeometricalParameters(**data)
+        InputParameters.model_validate(data)
+        obj = InputParameters(**data)
         print(obj)
 
 

cleedpy/config.py

@@ -1,6 +1,7 @@
 from pathlib import Path
 from typing import Any, Literal, NamedTuple
 
+import yaml
 from numpy.typing import ArrayLike
 from pydantic import BaseModel
 
@@ -16,8 +17,8 @@ class UnitCellParameters(BaseModel):
 class SuperstructureMatrix(BaseModel):
     """Superstructure matrix for bulk calculations"""
 
-    m1: tuple[int, int]
-    m2: tuple[int, int]
+    m1: tuple[float, float]
+    m2: tuple[float, float]
 
 
 class VibrationalDisplacementParameters(BaseModel):
@@ -30,7 +31,6 @@ class DebyeTemperatureMass(NamedTuple):
     version: Literal["dmt"]
     debye_temperature: float
     atomic_mass: float
-    temperature: float
 
 
 class RadialMeanSquareDisplacement(NamedTuple):
@@ -107,19 +107,28 @@ class EnergyRangeParameters(BaseModel):
     step: float
 
 
-class NonGeometricalParameters(BaseModel):
+class MinimumAtomRadius(NamedTuple):
+    """Minimum atom radius for search. This is used to avoid atoms overlapping."""
+
+    name: str
+    radius: float
+
+
+class InputParameters(BaseModel):
     """Non geometrical parameters for bulk calculations"""
 
     unit_cell: UnitCellParameters
     superstructure_matrix: SuperstructureMatrix
     overlayers: list[AtomParametersVariants]
     bulk_layers: list[AtomParametersVariants]
+    minimum_radius: list[MinimumAtomRadius]
     optical_potential: tuple[float, float] = (8, 4)
     energy_range: EnergyRangeParameters
     polar_incidence_angle: float = 0
     azimuthal_incidence_angle: float = 0
     epsilon: float = 1e-2
     maximum_angular_momentum: int = 8
+    sample_temperature: float = 300.0
 
 
 class SearchRadiusParameters(NamedTuple):
@@ -156,3 +165,11 @@ def model_post_init(self, __context: Any) -> None:
 
     def get_iv_curve(self) -> ArrayLike:
         """Get the IV curve from the experimental data"""
+
+
+def load_parameters(parameters_file: Path):
+    """Load the parameters file."""
+    with open(parameters_file) as f:
+        data = yaml.load(f, Loader=yaml.FullLoader)
+        InputParameters.model_validate(data)
+        return InputParameters(**data)

cleedpy/physics/constants.py

@@ -1,5 +1,6 @@
 HART = 27.2113962  # Hartree in eV
-BOHR = 0.529177249  # Bohr radius in Angstroms
+BOHR_TO_ANGSTROM = 0.529177249  # Bohr radius in Angstroms
+ANGSTROM_TO_BOHR = 1.0 / BOHR_TO_ANGSTROM
 MEL_U = 5.48579903e-4  # electron mass in amu (atomic mass units)
 U_MEL = 1822.88850636  # 1 amu in multiples of the electron mass
 KB = 3.16682941e-6  # Boltzmann constant in Hartree/K

cleedpy/preprocess.py

@@ -0,0 +1,84 @@
+import numpy as np
+
+from .interface.cleed import Crystal
+from .physics.constants import ANGSTROM_TO_BOHR, BOHR_TO_ANGSTROM
+
+GEO_TOLERANCE = 0.0001
+
+
+def transform_cells(parameters):
+    """Converto to right handed coordinate system."""
+    a1 = np.array(parameters.unit_cell.a1)
+    a2 = np.array(parameters.unit_cell.a2)
+    m_transformation = np.array([[1.0, 0], [0, 1.0]])
+
+    # Check the angle between a1 and a2 (1: must be >= pi/2, i.e. a1*a2 < 0).
+    # if not: take a1 and -a2 as basic vectors and modify transformation matrix.
+    if np.dot(a1, a2) > 0.0:
+        a2 *= -1.0
+        m_transformation[1] *= -1.0
+
+    # Check the angle between a1 and a2 (2: must be < pi, i.e. a1*a2 > 0)
+    # if not: exchange a1 and a2 and modify transformation matrix.
+    if np.cross(a1, a2)[2] < 0.0:
+        a1, a2 = a2, a1
+        m_transformation[0, 1] = m_transformation[1, 0]
+
+    parameters.unit_cell.a1, parameters.unit_cell.a2 = tuple(a1), tuple(a2)
+
+    return m_transformation
+
+
+def extract_bulk_parameters(parameters, transformation_matrix):
+    bulk = Crystal()
+
+    bulk.vr, bulk.vi = parameters.optical_potential
+    bulk.temp = parameters.sample_temperature
+
+    def to_matrix(v1, v2):
+        """Converts two vectors to a matrix."""
+        return np.array([v1, v2])
+
+    # Unit cell.
+    a = (
+        to_matrix(parameters.unit_cell.a1[:2], parameters.unit_cell.a2[:2])
+        * ANGSTROM_TO_BOHR
+    )
+    bulk.a = (0.0, a[0, 0], a[1, 0], a[0, 1], a[1, 1])
+
+    # Reciprocal lattice vectors.
+    a_recip = np.linalg.inv(a).T * 2.0 * np.pi
+    print(a_recip / BOHR_TO_ANGSTROM)
+    bulk.a_1 = (0.0, a_recip[0, 0], a_recip[0, 1], a_recip[1, 0], a_recip[1, 1])
+
+    # Superstructure matrix
+    superstructure_matrix = to_matrix(
+        parameters.superstructure_matrix.m1, parameters.superstructure_matrix.m2
+    )
+    superstructure_matrix @= transformation_matrix
+    bulk.m_super = (
+        0.0,
+        superstructure_matrix[0, 0],
+        superstructure_matrix[1, 0],
+        superstructure_matrix[0, 1],
+        superstructure_matrix[1, 1],
+    )
+    bulk.m_trans = (
+        0.0,
+        transformation_matrix[0, 0],
+        transformation_matrix[1, 0],
+        transformation_matrix[0, 1],
+        transformation_matrix[1, 1],
+    )
+
+    # Reciprocal superstructure matrix.
+    reciprocal_matrix = np.linalg.inv(superstructure_matrix).T
+    bulk.m_recip = (
+        0.0,
+        reciprocal_matrix[0, 0],
+        reciprocal_matrix[0, 1],
+        reciprocal_matrix[1, 0],
+        reciprocal_matrix[1, 1],
+    )
+
+    return bulk

examples/cleed_example_new/Ni111_Cu.yml

@@ -0,0 +1,46 @@
+unit_cell:
+  a1: [1.2450,  2.1564, 0.0]
+  a2: [1.2450, -2.1564, 0.0]
+  a3: [0.0000, 0.0000, -6.0990]
+superstructure_matrix:
+  m1: [1.0, 0]
+  m2: [0, 1.0]
+
+overlayers:
+  - phase_file: "Cu_BVH"
+    position: [0.0000, -0.0000, 6.0900]
+    vibrational_displacement: [dr3, 0.032, 0.032, 0.032]
+  - phase_file: "Ni_BVH"
+    position: [1.2450, -0.7188, 4.0600]
+    vibrational_displacement: [dr3, 0.025, 0.025, 0.025]
+  - phase_file: "Ni_BVH"
+    position: [1.2450, 0.7188, 2.0300]
+    vibrational_displacement: [dr3, 0.025, 0.025, 0.025]
+
+bulk_layers:
+  - phase_file: "Ni_BVH"
+    position: [0.0000, 0.0000,  0.0000]
+    vibrational_displacement: [dr3, 0.025, 0.025, 0.025]
+  - phase_file: "Ni_BVH"
+    position: [1.2450, -0.7188, -2.0300]
+    vibrational_displacement: [dr3, 0.025, 0.025, 0.025]
+  - phase_file: "Ni_BVH"
+    position: [1.2450, 0.7188, -4.0600]
+    vibrational_displacement: [dr3, 0.025, 0.025, 0.025]
+
+minimum_radius:
+  - name: Ni_BVH
+    radius: 0.9
+  - name: Cu_BVH
+    radius: 0.9
+
+optical_potential: [-8.00, 4.0]
+energy_range:
+  initial: 70.0
+  final: 498.0
+  step: 4.0
+
+polar_incidence_angle: 0.0
+azimuthal_incidence_angle: 0.0
+epsilon: 1.0e-2
+maximum_angular_momentum: 9