diff --git a/cleedpy/cleed/CMakeLists.txt b/cleedpy/cleed/CMakeLists.txt
index 1139401..a193b87 100644
--- a/cleedpy/cleed/CMakeLists.txt
+++ b/cleedpy/cleed/CMakeLists.txt
@@ -29,9 +29,3 @@ include_directories(include)
 
 # Add the src directory
 add_subdirectory(src)
-
-# Build the test
-add_executable(test_cleed test_cleed.c)
-
-# Link the test
-target_link_libraries(test_cleed cleed)
diff --git a/cleedpy/cleed/include/leed_func.h b/cleedpy/cleed/include/leed_func.h
index 2f7a809..a29088e 100644
--- a/cleedpy/cleed/include/leed_func.h
+++ b/cleedpy/cleed/include/leed_func.h
@@ -90,11 +90,8 @@ int pc_mk_ms(mat * , mat *, mat *, mat *, mat *, mat *, int );
 /*********************************************************************
  Output
 *********************************************************************/
-int out_head(struct cryst_str *, FILE *);
-int out_head_2(struct cryst_str *, const char *, const char *, FILE *);
-int out_bmlist(struct beam_str **, struct beam_str *, struct eng_str *, FILE *);
-int out_int(mat , struct beam_str *, struct beam_str *, struct var_str *, FILE * );
-int out_intsym(mat , struct beam_str *, struct beam_str *, struct var_str *, FILE * );
+int out_bmlist(struct beam_str **, struct beam_str *, struct eng_str *, real **, real **, int **);
+int out_int(mat , struct beam_str *, struct beam_str *, struct var_str *, real *);
 
     /* check cpu time */
 double cpu_time(FILE *, const char *);
diff --git a/cleedpy/cleed/original_test_cleed.c b/cleedpy/cleed/original_test_cleed.c
deleted file mode 100644
index ee0f6ad..0000000
--- a/cleedpy/cleed/original_test_cleed.c
+++ /dev/null
@@ -1,573 +0,0 @@
-/*********************************************************************
-GH/05.10.00
-  file contains functions:
-
-  main
-     Main program for LEED calculations (only for bravaislayer)
-
-Changes:
-
-GH/05.10.00 - early return option (-e).
-
-
-*********************************************************************/
-
-#include <stdio.h>
-#include <string.h>
-
-#include "gh_stddef.h"
-#include "leed.h"
-
-
-
-/*
-#define CONTROL_X
-*/
-#define CONTROL_FLOW
-
-
-#define CONTROL_IO
-#define CONTROL
-#define WARNING
-#define ERROR
-
-#define CTR_NORMAL       998
-#define CTR_EARLY_RETURN 999
-
-/*======================================================================*/
-
-main(int argc, char *argv[])
-
-/*********************************************************************
- Perform a LEED calculation for anisotropic vibrations a general case
-
- INPUT:
-
- DESIGN:
-
-*********************************************************************/
-{
-struct cryst_str *bulk;
-struct cryst_str *over;
-struct phs_str *phs_shifts;
-struct beam_str *beams_all;
-struct beam_str *beams_out;
-struct beam_str *beams_now;
-struct beam_str *beams_set;
-struct var_str *v_par;
-struct eng_str *eng;
-
-mat Tpp,   Tmm,   Rpm,   Rmp;
-mat Tpp_s, Tmm_s, Rpm_s, Rmp_s;
-mat R_bulk, R_tot;
-mat Amp;
-
-int ctr_flag;
-int i_c, i_arg;
-int n_beams_now, n_beams_set;
-int i_set, n_set, offset;
-int i_layer;
-
-real energy;
-real vec[4];
-
-char linebuffer[STRSZ];
-
-char bul_file[STRSZ];                 /* input/output files */
-char par_file[STRSZ];
-char pro_name[STRSZ];
-char res_file[STRSZ];
-
-FILE *pro_stream;
-FILE *res_stream;
-
-  Tpp   =  Tmm   =  Rpm   =  Rmp   = NULL;
-  Tpp_s =  Tmm_s =  Rpm_s =  Rmp_s = NULL;
-  R_bulk = R_tot = NULL;
-  Amp = NULL;
-
-  bulk = over = NULL;
-  phs_shifts  = NULL;
-  beams_all   = NULL;
-  beams_now   = NULL;
-  beams_set   = NULL;
-  beams_out   = NULL;
-  v_par = NULL;
-  eng   = NULL;
-
-
-/*********************************************************************
-  Preset parameters (file names) set by arguments
-*********************************************************************/
-
-  ctr_flag = CTR_NORMAL;
-
-  strncpy(bul_file,"---", STRSZ);
-
-  strncpy(par_file,"---", STRSZ);
-  strncpy(res_file,"leed.res", STRSZ);
-  strncpy(pro_name,"leed.pro", STRSZ);
-
-/*********************************************************************
-  Decode arguments:
-
-    -b <bul_file> - (optional input file) bulk and non-geometrical
-                    parameters.
-    -i <par_file> - (mandatory input file) overlayer parameters of all
-                    parameters (if bul_file does not exist).
-    -o <res_file> - (output file) IV output.
-*********************************************************************/
-
-  for (i_arg = 1; i_arg < argc; i_arg++)
-  {
-    if(*argv[i_arg] != '-')
-    {
-#ifdef ERROR
-      fprintf(STDERR,"*** error (LEED_TEMP):\tsyntax error:\n");
-      fprintf(STDERR,"\tusage: \tleed -i <par_file> -o <res_file>");
-      fprintf(STDERR," [-b <bul_file> -e]\n");
-#endif
-      exit(1);
-    }
-    else
-    {
-
-/* Read parameter input file */
-      if(strncmp(argv[i_arg], "-b", 2) == 0)
-      {
-        i_arg++;
-        strncpy(bul_file, argv[i_arg], STRSZ);
-      } /* -b */
-
-/* Read parameter input file */
-      if(strncmp(argv[i_arg], "-i", 2) == 0)
-      {
-        i_arg++;
-        strncpy(par_file, argv[i_arg], STRSZ);
-      } /* -i */
-
-/* Read and open results file */
-      if(strncmp(argv[i_arg], "-o", 2) == 0)
-      {
-        i_arg++;
-        strncpy(res_file, argv[i_arg], STRSZ);
-        if ((res_stream = fopen(res_file,"w")) == NULL)
-        {
-#ifdef ERROR
-          fprintf(STDERR,
-          "*** error (LEED_TEMP): could not open output file \"%s\"\n",
-          res_file);
-#endif
-          exit(1);
-        }
-      }  /* -o */
-
-/* Read parameter input file */
-      if(strncmp(argv[i_arg], "-e", 2) == 0)
-      {
-        ctr_flag = CTR_EARLY_RETURN;
-      } /* -e */
-
-
-    }  /* else */
-  }  /* for i_arg */
-
-/*********************************************************************
-  Check arguments:
-  - check existence of par_file.
-  - if bul_file is not specified, use par_file instead.
-  - check existence of res_file.
-*********************************************************************/
-
-  if(strncmp(par_file, "---", 3) == 0)
-  {
-#ifdef ERROR
-    fprintf(STDERR,
-    "*** error (LEED_TEMP): no parameter input file (option -i) specified\n");
-#endif
-    exit(1);
-  }
-
-  if(strncmp(bul_file, "---", 3) == 0)
-  {
-    strncpy(bul_file, par_file, STRSZ);
-  }
-
-  if(strncmp(res_file, "leed.res", 8) == 0)
-  {
-#ifdef WARNING
-    fprintf(STDWAR,
-            "* warning (LEED_TEMP): no output file (option -o) specified\n");
-    fprintf(STDWAR,"\toutput will be written to file \"%s\"\n", res_file);
-#endif
-
-    if ((res_stream = fopen(res_file,"w")) == NULL)
-    {
-#ifdef ERROR
-      fprintf(STDERR,
-      "*** error (LEED_TEMP): could not open output file \"%s\"\n",
-      res_file);
-#endif
-     exit(1);
-    }
-  }
-
-/*********************************************************************
-  Read input parameters
-*********************************************************************/
-
-  inp_rdbul_nd(&bulk, &phs_shifts, bul_file);
-  inp_rdpar(&v_par, &eng, bulk, bul_file);
-  inp_rdovl_nd(&over, &phs_shifts, bulk, par_file);
-  n_set = bm_gen(&beams_all, bulk, v_par, eng->fin);
-
-  inp_showbop(bulk, over, phs_shifts);
-
-  if( ctr_flag == CTR_EARLY_RETURN )
-  {
-    fprintf(STDCTR, "(LEED_TEMP): EARLY RETURN \n");
-    exit(0);
-  }
-
-  out_head (bulk, res_stream);
-  out_bmlist(&beams_out, beams_all, eng, res_stream);
-
-/*********************************************************************
- Prepare some often used parameters.
-*********************************************************************/
-
-  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): n_set = %d\n", n_set);
-#endif
-
-/*********************************************************************
- Energy Loop
-*********************************************************************/
-
-  for( energy = eng->ini; energy < eng->fin + E_TOLERANCE; energy += eng->stp)
-  {
-    pc_update(v_par, phs_shifts, energy);
-    n_beams_now = bm_select(&beams_now, beams_all, v_par, bulk->dmin);
-
-#ifdef CONTROL
-    fprintf(STDCTR, "(LEED_TEMP):\n\t => E = %.1f eV (%d beams used) <=\n\n",
-              v_par->eng_v*HART, n_beams_now);
-#endif
-
-
-/*********************************************************************
-  BULK:
-  Loop over beam sets
-
-  Create matrix R_bulk that will eventually contain the bulk
-  reflection matrix
-*********************************************************************/
-    R_bulk = matalloc(R_bulk, n_beams_now, n_beams_now, NUM_COMPLEX);
-    for(offset = 1, i_set = 0; i_set < n_set; i_set ++)
-    {
-      n_beams_set = bm_set(&beams_set, beams_now, i_set);
-
-/*********************************************************************
-    Loop over periodic bulk layers
-*********************************************************************/
-
-    /**********************************************************
-     Compute scattering matrices for bottom-most bulk layer:
-     - single Bravais layer or composite layer
-    **********************************************************/
-
-#ifdef CONTROL_FLOW
-      fprintf(STDCTR, "(LEED_TEMP periodic): bulk layer %d/%d, set %d/%d\n",
-                        0, bulk->nlayers - 1, i_set, n_set - 1);
-#endif
-
-      if( (bulk->layers + 0)->natoms == 1)
-      {
-        ms_bravl_nd( &Tpp, &Tmm, &Rpm, &Rmp,
-                     v_par, (bulk->layers + 0), beams_set);
-      }
-      else
-      {
-        ms_compl_nd( &Tpp, &Tmm, &Rpm, &Rmp,
-                     v_par, (bulk->layers + 0), beams_set);
-      }
-
-#ifdef CONTROL_X
-      fprintf(STDCTR, "(LEED_TEMP): after ms_bravl_nd: Tpp:");
-      matshow(Tpp);
-      fprintf(STDCTR, "(LEED_TEMP): after ms_bravl_nd: Tmm:");
-      matshow(Tmm);
-      fprintf(STDCTR, "(LEED_TEMP): after ms_bravl_nd: Rpm:");
-      matshow(Rpm);
-      fprintf(STDCTR, "(LEED_TEMP): after ms_bravl_nd: Rmp:");
-      matshow(Rmp);
-#endif
-
-    /**********************************************************
-      Loop over the other bulk layers
-    **********************************************************/
-
-      for(i_layer = 1;
-          ( (bulk->layers+i_layer)->periodic == 1) &&
-          (i_layer < bulk->nlayers);
-          i_layer ++)
-      {
-#ifdef CONTROL_FLOW
-        fprintf(STDCTR, "(LEED_TEMP periodic): bulk layer %d/%d, set %d/%d\n",
-                        i_layer, bulk->nlayers - 1, i_set, n_set - 1);
-#endif
-
-    /**************************************************************
-      Compute scattering matrices R/T_s for a single bulk layer
-       - single Bravais layer or composite layer
-    ***************************************************************/
-
-        if( (bulk->layers + i_layer)->natoms == 1)
-        {
-          ms_bravl_nd ( &Tpp_s, &Tmm_s, &Rpm_s, &Rmp_s,
-                        v_par, (bulk->layers + i_layer), beams_set);
-        }
-        else
-        {
-          ms_compl_nd( &Tpp_s, &Tmm_s, &Rpm_s, &Rmp_s,
-                       v_par, (bulk->layers + i_layer), beams_set);
-        }
-
-    /***************************************************************************
-       Add the single layer matrices to the rest by layer doubling
-       - inter layer vector is the vector between layers
-         (i_layer - 1) and (i_layer):
-         (bulk->layers + i_layer)->vec_from_last
-    ****************************************************************************/
-#ifdef CONTROL_FLOW
-        fprintf(STDCTR,
-                "(LEED_TEMP): before ld_2lay vec_from...(%.2f %.2f %.2f)\n",
-                       (bulk->layers + i_layer)->vec_from_last[1] * BOHR,
-                       (bulk->layers + i_layer)->vec_from_last[2] * BOHR,
-                       (bulk->layers + i_layer)->vec_from_last[3] * BOHR);
-#endif
-
-        ld_2lay( &Tpp,  &Tmm,  &Rpm,  &Rmp,
-                 Tpp,   Tmm,   Rpm,   Rmp,
-                 Tpp_s, Tmm_s, Rpm_s, Rmp_s,
-                 beams_set, (bulk->layers + i_layer)->vec_from_last);
-
-      } /* for i_layer (bulk) */
-
-   /*********************************************************************
-      Layer doubling for all periodic bulk layers until convergence is
-      reached:
-       - inter layer vector is (bulk->layers + 0)->vec_from_last
-   **********************************************************************/
-#ifdef CONTROL_FLOW
-      fprintf(STDCTR, "(LEED_TEMP): before ld_2n vec_from...(%.2f %.2f %.2f)\n",
-                      (bulk->layers + 0)->vec_from_last[1] * BOHR,
-                      (bulk->layers + 0)->vec_from_last[2] * BOHR,
-                      (bulk->layers + 0)->vec_from_last[3] * BOHR);
-#endif
-
-      Rpm = ld_2n( Rpm, Tpp, Tmm, Rpm, Rmp,
-                   beams_set, (bulk->layers + 0)->vec_from_last);
-
-   /*******************************************************************
-     Compute scattering matrices for top-most bulk layer if it is
-     not periodic.
-      - single Bravais layer or composite layer
-   **********************************************************************/
-
-      if( i_layer == bulk->nlayers - 1 )
-      {
-#ifdef CONTROL_FLOW
-        fprintf(STDCTR,
-                "(LEED_TEMP not periodic): bulk layer %d/%d, set %d/%d\n",
-                i_layer, bulk->nlayers - 1, i_set, n_set - 1);
-#endif
-
-        if( (bulk->layers + i_layer)->natoms == 1)
-        {
-          ms_bravl_nd( &Tpp_s, &Tmm_s, &Rpm_s, &Rmp_s,
-                    v_par, (bulk->layers + i_layer), beams_set);
-        }
-        else
-        {
-          ms_compl_nd( &Tpp_s, &Tmm_s, &Rpm_s, &Rmp_s,
-                       v_par, (bulk->layers + i_layer), beams_set);
-        }
-
-    /**************************************************************************
-       Add the single layer matrices of the top-most layer to the rest
-       by layer doubling:
-       - inter layer vector is the vector between layers
-         (i_layer - 1) and (i_layer):
-         (bulk->layers + i_layer)->vec_from_last
-    ***************************************************************************/
-
-        Rpm = ld_2lay_rpm(Rpm, Rpm, Tpp_s, Tmm_s, Rpm_s, Rmp_s,
-                          beams_set, (bulk->layers + i_layer)->vec_from_last);
-
-      }  /* if( i_layer == bulk->nlayers - 1 ) */
-
-   /*******************************************************
-     Insert reflection matrix for this beam set into R_bulk.
-   ********************************************************/
-
-      R_bulk = matins(R_bulk, Rpm, offset, offset);
-      offset += n_beams_set;
-
-   /*************************
-     Write cpu time to output
-   **************************/
-
-      sprintf(linebuffer,"(LEED_TEMP): bulk layers set %d, E = %.1f",
-              i_set, energy*HART);
-      cpu_time(STDCPU,linebuffer);
-    }  /* for i_set */
-
-/*********************************************************************
-  OVERLAYER
-  Loop over all overlayer layers
-*********************************************************************/
-
-    for(i_layer = 0; i_layer < over->nlayers; i_layer ++)
-    {
-#ifdef CONTROL_FLOW
-      fprintf(STDCTR, "(LEED_TEMP): overlayer %d/%d\n", i_layer, over->nlayers - 1);
-#endif
-   /***********************************************************
-     Calculate scattering matrices for a single overlayer layer
-      - only single Bravais layer
-   ************************************************************/
-
-      if( (over->layers + i_layer)->natoms == 1)
-      {
-        ms_bravl_nd( &Tpp_s, &Tmm_s, &Rpm_s, &Rmp_s,
-                     v_par, (over->layers + i_layer), beams_now);
-      }
-      else
-      {
-        ms_compl_nd( &Tpp_s, &Tmm_s, &Rpm_s, &Rmp_s,
-                     v_par, (over->layers + i_layer), beams_now);
-      }
-
-#ifdef CONTROL_X
- fprintf(STDCTR, "\n(LEED_TEMP):overlayer %d  ...\n",i_layer);
- fprintf(STDCTR, "\n(LEED_TEMP): Tpp:\n");
- matshowabs(Tpp_s);
- fprintf(STDCTR, "\n(LEED_TEMP): Tmm:\n");
- matshowabs(Tmm_s);
- fprintf(STDCTR, "\n(LEED_TEMP): Rpm:\n");
- matshowabs(Rpm_s);
- fprintf(STDCTR, "\n(LEED_TEMP): Rmp:\n");
- matshowabs(Rmp_s);
-#endif
-
-  /****************************************************************
-     Add the single layer matrices to the rest by layer doubling:
-     - if the current layer is the bottom-most (i_layer == 0),
-       the inter layer vector is calculated from the vectors between
-       top-most bulk layer and origin
-       ( (bulk->layers + nlayers)->vec_to_next )
-       and origin and bottom-most overlayer
-       (over->layers + 0)->vec_from_last.
-
-     - inter layer vector is the vector between layers
-       (i_layer - 1) and (i_layer): (over->layers + i_layer)->vec_from_last
-  **********************************************************************/
-
-      if (i_layer == 0)
-      {
-        for(i_c = 1; i_c <= 3; i_c ++)
-        {
-          vec[i_c] = (bulk->layers + bulk->nlayers - 1)->vec_to_next[i_c]
-                     + (over->layers + 0)->vec_from_last[i_c];
-        }
-
-#ifdef CONTROL_FLOW
-        fprintf(STDCTR,
-                "(LEED):over0 before ld_2lay_rpm vec..(%.2f %.2f %.2f)\n",
-                vec[1] * BOHR,vec[2] * BOHR, vec[3] * BOHR);
-#endif
-
-        R_tot = ld_2lay_rpm(R_tot, R_bulk, Tpp_s, Tmm_s, Rpm_s, Rmp_s,
-                            beams_now, vec);
-      }
-      else
-      {
-#ifdef CONTROL_FLOW
-        fprintf(STDCTR,
-                "(LEED):over%d  before ld_2lay_rpm vec..(%.2f %.2f %.2f)\n",
-                i_layer,(over->layers + i_layer)->vec_from_last[1] * BOHR,
-                        (over->layers + i_layer)->vec_from_last[2] * BOHR,
-                        (over->layers + i_layer)->vec_from_last[3] * BOHR);
-#endif
-
-        R_tot = ld_2lay_rpm(R_tot, R_tot, Tpp_s, Tmm_s, Rpm_s, Rmp_s,
-                          beams_now, (over->layers + i_layer)->vec_from_last);
-      }
-
-   /**************************
-     Write cpu time to output
-   **************************/
-
-      sprintf(linebuffer,"(LEED): overlayer %d, E = %.1f",
-              i_layer, energy * HART);
-      cpu_time(STDCPU,linebuffer);
-
-    }  /* for i_layer (overlayer) */
-
-/*********************************************
-   Add propagation towards the potential step.
-**********************************************/
-
-    vec[1] = vec[2] = 0.;
-    vec[3] = 1.25 / BOHR;
-
-/********************************************
-    No scattering at pot. step
-********************************************/
-
-    Amp = ld_potstep0(Amp, R_tot, beams_now, v_par->eng_v, vec);
-    out_int(Amp, beams_now, beams_out, v_par, res_stream);
-
-/********************************************
-    Write cpu time to output
-********************************************/
-
-    sprintf(linebuffer,"  %.1f   %d  ",energy * HART,n_beams_now);
-    cpu_time(STDWAR,linebuffer);
-
-
-  } /* end of energy loop */
-
-
-#ifdef CONTROL_IO
-  fprintf(STDCTR, "(LEED): end of energy loop: close files\n");
-#endif
-
-  fclose(res_stream);
-
-#ifdef CONTROL
-  fprintf(STDCTR, "\n\n(LEED):\tCORRECT TERMINATION");
-#endif
-
-/********************************************
-    Write cpu time to output
-********************************************/
-
-  cpu_time(STDCPU,"");
-
-/********************************************
-    set exit status explicitly
-********************************************/
-
-  printf("The original code is working fine.\n");
-
-  exit(0);
-
-} /* end of main */
diff --git a/cleedpy/cleed/src/leed.c b/cleedpy/cleed/src/leed.c
index ca95350..961ca95 100644
--- a/cleedpy/cleed/src/leed.c
+++ b/cleedpy/cleed/src/leed.c
@@ -1,5 +1,15 @@
 #include "leed.h"
 
+typedef struct {
+    int n_beams;
+    real * beam_index1;
+    real * beam_index2;
+    int * beam_set;
+    int n_energies;
+    real * energies;
+    real * iv_curves;
+} CleedResult;
+
 void print_phase_shift(struct phs_str phs_shift)
 {
     int i;
@@ -23,11 +33,9 @@ void print_phase_shift(struct phs_str phs_shift)
 
 }
 
-real ** leed(
-    char * par_file,
-    char * bul_file,
-    char *res_file
-    ){
+
+CleedResult leed(char * par_file, char * bul_file)
+{
     struct cryst_str *bulk=NULL;
     struct cryst_str *over=NULL;
     struct phs_str *phs_shifts=NULL;
@@ -41,17 +49,14 @@ real ** leed(
 
     int n_beams_now, n_beams_set;
 
+    CleedResult results;
+
     int i_c, i_set, offset, i;
     int i_layer;
     int energy_index;
     int n_set;
-    int energy_list_size;
-    real *energy_list;
     real energy;
     real vec[4];
-
-    real **iv_curves=NULL;
-
     mat R_bulk=NULL, R_tot=NULL;
     mat Amp=NULL;
 
@@ -60,8 +65,6 @@ real ** leed(
 
     struct eng_str *eng=NULL;
 
-    FILE *res_stream;
-
     // Read input parameters
     inp_rdbul_nd(&bulk, &phs_shifts, bul_file);
     inp_rdpar(&v_par, &eng, bulk, bul_file);
@@ -71,19 +74,18 @@ real ** leed(
         print_phase_shift(phs_shifts[i]);
 
 
-    // Construct energy list
-    energy_list_size = (eng->fin - eng->ini)/eng->stp + 1;
-    energy_list = (real *) malloc(energy_list_size * sizeof(real));
-    for (energy_index=0; energy_index<energy_list_size; energy_index++)
-    {
-        energy_list[energy_index] = eng->ini + energy_index * eng->stp;
-    }
 
+    // Construct energy list
+    results.n_energies = (eng->fin - eng->ini)/eng->stp + 1;
+    results.energies = (real *) malloc(results.n_energies * sizeof(real));
 
-    res_stream = fopen(res_file,"w");
+    for (energy_index=0; energy_index < results.n_energies; energy_index++)
+        results.energies[energy_index] = eng->ini + energy_index * eng->stp;
 
-    out_head(bulk, res_stream);
+    eng->fin = results.energies[results.n_energies - 1];
 
+    // Printing stuff
+    inp_showbop(bulk, over, phs_shifts);
 
     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.
@@ -93,13 +95,13 @@ real ** leed(
 
     /* Generate beams out */
     n_set = bm_gen(&beams_all, bulk, v_par, eng->fin);
-    out_bmlist(&beams_out, beams_all, eng, res_stream);
-
+    results.n_beams = out_bmlist(&beams_out, beams_all, eng, &results.beam_index1, &results.beam_index2, &results.beam_set);
+    results.iv_curves = (real *) calloc(results.n_energies * results.n_beams, sizeof(real));
 
     /* Main Energy Loop */
 
-    for(energy_index=0; energy_index<energy_list_size; energy_index++){
-        pc_update(v_par, phs_shifts, energy_list[energy_index]);
+    for(energy_index=0; energy_index < results.n_energies; energy_index++){
+        pc_update(v_par, phs_shifts, results.energies[energy_index]);
         n_beams_now = bm_select(&beams_now, beams_all, v_par, bulk->dmin);
 
         /*********************************************************************
@@ -280,23 +282,9 @@ real ** leed(
         ********************************************/
 
         Amp = ld_potstep0(Amp, R_tot, beams_now, v_par->eng_v, vec);
-        out_int(Amp, beams_now, beams_out, v_par, res_stream);
+        out_int(Amp, beams_now, beams_out, v_par, &results.iv_curves[energy_index * results.n_beams]);
 
     }  /* end of energy loop */
 
-    fclose(res_stream);
-
-    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;
+    return results;
 }
diff --git a/cleedpy/cleed/src/loutbmlist.c b/cleedpy/cleed/src/loutbmlist.c
index 2030aba..53db5c2 100644
--- a/cleedpy/cleed/src/loutbmlist.c
+++ b/cleedpy/cleed/src/loutbmlist.c
@@ -26,20 +26,14 @@
 
 #include "leed.h"
 
-/*
-#define WARNING
-*/
-
-#define CONTROL
-#define ERROR
-
-#define EXIT_ON_ERROR
-
 
 int out_bmlist(struct beam_str ** p_beams_out,
                struct beam_str *beams_all,
                struct eng_str *eng,
-               FILE * outfile)
+               real ** index1,
+               real ** index2,
+               int ** beam_set
+               )
 
 /************************************************************************
 
@@ -103,121 +97,92 @@ int out_bmlist(struct beam_str ** p_beams_out,
 *************************************************************************/
 {
 
-int n_beams, n_eng;
-int i_bm_all, i_bm_out;
+    int n_beams, n_eng;
+    int i_bm_all, i_bm_out;
 
-real k_max;
-real faux;
+    real k_max;
+    real faux;
 
-struct beam_str *beams_out;
+    struct beam_str *beams_out;
 
-/************************************************************************
-  Find number of beams in list beams_all (n_beams) and
-  allocate *p_beams_out = beams_out of the same size.
-************************************************************************/
+    /************************************************************************
+     Find number of beams in list beams_all (n_beams) and
+    allocate *p_beams_out = beams_out of the same size.
+    ************************************************************************/
 
- for(n_beams = 0; (beams_all + n_beams)->k_par != F_END_OF_LIST; n_beams ++)
- { ; }
+    for(n_beams = 0; (beams_all + n_beams)->k_par != F_END_OF_LIST; n_beams ++);
 
-#ifdef CONTROL
- fprintf(STDCTR,"(out_bmlist): \tTotal No. of beams %.f eV = %d\n",
-                 eng->fin*HART, n_beams);
-#endif
 
- if (*p_beams_out == NULL)
-   *p_beams_out = beams_out = (struct beam_str *)
-                 calloc(n_beams + 1, sizeof(struct beam_str));
- else
-   *p_beams_out = beams_out = (struct beam_str *)
-                 realloc(*p_beams_out, (n_beams+1) * sizeof(struct beam_str));
-
- if(beams_out == NULL)
- {
-#ifdef ERROR
-   fprintf(STDERR," *** error (out_bmlist): allocation error.\n");
-#endif
-#ifdef EXIT_ON_ERROR
-   exit(1);
-#else
-   return(-1)
-#endif
- }
+    if (*p_beams_out == NULL)
+        *p_beams_out = beams_out = (struct beam_str *) calloc(n_beams + 1, sizeof(struct beam_str));
+    else
+        *p_beams_out = beams_out = (struct beam_str *) realloc(*p_beams_out, (n_beams+1) * sizeof(struct beam_str));
 
-/************************************************************************
-  Write appropriate beams to beams_out
-
-  k_max is the square of the maximum pareallel vector component of
-        non-evanescent wave vectors at eng->fin.
-  n_beams is set to number of output beams afterwards.
-************************************************************************/
-
- k_max =  2. * eng->fin;
-
- for(i_bm_all = 0, i_bm_out = 0; i_bm_all < n_beams; i_bm_all ++)
- {
-   if( (beams_all + i_bm_all)->k_par <= k_max )
-   {
-     memcpy( beams_out + i_bm_out,
-             beams_all + i_bm_all,
-             sizeof(struct beam_str) );
-     i_bm_out ++;
-   }
- }
-
-/* terminate list beams_out */
- (beams_out + i_bm_out)->k_par = F_END_OF_LIST;
- n_beams = i_bm_out;
-
-#ifdef CONTROL
- fprintf(STDCTR,"(out_bmlist): \tNo. of non_evanescent beams = %d\n",
-                  n_beams);
-#endif
+    if(beams_out == NULL)
+    {
+        fprintf(STDERR," *** error (out_bmlist): allocation error.\n");
+        exit(1);
+    }
 
+    /************************************************************************
+     Write appropriate beams to beams_out
 
-/************************************************************************
-  Write energies, number of beams, and beams to output
+    k_max is the square of the maximum pareallel vector component of
+            non-evanescent wave vectors at eng->fin.
+    n_beams is set to number of output beams afterwards.
+    ************************************************************************/
 
-  "#en" energies:
-        number of energies, initial energy, final energy, energy step.
-  "#bn" number of beams
-  "#bi" beam indices:
-        number (starting from zero), 1st index, 2nd index, beam set.
+    k_max =  2. * eng->fin;
+
+    for(i_bm_all = 0, i_bm_out = 0; i_bm_all < n_beams; i_bm_all ++)
+    {
+        if( (beams_all + i_bm_all)->k_par <= k_max )
+        {
+            memcpy(beams_out + i_bm_out, beams_all + i_bm_all, sizeof(struct beam_str) );
+            i_bm_out ++;
+        }
+    }
+
+    /* terminate list beams_out */
+    (beams_out + i_bm_out)->k_par = F_END_OF_LIST;
+    n_beams = i_bm_out;
+
+
+    /************************************************************************
+     Write energies, number of beams, and beams to output
+
+    "#en" energies:
+            number of energies, initial energy, final energy, energy step.
+    "#bn" number of beams
+    "#bi" beam indices:
+            number (starting from zero), 1st index, 2nd index, beam set.
+
+    - Count energy points (Running through the same loop0 as the main
+        program does avoids inconsistencies due to rounding errors)
+    -
+    ************************************************************************/
+
+    /* energies */
+    for(faux = eng->ini, n_eng = 0; faux <= eng->fin; faux += eng->stp, n_eng++);
+
+
+    // Allocate memory for index1 and index2
+    (* index1) = (real *) calloc(n_beams, sizeof(real));
+    (* index2) = (real *) calloc(n_beams, sizeof(real));
+    (* beam_set) = (int *) calloc(n_beams, sizeof(int));
+
+
+    /* beam indices */
+    for(i_bm_out = 0; i_bm_out < n_beams; i_bm_out ++)
+    {
+        (* index1)[i_bm_out] = (beams_out+i_bm_out)->ind_1;
+        (* index2)[i_bm_out] = (beams_out+i_bm_out)->ind_2;
+        (* beam_set)[i_bm_out] = (beams_out+i_bm_out)->set;
+    }
+
+
+    /* write beams_out back to pointer */
+    *p_beams_out = beams_out;
 
-  - Count energy points (Running through the same loop0 as the main
-    program does avoids inconsistencies due to rounding errors)
-  -
-************************************************************************/
-
-/* energies */
- for(faux = eng->ini, n_eng = 0;
-     faux <= eng->fin;
-     faux += eng->stp, n_eng++)
- { ; }
-
- fprintf(outfile, "#en %d %f %f %f\n",
-                   n_eng,
-                   eng->ini*HART,
-                   eng->fin*HART,
-                   eng->stp*HART);
-
-/* number of beams */
- fprintf(outfile, "#bn %d\n", n_beams);
-
-/* beam indices */
- for(i_bm_out = 0; i_bm_out < n_beams; i_bm_out ++)
- {
-   fprintf(outfile, "#bi %d %f %f %d\n",
-                     i_bm_out,
-                     (beams_out+i_bm_out)->ind_1,
-                     (beams_out+i_bm_out)->ind_2,
-                     (beams_out+i_bm_out)->set);
- }
-/* flush output file */
- fflush(outfile);
-
-/* write beams_out back to pointer */
- *p_beams_out = beams_out;
-
- return(n_beams);
-}  /* end of function out_bmlist */
-/************************************************************************/
+    return n_beams;
+}
diff --git a/cleedpy/cleed/src/louthead.c b/cleedpy/cleed/src/louthead.c
deleted file mode 100644
index 27f70dd..0000000
--- a/cleedpy/cleed/src/louthead.c
+++ /dev/null
@@ -1,85 +0,0 @@
-/*********************************************************************
-GH/11.08.95
-  file contains function:
-
-  int out_head( FILE * outfile)
-
- Write header information to output file.
-
- Changes:
-
- GH/11.08.95 - Creation
-
-*********************************************************************/
-
-#include <stdio.h>
-#include <time.h>
-
-#include "leed.h"
-
-/*
-#define WARNING
-*/
-
-#define CONTROL
-#define ERROR
-
-#define EXIT_ON_ERROR
-
-#ifndef LEED_VERSION               /* should be defined in leed_def.h */
-#define LEED_VERSION "0.0 (test version GH/11.08.95)"
-#endif
-
-int out_head(struct cryst_str *cryst_par,
-             FILE * outfile)
-
-/************************************************************************
-
- write header information to output file.
-
- INPUT:
-
-  struct cryst_str *bulk_par - (input)
-  FILE * outfile - (input) pointer to the output file were the output
-            is written to.
-
- DESIGN:
-
-  "#vn" No of version.
-  "#ts" start time and date.
-
- RETURN VALUES:
-
-   1   if o.k.
-  -1   if failed
-
-*************************************************************************/
-{
-
-struct tm *l_time;
-time_t t_time;
-
- fprintf(outfile, "# ####################################### #\n");
- fprintf(outfile, "#            output from CLEED            #\n");
- fprintf(outfile, "# ####################################### #\n");
-
-/************************************************************************
-  Write version number and start time to output file
-************************************************************************/
-
- fprintf(outfile, "#vn %s\n", LEED_VERSION);
-
- t_time = time(NULL);
- l_time = localtime(&t_time);
-
- fprintf(outfile, "#ts %s", asctime(l_time) );
- fprintf(outfile, "#\n");
-
-#ifdef CONTROL
- fprintf(STDCTR,"(out_head): Start date: %s", asctime(l_time) );
-#endif
-
-
- return(1);
-}  /* end of function out_head */
-/************************************************************************/
diff --git a/cleedpy/cleed/src/loutint.c b/cleedpy/cleed/src/loutint.c
index 66be862..0995280 100644
--- a/cleedpy/cleed/src/loutint.c
+++ b/cleedpy/cleed/src/loutint.c
@@ -22,132 +22,54 @@
 
 #include "leed.h"
 
-/*
-#define CONTROL_ALL
-#define WARNING
-*/
-#define CONTROL
-#define ERROR
-
-#define EXIT_ON_ERROR
-
-#ifndef INT_TOLERANCE            /* should be defined in "leed_def.h" */
-#define INT_TOLERANCE 1.e-10               /* min intensity != 0. */
-#endif
-
-int out_int(mat Amp, struct beam_str *beams_now, struct beam_str *beams_all,
-            struct var_str *par, FILE * outfile)
 
+int out_int(mat Amp, struct beam_str *beams_now, struct beam_str *beams_all, struct var_str *par, real * beam_intensities)
 /************************************************************************
+    output of beam intensities
 
- output of beam intensities
-
- INPUT:
-
-  mat Amp - (input) vector containing the beam amplitudes of all beams
-           included at the current energy.
-  struct beam_str *beams_now  -  all beams included at the current energy.
-  struct beam_str *beams_all  -  all beams included at the highest energy.
-  FILE * outfile - pointer to the output file were the intensities are
-           written to.
+    INPUT:
 
- RETURN VALUES:
+    mat Amp - (input) vector containing the beam amplitudes of all beams
+            included at the current energy.
+    struct beam_str *beams_now  -  all beams included at the current energy.
+    struct beam_str *beams_all  -  all beams included at the highest energy.
+    FILE * outfile - pointer to the output file were the intensities are
+            written to.
 
-  number of intensities written to outfile.
-  -1 if failed.
+    RETURN VALUES:
 
- DESIGN:
+    number of intensities written to outfile.
+    -1 if failed.
 
+    DESIGN:
 *************************************************************************/
 {
 
-int i_beams_now, i_beams_all, i_out;
-
-mat Int;
-real k_r;
-
- Int = NULL;
-/*********************************************************
-   Calculate intensitied as the square of the moduli of the
-   amplitudes.
-*********************************************************/
-
- Int = matsqmod(Int, Amp);
-
- k_r = R_sqrt(2*par->eng_v);
-
-/*********************************************************
-   Print intensities for non-evanescent beams
-*********************************************************/
-
-#ifdef CONTROL
- fprintf(STDCTR,"(out_int):\t     beam\t  intensity\n\t\t\t== %.2f eV ==\n",
-                par->eng_v*HART);
- for (i_beams_now = 0, i_out = 0; i_beams_now < Int->rows; i_beams_now ++)
- {
-   if((beams_now + i_beams_now)->k_par <= k_r)
-   {
-     if(Int->rel[i_beams_now + 1] > INT_TOLERANCE)
-       fprintf(STDCTR,"\t\t(%5.2f,%5.2f):\t  %.4e\n",
-        (beams_now + i_beams_now)->ind_1, (beams_now + i_beams_now)->ind_2,
-         Int->rel[i_beams_now + 1]);
-     else
-       fprintf(STDCTR,"\t\t(%5.2f,%5.2f):\t   < %.0e\n",
-        (beams_now + i_beams_now)->ind_1, (beams_now + i_beams_now)->ind_2,
-         INT_TOLERANCE);
-     i_out ++;
-   }
-/*  print also evnescent beams
-   else
-     fprintf(STDCTR,"\t\t(%5.2f,%5.2f):\t (evanescent)\n",
-             (beams_now + i_beams_now)->ind_1, (beams_now + i_beams_now)->ind_2,
-             Int->rel[i_beams_now + 1]);
-*/
- }
- fprintf(STDCTR,"\t\t(%d/%d)\n",i_beams_now, i_out);
-#endif
-
-#ifdef CONTROL_ALL
- fprintf(STDCTR,"\n(out_int): all beams: \n");
-
- for (i_beams_all = 0; (beams_all + i_beams_all)->k_par != F_END_OF_LIST;
-      i_beams_all ++)
- {
-   fprintf(STDCTR,"\t\t(%5.2f,%5.2f)\n",
-   (beams_all + i_beams_all)->ind_1, (beams_all + i_beams_all)->ind_2);
- }
- fprintf(STDCTR,"\t\t(%d)\n",i_beams_all);
-#endif
-
-
- fprintf(outfile,"%.2f ", par->eng_v*HART);
-
- for(i_beams_all = 0; (beams_all + i_beams_all)->k_par != F_END_OF_LIST;
-     i_beams_all ++)
- {
-   for(i_beams_now = 0; i_beams_now < Int->rows; i_beams_now ++)
-   {
-     if( ((beams_all+i_beams_all)->ind_1 == (beams_now+i_beams_now)->ind_1) &&
-         ((beams_all+i_beams_all)->ind_2 == (beams_now+i_beams_now)->ind_2)  )
-     {
-       if((beams_now + i_beams_now)->k_par <= k_r)
-       {
-         if(Int->rel[i_beams_now + 1] > INT_TOLERANCE)
-           fprintf(outfile,"%.6e ", Int->rel[i_beams_now + 1]);
-         else
-           fprintf(outfile,"%.6e ",0.);
-       }
-       else
-         fprintf(outfile,"%.6e ",0.);
-       break;
-     }  /* if index = index */
-   }  /* for i_beams_now */
-
-   if( i_beams_now == Int->rows) fprintf(outfile,"%.6e ",0.);
- } /* for i_beams_all */
-
- fprintf(outfile,"\n");
-
- fflush(outfile);
- return(i_beams_now);
-}  /* end of function out_int */
+    int i_beams_now, i_beams_all, i_out;
+    mat Int=NULL;
+    real k_r;
+
+    // Calculate intensitied as the square of the moduli of the amplitudes.
+    Int = matsqmod(Int, Amp);
+    k_r = R_sqrt(2*par->eng_v);
+
+    for(i_beams_all = 0; (beams_all + i_beams_all)->k_par != F_END_OF_LIST; i_beams_all ++)
+    {
+        for(i_beams_now = 0; i_beams_now < Int->rows; i_beams_now ++)
+        {
+            if( ((beams_all+i_beams_all)->ind_1 == (beams_now+i_beams_now)->ind_1) && ((beams_all+i_beams_all)->ind_2 == (beams_now+i_beams_now)->ind_2)  )
+            {
+                if (((beams_now+i_beams_now)->k_par <= k_r) && (Int->rel[i_beams_now + 1] > INT_TOLERANCE))
+                    beam_intensities[i_beams_all] = Int->rel[i_beams_now + 1];
+                else
+                    beam_intensities[i_beams_all] = 0.;
+                break;
+            }
+        }  /* for i_beams_now */
+
+        if( i_beams_now == Int->rows)
+            beam_intensities[i_beams_all] = 0.;
+    } /* for i_beams_all */
+
+    return(i_beams_now);
+}
diff --git a/cleedpy/cleed/test_cleed.c b/cleedpy/cleed/test_cleed.c
deleted file mode 100644
index 575078a..0000000
--- a/cleedpy/cleed/test_cleed.c
+++ /dev/null
@@ -1,17 +0,0 @@
-#include <stdio.h>
-#include <string.h>
-
-#include "gh_stddef.h"
-#include "leed.h"
-
-main(){
-  // File names
-  char *bul_file="Ni111_Cu.inp"; // bulk and non-geometrical parameters.
-  char *par_file="Ni111_Cu.inp"; // overlayer parameters of all parameters (if bul_file does not exist)
-  char *res_file="leed.res"; // output file containing the LEED intensities
-
-  // Calculate LEED intensities
-  leed(bul_file, par_file, res_file);
-
-  exit(0);
-}
diff --git a/cleedpy/cli/leed.py b/cleedpy/cli/leed.py
index 83e8e45..8f99a0e 100644
--- a/cleedpy/cli/leed.py
+++ b/cleedpy/cli/leed.py
@@ -1,11 +1,13 @@
+import datetime
 from pathlib import Path
 
 import typer
 
 from ..config import load_parameters
 from ..interface.cleed import call_cleed
-from ..physics.constants import BOHR_TO_ANGSTROM
+from ..physics.constants import BOHR_TO_ANGSTROM, HART
 from ..preprocess import extract_bulk_parameters, transform_cells
+from ..version import __version__
 
 
 def leed_bak(inptut_file: Path):
@@ -25,6 +27,43 @@ def leed_bak(inptut_file: Path):
     )
 
 
+CLEED_OUT_HEADER = f"""# ####################################### #
+#            output from CLEED            #
+# ####################################### #
+# vn {__version__}
+# ts {datetime.datetime.now().ctime()}
+#
+"""
+
+
+def print_cleed_results(result, output_file):
+    """Function that prints the results of cleed."""
+    with open(output_file, "w") as f:
+        f.write(CLEED_OUT_HEADER)
+        energy_step = (result.energies[result.n_energies - 1] - result.energies[0]) / (
+            result.n_energies - 1
+        )
+        f.write(
+            f"#en {result.n_energies} {result.energies[0]*HART:.6f} {result.energies[result.n_energies-1]*HART:.6f} {energy_step*HART:.6f}\n"
+        )
+        f.write(f"#bn {result.n_beams}\n")
+        for i in range(result.n_beams):
+            f.write(
+                f"#bi {i} {result.beam_index1[i]:.6f} {result.beam_index2[i]:.6f} {result.beam_set[i]}\n"
+            )
+        for i in range(result.n_energies):
+            f.write(
+                " ".join(
+                    [f"{result.energies[i]*HART:.2f}"]
+                    + [
+                        f"{result.iv_curves[i*result.n_beams+j]:.6e}"
+                        for j in range(result.n_beams)
+                    ]
+                )
+            )
+            f.write("\n")
+
+
 def leed(
     parameters_file: str = typer.Option(  # noqa: B008
         "leed.inp", "--input", "-i", help="Input file with parameters"
@@ -35,8 +74,9 @@ def leed(
     ),  # noqa: B008
 ):
     """Leed CLI."""
-    result = call_cleed(parameters_file, bulk_file, output_file)
-    print(result)
+    result = call_cleed(parameters_file, bulk_file)
+
+    print_cleed_results(result, output_file)
 
 
 def cli():
diff --git a/cleedpy/interface/cleed.py b/cleedpy/interface/cleed.py
index 5c55afd..5336c7a 100644
--- a/cleedpy/interface/cleed.py
+++ b/cleedpy/interface/cleed.py
@@ -107,6 +107,30 @@ class EnergyLoopVariables(Structure):
     ]
 
 
+class CleedResult(Structure):
+    """Parses C structure to python class:
+    struct leed_results {
+        int n_beams;
+        real * beam_index1;
+        real * beam_index2;
+        int * beam_set;
+        int n_energies;
+        real * energies;
+        real * iv_curves;
+    };
+    """
+
+    _fields_ = [
+        ("n_beams", c_int),
+        ("beam_index1", POINTER(c_double)),
+        ("beam_index2", POINTER(c_double)),
+        ("beam_set", POINTER(c_int)),
+        ("n_energies", c_int),
+        ("energies", POINTER(c_double)),
+        ("iv_curves", POINTER(c_double)),
+    ]
+
+
 @dataclass
 class CleedInputs:
     bulk: Crystal
@@ -210,15 +234,13 @@ def get_cleed_lib() -> CDLL:
     return cdll.LoadLibrary(cleed_lib)
 
 
-def call_cleed(parameters_file, bulk_file, output_file):
+def call_cleed(parameters_file, bulk_file):
     lib = get_cleed_lib()
 
-    lib.leed.argtypes = [c_char_p, c_char_p, c_char_p]
-    lib.leed.restype = POINTER(POINTER(c_double))
+    lib.leed.argtypes = [c_char_p, c_char_p]
+    lib.leed.restype = CleedResult
 
-    result = lib.leed(
-        parameters_file.encode(), bulk_file.encode(), output_file.encode()
-    )
+    result = lib.leed(parameters_file.encode(), bulk_file.encode())
 
     return result
 
diff --git a/examples/cleed_example_old/leed.res_new b/examples/cleed_example_old/leed.res_new
index 640419c..0f27747 100644
--- a/examples/cleed_example_old/leed.res_new
+++ b/examples/cleed_example_old/leed.res_new
@@ -4,7 +4,7 @@
 #vn 0.0 (test version GH/11.08.95)
 #ts Fri Jan 12 09:25:30 2024
 #
-#en 107 70.000000 498.000000 4.000000
+#en 108 70.000000 498.000000 4.000000
 #bn 55
 #bi 0 0.000000 0.000000 0
 #bi 1 -1.000000 0.000000 0
diff --git a/examples/cleedpy_example_intermediate/run.sh b/examples/cleedpy_example_intermediate/run.sh
index e4f5f8a..3ca2620 100755
--- a/examples/cleedpy_example_intermediate/run.sh
+++ b/examples/cleedpy_example_intermediate/run.sh
@@ -1 +1,3 @@
+export CLEED_PHASE=./PHASE
+
 cleedpy-leed -i leed.inp -b leed.inp -o leed.res
diff --git a/examples/ni111_2x2O/Ni111_2x2O.inp b/examples/ni111_2x2O/Ni111_2x2O.inp
new file mode 100644
index 0000000..558e59a
--- /dev/null
+++ b/examples/ni111_2x2O/Ni111_2x2O.inp
@@ -0,0 +1,47 @@
+# input file for SEARCH
+# Ni(111) + 2x2 O center on fcc site
+# 25 Mai 2001
+# lattice parameters
+a1:       1.2450  -2.1564   0.0000
+a2:       1.2450   2.1564   0.0000
+a3:       0.0000   0.0000  -6.0990
+#
+m1:  2.  0.
+m2:  0.  2.
+#
+vr:   -8.00
+vi:     4.00
+#
+#surface:
+po: O_CO_Pendry_cs    0.0000  0.0000  5.2000       dr3  0.061  0.061  0.061
+#
+po: Ni_Wakoh_cs       1.2450 -0.7188  4.1000       dr3  0.025  0.025  0.025
+po: Ni_Wakoh_cs      -1.2450 -0.7188  4.1000       dr3  0.025  0.025  0.025
+po: Ni_Wakoh_cs       2.4900  1.4376  4.1000       dr3  0.025  0.025  0.025
+po: Ni_Wakoh_cs       0.0000  1.4376  4.1000       dr3  0.025  0.025  0.025
+#
+po: Ni_Wakoh_cs       1.2450  0.7188  2.0000       dr3  0.025  0.025  0.025
+po: Ni_Wakoh_cs      -1.2450  0.7188  2.0000       dr3  0.025  0.025  0.025
+po: Ni_Wakoh_cs       0.0000 -1.4376  2.0000       dr3  0.025  0.025  0.025
+po: Ni_Wakoh_cs      -2.4900 -1.4376  2.0000       dr3  0.025  0.025  0.025
+#
+# bulk:
+pb: Ni_Wakoh_cs  0.0000    +0.0000   0.0000  dr3 0.025 0.025 0.025
+pb: Ni_Wakoh_cs  1.2450    -0.7188  -2.0330  dr3 0.025 0.025 0.025
+pb: Ni_Wakoh_cs  1.2450    +0.7188  -4.0660  dr3 0.025 0.025 0.025
+# minimum radii:
+rm: Ni_Wakoh_cs  1.10
+rm: O_CO_Pendry_cs   0.40
+# z range
+zr: 1.60  6.00
+# sz: 0 (xyz search), 1 (z only)
+sz: 0
+# sr: rotational axis
+sr: 3  0.000 0.00
+ei: 70.
+ef: 300.1
+es: 4.
+it: 0.
+ip: 0.
+ep: 1.e-2
+lm: 7
diff --git a/examples/ni111_2x2O/PHASE/Ni_Wakoh_cs.phs b/examples/ni111_2x2O/PHASE/Ni_Wakoh_cs.phs
new file mode 100644
index 0000000..183b53a
--- /dev/null
+++ b/examples/ni111_2x2O/PHASE/Ni_Wakoh_cs.phs
@@ -0,0 +1,201 @@
+100 11 neng lmax (Ni, 58.70 amu, potential (WAKOH) from potpsb., zero at MT radius)
+  .2000
+ -.2548 -.0144  .0765  .0002  .0000  .0000  .0000  .0000  .0000  .0000 .0000 .0000
+  .4000
+ -.4081 -.0365 -.3223  .0019  .0001  .0000  .0000  .0000  .0000  .0000 .0000 .0000
+  .6000
+ -.5439 -.0638 -.2437  .0071  .0003  .0000  .0000  .0000  .0000  .0000 .0000 .0000
+  .8000
+ -.6675 -.0965 -.2424  .0175  .0010  .0000  .0000  .0000  .0000  .0000 .0000 .0000
+ 1.0000
+ -.7813 -.1340 -.2438  .0344  .0024  .0001  .0000  .0000  .0000  .0000 .0000 .0000
+ 1.2000
+ -.8867 -.1753 -.2407  .0588  .0048  .0003  .0000  .0000  .0000  .0000 .0000 .0000
+ 1.4000
+ -.9847 -.2193 -.2328  .0914  .0085  .0007  .0000  .0000  .0000  .0000 .0000 .0000
+ 1.6000
+-1.0761 -.2650 -.2211  .1319  .0138  .0013  .0001  .0000  .0000  .0000 .0000 .0000
+ 1.8000
+-1.1616 -.3115 -.2073  .1799  .0208  .0022  .0002  .0000  .0000  .0000 .0000 .0000
+ 2.0000
+-1.2418 -.3581 -.1928  .2342  .0299  .0035  .0003  .0000  .0000  .0000 .0000 .0000
+ 2.2000
+-1.3171 -.4043 -.1789  .2931  .0409  .0053  .0006  .0000  .0000  .0000 .0000 .0000
+ 2.4000
+-1.3880 -.4497 -.1665  .3547  .0540  .0077  .0009  .0001  .0000  .0000 .0000 .0000
+ 2.6000
+-1.4548 -.4941 -.1561  .4170  .0691  .0107  .0014  .0001  .0000  .0000 .0000 .0000
+ 2.8000
+-1.5179 -.5372 -.1478  .4783  .0860  .0144  .0020  .0002  .0000  .0000 .0000 .0000
+ 3.0000
+-1.5776 -.5789 -.1417  .5372  .1047  .0188  .0028  .0003  .0000  .0000 .0000 .0000
+ 3.2000
+-1.6343 -.6193 -.1377  .5928  .1248  .0239  .0038  .0005  .0001  .0000 .0000 .0000
+ 3.4000
+-1.6880 -.6581 -.1355  .6445  .1460  .0299  .0051  .0007  .0001  .0000 .0000 .0000
+ 3.6000
+-1.7392 -.6956 -.1348  .6923  .1682  .0366  .0066  .0010  .0001  .0000 .0000 .0000
+ 3.8000
+-1.7880 -.7316 -.1354  .7363  .1909  .0441  .0085  .0013  .0002  .0000 .0000 .0000
+ 4.0000
+-1.8346 -.7662 -.1370  .7767  .2140  .0523  .0106  .0017  .0002  .0000 .0000 .0000
+ 4.2000
+-1.8791 -.7995 -.1394  .8139  .2371  .0612  .0131  .0023  .0003  .0000 .0000 .0000
+ 4.4000
+-1.9218 -.8316 -.1425  .8484  .2600  .0707  .0159  .0029  .0004  .0001 .0000 .0000
+ 4.6000
+-1.9629 -.8625 -.1460  .8804  .2824  .0807  .0191  .0037  .0006  .0001 .0000 .0000
+ 4.8000
+-2.0024 -.8923 -.1499  .9103  .3043  .0913  .0227  .0046  .0008  .0001 .0000 .0000
+ 5.0000
+-2.0404 -.9211 -.1539  .9384  .3256  .1022  .0266  .0056  .0010  .0001 .0000 .0000
+ 5.2000
+-2.0772 -.9489 -.1582  .9648  .3461  .1135  .0308  .0068  .0012  .0002 .0000 .0000
+ 5.4000
+-2.1127 -.9758 -.1626  .9900  .3658  .1250  .0355  .0082  .0015  .0002 .0000 .0000
+ 5.6000
+-2.1471-1.0019 -.1670 1.0138  .3846  .1366  .0404  .0098  .0019  .0003 .0000 .0000
+ 5.8000
+-2.1804-1.0272 -.1715 1.0366  .4027  .1484  .0456  .0115  .0024  .0004 .0001 .0000
+ 6.0000
+-2.2128-1.0518 -.1760 1.0582  .4200  .1601  .0512  .0134  .0029  .0005 .0001 .0000
+ 6.2000
+-2.2443-1.0757 -.1805 1.0790  .4366  .1718  .0570  .0155  .0034  .0006 .0001 .0000
+ 6.4000
+-2.2749-1.0990 -.1851 1.0988  .4525  .1833  .0630  .0178  .0041  .0008 .0001 .0000
+ 6.6000
+-2.3047-1.1218 -.1897 1.1177  .4677  .1947  .0693  .0202  .0048  .0009 .0002 .0000
+ 6.8000
+-2.3338-1.1439 -.1943 1.1359  .4823  .2059  .0757  .0229  .0057  .0011 .0002 .0000
+ 7.0000
+-2.3622-1.1655 -.1989 1.1532  .4964  .2169  .0823  .0257  .0066  .0014 .0002 .0000
+ 7.2000
+-2.3898-1.1867 -.2035 1.1698  .5100  .2276  .0889  .0287  .0076  .0016 .0003 .0000
+ 7.4000
+-2.4169-1.2074 -.2082 1.1856  .5231  .2381  .0957  .0319  .0087  .0020 .0004 .0001
+ 7.6000
+-2.4433-1.2276 -.2129 1.2008  .5358  .2483  .1025  .0353  .0100  .0023 .0004 .0001
+ 7.8000
+-2.4691-1.2474 -.2176 1.2154  .5481  .2582  .1094  .0388  .0113  .0027 .0005 .0001
+ 8.0000
+-2.4943-1.2668 -.2223 1.2293  .5600  .2678  .1163  .0424  .0127  .0031 .0006 .0001
+ 8.2000
+-2.5190-1.2858 -.2270 1.2426  .5716  .2772  .1231  .0462  .0143  .0036 .0008 .0001
+ 8.4000
+-2.5432-1.3044 -.2317 1.2554  .5828  .2863  .1299  .0501  .0159  .0042 .0009 .0002
+ 8.6000
+-2.5668-1.3226 -.2364 1.2677  .5938  .2952  .1367  .0541  .0177  .0048 .0011 .0002
+ 8.8000
+-2.5900-1.3405 -.2411 1.2795  .6044  .3039  .1434  .0582  .0196  .0054 .0012 .0002
+ 9.0000
+-2.6127-1.3580 -.2458 1.2909  .6148  .3123  .1500  .0624  .0215  .0061 .0014 .0003
+ 9.2000
+-2.6349-1.3752 -.2504 1.3018  .6249  .3205  .1565  .0666  .0236  .0069 .0017 .0003
+ 9.4000
+-2.6567-1.3920 -.2551 1.3123  .6348  .3285  .1629  .0709  .0258  .0077 .0019 .0004
+ 9.6000
+-2.6781-1.4085 -.2597 1.3225  .6443  .3364  .1692  .0752  .0281  .0086 .0022 .0005
+ 9.8000
+-2.6991-1.4247 -.2642 1.3323  .6536  .3440  .1754  .0796  .0304  .0096 .0025 .0005
+10.0000
+-2.7196-1.4407 -.2688 1.3418  .6627  .3515  .1815  .0840  .0328  .0106 .0029 .0006
+10.2000
+-2.7398-1.4563 -.2732 1.3509  .6715  .3588  .1875  .0883  .0354  .0117 .0032 .0007
+10.4000
+-2.7597-1.4716 -.2777 1.3598  .6801  .3660  .1934  .0927  .0379  .0129 .0036 .0009
+10.6000
+-2.7791-1.4867 -.2820 1.3684  .6885  .3730  .1991  .0971  .0406  .0141 .0041 .0010
+10.8000
+-2.7983-1.5015 -.2864 1.3767  .6967  .3799  .2048  .1014  .0433  .0154 .0046 .0011
+11.0000
+-2.8171-1.5160 -.2907 1.3847  .7046  .3867  .2103  .1057  .0461  .0168 .0051 .0013
+11.2000
+-2.8355-1.5303 -.2949 1.3924  .7124  .3933  .2158  .1101  .0489  .0182 .0056 .0015
+11.4000
+-2.8537-1.5444 -.2991 1.3999  .7199  .3998  .2211  .1143  .0517  .0197 .0062 .0016
+11.6000
+-2.8716-1.5582 -.3032 1.4072  .7273  .4062  .2263  .1185  .0546  .0212 .0069 .0019
+11.8000
+-2.8892-1.5718 -.3073 1.4142  .7345  .4124  .2316  .1227  .0575  .0228 .0075 .0021
+12.0000
+-2.9065-1.5852 -.3114 1.4210  .7415  .4186  .2367  .1268  .0605  .0245 .0083 .0023
+12.2000
+-2.9236-1.5984 -.3154 1.4276  .7484  .4246  .2416  .1308  .0634  .0262 .0090 .0026
+12.4000
+-2.9404-1.6113 -.3194 1.4340  .7551  .4305  .2466  .1349  .0663  .0279 .0098 .0029
+12.6000
+-2.9569-1.6241 -.3233 1.4401  .7617  .4362  .2514  .1388  .0693  .0297 .0107 .0032
+12.8000
+-2.9733-1.6367 -.3272 1.4460  .7682  .4419  .2562  .1427  .0723  .0316 .0116 .0035
+13.0000
+-2.9893-1.6492 -.3311 1.4518  .7745  .4474  .2609  .1466  .0753  .0334 .0125 .0039
+13.2000
+-3.0052-1.6614 -.3349 1.4573  .7807  .4529  .2655  .1504  .0782  .0353 .0135 .0043
+13.4000
+-3.0209-1.6735 -.3387 1.4627  .7868  .4582  .2702  .1542  .0812  .0373 .0145 .0047
+13.6000
+-3.0363-1.6854 -.3425 1.4679  .7928  .4634  .2746  .1579  .0841  .0393 .0155 .0051
+13.8000
+-3.0515-1.6972 -.3463 1.4729  .7987  .4685  .2791  .1615  .0870  .0413 .0166 .0056
+14.0000
+-3.0666-1.7088 -.3500 1.4777  .8044  .4736  .2835  .1651  .0898  .0433 .0177 .0061
+14.2000
+-3.0814-1.7203 -.3537 1.4824  .8101  .4785  .2878  .1687  .0928  .0454 .0189 .0066
+14.4000
+-3.0961-1.7316 -.3574 1.4869  .8157  .4833  .2920  .1722  .0956  .0474 .0201 .0072
+14.6000
+-3.1106-1.7428 -.3611 1.4913  .8212  .4881  .2962  .1758  .0985  .0495 .0214 .0078
+14.8000
+-3.1249-1.7538 -.3648 1.4955  .8265  .4927  .3004  .1792  .1013  .0515 .0226 .0084
+15.0000
+-3.1390-1.7647 -.3684 1.4996  .8318  .4973  .3045  .1825  .1040  .0537 .0239 .0090
+15.2000
+-3.1529-1.7755 -.3720 1.5036  .8371  .5018  .3085  .1859  .1068  .0558 .0253 .0097
+15.4000
+-3.1667-1.7862 -.3756 1.5074  .8422  .5062  .3124  .1893  .1096  .0579 .0266 .0104
+15.6000
+-3.1804-1.7967 -.3792 1.5112  .8472  .5105  .3163  .1926  .1123  .0600 .0280 .0111
+15.8000
+-3.1938-1.8071 -.3827 1.5148  .8522  .5148  .3202  .1959  .1149  .0621 .0293 .0118
+16.0000
+-3.2071-1.8174 -.3862 1.5183  .8570  .5190  .3239  .1991  .1176  .0642 .0307 .0126
+16.2000
+-3.2203-1.8276 -.3897 1.5217  .8618  .5231  .3276  .2023  .1202  .0661 .0323 .0134
+16.4000
+-3.2333-1.8376 -.3932 1.5251  .8665  .5272  .3313  .2055  .1228  .0683 .0336 .0142
+16.6000
+-3.2461-1.8476 -.3966 1.5283  .8711  .5312  .3349  .2086  .1255  .0705 .0352 .0151
+16.8000
+-3.2589-1.8574 -.4001 1.5314  .8757  .5351  .3384  .2118  .1279  .0723 .0367 .0159
+17.0000
+-3.2714-1.8671 -.4034 1.5345  .8801  .5390  .3419  .2148  .1306  .0746 .0381 .0168
+17.2000
+-3.2839-1.8767 -.4068 1.5375  .8845  .5429  .3454  .2179  .1330  .0766 .0396 .0177
+17.4000
+-3.2962-1.8862 -.4102 1.5404  .8888  .5467  .3487  .2209  .1355  .0786 .0412 .0186
+17.6000
+-3.3083-1.8956 -.4135 1.5432  .8930  .5504  .3521  .2238  .1380  .0806 .0428 .0197
+17.8000
+-3.3203-1.9049 -.4168 1.5459  .8972  .5542  .3553  .2268  .1404  .0826 .0441 .0207
+18.0000
+-3.3322-1.9141 -.4200 1.5486  .9012  .5578  .3585  .2297  .1429  .0846 .0460 .0215
+18.2000
+-3.3440-1.9232 -.4233 1.5512  .9052  .5614  .3617  .2326  .1453  .0866 .0474 .0226
+18.4000
+-3.3557-1.9322 -.4265 1.5537  .9092  .5650  .3649  .2354  .1477  .0885 .0488 .0237
+18.6000
+-3.3672-1.9411 -.4297 1.5562  .9130  .5685  .3679  .2382  .1501  .0905 .0503 .0248
+18.8000
+-3.3786-1.9499 -.4328 1.5586  .9168  .5720  .3710  .2410  .1524  .0923 .0517 .0258
+19.0000
+-3.3899-1.9586 -.4360 1.5609  .9206  .5755  .3740  .2437  .1548  .0943 .0534 .0269
+19.2000
+-3.4011-1.9672 -.4391 1.5632  .9242  .5789  .3769  .2464  .1571  .0963 .0550 .0280
+19.4000
+-3.4121-1.9757 -.4422 1.5654  .9278  .5823  .3799  .2491  .1594  .0982 .0565 .0289
+19.6000
+-3.4231-1.9842 -.4452 1.5675  .9314  .5856  .3827  .2518  .1617  .1000 .0581 .0300
+19.8000
+-3.4339-1.9925 -.4483 1.5696  .9348  .5889  .3856  .2544  .1640  .1019 .0595 .0314
+20.0000
+-3.4446-2.0008 -.4513 1.5716  .9383  .5921  .3884  .2569  .1663  .1038 .0611 .0325
diff --git a/examples/ni111_2x2O/PHASE/O_CO_Pendry_cs.phs b/examples/ni111_2x2O/PHASE/O_CO_Pendry_cs.phs
new file mode 100644
index 0000000..5ba9e2e
--- /dev/null
+++ b/examples/ni111_2x2O/PHASE/O_CO_Pendry_cs.phs
@@ -0,0 +1,151 @@
+75  9   neng lmax (O from CO, Z = 8, 16 amu, potential (PENDRY), CS = 5.194)
+  .2000
+  .5302  .0417  .0001  .0000  .0000  .0000  .0000  .0000  .0000  .0000
+  .4000
+  .5618  .1284  .0008  .0000  .0000  .0000  .0000  .0000  .0000  .0000
+  .6000
+  .5469  .2535  .0021  .0000  .0000  .0000  .0000  .0000  .0000  .0000
+  .8000
+  .5178  .4102  .0042  .0001  .0000  .0000  .0000  .0000  .0000  .0000
+ 1.0000
+  .4835  .5835  .0071  .0002  .0000  .0000  .0000  .0000  .0000  .0000
+ 1.2000
+  .4475  .7537  .0110  .0003  .0000  .0000  .0000  .0000  .0000  .0000
+ 1.4000
+  .4112  .9040  .0157  .0005  .0000  .0000  .0000  .0000  .0000  .0000
+ 1.6000
+  .3753 1.0271  .0214  .0008  .0000  .0000  .0000  .0000  .0000  .0000
+ 1.8000
+  .3402 1.1231  .0280  .0011  .0000  .0000  .0000  .0000  .0000  .0000
+ 2.0000
+  .3060 1.1961  .0356  .0016  .0001  .0000  .0000  .0000  .0000  .0000
+ 2.2000
+  .2728 1.2508  .0440  .0022  .0001  .0000  .0000  .0000  .0000  .0000
+ 2.4000
+  .2406 1.2916  .0533  .0029  .0001  .0000  .0000  .0000  .0000  .0000
+ 2.6000
+  .2094 1.3217  .0634  .0037  .0002  .0000  .0000  .0000  .0000  .0000
+ 2.8000
+  .1792 1.3438  .0743  .0046  .0002  .0000  .0000  .0000  .0000  .0000
+ 3.0000
+  .1500 1.3598  .0859  .0057  .0003  .0000  .0000  .0000  .0000  .0000
+ 3.2000
+  .1216 1.3711  .0982  .0070  .0004  .0000  .0000  .0000  .0000  .0000
+ 3.4000
+  .0942 1.3788  .1111  .0084  .0005  .0000  .0000  .0000  .0000  .0000
+ 3.6000
+  .0676 1.3838  .1245  .0099  .0006  .0000  .0000  .0000  .0000  .0000
+ 3.8000
+  .0419 1.3866  .1385  .0116  .0008  .0000  .0000  .0000  .0000  .0000
+ 4.0000
+  .0169 1.3877  .1528  .0135  .0010  .0001  .0000  .0000  .0000  .0000
+ 4.2000
+ -.0074 1.3875  .1674  .0156  .0012  .0001  .0000  .0000  .0000  .0000
+ 4.4000
+ -.0309 1.3863  .1824  .0178  .0014  .0001  .0000  .0000  .0000  .0000
+ 4.6000
+ -.0537 1.3844  .1974  .0202  .0017  .0001  .0000  .0000  .0000  .0000
+ 4.8000
+ -.0759 1.3818  .2126  .0228  .0020  .0001  .0000  .0000  .0000  .0000
+ 5.0000
+ -.0975 1.3788  .2279  .0255  .0023  .0002  .0000  .0000  .0000  .0000
+ 5.2000
+ -.1184 1.3755  .2431  .0284  .0027  .0002  .0000  .0000  .0000  .0000
+ 5.4000
+ -.1388 1.3720  .2582  .0315  .0031  .0002  .0000  .0000  .0000  .0000
+ 5.6000
+ -.1586 1.3682  .2731  .0347  .0035  .0003  .0000  .0000  .0000  .0000
+ 5.8000
+ -.1779 1.3644  .2878  .0381  .0040  .0003  .0000  .0000  .0000  .0000
+ 6.0000
+ -.1966 1.3606  .3023  .0417  .0046  .0004  .0000  .0000  .0000  .0000
+ 6.2000
+ -.2148 1.3567  .3165  .0454  .0051  .0004  .0000  .0000  .0000  .0000
+ 6.4000
+ -.2326 1.3528  .3303  .0493  .0058  .0005  .0000  .0000  .0000  .0000
+ 6.6000
+ -.2499 1.3490  .3438  .0533  .0064  .0006  .0000  .0000  .0000  .0000
+ 6.8000
+ -.2667 1.3452  .3569  .0574  .0072  .0007  .0000  .0000  .0000  .0000
+ 7.0000
+ -.2832 1.3415  .3696  .0617  .0079  .0008  .0001  .0000  .0000  .0000
+ 7.2000
+ -.2992 1.3380  .3818  .0661  .0088  .0009  .0001  .0000  .0000  .0000
+ 7.4000
+ -.3148 1.3345  .3936  .0706  .0096  .0010  .0001  .0000  .0000  .0000
+ 7.6000
+ -.3300 1.3311  .4049  .0752  .0106  .0011  .0001  .0000  .0000  .0000
+ 7.8000
+ -.3448 1.3278  .4158  .0799  .0115  .0013  .0001  .0000  .0000  .0000
+ 8.0000
+ -.3593 1.3246  .4262  .0847  .0126  .0014  .0001  .0000  .0000  .0000
+ 8.2000
+ -.3734 1.3215  .4362  .0896  .0137  .0016  .0001  .0000  .0000  .0000
+ 8.4000
+ -.3872 1.3186  .4458  .0946  .0148  .0018  .0002  .0000  .0000  .0000
+ 8.6000
+ -.4007 1.3157  .4549  .0996  .0160  .0020  .0002  .0000  .0000  .0000
+ 8.8000
+ -.4138 1.3130  .4635  .1047  .0173  .0022  .0002  .0000  .0000  .0000
+ 9.0000
+ -.4266 1.3103  .4718  .1098  .0186  .0024  .0002  .0000  .0000  .0000
+ 9.2000
+ -.4392 1.3078  .4797  .1150  .0200  .0026  .0003  .0000  .0000  .0000
+ 9.4000
+ -.4515 1.3053  .4871  .1202  .0214  .0029  .0003  .0000  .0000  .0000
+ 9.6000
+ -.4634 1.3029  .4942  .1254  .0229  .0032  .0003  .0000  .0000  .0000
+ 9.8000
+ -.4752 1.3006  .5009  .1307  .0244  .0035  .0004  .0000  .0000  .0000
+10.0000
+ -.4866 1.2984  .5073  .1360  .0260  .0038  .0004  .0000  .0000  .0000
+10.2000
+ -.4979 1.2963  .5134  .1412  .0276  .0041  .0005  .0000  .0000  .0000
+10.4000
+ -.5088 1.2942  .5191  .1465  .0293  .0044  .0005  .0000  .0000  .0000
+10.6000
+ -.5196 1.2922  .5245  .1517  .0310  .0048  .0006  .0001  .0000  .0000
+10.8000
+ -.5301 1.2902  .5296  .1570  .0328  .0052  .0006  .0001  .0000  .0000
+11.0000
+ -.5405 1.2883  .5345  .1622  .0347  .0056  .0007  .0001  .0000  .0000
+11.2000
+ -.5506 1.2864  .5391  .1673  .0365  .0060  .0008  .0001  .0000  .0000
+11.4000
+ -.5605 1.2846  .5435  .1725  .0385  .0064  .0008  .0001  .0000  .0000
+11.6000
+ -.5702 1.2828  .5476  .1775  .0405  .0069  .0009  .0001  .0000  .0000
+11.8000
+ -.5797 1.2811  .5516  .1826  .0425  .0074  .0010  .0001  .0000  .0000
+12.0000
+ -.5891 1.2794  .5553  .1876  .0445  .0079  .0011  .0001  .0000  .0000
+12.2000
+ -.5982 1.2777  .5588  .1925  .0466  .0084  .0012  .0001  .0000  .0000
+12.4000
+ -.6073 1.2760  .5622  .1973  .0488  .0090  .0013  .0001  .0000  .0000
+12.6000
+ -.6161 1.2744  .5653  .2021  .0510  .0096  .0014  .0002  .0000  .0000
+12.8000
+ -.6248 1.2727  .5684  .2068  .0532  .0102  .0015  .0002  .0000  .0000
+13.0000
+ -.6333 1.2711  .5712  .2115  .0554  .0108  .0016  .0002  .0000  .0000
+13.2000
+ -.6417 1.2695  .5740  .2161  .0577  .0114  .0017  .0002  .0000  .0000
+13.4000
+ -.6500 1.2679  .5766  .2205  .0600  .0121  .0019  .0002  .0000  .0000
+13.6000
+ -.6581 1.2663  .5791  .2250  .0623  .0128  .0020  .0002  .0000  .0000
+13.8000
+ -.6661 1.2647  .5814  .2293  .0647  .0135  .0021  .0003  .0000  .0000
+14.0000
+ -.6740 1.2631  .5837  .2335  .0671  .0142  .0023  .0003  .0000  .0000
+14.2000
+ -.6817 1.2616  .5859  .2377  .0695  .0150  .0025  .0003  .0000  .0000
+14.4000
+ -.6893 1.2600  .5880  .2417  .0719  .0158  .0026  .0003  .0000  .0000
+14.6000
+ -.6968 1.2584  .5900  .2457  .0744  .0166  .0028  .0004  .0000  .0000
+14.8000
+ -.7042 1.2568  .5919  .2496  .0769  .0174  .0030  .0004  .0000  .0000
+15.0000
+ -.7115 1.2552  .5937  .2534  .0793  .0182  .0032  .0004  .0000  .0000
diff --git a/examples/ni111_2x2O/run_cleedpy.sh b/examples/ni111_2x2O/run_cleedpy.sh
new file mode 100755
index 0000000..9a9ec7b
--- /dev/null
+++ b/examples/ni111_2x2O/run_cleedpy.sh
@@ -0,0 +1,3 @@
+export CLEED_PHASE=./PHASE
+
+cleedpy-leed -i Ni111_2x2O.inp -b Ni111_2x2O.inp -o cleedpy.res