diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90
index c262105f..32402c74 100644
--- a/src/Integrals_Monoelec/pseudopot.f90
+++ b/src/Integrals_Monoelec/pseudopot.f90
@@ -1585,148 +1585,8 @@ end
        bessel_mod_exp=x**n*accu
        end
 
-!        double precision function bessel_mod(x,n)
-!        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-!        parameter(NBESSMAX=100)
-!        dimension SI(0:NBESSMAX),DI(0:NBESSMAX)
-!        if(n.lt.0.or.n.gt.NBESSMAX)stop 'pb with argument of bessel_mod'
-!        CALL SPHI(N,X,NBESSMAX,SI,DI)
-!        bessel_mod=si(n)
-!        end
-
-        SUBROUTINE SPHI(N,X,NMAX,SI,DI)
-!
-!       ========================================================
-!       Purpose: Compute modified spherical Bessel functions
-!                of the first kind, in(x) and in'(x)
-!       Input :  x --- Argument of in(x)
-!                n --- Order of in(x) ( n = 0,1,2,... )
-!       Output:  SI(n) --- in(x)
-!                DI(n) --- in'(x)
-!                NM --- Highest order computed
-!       Routines called:
-!                MSTA1 and MSTA2 for computing the starting
-!                point for backward recurrence
-!       ========================================================
-!
-        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        DIMENSION SI(0:NMAX),DI(0:NMAX)
-        NM=N
-        IF (DABS(X).LT.1.0D-100) THEN
-           DO 10 K=0,N
-              SI(K)=0.0D0
-10            DI(K)=0.0D0
-           SI(0)=1.0D0
-           DI(1)=0.333333333333333D0
-           RETURN
-        ENDIF
-        SI(0)=DSINH(X)/X
-        SI(1)=-(DSINH(X)/X-DCOSH(X))/X
-        SI0=SI(0)
-        IF (N.GE.2) THEN
-           M=MSTA1(X,200)
-
-           write(34,*)'m=',m
-
-           IF (M.LT.N) THEN
-              NM=M
-           ELSE
-              M=MSTA2(X,N,15)
-           write(34,*)'m=',m
-           ENDIF
-           F0=0.0D0
-           F1=1.0D0-100
-           DO 15 K=M,0,-1
-              F=(2.0D0*K+3.0D0)*F1/X+F0
-              IF (K.LE.NM) SI(K)=F
-              F0=F1
-15            F1=F
-           CS=SI0/F
-           write(34,*)'cs=',cs
-           DO 20 K=0,NM
-20            SI(K)=CS*SI(K)
-        ENDIF
-        DI(0)=SI(1)
-        DO 25 K=1,NM
-25         DI(K)=SI(K-1)-(K+1.0D0)/X*SI(K)
-        RETURN
-        END
 
 
-        INTEGER FUNCTION MSTA1(X,MP)
-!
-!       ===================================================
-!       Purpose: Determine the starting point for backward  
-!                recurrence such that the magnitude of    
-!                Jn(x) at that point is about 10^(-MP)
-!       Input :  x     --- Argument of Jn(x)
-!                MP    --- Value of magnitude
-!       Output:  MSTA1 --- Starting point   
-!       ===================================================
-!
-        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        A0=DABS(X)
-        N0=INT(1.1*A0)+1
-        F0=ENVJ(N0,A0)-MP
-        N1=N0+5
-        F1=ENVJ(N1,A0)-MP
-        DO 10 IT=1,20             
-           NN=N1-(N1-N0)/(1.0D0-F0/F1)                  
-           F=ENVJ(NN,A0)-MP
-           IF(ABS(NN-N1).LT.1) GO TO 20
-           N0=N1
-           F0=F1
-           N1=NN
- 10        F1=F
- 20     MSTA1=NN
-        RETURN
-        END
-
-
-        INTEGER FUNCTION MSTA2(X,N,MP)
-!
-!       ===================================================
-!       Purpose: Determine the starting point for backward
-!                recurrence such that all Jn(x) has MP
-!                significant digits
-!       Input :  x  --- Argument of Jn(x)
-!                n  --- Order of Jn(x)
-!                MP --- Significant digit
-!       Output:  MSTA2 --- Starting point
-!       ===================================================
-!
-        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        A0=DABS(X)
-        HMP=0.5D0*MP
-        EJN=ENVJ(N,A0)
-        IF (EJN.LE.HMP) THEN
-           OBJ=MP
-           N0=INT(1.1*A0)
-        ELSE
-           OBJ=HMP+EJN
-           N0=N
-        ENDIF
-        F0=ENVJ(N0,A0)-OBJ
-        N1=N0+5
-        F1=ENVJ(N1,A0)-OBJ
-        DO 10 IT=1,20
-           NN=N1-(N1-N0)/(1.0D0-F0/F1)
-           F=ENVJ(NN,A0)-OBJ
-           IF (iABS(NN-N1).LT.1) GO TO 20
-           N0=N1
-           F0=F1
-           N1=NN
-10         F1=F
-20      MSTA2=NN+10
-        RETURN
-        END
-
-        double precision FUNCTION ENVJ(N,X)
-        DOUBLE PRECISION X
-        integer N
-        ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
-        RETURN
-        END
 
 !c  Computation of real spherical harmonics Ylm(theta,phi)
 !c