diff --git a/plugins/local/spher_harm/.gitignore b/plugins/local/spher_harm/.gitignore
new file mode 100644
index 00000000..1561915b
--- /dev/null
+++ b/plugins/local/spher_harm/.gitignore
@@ -0,0 +1,59 @@
+IRPF90_temp/
+IRPF90_man/
+build.ninja
+irpf90.make
+ezfio_interface.irp.f
+irpf90_entities
+tags
+Makefile
+ao_basis
+ao_one_e_ints
+ao_two_e_erf_ints
+ao_two_e_ints
+aux_quantities
+becke_numerical_grid
+bitmask
+cis
+cisd
+cipsi
+davidson
+davidson_dressed
+davidson_undressed
+density_for_dft
+determinants
+dft_keywords
+dft_utils_in_r
+dft_utils_one_e
+dft_utils_two_body
+dressing
+dummy
+electrons
+ezfio_files
+fci
+generators_cas
+generators_full
+hartree_fock
+iterations
+kohn_sham
+kohn_sham_rs
+mo_basis
+mo_guess
+mo_one_e_ints
+mo_two_e_erf_ints
+mo_two_e_ints
+mpi
+mrpt_utils
+nuclei
+perturbation
+pseudo
+psiref_cas
+psiref_utils
+scf_utils
+selectors_cassd
+selectors_full
+selectors_utils
+single_ref_method
+slave
+tools
+utils
+zmq
diff --git a/plugins/local/spher_harm/NEED b/plugins/local/spher_harm/NEED
new file mode 100644
index 00000000..92df7f12
--- /dev/null
+++ b/plugins/local/spher_harm/NEED
@@ -0,0 +1 @@
+dft_utils_in_r
diff --git a/plugins/local/spher_harm/README.rst b/plugins/local/spher_harm/README.rst
new file mode 100644
index 00000000..bf897f73
--- /dev/null
+++ b/plugins/local/spher_harm/README.rst
@@ -0,0 +1,4 @@
+==========
+spher_harm
+==========
+
diff --git a/plugins/local/spher_harm/assoc_gaus_pol.irp.f b/plugins/local/spher_harm/assoc_gaus_pol.irp.f
new file mode 100644
index 00000000..fa790307
--- /dev/null
+++ b/plugins/local/spher_harm/assoc_gaus_pol.irp.f
@@ -0,0 +1,50 @@
+double precision function plgndr(l,m,x)
+ integer, intent(in) :: l,m
+ double precision, intent(in) :: x
+ BEGIN_DOC
+ ! associated Legenre polynom P_l,m(x). Used for the Y_lm(theta,phi)
+ ! Taken from https://iate.oac.uncor.edu/~mario/materia/nr/numrec/f6-8.pdf
+ END_DOC
+ integer :: i,ll
+ double precision :: fact,pll,pmm,pmmp1,somx2
+ if(m.lt.0.or.m.gt.l.or.dabs(x).gt.1.d0)then 
+  print*,'bad arguments in plgndr'
+  pause
+ endif
+ pmm=1.d0
+ if(m.gt.0) then
+  somx2=dsqrt((1.d0-x)*(1.d0+x))
+  fact=1.d0
+  do i=1,m
+   pmm=-pmm*fact*somx2
+   fact=fact+2.d0
+  enddo 
+ endif ! m > 0
+ if(l.eq.m) then
+  plgndr=pmm 
+ else
+  pmmp1=x*(2*m+1)*pmm ! Compute P_m+1^m
+   if(l.eq.m+1) then
+    plgndr=pmmp1
+   else ! Compute P_l^m, l> m+1
+    do ll=m+2,l
+     pll=(x*dble(2*ll-1)*pmmp1-dble(ll+m-1)*pmm)/(ll-m)
+     pmm=pmmp1
+     pmmp1=pll
+    enddo
+    plgndr=pll
+   endif ! l.eq.m+1
+ endif ! l.eq.m
+ return
+end
+
+double precision function ortho_assoc_gaus_pol(l1,m1,l2)
+ implicit none
+ integer, intent(in) :: l1,m1,l2
+ double precision :: fact
+ if(l1.ne.l2)then
+  ortho_assoc_gaus_pol= 0.d0
+ else
+  ortho_assoc_gaus_pol = 2.d0*fact(l1+m1) / (dble(2*l1+1)*fact(l1-m1))
+ endif
+end
diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f
new file mode 100644
index 00000000..40661db1
--- /dev/null
+++ b/plugins/local/spher_harm/spher_harm.irp.f
@@ -0,0 +1,217 @@
+program spher_harm
+ implicit none
+ call test_spher_harm
+! call test_cart
+! call test_brutal_spheric
+end
+
+subroutine test_cart
+ implicit none
+ include 'constants.include.F'
+ double precision :: r(3),theta,phi,r_abs 
+ print*,''
+ r = 0.d0
+ r(1) = 1.d0
+ r(2) = 1.d0
+ call cartesian_to_spherical(r,theta,phi,r_abs)
+ print*,r
+ print*,phi/pi
+ print*,''
+ r = 0.d0
+ r(1) =-1.d0
+ r(2) = 1.d0
+ call cartesian_to_spherical(r,theta,phi,r_abs)
+ print*,r
+ print*,phi/pi
+ print*,''
+ r = 0.d0
+ r(1) =-1.d0
+ r(2) =-1.d0
+ call cartesian_to_spherical(r,theta,phi,r_abs)
+ print*,r
+ print*,phi/pi
+ print*,''
+ r = 0.d0
+ r(1) = 1.d0
+ r(2) =-1.d0
+ call cartesian_to_spherical(r,theta,phi,r_abs)
+ print*,r
+ print*,phi/pi
+end
+
+subroutine test_spher_harm
+ implicit none
+  include 'constants.include.F'
+ integer :: l1,m1,i,l2,m2,lmax
+ double precision :: r(3),weight,accu_re, accu_im,accu
+ double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2
+ l1 = 0
+ m1 = 0
+ l2 = 0
+ m2 = 0
+ lmax = 5
+ do l1 = 0,lmax
+  do m1 =  -l1 ,l1
+   do l2 = 0,lmax
+    do m2 =  -l2 ,l2
+     accu_re = 0.d0
+     accu_im = 0.d0
+     ! <l1,m1|l2,m2> = \int dOmega Y_l1,m1^* Y_l2,m2 
+     !               = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2)
+     !               = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2)
+     accu = 0.d0
+     do i = 1, n_points_integration_angular 
+      double precision :: theta,phi,r_abs
+      r(1:3) = angular_quadrature_points(i,1:3)
+      weight = weights_angular_points(i)
+      call cartesian_to_spherical(r,theta,phi,r_abs)
+      if(theta.gt.pi.or.theta.lt.0.d0)then
+       print*,'pb with theta',theta
+       print*,r
+      endif
+      if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then
+       print*,'pb with phi',phi/pi
+       print*,r
+      endif
+      call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1)
+      call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2)
+!      call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
+!      call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
+!      call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
+!      call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
+!      call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
+!      call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
+      accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2)
+      accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2)
+      accu += weight
+      write(33,'(10(F16.10,X))')phi/pi
+     enddo
+     if(l1.ne.l2.or.m1.ne.m2)then
+      if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then
+       print*,'pb OFF DIAG !!!!! '
+       print*,'l1,m1,l2,m2',l1,m1,l2,m2
+       print*,'accu_re = ',accu_re
+       print*,'accu_im = ',accu_im
+      endif
+     endif
+     if(l1==l2.and.m1==m2)then
+      if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then
+       print*,'pb DIAG !!!!! '
+       print*,'l1,m1,l2,m2',l1,m1,l2,m2
+       print*,'accu_re = ',accu_re
+       print*,'accu_im = ',accu_im
+      endif
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ double precision :: x,dx,xmax,xmin
+ integer:: nx
+ nx = 10000
+ xmin = -5.d0
+ xmax =  5.d0
+ dx = (xmax - xmin)/dble(nx)
+ x = xmin
+ do i = 1, nx
+  write(34,*)x,datan(x),dacos(x)
+  x += dx
+ enddo
+end
+
+subroutine test_brutal_spheric
+ implicit none
+  include 'constants.include.F'
+ integer :: itheta, iphi,ntheta,nphi
+ double precision :: theta_min, theta_max, dtheta,theta
+ double precision :: phi_min, phi_max, dphi,phi
+ double precision :: accu_re, accu_im,weight 
+ double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu
+ integer :: l1,m1,i,l2,m2,lmax
+ phi_min = 0.d0
+ phi_max = 2.D0 * pi
+ theta_min = 0.d0
+ theta_max = 1.D0 * pi
+ ntheta = 1000
+ nphi = 1000
+ dphi = (phi_max - phi_min)/dble(nphi)
+ dtheta = (theta_max - theta_min)/dble(ntheta)
+
+ lmax = 3
+ do l1 = 0,lmax
+  do m1 =  0 ,l1
+   do l2 = 0,lmax
+    do m2 =  0 ,l2
+     accu_re = 0.d0
+     accu_im = 0.d0
+     accu  = 0.d0
+     theta = theta_min
+     do itheta = 1, ntheta
+      phi = phi_min
+      do iphi = 1, nphi
+!      call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
+!      call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
+       call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1)
+       call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2)
+       weight = dtheta * dphi * dsin(theta) 
+       accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2)
+       accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2)
+       accu += weight
+       phi += dphi
+      enddo
+      theta += dtheta
+     enddo
+     print*,'l1,m1,l2,m2',l1,m1,l2,m2
+     print*,'accu_re = ',accu_re
+     print*,'accu_im = ',accu_im
+     print*,'accu    = ',accu
+     if(l1.ne.l2.or.m1.ne.m2)then
+      if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then
+       print*,'pb OFF DIAG !!!!! '
+      endif
+     endif
+     if(l1==l2.and.m1==m2)then
+      if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then
+       print*,'pb DIAG !!!!! '
+      endif
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+ 
+
+end
+
+subroutine test_assoc_leg_pol
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  print *, 'Hello world'
+ integer :: l1,m1,ngrid,i,l2,m2
+ l1 = 0
+ m1 = 0
+ l2 = 2
+ m2 = 0
+ double precision :: x, dx,xmax,accu,xmin
+ double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol
+ ngrid = 100000
+ xmax = 1.d0
+ xmin = -1.d0
+ dx = (xmax-xmin)/dble(ngrid)
+ do l2 = 0,10
+  x = xmin
+  accu = 0.d0
+  do i = 1, ngrid
+   func_1 = plgndr(l1,m1,x)
+   func_2 = plgndr(l2,m2,x)
+   write(33,*)x, func_1,func_2
+   accu += func_1 * func_2 * dx
+   x += dx
+  enddo
+  print*,'l2 = ',l2
+  print*,'accu = ',accu
+  print*,ortho_assoc_gaus_pol(l1,m1,l2)
+ enddo
+end
diff --git a/plugins/local/spher_harm/spher_harm_func.irp.f b/plugins/local/spher_harm/spher_harm_func.irp.f
new file mode 100644
index 00000000..825bd8ac
--- /dev/null
+++ b/plugins/local/spher_harm/spher_harm_func.irp.f
@@ -0,0 +1,151 @@
+subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm)
+ implicit none
+ integer, intent(in) :: l,m
+ double precision, intent(in)  :: r(3)
+ double precision, intent(out) :: re_ylm, im_ylm
+
+ double precision :: theta, phi,r_abs
+ call cartesian_to_spherical(r,theta,phi,r_abs)
+ call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
+end
+
+
+subroutine spher_harm_func_m_pos(l,m,theta,phi,re_ylm, im_ylm)
+  include 'constants.include.F'
+ implicit none
+ BEGIN_DOC
+! Y_lm(theta,phi) with m >0
+!
+ END_DOC
+ double precision, intent(in) :: theta, phi
+ integer, intent(in)          :: l,m
+ double precision, intent(out):: re_ylm,im_ylm
+ double precision :: prefact,fact,cos_theta,plgndr,p_lm
+ double precision :: tmp
+ prefact = dble(2*l+1)*fact(l-m)/(dfour_pi * fact(l+m))
+ prefact = dsqrt(prefact)
+ cos_theta = dcos(theta)
+ p_lm = plgndr(l,m,cos_theta)
+ tmp  = prefact * p_lm
+ re_ylm = dcos(dble(m)*phi) * tmp
+ im_ylm = dsin(dble(m)*phi) * tmp
+end
+
+subroutine spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
+ implicit none 
+ BEGIN_DOC
+ ! Y_lm(theta,phi) with -l<m<+l
+ !
+ END_DOC
+ double precision, intent(in) :: theta, phi
+ integer, intent(in)          :: l,m
+ double precision, intent(out):: re_ylm,im_ylm
+ double precision :: re_ylm_pos,im_ylm_pos,tmp
+ integer :: minus_m
+ if(abs(m).gt.l)then
+  print*,'|m| > l in spher_harm_func !! stopping ...'
+  stop
+ endif
+ if(m.ge.0)then
+  call spher_harm_func_m_pos(l,m,theta,phi,re_ylm_pos, im_ylm_pos)
+  re_ylm = re_ylm_pos
+  im_ylm = im_ylm_pos
+ else
+  minus_m = -m !> 0
+  call spher_harm_func_m_pos(l,minus_m,theta,phi,re_ylm_pos, im_ylm_pos)
+  tmp = (-1)**minus_m
+  re_ylm =   tmp  * re_ylm_pos
+  im_ylm =  -tmp  * im_ylm_pos ! complex conjugate 
+ endif
+end
+
+subroutine cartesian_to_spherical(r,theta,phi,r_abs)
+ implicit none 
+ double precision, intent(in) :: r(3)
+ double precision, intent(out):: theta, phi,r_abs
+ double precision :: r_2,x_2_y_2,tmp
+ include 'constants.include.F'
+ x_2_y_2 = r(1)*r(1) + r(2)*r(2)
+ r_2 = x_2_y_2 + r(3)*r(3)
+ r_abs = dsqrt(r_2)
+
+ if(r_abs.gt.1.d-20)then
+  theta = dacos(r(3)/r_abs)
+ else
+  theta = 0.d0
+ endif
+
+ if(.true.)then
+  if(dabs(r(1)).gt.0.d0)then
+   tmp = datan(r(2)/r(1))
+!   phi = datan2(r(2),r(1))
+  endif
+  ! From Wikipedia on Spherical Harmonics
+  if(r(1).gt.0.d0)then
+   phi = tmp
+  else if(r(1).lt.0.d0.and.r(2).ge.0.d0)then
+   phi = tmp + pi
+  else if(r(1).lt.0.d0.and.r(2).lt.0.d0)then
+   phi = tmp - pi
+  else if(r(1)==0.d0.and.r(2).gt.0.d0)then
+   phi = 0.5d0*pi
+  else if(r(1)==0.d0.and.r(2).lt.0.d0)then
+   phi =-0.5d0*pi
+  else if(r(1)==0.d0.and.r(2)==0.d0)then
+   phi = 0.d0
+  endif
+  if(r(2).lt.0.d0.and.r(1).le.0.d0)then
+   tmp = pi - dabs(phi)
+   phi = pi + tmp
+  else if(r(2).lt.0.d0.and.r(1).gt.0.d0)then
+   phi = dtwo_pi + phi
+  endif
+ endif
+
+ if(.false.)then
+  x_2_y_2 = dsqrt(x_2_y_2)  
+  if(dabs(x_2_y_2).gt.1.d-20.and.dabs(r(2)).gt.1.d-20)then
+   phi = dabs(r(2))/r(2) * dacos(r(1)/x_2_y_2) 
+  else 
+   phi = 0.d0
+  endif
+ endif
+end
+
+
+subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
+ implicit none 
+ BEGIN_DOC
+ ! Y_lm(theta,phi) with -l<m<+l and 0<= l <=2
+ !
+ END_DOC
+ double precision, intent(in) :: theta, phi
+ integer, intent(in)          :: l,m
+ double precision, intent(out):: re_ylm,im_ylm
+ double precision :: tmp
+ include 'constants.include.F'
+ if(l==0.and.m==0)then
+  re_ylm = 0.5d0 * inv_sq_pi 
+  im_ylm = 0.d0
+ else if(l==1.and.m==1)then
+  tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta) 
+  re_ylm = tmp * dcos(phi)
+  im_ylm = tmp * dsin(phi)
+ else if(l==1.and.m==0)then
+  tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta) 
+  re_ylm = tmp 
+  im_ylm = 0.d0
+ else if(l==2.and.m==2)then
+  tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
+  re_ylm = tmp * dcos(2.d0*phi)
+  im_ylm = tmp * dsin(2.d0*phi)
+ else if(l==2.and.m==1)then
+  tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
+  re_ylm = tmp * dcos(phi)
+  im_ylm = tmp * dsin(phi)
+ else if(l==2.and.m==0)then
+  tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
+  re_ylm = tmp
+  im_ylm = 0.d0
+ endif
+end