Define some data structures used inside C code (#41)

Co-authored-by: Simone Baffelli <[email protected]>

.flake8

@@ -6,3 +6,6 @@ extend-ignore =
     E203
 exclude =
     docs/source/conf.py
+
+per-file-ignores =
+    cleedpy/interface/cleed.py: F841, TRY003

.gitignore

@@ -213,3 +213,6 @@ dkms.conf
 
 # CMake builds
 cmake-build*/
+
+# Vscode
+.vscode/

cleedpy/cleed/src/leed.c

@@ -35,6 +35,10 @@ real ** leed(
 
     struct eng_str eng;
 
+    // Printing stuff
+    inp_showbop(bulk, over, phs_shifts);
+    return 0;
+
     eng.ini = energy_list[0];
     eng.stp = energy_list[1] - energy_list[0];
     eng.fin = energy_list[energy_list_size-1];
@@ -235,3 +239,15 @@ real ** leed(
 
     return iv_curves;
 }
+
+
+int my_test_function(int a, int b, struct cryst_str * bulk){
+    printf("vr: %lf\n", bulk->vr);
+    printf("vi: %lf\n", bulk->vi);
+
+    printf("\nbulk 2-dim. unit cell:\n");
+    printf("a1:  (%7.4lf  %7.4lf)\n", bulk->a[1]*BOHR, bulk->a[3]*BOHR);
+    printf("a2:  (%7.4lf  %7.4lf)\n", bulk->a[2]*BOHR, bulk->a[4]*BOHR);
+
+    return a+b;
+}

cleedpy/cleed/test_cleed.c

@@ -240,14 +240,6 @@ FILE *res_stream;
   mk_cg_coef (2*v_par->l_max);
   mk_ylm_coef(2*v_par->l_max);
 
-#ifdef CONTROL
-  fprintf(STDCTR, "(LEED_TEMP): E_ini = %.1f, E_fin = %.1f, E_stp %.1f\n",
-          eng->ini*HART, eng->fin*HART, eng->stp*HART);
-
-  //fprintf(STDCTR, "(LEED_TEMP): lset = %d\n", n_set);
-#endif
-
-
   // Construct energy list
   energy_list_size = (eng->fin - eng->ini)/eng->stp + 1;
   energy_list = (real *) malloc(energy_list_size * sizeof(real));

cleedpy/interface/cleed.py

@@ -0,0 +1,209 @@
+import math
+from ctypes import POINTER, Structure, c_char_p, c_double, c_int, cdll
+from dataclasses import dataclass
+from pathlib import Path
+
+from ..config import (
+    AtomParametersStructured,
+    AtomParametersVariants,
+    Position,
+    PositionOptimizationParameters,
+    SearchParameters,
+)
+from ..physics import constants
+from .matrix import MatPtr
+
+
+class Atom(Structure):
+    _fields_ = [
+        ("layer", c_int),
+        ("type", c_int),
+        ("t_type", c_int),
+        ("pos", c_double * 4),
+        ("dwf", c_double),
+    ]
+
+
+class Layer(Structure):
+    _fields_ = [
+        ("no_of_layer", c_int),
+        ("bulk_over", c_int),
+        ("periodic", c_int),
+        ("natoms", c_int),
+        ("a_lat", c_double * 5),
+        ("rel_area", c_double),
+        ("reg_shift", c_double * 4),
+        ("vec_from_last", c_double * 4),
+        ("vec_to_next", c_double * 4),
+        ("atoms", POINTER(Atom)),
+    ]
+
+
+class Crystal(Structure):
+    _fields_ = [
+        ("vr", c_double),
+        ("vi", c_double),
+        ("temp", c_double),
+        # Not used anymore
+        # From here
+        ("n_rot", c_int),
+        ("rot_axis", c_double * 4),
+        ("n_mir", c_int),
+        ("m_plane", POINTER(c_double)),
+        ("alpha", POINTER(c_double)),
+        ("symmetry", c_int),
+        # To here
+        ("a", c_double * 5),
+        ("a_1", c_double * 5),
+        ("area", c_double),
+        ("m_trans", c_double * 5),
+        ("m_super", c_double * 5),
+        ("m_recip", c_double * 5),
+        ("b", c_double * 5),
+        ("b_1", c_double * 5),
+        ("rel_area_sup", c_double),
+        ("nlayers", c_int),
+        ("layers", POINTER(Layer)),
+        ("dmin", c_double),
+        ("natoms", c_int),
+        ("ntypes", c_int),
+        ("comments", POINTER(c_char_p)),
+    ]
+
+
+class PhaseShifts(Structure):
+    _fields_ = [
+        ("lmax", c_int),
+        ("neng", c_int),
+        ("t_type", c_int),
+        ("eng_max", c_double),
+        ("eng_min", c_double),
+        ("energy", POINTER(c_double)),
+        ("pshift", POINTER(c_double)),
+        ("dr", c_double * 4),
+        ("input_file", c_char_p),
+    ]
+
+
+class LeedResult(Structure):
+    _fields_ = [("x", c_double), ("y", c_double), ("result", c_double)]
+
+
+class EnergyLoopVariables(Structure):
+    _fields_ = [
+        ("eng_r", c_double),
+        ("eng_i", c_double),
+        ("eng_v", c_double),
+        ("vr", c_double),
+        ("vi_prev", c_double),
+        ("vi_exp", c_double),
+        ("theta", c_double),
+        ("phi", c_double),
+        ("k_in", c_double * 4),
+        ("epsilon", c_double),
+        ("l_max", c_int),
+        ("p_t1", POINTER(MatPtr)),
+    ]
+
+
+@dataclass
+class CleedInputs:
+    bulk: Crystal
+    overlayers: Crystal
+    energies: list[float]
+    phase_shifts: PhaseShifts
+    params: EnergyLoopVariables
+
+
+def generate_energies(inp: list[float]) -> list[c_double]:
+    return [c_double(i) for i in inp]
+
+
+def convert_energy_loop_variables(inp: SearchParameters) -> EnergyLoopVariables:
+    """
+    This corresponds to the following c function: inp_rdpar
+    """
+    return EnergyLoopVariables(
+        eng_r=c_double(inp.energy_range.initial / constants.HART),
+        eng_i=c_double(inp.energy_range.final / constants.HART),
+        eng_v=c_double(inp.energy_range.step / constants.HART),
+        vr=c_double(inp.optical_potential[0]),
+        vi_prev=c_double(inp.optical_potential[1]),
+        vi_exp=c_double(inp.optical_potential[1]),
+        theta=c_double(math.radians(inp.polar_incidence_angle)),
+        phi=c_double(math.radians(inp.azimuthal_incidence_angle)),
+        k_in=(c_double * 4)(0, 0, 0, 0),
+        epsilon=c_double(inp.epsilon),
+        l_max=c_int(inp.maximum_angular_momentum),
+        p_t1=POINTER(MatPtr),
+    )
+
+
+def generate_single_atom_structure(inp: AtomParametersVariants) -> Atom:
+    """
+    Generate a single atom for cleed from the inp atom parameters
+    """
+    match inp:
+        case AtomParametersStructured(
+            phase_file, Position(x, y, z), vibrational_displacement
+        ):
+            a = Atom()
+            a.layer = 0
+            a.type = 0
+            a.t_type = 0
+            a.pos = (x, y, z, 0)
+            a.dwf
+
+        case AtomParametersStructured(
+            phase_file, PositionOptimizationParameters(), vibrational_displacement
+        ):
+            raise ValueError(
+                "Cleed needs a position for each atom, not an optimization range"
+            )
+        case _:
+            raise ValueError("Not implemented")
+
+
+def generate_atom_structure(inp: list[AtomParametersVariants]) -> list[Atom]:
+    """
+    Generate a list of atoms for cleed from the inp list of atom parameters
+    It also sorts the atoms by the z-coordinate
+    """
+    return [Atom() for i in inp]
+
+
+def generate_crystal_structure(inp: AtomParametersVariants) -> Crystal:
+    match inp:
+        case AtomParametersStructured(
+            phase_file, Position(x, y, z), vibrational_displacement
+        ):
+            pass
+        case AtomParametersStructured(
+            phase_file, PositionOptimizationParameters(), vibrational_displacement
+        ):
+            raise ValueError(
+                "Cleed needs a position for each atom, not an optimization range"
+            )
+        case _:
+            raise ValueError("Not implemented")
+
+
+def call_cleed():
+    path = Path(__file__).parent.parent / "cleed" / "build" / "lib" / "libcleed.dylib"
+    lib = cdll.LoadLibrary(str(path))
+
+    lib.my_test_function.argtypes = [c_int, c_int, POINTER(Crystal)]
+    lib.my_test_function.restype = c_int
+
+    bulk = Crystal()
+    bulk.vr = 1.0
+    bulk.vi = 2.0
+    bulk.temp = 300.0
+
+    bulk.a = (0.0, 1.0, 2.0, 3.0, 4.0)
+
+    print("Restul is ", lib.my_test_function(1, 2, bulk))
+
+
+if __name__ == "__main__":
+    call_cleed()

cleedpy/interface/matrix.py

@@ -0,0 +1,17 @@
+from ctypes import POINTER, Structure, c_double, c_int
+
+
+class Matrix(Structure):
+    _fields_ = [
+        ("mag_no", c_int),
+        ("blk_type", c_int),
+        ("mat_type", c_int),
+        ("num_type", c_int),
+        ("rows", c_int),
+        ("cols", c_int),
+        ("rel", POINTER(c_double)),
+        ("img", POINTER(c_double)),
+    ]
+
+
+MatPtr = POINTER(Matrix)

cleedpy/physics/constants.py

@@ -0,0 +1,5 @@
+HART = 27.2113962  # Hartree in eV
+BOHR = 0.529177249  # Bohr radius in Angstroms
+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

examples/cleed_example_old/run.sh

@@ -1,4 +1,4 @@
 export CLEED_PHASE=./PHASE
-EXECUTABLE=../../cleedpy/cleed/build/test_cleed
+EXECUTABLE=../../cleedpy/cleed/build/bin/test_cleed
 
 ${EXECUTABLE} -i Ni111_Cu.inp