Initial python api support

Global variables and read_input function

meson.build

@@ -0,0 +1,22 @@
+project('edipy2', 'fortran')
+
+python = import('python').find_installation(pure: false)
+
+scifor_dep= dependency('scifor', required:true)
+edipack_dep= dependency('edipack', required:true)
+
+fortran_src = ['src/ED_INPUT_VARS.f90','src/edi2py/edi2py.f90']
+python_src  = ['python/edi2py.py', 'python/__init__.py']
+
+library('edi2py', 
+                    fortran_src,
+                    fortran_args: ['-ffree-line-length-none', '-cpp', '-D_MPI'],
+                    dependencies: [scifor_dep,edipack_dep],
+                    install: true,
+                    install_dir: python.get_install_dir() / 'edipy2'
+)
+
+python.install_sources(
+    python_src,
+    subdir: 'edipy2'
+)

pyproject.toml

@@ -0,0 +1,16 @@
+[build-system]
+build-backend = 'mesonpy'
+requires = ['meson-python']
+
+
+[project]
+name = 'edipy'
+version = '4.0.0'
+description = 'Edipack python API'
+authors = [
+  {name = 'Lorenzo Crippa', email = 'crippa.lorenzo@gmail.com'},
+  {name = 'Adriano Amaricci', email = 'amaricci@sissa.it'},
+]
+
+[project.urls]
+Repository = "https://github.com/aamaricci/EDIpack.git"

python/__init__.py

@@ -0,0 +1,7 @@
+from edipy2 import edi2py
+from ctypes import *
+import numpy as np
+
+#this ed class contains all the global variables and methods
+
+global_env = edi2py.global_env

python/edi2py.py

@@ -0,0 +1,112 @@
+from ctypes import *
+import numpy as np
+import os,sys
+import types
+
+#################################
+#link to access global variables
+#################################
+
+#dummy class, to be filled
+class Link:
+    pass
+
+#function that will add a variable to the dummy class, will be called in variable definition
+def add_global_variable(obj, dynamic_name, target_object, target_attribute):
+    @property
+    def getter(self):
+        try:
+            attrib = getattr(target_object, target_attribute)
+            try: #this is for strings
+                attrib = attrib.decode()
+            except:
+                pass
+        except: #this is for arrays
+            if(len(target_object)>1):
+                return [target_object[x] for x in range(len(target_object))]
+        return attrib
+
+    @getter.setter
+    def setter(self, new_value):
+        try: #this is for arrays
+            if(len(target_object)>1):
+                minlength=min(len(target_object),len(new_value))
+                target_object[0:minlength]=new_value[0:minlength]
+        except:
+            try:
+                new_value = new_value.encode()
+            except:
+                pass
+            setattr(target_object, target_attribute, new_value)
+
+    # Dynamically add the property to the class
+    setattr(obj.__class__, dynamic_name, getter)
+    setattr(obj.__class__, dynamic_name, setter)
+
+
+######################################
+# Load shared library with C-bindings
+######################################
+
+system = sys.platform
+libext = '.so' 
+if(system=='darwin'):
+    libext = '.dylib'
+libpath = os.path.dirname(os.path.realpath(__file__))
+libfile = os.path.join(libpath, 'libedi2py'+libext)
+libedi2py = CDLL(libfile)
+
+
+######################################
+# READ_INPUT
+######################################
+
+read_input_wrap = libedi2py.read_input
+read_input_wrap.argtypes = [c_char_p]  
+read_input_wrap.restype = None
+
+def read_input(self,input_string):
+    c_string = c_char_p(input_string.encode())
+    read_input_wrap(c_string)
+
+
+
+####################################################################
+# Create the global_env class (this is what the python module sees)
+####################################################################
+
+global_env=Link()
+
+#variables
+add_global_variable(global_env, "Nbath", c_int.in_dll(libedi2py, "Nbath"), "value")
+add_global_variable(global_env, "Norb", c_int.in_dll(libedi2py, "Norb"), "value")
+add_global_variable(global_env, "Nspin", c_int.in_dll(libedi2py, "Nspin"), "value")
+add_global_variable(global_env, "Nloop", c_int.in_dll(libedi2py, "Nloop"), "value")
+add_global_variable(global_env, "Nph", c_int.in_dll(libedi2py, "Nph"), "value")
+add_global_variable(global_env, "Nsuccess", c_int.in_dll(libedi2py, "Nsuccess"), "value")
+add_global_variable(global_env, "Lmats", c_int.in_dll(libedi2py, "Lmats"), "value")
+add_global_variable(global_env, "Lreal", c_int.in_dll(libedi2py, "Lreal"), "value")
+add_global_variable(global_env, "Ltau", c_int.in_dll(libedi2py, "Ltau"), "value")
+add_global_variable(global_env, "Lpos", c_int.in_dll(libedi2py, "Lpos"), "value")
+add_global_variable(global_env, "LOGfile", c_int.in_dll(libedi2py, "LOGfile"), "value")
+
+add_global_variable(global_env, "Uloc", ARRAY(c_double, 5).in_dll(libedi2py, "Uloc"), "value")
+add_global_variable(global_env, "Ust", c_double.in_dll(libedi2py, "Ust"), "value")
+add_global_variable(global_env, "Jh", c_double.in_dll(libedi2py, "Jh"), "value")
+add_global_variable(global_env, "Jx", c_double.in_dll(libedi2py, "Jx"), "value")
+add_global_variable(global_env, "Jp", c_double.in_dll(libedi2py, "Jp"), "value")
+add_global_variable(global_env, "xmu", c_double.in_dll(libedi2py, "xmu"), "value")
+add_global_variable(global_env, "beta", c_double.in_dll(libedi2py, "beta"), "value")
+add_global_variable(global_env, "dmft_error", c_double.in_dll(libedi2py, "dmft_error"), "value")
+add_global_variable(global_env, "eps", c_double.in_dll(libedi2py, "eps"), "value")
+add_global_variable(global_env, "wini", c_double.in_dll(libedi2py, "wini"), "value")
+add_global_variable(global_env, "wfin", c_double.in_dll(libedi2py, "wfin"), "value")
+add_global_variable(global_env, "xmin", c_double.in_dll(libedi2py, "xmin"), "value")
+add_global_variable(global_env, "xmax", c_double.in_dll(libedi2py, "xmax"), "value")
+add_global_variable(global_env, "sb_field", c_double.in_dll(libedi2py, "sb_field"), "value")
+add_global_variable(global_env, "nread", c_double.in_dll(libedi2py, "nread"), "value")
+
+
+#functions
+global_env.read_input = types.MethodType(read_input, global_env)
+

src/ED_INPUT_VARS.f90

@@ -3,40 +3,47 @@ MODULE ED_INPUT_VARS
   USE SF_PARSE_INPUT
   USE SF_IOTOOLS, only:str,free_unit
   USE ED_VERSION
+  use iso_c_binding
   implicit none
 
 
   !input variables
   !=========================================================
-  integer              :: Nbath               !Nbath=# of bath sites (per orbital or not depending on bath_type)
-  integer              :: Norb                !Norb =# of impurity orbitals
-  integer              :: Nspin               !Nspin=# spin degeneracy (max 2)
-  integer              :: nloop               !max dmft loop variables
-  integer              :: Nph                 !max number of phonons allowed (cut off)
-  real(8),allocatable  :: Uloc(:)             !local interactions
-  real(8)              :: Ust                 !intra-orbitals interactions
-  real(8)              :: Jh                  !J_Hund: Hunds' coupling constant 
-  real(8)              :: Jx                  !J_X: coupling constant for the spin-eXchange interaction term
-  real(8)              :: Jp                  !J_P: coupling constant for the Pair-hopping interaction term 
-  real(8)              :: xmu                 !chemical potential
-  real(8)              :: beta                !inverse temperature
+  integer(c_int), bind(c, name="Nbath")                              :: Nbath               !Nbath=# of bath sites (per orbital or not depending on bath_type)
+  integer(c_int), bind(c, name="Norb")                               :: Norb                !Norb =# of impurity orbitals
+  integer(c_int), bind(c, name="Nspin")                              :: Nspin               !Nspin=# spin degeneracy (max 2)
+
+  integer(c_int), bind(c, name="Nloop")                              :: Nloop               !max dmft loop variables
+  integer(c_int), bind(c, name="Nph")                                :: Nph                 !max number of phonons allowed (cut off)
+  real(c_double),dimension(5),bind(c, name="Uloc")                   :: Uloc                !local interactions
+  real(c_double),bind(c, name="Ust")                                 :: Ust                 !intra-orbitals interactions
+  real(c_double),bind(c, name="Jh")                                  :: Jh                  !J_Hund: Hunds' coupling constant 
+  real(c_double),bind(c, name="Jx")                                  :: Jx                  !J_X: coupling constant for the spin-eXchange interaction term
+  real(c_double),bind(c, name="Jp")                                  :: Jp                  !J_P: coupling constant for the Pair-hopping interaction term 
+  real(c_double),bind(c, name="xmu")                                 :: xmu                 !chemical potential
+  real(c_double),bind(c, name="beta")                                :: beta                !inverse temperature
   !
-  integer              :: ph_type             !shape of the e part of the e-ph interaction: 1=orbital occupation, 2=orbital hybridization
-  complex(8),allocatable  :: g_ph(:,:)        !g_ph: electron-phonon coupling constant all
-  real(8),allocatable  :: g_ph_diag(:)        !g_ph: electron-phonon coupling constant diagonal (density)
-  real(8)              :: w0_ph               !w0_ph: phonon frequency (constant)
-  real(8)              :: A_ph                !A_ph: phonon field coupled to displacement operator (constant)
   !
-  integer              :: Nsuccess            !# of repeated success to fall below convergence threshold  
-  real(8)              :: dmft_error          !dmft convergence threshold
-  real(8)              :: eps                 !broadening
-  real(8)              :: wini,wfin           !frequency range
-  real(8)              :: xmin,xmax           !x-range for the local lattice probability distribution function (phonons)
+  integer(c_int), bind(c, name="Nsuccess")                           :: Nsuccess            !# of repeated success to fall below convergence threshold  
+  real(c_double),bind(c, name="dmft_error")                          :: dmft_error          !dmft convergence threshold
+  real(c_double),bind(c, name="eps")                                 :: eps                 !broadening
+  real(c_double),bind(c, name="wini")                                :: wini           !frequency range min
+  real(c_double),bind(c, name="wfin")                                :: wfin                !frequency range max
+  real(c_double),bind(c, name="xmin")                                :: xmin                !x-range for the local lattice probability distribution function (phonons)
+  real(c_double),bind(c, name="xmax")                                :: xmax                !x-range for the local lattice probability distribution function (phonons)
+  real(c_double),bind(c, name="sb_field")                            :: sb_field            !symmetry breaking field
+  real(c_double),bind(c, name="nread")                               :: nread               !fixed density. if 0.d0 fixed chemical potential calculation.
+  !
   logical              :: HFmode              !flag for HF interaction form U(n-1/2)(n-1/2) VS Unn
   real(8)              :: cutoff              !cutoff for spectral summation
   real(8)              :: gs_threshold        !Energy threshold for ground state degeneracy loop up
   real(8)              :: deltasc             !breaking symmetry field
-  real(8)              :: sb_field            !symmetry breaking field
+  !
+  integer              :: ph_type             !shape of the e part of the e-ph interaction: 1=orbital occupation, 2=orbital hybridization
+  real(8)              :: A_ph                !A_ph: phonon field coupled to displacement operator (constant)
+  complex(8),allocatable  :: g_ph(:,:)        !g_ph: electron-phonon coupling constant all
+  real(8)              :: w0_ph               !w0_ph: phonon frequency (constant)
+  real(8),allocatable  :: g_ph_diag(:)        !g_ph: electron-phonon coupling constant diagonal (density)
   !
   real(8),allocatable  :: spin_field_x(:)        !magnetic field per orbital coupling to X-spin component
   real(8),allocatable  :: spin_field_y(:)        !magnetic field per orbital coupling to Y-spin component
@@ -92,7 +99,6 @@ MODULE ED_INPUT_VARS
   logical              :: finiteT             !flag for finite temperature calculation
   character(len=7)     :: bath_type           !flag to set bath type: normal (1bath/imp), hybrid(1bath)
   !
-  real(8)              :: nread               !fixed density. if 0.d0 fixed chemical potential calculation.
   real(8)              :: nerr                !fix density threshold. a loop over from 1.d-1 to required nerr is performed
   real(8)              :: ndelta              !initial chemical potential step
   real(8)              :: ncoeff              !multiplier for the initial ndelta read from a file (ndelta-->ndelta*ncoeff)
@@ -102,17 +108,17 @@ MODULE ED_INPUT_VARS
   real(8)              :: Jz_max_value
 
   !Some parameters for function dimension:
-  !=========================================================
-  integer              :: Lmats
-  integer              :: Lreal
-  integer              :: Lfit
-  integer              :: Ltau
-  integer              :: Lpos
+  integer(c_int),bind(c, name="Lmats")             :: Lmats
+  integer(c_int),bind(c, name="Lreal")             :: Lreal
+  integer(c_int),bind(c, name="Lfit")              :: Lfit
+
+  integer(c_int),bind(c, name="Ltau")              :: Ltau
+  integer(c_int),bind(c, name="Lpos")              :: Lpos
 
   !LOG AND Hamiltonian UNITS
   !=========================================================
   character(len=100)   :: Hfile,HLOCfile,SectorFile,GPHfile
-  integer,save         :: LOGfile
+  integer(c_int),bind(c, name="LOGfile"),save             :: LOGfile
 
   !THIS IS JUST A RELOCATED GLOBAL VARIABLE
   character(len=200)                                 :: ed_input_file=""
@@ -153,8 +159,9 @@ subroutine ed_read_input(INPUTunit)
     call parse_input_variable(Nph,"NPH",INPUTunit,default=0,comment="Max number of phonons allowed (cut off)")
     call parse_input_variable(bath_type,"BATH_TYPE",INPUTunit,default='normal',comment="flag to set bath type: normal (1bath/imp), hybrid(1bath), replica(1replica/imp), general(replica++)")
     !
-    allocate(Uloc(Norb))
-    call parse_input_variable(uloc,"ULOC",INPUTunit,default=(/( 2d0,i=1,size(Uloc) )/),comment="Values of the local interaction per orbital")
+    !allocate(Uloc(Norb)) #TODO: put me back!
+    !call parse_input_variable(uloc,"ULOC",INPUTunit,default=(/( 2d0,i=1,size(Uloc) )/),comment="Values of the local interaction per orbital")
+    call parse_input_variable(uloc,"ULOC",INPUTunit,default=[2d0,0d0,0d0,0d0,0d0],comment="Values of the local interaction per orbital (max 5)")
     call parse_input_variable(ust,"UST",INPUTunit,default=0.d0,comment="Value of the inter-orbital interaction term")
     call parse_input_variable(Jh,"JH",INPUTunit,default=0.d0,comment="Hunds coupling")
     call parse_input_variable(Jx,"JX",INPUTunit,default=0.d0,comment="S-E coupling")

src/edi2py/edi2py.f90

@@ -0,0 +1,61 @@
+module edi2py_bindings
+    use edipack
+    use scifor
+    use iso_c_binding
+    implicit none
+
+
+contains
+
+    !integer to logical
+    function i2l(var_integer) result (var_logical)
+      integer    :: var_integer
+      logical    :: var_logical   
+
+      if (var_integer == 1) then
+        var_logical = .true.
+      else
+        var_logical = .false.
+      endif
+    end function i2l
+
+    !logical to integer
+    function l2i(var_logical) result (var_integer)
+      integer    :: var_integer
+      logical    :: var_logical   
+
+      if (var_logical) then
+        var_integer = 1
+      else
+        var_integer = 0
+      endif
+    end function l2i
+
+    !c string to fortran string
+    subroutine c2f(c_str)
+        character(kind=c_char), dimension(*),intent(IN) :: c_str
+        character(len=120), allocatable                 :: f_str
+        integer                                         :: length
+        integer                                         :: i
+
+        length=0
+        f_str=" "
+        do
+           if (c_str(length+1) == C_NULL_CHAR) exit
+           length = length + 1
+        end do
+        do i = 1, length
+          f_str(i:i) = c_str(i)
+        enddo
+        f_str=trim(f_str)
+    end subroutine c2f
+
+    !include library functions
+    include "edi2py_read_input.f90"
+    !include "edi2py_main.f90"
+    !include "edi2py_bath.f90"
+    !include "edi2py_io.f90"
+    !include "edi2py_bath_fit.f90"
+    !include "edi2py_aux_funx.f90"
+
+end module edi2py_bindings