diff --git a/slateratom/lib/core_overlap.f90 b/slateratom/lib/core_overlap.f90
index 22923ed6..ae484a03 100644
--- a/slateratom/lib/core_overlap.f90
+++ b/slateratom/lib/core_overlap.f90
@@ -167,7 +167,7 @@ end subroutine nuclear
   !            oo=oo+1
   !            nlq=mm+ii
   !
-  !            normalization=float(2**(nlp+nlq+1))/&
+  !            normalization=real(2**(nlp+nlq+1),dp)/&
   !                sqrt(v(alpha(ii,jj),2*nlp)*v(alpha(ii,kk),2*nlq))
   !
   !            part1=exp_int(alpha2,nlp+nlq-1,r0)-exp_int(alpha2,nlp+nlq-1,0.0d0)
@@ -387,7 +387,7 @@ function w(x,i,j) ! W_{ij}(x)
     integer, intent(in) :: i,j
     real(dp) :: w
 
-    w=2.0d0*float((j-i-1))/x
+    w=2.0d0*real((j-i-1),dp)/x
 
     return
   end function w
diff --git a/slateratom/lib/coulomb_hfex.f90 b/slateratom/lib/coulomb_hfex.f90
index 04631eec..36d64272 100644
--- a/slateratom/lib/coulomb_hfex.f90
+++ b/slateratom/lib/coulomb_hfex.f90
@@ -85,7 +85,7 @@ subroutine coulomb(j,max_l,num_alpha,alpha,poly_order,u,s)
     !  do ii=0,max_l
     !    do jj=0,max_l
     !      j(ii,:,:,jj,:,:)=j(ii,:,:,jj,:,:)/&
-    !         &((2.0d0*float(ii)+1.0d0)*(2.0d0*float(jj)+1.0d0))
+    !         &((2.0d0*real(ii,dp)+1.0d0)*(2.0d0*real(jj,dp)+1.0d0))
     !    end do
     !  end do
 
@@ -222,7 +222,7 @@ subroutine hfex(k,max_l,num_alpha,alpha,poly_order,problemsize)
     !  do ii=0,max_l
     !    do jj=0,max_l
     !      k(ii,:,:,jj,:,:)=k(ii,:,:,jj,:,:)/&
-    !         &((2.0d0*float(ii)+1.0d0)*(2.0d0*float(jj)+1.0d0))
+    !         &((2.0d0*real(ii,dp)+1.0d0)*(2.0d0*real(jj,dp)+1.0d0))
     !    end do
     !  end do
 
@@ -304,7 +304,7 @@ function almn(lambda,mu,nu)
     real(dp) :: almn
 
     almn=a(lambda+mu-nu)*a(lambda-mu+nu)*a(mu-lambda+nu)/&
-        &(float(lambda+mu+nu+1)*a(lambda+mu+nu))
+        &(real(lambda+mu+nu+1,dp)*a(lambda+mu+nu))
 
   end function almn
 
diff --git a/slateratom/lib/coulomb_potential.f90 b/slateratom/lib/coulomb_potential.f90
index 13aa5d73..5f183884 100644
--- a/slateratom/lib/coulomb_potential.f90
+++ b/slateratom/lib/coulomb_potential.f90
@@ -71,9 +71,9 @@ subroutine cou_pot(p,max_l,num_alpha,poly_order,alpha,problemsize,&
 
               ! add normalization of basis functions
               ! watch out for 2**(nlp+nlq+1) needed because variable integration ranges
-              help1(:,ii,ll,oo)=help1(:,ii,ll,oo)*float(2**(nlp+nlq+1))/&
+              help1(:,ii,ll,oo)=help1(:,ii,ll,oo)*real(2**(nlp+nlq+1),dp)/&
                   &sqrt(v(alpha(ii,jj),2*nlp)*v(alpha(ii,mm),2*nlq))
-              help2(:,ii,ll,oo)=help2(:,ii,ll,oo)*float(2**(nlp+nlq+1))/&
+              help2(:,ii,ll,oo)=help2(:,ii,ll,oo)*real(2**(nlp+nlq+1),dp)/&
                   &sqrt(v(alpha(ii,jj),2*nlp)*v(alpha(ii,mm),2*nlq))
 
             end do
diff --git a/slateratom/lib/density.f90 b/slateratom/lib/density.f90
index 5a662524..89ef944c 100644
--- a/slateratom/lib/density.f90
+++ b/slateratom/lib/density.f90
@@ -315,10 +315,10 @@ function basis_1st(alpha,poly_order,l,r)
     if ((r==0.0d0).and.((poly_order+l-1)==0)) then
       basis_1st=normalization*(-alpha*exp(-alpha*r))
     else if ((r==0.0d0).and.((poly_order+l-2)==0)) then
-      basis_1st=normalization*(float(poly_order+l-1)*&
+      basis_1st=normalization*(real(poly_order+l-1,dp)*&
           &exp(-alpha*r)-alpha*r**(poly_order+l-1)*exp(-alpha*r))
     else
-      basis_1st=normalization*(float(poly_order+l-1)*r**(poly_order+l-2)*&
+      basis_1st=normalization*(real(poly_order+l-1,dp)*r**(poly_order+l-2)*&
           &exp(-alpha*r)-alpha*r**(poly_order+l-1)*exp(-alpha*r))
     end if
 
@@ -339,16 +339,16 @@ function basis_2nd(alpha,poly_order,l,r)
 
     ! catch 0^0
     if ((r==0.0d0).and.((poly_order+l-3)==0)) then
-      basis_2nd=normalization*(float(poly_order+l-1)*float(poly_order+l-2)*&
+      basis_2nd=normalization*(real(poly_order+l-1,dp)*real(poly_order+l-2,dp)*&
           &exp(-alpha*r))
     else if ((r==0.0d0).and.((poly_order+l-2)==0)) then
-      basis_2nd=normalization*(-2.0d0*alpha*float(poly_order+l-1)*&
+      basis_2nd=normalization*(-2.0d0*alpha*real(poly_order+l-1,dp)*&
           &exp(-alpha*r))
     else if ((r==0.0d0).and.((poly_order+l-1)==0)) then
       basis_2nd=normalization*(alpha**2*exp(-alpha*r))
     else
-      basis_2nd=normalization*(float(poly_order+l-1)*float(poly_order+l-2)*&
-          &r**(poly_order+l-3)*exp(-alpha*r)-2.0d0*alpha*float(poly_order+l-1)*&
+      basis_2nd=normalization*(real(poly_order+l-1,dp)*real(poly_order+l-2,dp)*&
+          &r**(poly_order+l-3)*exp(-alpha*r)-2.0d0*alpha*real(poly_order+l-1,dp)*&
           &r**(poly_order+l-2)*exp(-alpha*r)+alpha**2*r**(poly_order+l-1)*&
           &exp(-alpha*r))
     end if
@@ -419,10 +419,10 @@ function basis_times_basis_1st(alpha,poly1,beta,poly2,l,r)
           &(-beta)*exp(ab*r)
     else if ((r==0.0d0).and.((m+n-3)==0)) then
       basis_times_basis_1st=normalization1*normalization2*&
-          &(float(n-1))*exp(ab*r)
+          &(real(n-1,dp))*exp(ab*r)
     else
       basis_times_basis_1st=normalization1*normalization2*&
-          &(float(n-1)*r**(m+n-3)-beta*r**(n+m-2))*exp(ab*r)
+          &(real(n-1,dp)*r**(m+n-3)-beta*r**(n+m-2))*exp(ab*r)
     end if
 
     if (abs(basis_times_basis_1st)<1.0d-20) basis_times_basis_1st=0.0d0
@@ -452,8 +452,8 @@ function basis_times_basis_2nd(alpha,poly1,beta,poly2,l,r)
 
     ! WARNING: without summing negative and positive contributions independently
     ! zora becomes completely unstable !
-    positive=float((n-1)*(n-2))*r**(m+n-4)+beta**2*r**(m+n-2)
-    negative=float(2*(n-1))*beta*r**(n+m-3)
+    positive=real((n-1)*(n-2),dp)*r**(m+n-4)+beta**2*r**(m+n-2)
+    negative=real(2*(n-1),dp)*beta*r**(n+m-3)
 
     basis_times_basis_2nd=normalization1*normalization2*&
         &(positive-negative)*exp(ab*r)
@@ -488,16 +488,16 @@ function basis_1st_times_basis_1st(alpha,poly1,beta,poly2,l,r)
     if ((r==0.0d0).and.((m+n-2)==0)) then
       positive=alpha*beta
     else if ((r==0.0d0).and.((m+n-4)==0)) then
-      positive=float((m-1)*(n-1))
+      positive=real((m-1)*(n-1),dp)
     else
-      positive=float((m-1)*(n-1))*r**(m+n-4)+&
+      positive=real((m-1)*(n-1),dp)*r**(m+n-4)+&
           &alpha*beta*r**(m+n-2)
     end if
 
     if ((r==0.0d0).and.((m+n-3)==0)) then
-      negative=(alpha*float(n-1)+beta*float(m-1))
+      negative=(alpha*real(n-1,dp)+beta*real(m-1,dp))
     else
-      negative=(alpha*float(n-1)+beta*float(m-1))*r**(m+n-3)
+      negative=(alpha*real(n-1,dp)+beta*real(m-1,dp))*r**(m+n-3)
     end if
 
     basis_1st_times_basis_1st=normalization1*normalization2*&
@@ -528,15 +528,15 @@ function basis_2nd_times_basis_2nd(alpha,poly1,beta,poly2,l,r)
 
     ! WARNING: without summing negative and positive contributions independently
     ! zora becomes completely unstable !
-    positive=float((m-1)*(m-2)*(n-1)*(n-2))*r**(n+m-6)+&
-        &r**(m+n-4)*(beta**2*float((m-1)*(m-2))+alpha**2*float((n-1)*(n-2))+&
-        &alpha*beta*float(4*(m-1)*(n-1)))+&
+    positive=real((m-1)*(m-2)*(n-1)*(n-2),dp)*r**(n+m-6)+&
+        &r**(m+n-4)*(beta**2*real((m-1)*(m-2),dp)+alpha**2*real((n-1)*(n-2),dp)+&
+        &alpha*beta*real(4*(m-1)*(n-1),dp))+&
         &alpha**2*beta**2*r**(m+n-2)
 
-    negative=r**(m+n-5)*(beta*float(2*(n-1)*(m-1)*(m-2))+&
-        &alpha*float(2*(m-1)*(n-1)*(n-2)))+&
-        &r**(m+n-3)*(alpha*beta**2*float(2*(m-1))+&
-        &beta*alpha**2*float(2*(n-1)))
+    negative=r**(m+n-5)*(beta*real(2*(n-1)*(m-1)*(m-2),dp)+&
+        &alpha*real(2*(m-1)*(n-1)*(n-2),dp))+&
+        &r**(m+n-3)*(alpha*beta**2*real(2*(m-1),dp)+&
+        &beta*alpha**2*real(2*(n-1),dp))
 
     basis_2nd_times_basis_2nd=normalization1*normalization2*&
         &(positive-negative)*exp(ab*r)
@@ -596,7 +596,7 @@ function basis_times_basis_1st_times_r2(alpha,poly1,beta,poly2,l,r)
     ! WARNING: without summing negative and positive contributions independently
     ! zora becomes completely unstable !
     basis_times_basis_1st_times_r2=normalization1*normalization2*&
-        &(float(n-1)*r**(m+n-1)-beta*r**(n+m))*exp(ab*r)
+        &(real(n-1,dp)*r**(m+n-1)-beta*r**(n+m))*exp(ab*r)
 
     if (abs(basis_times_basis_1st_times_r2)<1.0d-20) &
         &basis_times_basis_1st_times_r2=0.0d0
@@ -626,8 +626,8 @@ function basis_times_basis_2nd_times_r2(alpha,poly1,beta,poly2,l,r)
 
     ! WARNING: without summing negative and positive contributions independently
     ! zora becomes completely unstable !
-    positive=float((n-1)*(n-2))*r**(m+n-2)+beta**2*r**(m+n)
-    negative=float(2*(n-1))*beta*r**(n+m-1)
+    positive=real((n-1)*(n-2),dp)*r**(m+n-2)+beta**2*r**(m+n)
+    negative=real(2*(n-1),dp)*beta*r**(n+m-1)
 
     basis_times_basis_2nd_times_r2=normalization1*normalization2*&
         &(positive-negative)*exp(ab*r)
@@ -661,7 +661,7 @@ function basis_times_basis_1st_times_r(alpha,poly1,beta,poly2,l,r)
     ! WARNING: without summing negative and positive contributions independently
     ! zora becomes completely unstable !
     basis_times_basis_1st_times_r=normalization1*normalization2*&
-        &(float(n-1)*r**(m+n-2)-beta*r**(n+m-1))*exp(ab*r)
+        &(real(n-1,dp)*r**(m+n-2)-beta*r**(n+m-1))*exp(ab*r)
 
     if (abs(basis_times_basis_1st_times_r)<1.0d-20) &
         &basis_times_basis_1st_times_r=0.0d0
@@ -689,9 +689,9 @@ function basis_1st_times_basis_1st_times_r2(alpha,poly1,beta,poly2,l,r)
 
     ! WARNING: without summing negative and positive contributions independently
     ! zora becomes completely unstable !
-    positive=float((m-1)*(n-1))*r**(m+n-2)+&
+    positive=real((m-1)*(n-1),dp)*r**(m+n-2)+&
         &alpha*beta*r**(m+n)
-    negative=(alpha*float(n-1)+beta*float(m-1))*r**(m+n-1)
+    negative=(alpha*real(n-1,dp)+beta*real(m-1,dp))*r**(m+n-1)
 
     basis_1st_times_basis_1st_times_r2=normalization1*normalization2*&
         &(positive-negative)*exp(ab*r)
diff --git a/slateratom/lib/dft.f90 b/slateratom/lib/dft.f90
index daeb9840..8719d8a6 100644
--- a/slateratom/lib/dft.f90
+++ b/slateratom/lib/dft.f90
@@ -33,14 +33,14 @@ subroutine dft_start_pot(abcissa,num_mesh_points,nuc,vxc)
     real(dp) :: b,t,x,rtx
     integer :: ii
 
-    b= (0.69395656d0/float(nuc))**(1.0d0/3.0d0)
+    b= (0.69395656d0/real(nuc,dp))**(1.0d0/3.0d0)
 
     do ii=1,num_mesh_points
 
       x= abcissa(ii)/b
       rtx= sqrt(x)
 
-      t= float(nuc)/(1.0d0+rtx*(0.02747d0-x*(0.1486d0-0.007298d0*x))&
+      t= real(nuc,dp)/(1.0d0+rtx*(0.02747d0-x*(0.1486d0-0.007298d0*x))&
           &+x*(1.243d0+x*(0.2302d0+0.006944d0*x)));
       if (t < 1.0d0) t= 1.0d0
 
diff --git a/slateratom/lib/hamiltonian.f90 b/slateratom/lib/hamiltonian.f90
index e20fcb27..865921e3 100644
--- a/slateratom/lib/hamiltonian.f90
+++ b/slateratom/lib/hamiltonian.f90
@@ -65,8 +65,8 @@ subroutine build_fock(iter,t,u,nuc,vconf,j,k,p,max_l,num_alpha,poly_order,&
 
     ! build mixer input 
 
-      pot_new(1,:,:,:)=-float(nuc)*u(:,:,:)+j_matrix(:,:,:)-k_matrix(1,:,:,:)
-      pot_new(2,:,:,:)=-float(nuc)*u(:,:,:)+j_matrix(:,:,:)-k_matrix(2,:,:,:)
+      pot_new(1,:,:,:)=-real(nuc,dp)*u(:,:,:)+j_matrix(:,:,:)-k_matrix(1,:,:,:)
+      pot_new(2,:,:,:)=-real(nuc,dp)*u(:,:,:)+j_matrix(:,:,:)-k_matrix(2,:,:,:)
 
 
     ! mixer
diff --git a/slateratom/lib/integration.f90 b/slateratom/lib/integration.f90
index 75647b12..45a51e5b 100644
--- a/slateratom/lib/integration.f90
+++ b/slateratom/lib/integration.f90
@@ -34,7 +34,7 @@ subroutine gauss_chebyshev_becke_mesh(N,nuc,w,r, dzdr, d2zdr2, dz)
     allocate(x(N)) 
     allocate(fak(N)) 
     !
-    temp=pi/float(N+1)
+    temp=pi/real(N+1,dp)
     dz = temp
     !    
     do ii=1,N
@@ -52,7 +52,7 @@ subroutine gauss_chebyshev_becke_mesh(N,nuc,w,r, dzdr, d2zdr2, dz)
           &/ (4.0_dp * bragg(nuc)**2 * (-1.0_dp + cosz) * sinz)
 
       ! r**2 times first derivative of x -> r mapping function
-      w(ii)=temp*(sin(float(ii)*temp))
+      w(ii)=temp*(sin(real(ii,dp)*temp))
       !      fak(ii)=2.0_dp*r(ii)**2*bragg(nuc)/(1.0_dp-x(ii))**2
       fak(ii)=2.0_dp*bragg(nuc)/(1.0_dp-x(ii))**2
 
@@ -79,12 +79,12 @@ subroutine get_abcissas(N,nuc,r,step)
 
     allocate(x(N)) 
 
-    step=pi/float(N+1)
+    step=pi/real(N+1,dp)
 
     do ii=1,N
 
       ! NOTE prefactor 
-      x(ii)=(-1.0_dp)*cos(step*float(ii)) ! gauss-chebyshev abcissas
+      x(ii)=(-1.0_dp)*cos(step*real(ii,dp)) ! gauss-chebyshev abcissas
       r(ii)=(1.0_dp+x(ii))/(1.0_dp-x(ii))*bragg(nuc)
 
     end do
@@ -103,12 +103,12 @@ subroutine get_abcissas_z_1st(N,nuc,dr,step)
     integer, intent(out) :: step ! generator step size
     integer :: ii
 
-    step=pi/float(N+1)
+    step=pi/real(N+1,dp)
 
     do ii=1,N
 
-      dr(ii)=2.0d0*bragg(nuc)*pi*sin(step*float(ii))/&
-          &(1.0d0+2.0d0*cos(step*float(ii))+cos(step*float(ii))**2)
+      dr(ii)=2.0d0*bragg(nuc)*pi*sin(step*real(ii,dp))/&
+          &(1.0d0+2.0d0*cos(step*real(ii,dp))+cos(step*real(ii,dp))**2)
 
     end do
 
@@ -124,12 +124,12 @@ subroutine get_abcissas_z_2nd(N,nuc,ddr,step)
     integer, intent(out) :: step ! generator step size
     integer :: ii
 
-    step=pi/float(N+1)
+    step=pi/real(N+1,dp)
 
     do ii=1,N
 
-      ddr(ii)=(-2.0d0*bragg(nuc)*pi**2)*(cos(step*float(ii))-2.0d0)/&
-          &(1.0d0+2.0d0*cos(step*float(ii))+cos(step*float(ii))**2)
+      ddr(ii)=(-2.0d0*bragg(nuc)*pi**2)*(cos(step*real(ii,dp))-2.0d0)/&
+          &(1.0d0+2.0d0*cos(step*real(ii,dp))+cos(step*real(ii,dp))**2)
 
     end do
 
@@ -240,7 +240,7 @@ function exp_int(alpha,power,r)
     exp_int=1.0d0/alpha*exp(alpha*r)
 
     do ii=1,power
-      exp_int=1.0d0/alpha*r**ii*exp(alpha*r)-float(ii)/alpha*exp_int
+      exp_int=1.0d0/alpha*r**ii*exp(alpha*r)-real(ii,dp)/alpha*exp_int
     end do
 
   end function exp_int
diff --git a/slateratom/lib/output.f90 b/slateratom/lib/output.f90
index 78ced5b6..f0cb7a2d 100644
--- a/slateratom/lib/output.f90
+++ b/slateratom/lib/output.f90
@@ -276,13 +276,13 @@ subroutine write_potentials_file_standard(num_mesh_points,abcissa,weight,&
 
     do ii=1,num_mesh_points
       write(95,'(6ES21.12E3)') abcissa(ii), weight(ii), &
-          &float(-nuc) / abcissa(ii), cpot(ii), vxc(ii,1), vxc(ii,2)
+          &real(-nuc,dp) / abcissa(ii), cpot(ii), vxc(ii,1), vxc(ii,2)
     end do
     close(95)
 
     do ii=1,num_mesh_points
       ecou=ecou+weight(ii)*rhotot(ii)*cpot(ii)*abcissa(ii)**2
-      enuc=enuc-weight(ii)*rhotot(ii)*float(nuc)*abcissa(ii)
+      enuc=enuc-weight(ii)*rhotot(ii)*real(nuc,dp)*abcissa(ii)
       vxcint(1)=vxcint(1)+weight(ii)*rho(ii,1)*vxc(ii,1)*abcissa(ii)**2
       vxcint(2)=vxcint(2)+weight(ii)*rho(ii,2)*vxc(ii,2)*abcissa(ii)**2
     end do
diff --git a/slateratom/lib/total_energy.f90 b/slateratom/lib/total_energy.f90
index 32b161f6..3481f97c 100644
--- a/slateratom/lib/total_energy.f90
+++ b/slateratom/lib/total_energy.f90
@@ -229,7 +229,7 @@ subroutine core_hamiltonian_energies(t,u,vconf,p_total,max_l,num_alpha,&
             do mm=1,poly_order(ii)
               tt=tt+1
               kinetic=kinetic+t(ii,ss,tt)*p_total(ii,ss,tt)
-              nuclear=nuclear-float(nuc)*u(ii,ss,tt)*p_total(ii,ss,tt)
+              nuclear=nuclear-real(nuc,dp)*u(ii,ss,tt)*p_total(ii,ss,tt)
               confinement=confinement+vconf(ii,ss,tt)*p_total(ii,ss,tt)
             end do
           end do
diff --git a/slateratom/lib/zora_routines.f90 b/slateratom/lib/zora_routines.f90
index a614c38a..3033bfbf 100644
--- a/slateratom/lib/zora_routines.f90
+++ b/slateratom/lib/zora_routines.f90
@@ -62,13 +62,13 @@ subroutine zora_t_correction(mode,t,max_l,num_alpha,alpha,poly_order,&
                     &kappa(1,:),alpha(ii,jj),ll,alpha(ii,kk),&
                     &mm,ii)+kinetic_part_2(num_mesh_points,weight,abcissa,&
                     &kappa(1,:),alpha(ii,jj),ll,alpha(ii,kk),&
-                    &mm,ii)*dfloat(ii*(ii+1))
+                    &mm,ii)*real(ii*(ii+1),dp)
 
                 t(2,ii,nn,oo)=kinetic_part_1(num_mesh_points,weight,abcissa,&
                     &kappa(2,:),alpha(ii,jj),ll,alpha(ii,kk),&
                     &mm,ii)+kinetic_part_2(num_mesh_points,weight,abcissa,&
                     &kappa(2,:),alpha(ii,jj),ll,alpha(ii,kk),&
-                    &mm,ii)*dfloat(ii*(ii+1))
+                    &mm,ii)*real(ii*(ii+1),dp)
 
               end if
 
@@ -82,13 +82,13 @@ subroutine zora_t_correction(mode,t,max_l,num_alpha,alpha,poly_order,&
                     &kappa2(1,:),alpha(ii,jj),ll,alpha(ii,kk),&
                     &mm,ii)+kinetic_part_2(num_mesh_points,weight,abcissa,&
                     &kappa2(1,:),alpha(ii,jj),ll,alpha(ii,kk),&
-                    &mm,ii)*dfloat(ii*(ii+1))
+                    &mm,ii)*real(ii*(ii+1),dp)
 
                 t(2,ii,nn,oo)=kinetic_part_1(num_mesh_points,weight,abcissa,&
                     &kappa2(2,:),alpha(ii,jj),ll,alpha(ii,kk),&
                     &mm,ii)+kinetic_part_2(num_mesh_points,weight,abcissa,&
                     &kappa2(2,:),alpha(ii,jj),ll,alpha(ii,kk),&
-                    &mm,ii)*dfloat(ii*(ii+1))
+                    &mm,ii)*real(ii*(ii+1),dp)
 
               end if
 
@@ -325,8 +325,8 @@ subroutine potential_to_mesh(num_mesh_points,abcissa,&
 
     do ii=1,num_mesh_points
 
-      vtot(1,ii)=-float(nuc)/abcissa(ii)+cpot(ii)+vxc(ii,1)
-      vtot(2,ii)=-float(nuc)/abcissa(ii)+cpot(ii)+vxc(ii,2)
+      vtot(1,ii)=-real(nuc,dp)/abcissa(ii)+cpot(ii)+vxc(ii,1)
+      vtot(2,ii)=-real(nuc,dp)/abcissa(ii)+cpot(ii)+vxc(ii,2)
 
     end do