diff --git a/meson.build b/meson.build new file mode 100644 index 00000000..e4b05b57 --- /dev/null +++ b/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' +) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..8f6d462a --- /dev/null +++ b/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" \ No newline at end of file diff --git a/python/__init__.py b/python/__init__.py new file mode 100644 index 00000000..9a01437d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/python/edi2py.py b/python/edi2py.py new file mode 100644 index 00000000..27481ef5 --- /dev/null +++ b/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) + diff --git a/src/ED_INPUT_VARS.f90 b/src/ED_INPUT_VARS.f90 index 6ff083d5..d6442f4b 100644 --- a/src/ED_INPUT_VARS.f90 +++ b/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") diff --git a/src/edi2py/edi2py.f90 b/src/edi2py/edi2py.f90 new file mode 100644 index 00000000..1fcaaec3 --- /dev/null +++ b/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 diff --git a/src/edi2py/edi2py_aux_funx.f90 b/src/edi2py/edi2py_aux_funx.f90 new file mode 100644 index 00000000..406fb1d5 --- /dev/null +++ b/src/edi2py/edi2py_aux_funx.f90 @@ -0,0 +1,74 @@ +!ED_AUX_FUNX: +subroutine search_variable(var,ntmp,converged) bind(c, name='search_variable') + real(c_double),dimension(1) :: var(1) + real(c_double),dimension(1) :: ntmp(1) + integer(c_int),dimension(1) :: converged(1) + logical :: bool + converged(1)=0 + call ed_search_variable(var(1),ntmp(1),bool) + if (bool) converged(1)=1 +end subroutine search_variable + + + + +!AUX: +subroutine check_convergence(Xnew,dim_xnew,eps,N1,N2,oerr,convergence) bind(c, name='check_convergence') + complex(c_double_complex) :: Xnew(dim_xnew) + integer(c_int),value :: dim_xnew + real(c_double),value :: eps + integer(c_int),value :: N1,N2 + real(c_double),dimension(1) :: oerr(1) + integer(c_int),dimension(1) :: convergence(1) + integer :: i,Msum + real(c_double) :: err + real(c_double) :: M,S + complex(8),save,allocatable :: Xold(:) + integer,save :: success=0,check=1 + character(len=100) :: file_ + file_='error.err' + Msum=size(Xnew) + if(.not.allocated(Xold))then + allocate(Xold(Msum)) + Xold=0.d0 + endif + S=0.d0 ; M=0.d0 + do i=1,Msum + M=M + abs(Xnew(i)-Xold(i)) + S=S + abs(Xnew(i)) + enddo + err= M/S + Xold=Xnew + open(10,file=reg(file_),position="append") + write(10,"(I5,ES15.7)")check,err + close(10) + if(err < eps)then + success=success+1 + else + success=0 + endif + convergence(1)=0 + if(success > N1)convergence(1)=1 + if(check>=N2)then + open(10,file="ERROR.README") + write(10,*)"" + close(10) + write(*,"(A,I4,A)")"Not converged after",N2," iterations." + endif + if(convergence(1)==1)then + write(*,"(A,ES15.7)")bold_green("error="),err + else + if(err < eps)then + write(*,"(A,ES15.7)")bold_yellow("error="),err + else + write(*,"(A,ES15.7)")bold_red("error="),err + endif + endif + oerr(1)=err + check=check+1 +end subroutine check_convergence + + + + + diff --git a/src/edi2py/edi2py_bath.f90 b/src/edi2py/edi2py_bath.f90 new file mode 100644 index 00000000..af74eab8 --- /dev/null +++ b/src/edi2py/edi2py_bath.f90 @@ -0,0 +1,178 @@ +!ED_BATH: +integer(c_int) function get_bath_dimension_c() result(Nb) bind(c, name='get_bath_dimension') + Nb=ed_get_bath_dimension() +end function get_bath_dimension_c + +!H_REPLICA SETUP +subroutine init_Hreplica_direct_nn_c(Hloc) bind(c, name='init_Hreplica_direct_nn') + real(c_double),dimension(Nspin,Nspin,Norb,Norb) :: Hloc + call ed_set_Hreplica(hloc) +end subroutine init_Hreplica_direct_nn_c + +subroutine init_Hreplica_direct_so_c(Hloc) bind(c, name='init_Hreplica_direct_so') + real(c_double),dimension(Nspin*Norb,Nspin*Norb) :: Hloc + call ed_set_Hreplica(hloc) +end subroutine init_Hreplica_direct_so_c + +subroutine init_Hreplica_symmetries_site_c(Hvec,lambdavec,Nsym) bind(c, name='init_Hreplica_symmetries_site') + integer(c_int),value :: Nsym + real(c_double),dimension(Nspin,Nspin,Norb,Norb,Nsym) :: Hvec + real(c_double),dimension(Nsym) :: lambdavec + call ed_set_Hreplica(Hvec,lambdavec) +end subroutine init_Hreplica_symmetries_site_C + +subroutine init_Hreplica_symmetries_ineq_c(Hvec,lambdavec,Nlat,Nsym) bind(c, name='init_Hreplica_symmetries_ineq') + integer(c_int),value :: Nlat, Nsym + real(c_double),dimension(Nspin,Nspin,Norb,Norb,Nsym) :: Hvec + real(c_double),dimension(Nlat,Nsym) :: lambdavec + call ed_set_Hreplica(Hvec,lambdavec) +end subroutine init_Hreplica_symmetries_ineq_c + + + + +!BREAK_SYMMETRY_BATH +subroutine break_symmetry_bath_site_c(bath,dim_bath,field,sgn,sav) bind(c, name='break_symmetry_bath_site') + integer(c_int64_t) :: dim_bath(1) + real(c_double),dimension(dim_bath(1)) :: bath + real(c_double),value :: field + real(c_double),value :: sgn + integer(c_int),value :: sav + call ed_break_symmetry_bath(bath,field,sgn,i2l(sav)) +end subroutine break_symmetry_bath_site_c +! +subroutine break_symmetry_bath_ineq_c(bath,dim_bath,field,sgn,sav) bind(c, name='break_symmetry_bath_ineq') + integer(c_int64_t) :: dim_bath(2) + real(c_double),dimension(dim_bath(1),dim_bath(2)) :: bath + real(c_double),value :: field + real(c_double),value :: sgn + integer(c_int),value :: sav + call ed_break_symmetry_bath(bath,field,sgn,i2l(sav)) +end subroutine break_symmetry_bath_ineq_c + + + +!SPIN_SYMMETRIZE_BATH +subroutine spin_symmetrize_bath_site_c(bath,dim_bath,sav) bind(c, name='spin_symmetrize_bath_site') + integer(c_int64_t) :: dim_bath(1) + integer(c_int),value :: sav + real(c_double),dimension(dim_bath(1)) :: bath + call ed_spin_symmetrize_bath(bath,i2l(sav)) +end subroutine spin_symmetrize_bath_site_c +! +subroutine spin_symmetrize_bath_ineq_c(bath,dim_bath,sav) bind(c, name='spin_symmetrize_bath_ineq') + integer(c_int64_t) :: dim_bath(2) + integer(c_int),value :: sav + real(c_double),dimension(dim_bath(1),dim_bath(2)) :: bath + call ed_spin_symmetrize_bath(bath,i2l(sav)) +end subroutine spin_symmetrize_bath_ineq_c + + + +!ORB_SYMMETRIZE_BATH +subroutine orb_symmetrize_bath_site_c(bath,dim_bath,sav) bind(c, name='orb_symmetrize_bath_site') + integer(c_int64_t) :: dim_bath(1) + integer(c_int),value :: sav + real(c_double),dimension(dim_bath(1)) :: bath + call ed_orb_symmetrize_bath(bath,i2l(sav)) +end subroutine orb_symmetrize_bath_site_c +! +subroutine orb_symmetrize_bath_ineq_c(bath,dim_bath,sav) bind(c, name='orb_symmetrize_bath_ineq') + integer(c_int64_t) :: dim_bath(2) + integer(c_int),value :: sav + real(c_double),dimension(dim_bath(1),dim_bath(2)) :: bath + call ed_orb_symmetrize_bath(bath,i2l(sav)) +end subroutine orb_symmetrize_bath_ineq_c + + +!ORB_EQUALITY_BATH +subroutine orb_equality_bath_site_c(bath,dim_bath,indx,sav) bind(c, name='orb_equality_bath_site') + integer(c_int64_t) :: dim_bath(1) + real(c_double),dimension(dim_bath(1)) :: bath + integer(c_int),value :: indx,sav + call ed_orb_equality_bath(bath,indx,i2l(sav)) +end subroutine orb_equality_bath_site_c +! +subroutine orb_equality_bath_ineq_c(bath,dim_bath,indx,sav) bind(c, name='orb_equality_bath_ineq') + integer(c_int64_t) :: dim_bath(2) + real(c_double),dimension(dim_bath(1),dim_bath(2)) :: bath + integer(c_int),value :: indx,sav + call ed_orb_equality_bath(bath,indx,i2l(sav)) +end subroutine orb_equality_bath_ineq_c + + + + + +!PH_SYMMETRIZE_BATH +subroutine ph_symmetrize_bath_site_c(bath,dim_bath,sav) bind(c, name='ph_symmetrize_bath_site') + integer(c_int64_t) :: dim_bath(1) + real(c_double),dimension(dim_bath(1)) :: bath + integer(c_int),value :: sav + call ed_ph_symmetrize_bath(bath,i2l(sav)) +end subroutine ph_symmetrize_bath_site_c +! +subroutine ph_symmetrize_bath_ineq_c(bath,dim_bath,sav) bind(c, name='ph_symmetrize_bath_ineq') + integer(c_int64_t) :: dim_bath(2) + real(c_double),dimension(dim_bath(1),dim_bath(2)) :: bath + integer(c_int),value :: sav + call ed_ph_symmetrize_bath(bath,i2l(sav)) +end subroutine ph_symmetrize_bath_ineq_c + + + +!PH_TRANS_BATH +subroutine ph_trans_bath_site_c(bath,dim_bath,sav) bind(c, name='ph_trans_bath_site') + integer(c_int64_t) :: dim_bath(1) + real(c_double),dimension(dim_bath(1)) :: bath + integer(c_int),value :: sav + call ed_ph_trans_bath(bath,i2l(sav)) +end subroutine ph_trans_bath_site_c +! +subroutine ph_trans_bath_ineq_c(bath,dim_bath,sav) bind(c, name='ph_trans_bath_ineq') + integer(c_int64_t) :: dim_bath(2) + real(c_double),dimension(dim_bath(1),dim_bath(1)) :: bath + integer(c_int),value :: sav + call ed_ph_trans_bath(bath,i2l(sav)) +end subroutine ph_trans_bath_ineq_c + + + +!BATH COMPONENT ROUTINES +subroutine get_bath_component_dimension_c(instr,Ndim) bind(c, name='get_bath_component_dimension') + character(kind=c_char), dimension(1) :: instr + character(len=1) :: typ + integer(c_int64_t) :: Ndim(3) + typ(1:1)=instr(1) + Ndim = ed_get_bath_component_dimension(typ) +end subroutine get_bath_component_dimension_c + +subroutine get_bath_component(array,dim_array,bath,dim_bath,instr) bind(c, name='get_bath_component') + integer(c_int64_t) :: dim_bath(1),dim_array(3) + real(c_double),dimension(dim_array(1),dim_array(2),dim_array(3)) :: array + real(c_double),dimension(dim_bath(1)) :: bath + character(kind=c_char), dimension(1) :: instr + character(len=1) :: typ + typ(1:1)=instr(1) + call ed_get_bath_component(array,bath,typ) +end subroutine get_bath_component + +subroutine set_bath_component(array,dim_array,bath,dim_bath,instr) bind(c, name='set_bath_component') + integer(c_int64_t) :: dim_bath(1),dim_array(3) + real(c_double),dimension(dim_array(1),dim_array(2),dim_array(3)) :: array + real(c_double),dimension(dim_bath(1)) :: bath + character(kind=c_char), dimension(1) :: instr + character(len=1) :: typ + typ(1:1)=instr(1) + call ed_set_bath_component(array,bath,typ) +end subroutine set_bath_component + +subroutine copy_bath_component(bathIN,dim_bathin,bathOUT,dim_bathout,instr) bind(c, name='copy_bath_component') + real(c_double),dimension(dim_bathin) :: bathIN + real(c_double),dimension(dim_bathout) :: bathOUT + character(kind=c_char), dimension(1) :: instr + character(len=1) :: typ + integer(c_int),value :: dim_bathin, dim_bathout + typ(1:1)=instr(1) + call ed_copy_bath_component(bathIN,bathOUT,typ) +end subroutine copy_bath_component diff --git a/src/edi2py/edi2py_bath_fit.f90 b/src/edi2py/edi2py_bath_fit.f90 new file mode 100644 index 00000000..05105cfc --- /dev/null +++ b/src/edi2py/edi2py_bath_fit.f90 @@ -0,0 +1,27 @@ +!USE ED_BATH_FIT: +subroutine chi2_fitgf_site_c(func,dim_func,bath,dim_bath,ispin,iorb) bind(c, name='chi2_fitgf_site') + integer(c_int64_t) :: dim_func(5),dim_bath(1) + real(c_double),dimension(dim_bath(1)),intent(inout) :: bath + complex(c_double_complex),dimension(dim_func(1),dim_func(2),dim_func(3),dim_func(4),dim_func(5)),intent(in) :: func + integer(c_int),value :: ispin + integer(c_int),value :: iorb + call assert_shape(func,[Nspin,Nspin,Norb,Norb,Lmats],"chi2_fitgf","func") + if(iorb>0)then + call ed_chi2_fitgf(func,bath,ispin,iorb) + else + call ed_chi2_fitgf(func,bath,ispin) + endif +end subroutine chi2_fitgf_site_c + +subroutine chi2_fitgf_ineq_c(func,dim_func,bath,dim_bath,hloc,dim_hloc,ispin) bind(c, name='chi2_fitgf_ineq') + integer(c_int64_t) :: dim_func(6),dim_bath(2),dim_hloc(5) + real(c_double),dimension(dim_bath(1),dim_bath(2)),intent(inout) :: bath + complex(c_double_complex),dimension(dim_func(1),dim_func(2),dim_func(3),dim_func(4),dim_func(5),dim_func(6)),intent(inout) :: func + real(c_double),dimension(dim_hloc(1),dim_hloc(2),dim_hloc(3),dim_hloc(4),dim_hloc(5)),intent(in) :: hloc + integer(c_int),intent(in) :: ispin + integer(c_int) :: Nineq + Nineq = size(bath,1) + call assert_shape(func,[Nineq,Nspin,Nspin,Norb,Norb,Lmats],"chi2_fitgf","func") + call assert_shape(Hloc,[Nineq,Nspin,Nspin,Norb,Norb],"chi2_fitgf","hloc") + call ed_chi2_fitgf(func,bath,hloc,ispin) +end subroutine chi2_fitgf_ineq_c diff --git a/src/edi2py/edi2py_io.f90 b/src/edi2py/edi2py_io.f90 new file mode 100644 index 00000000..8ea21e58 --- /dev/null +++ b/src/edi2py/edi2py_io.f90 @@ -0,0 +1,154 @@ +!ED_IO: +subroutine get_sigma_matsubara_site_c(sigma,d) bind(c, name='get_sigma_matsubara_site') + integer(c_int64_t) :: d(5) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5)),intent(inout) :: sigma + call assert_shape(sigma,[Nspin,Nspin,Norb,Norb,Lmats],"get_sigma_matsubara","Sigma") + call ed_get_sigma_matsubara(sigma) +end subroutine get_sigma_matsubara_site_c +! +subroutine get_sigma_matsubara_ineq_c(sigma,d) bind(c, name='get_sigma_matsubara_ineq') + integer(c_int64_t) :: d(6) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5),d(6)),intent(inout) :: sigma + integer(c_int) :: Nsites + Nsites=size(sigma,1) + call assert_shape(sigma,[Nsites,Nspin,Nspin,Norb,Norb,Lmats],"get_sigma_matsubara","Sigma") + call ed_get_sigma_matsubara(sigma,Nsites) +end subroutine get_sigma_matsubara_ineq_c + + + +subroutine get_sigma_realaxis_site_c(sigma,d) bind(c, name='get_sigma_realaxis_site') + integer(c_int64_t) :: d(5) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5)),intent(inout) :: sigma + call assert_shape(sigma,[Nspin,Nspin,Norb,Norb,Lreal],"get_sigma_realaxis","sigma") + call ed_get_sigma_realaxis(sigma) +end subroutine get_sigma_realaxis_site_c +! +subroutine get_sigma_realaxis_ineq_c(sigma,d) bind(c, name='get_sigma_realaxis_ineq') + integer(c_int64_t) :: d(6) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5),d(6)),intent(inout) :: sigma + integer(c_int) :: Nsites + Nsites=size(sigma,1) + call assert_shape(sigma,[Nsites,Nspin,Nspin,Norb,Norb,Lreal],"get_sigma_realaxis","sigma") + call ed_get_sigma_realaxis(sigma,Nsites) +end subroutine get_sigma_realaxis_ineq_c + + + +subroutine get_gimp_matsubara_site_c(gimp,d) bind(c, name='get_gimp_matsubara_site') + integer(c_int64_t) :: d(5) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5)),intent(inout) :: gimp + call assert_shape(gimp,[Nspin,Nspin,Norb,Norb,Lmats],"get_gimp_matsubara","gimp") + call ed_get_gimp_matsubara(gimp) +end subroutine get_gimp_matsubara_site_c +! +subroutine get_gimp_matsubara_ineq_c(gimp,d) bind(c, name='get_gimp_matsubara_ineq') + integer(c_int64_t) :: d(6) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5),d(6)),intent(inout) :: gimp + integer(c_int) :: Nsites + Nsites=size(gimp,1) + call assert_shape(gimp,[Nsites,Nspin,Nspin,Norb,Norb,Lmats],"get_gimp_matsubara","gimp") + call ed_get_gimp_matsubara(gimp,Nsites) +end subroutine get_gimp_matsubara_ineq_c + + + +subroutine get_gimp_realaxis_site_c(gimp,d) bind(c, name='get_gimp_realaxis_site') + integer(c_int64_t) :: d(5) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5)),intent(inout) :: gimp + call assert_shape(gimp,[Nspin,Nspin,Norb,Norb,Lreal],"get_gimp_realaxis","gimp") + call ed_get_gimp_realaxis(gimp) +end subroutine get_gimp_realaxis_site_c +! +subroutine get_gimp_realaxis_ineq_c(gimp,d) bind(c, name='get_gimp_realaxis_ineq') + integer(c_int64_t) :: d(6) + complex(c_double_complex),dimension(d(1),d(2),d(3),d(4),d(5),d(6)),intent(inout) :: gimp + integer(c_int) :: Nsites + Nsites=size(gimp,1) + call assert_shape(gimp,[Nsites,Nspin,Nspin,Norb,Norb,Lreal],"get_gimp_realaxis","gimp") + call ed_get_gimp_realaxis(gimp,Nsites) +end subroutine get_gimp_realaxis_ineq_c + + + +subroutine get_dens_site_c(arg,arg_dim1) bind(c, name='get_dens_site') + integer(c_int),value :: arg_dim1 + real(c_double),dimension(arg_dim1) :: arg + call assert_shape(arg,[Norb],"get_dens_site","arg") + call ed_get_dens(arg) +end subroutine get_dens_site_c +! +subroutine get_dens_ineq_c(arg,arg_dim1,arg_dim2,nlat) bind(c, name='get_dens_ineq') + integer(c_int),value :: arg_dim1,arg_dim2,nlat + real(c_double),dimension(arg_dim1,arg_dim2) :: arg + call assert_shape(arg,[Nlat,Norb],"get_dens_ineq","arg") + call ed_get_dens(arg,nlat) +end subroutine get_dens_ineq_c + + + +subroutine get_mag_site_c(arg,arg_dim1) bind(c, name='get_mag_site') + integer(c_int),value :: arg_dim1 + real(c_double),dimension(arg_dim1) :: arg + call assert_shape(arg,[Norb],"get_mag_site","arg") + call ed_get_mag(arg) +end subroutine get_mag_site_c +! +subroutine get_mag_ineq_c(arg,arg_dim1,arg_dim2,nlat) bind(c, name='get_mag_ineq') + integer(c_int),value :: arg_dim1,arg_dim2,nlat + real(c_double),dimension(arg_dim1,arg_dim2) :: arg + call assert_shape(arg,[Nlat,Norb],"get_mag_ineq","arg") + call ed_get_mag(arg,nlat) +end subroutine get_mag_ineq_c + + +subroutine get_docc_site_c(arg,arg_dim1) bind(c, name='get_docc_site') + integer(c_int),value :: arg_dim1 + real(c_double),dimension(arg_dim1) :: arg + call assert_shape(arg,[Norb],"get_doc_site","arg") + call ed_get_docc(arg) +end subroutine get_docc_site_c +! +subroutine get_docc_ineq_c(arg,arg_dim1,arg_dim2,nlat) bind(c, name='get_docc_ineq') + integer(c_int),value :: arg_dim1,arg_dim2,nlat + real(c_double),dimension(arg_dim1,arg_dim2) :: arg + call assert_shape(arg,[Nlat,Norb],"get_docc_ineq","arg") + call ed_get_docc(arg,nlat) +end subroutine get_docc_ineq_c + + + +subroutine get_eimp_site_c(arg) bind(c, name='get_eimp_site') + real(c_double),dimension(4) :: arg + call assert_shape(arg,[4],"get_eimp_site","arg") + call ed_get_eimp(arg) +end subroutine get_eimp_site_c +! +subroutine get_eimp_ineq_c(arg,nlat) bind(c, name='get_eimp_ineq') + integer(c_int),value :: nlat + real(c_double),dimension(Nlat,4) :: arg + call assert_shape(arg,[Nlat,4],"get_eimp_ineq","arg") + call ed_get_eimp(arg,nlat) +end subroutine get_eimp_ineq_c + + + +subroutine get_doubles_site_c(arg) bind(c, name='get_doubles_site') + real(c_double),dimension(4) :: arg + call assert_shape(arg,[4],"get_doubles_site","arg") + call ed_get_doubles(arg) +end subroutine get_doubles_site_c +! +subroutine get_doubles_ineq_c(arg,nlat) bind(c, name='get_doubles_ineq') + integer(c_int),value :: nlat + real(c_double),dimension(Nlat,4) :: arg + call assert_shape(arg,[Nlat,4],"get_doubles_ineq","arg") + call ed_get_doubles(arg,nlat) +end subroutine get_doubles_ineq_c + + + + + + + diff --git a/src/edi2py/edi2py_main.f90 b/src/edi2py/edi2py_main.f90 new file mode 100644 index 00000000..5884b7dc --- /dev/null +++ b/src/edi2py/edi2py_main.f90 @@ -0,0 +1,35 @@ +!ED_MAIN: +subroutine init_solver_site_c(bath,dim_bath) bind(c, name='init_solver_site') + integer(c_int64_t),dimension(1),intent(in) :: dim_bath + real(c_double),dimension(dim_bath(1)),intent(inout) :: bath + call ed_init_solver(bath) +end subroutine init_solver_site_c +! +subroutine init_solver_ineq_c(bath,dim_bath) bind(c, name='init_solver_ineq') + integer(c_int64_t),dimension(2),intent(in) :: dim_bath + real(c_double),dimension(dim_bath(1),dim_bath(2)),intent(inout) :: bath + call ed_init_solver(bath) +end subroutine init_solver_ineq_c + + +subroutine solve_site_c(bath,dim_bath,hloc,dim_hloc,sflag) bind(c, name='solve_site') + integer(c_int64_t),dimension(1),intent(in) :: dim_bath + integer(c_int64_t),dimension(4),intent(in) :: dim_hloc + integer(c_int),value :: sflag + real(c_double),dimension(dim_bath(1)),intent(in) :: bath + real(c_double),dimension(dim_hloc(1),dim_hloc(2),dim_hloc(3),dim_hloc(4)),intent(in) :: hloc + call assert_shape(hloc,[Nspin,Nspin,Norb,Norb],"solve","hloc") + call ed_solve(bath,hloc,sflag=i2l(sflag)) +end subroutine solve_site_c +! +subroutine solve_ineq_c(bath,dim_bath,hloc,dim_hloc,mpi_lanc) bind(c, name='solve_ineq') + integer(c_int64_t),dimension(2),intent(in) :: dim_bath + integer(c_int64_t),dimension(5),intent(in) :: dim_hloc + integer(c_int),value :: mpi_lanc + real(c_double),dimension(dim_bath(1),dim_bath(2)),intent(in) :: bath + real(c_double),dimension(dim_hloc(1),dim_hloc(2),dim_hloc(3),dim_hloc(4),dim_hloc(5)),intent(in) :: hloc + integer :: Nineq + Nineq = size(bath,1) + call assert_shape(Hloc,[Nineq,Nspin,Nspin,Norb,Norb],"solve","hloc") + call ed_solve(bath,hloc,i2l(mpi_lanc)) +end subroutine solve_ineq_c diff --git a/src/edi2py/edi2py_read_input.f90 b/src/edi2py/edi2py_read_input.f90 new file mode 100644 index 00000000..4bc7e77b --- /dev/null +++ b/src/edi2py/edi2py_read_input.f90 @@ -0,0 +1,18 @@ +!ED_MAIN: +subroutine ed_read_input_c(instr) bind(c, name='read_input') + character(kind=c_char), dimension(*), intent(IN) :: instr + character(len=20), allocatable :: INPUTunit + integer :: length + integer :: i + length=0 + INPUTunit=" " + do + if (instr(length+1) == C_NULL_CHAR) exit + length = length + 1 + end do + do i = 1, length + INPUTunit(i:i) = instr(i) + enddo + INPUTunit=trim(INPUTunit) + call ed_read_input(INPUTunit) +end subroutine ed_read_input_c