diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt
index 9dea55e3..c23d9ad2 100644
--- a/slateratom/lib/CMakeLists.txt
+++ b/slateratom/lib/CMakeLists.txt
@@ -9,7 +9,6 @@ list(APPEND fypp_flags -I${projectdir}/common/include -DRELEASE="'${RELEASE}'")
 set(sources-f90
   avgpot.f90
   blas.f90
-  confinement.f90
   core_overlap.f90
   coulomb_hfex.f90
   coulomb_potential.f90
@@ -33,6 +32,7 @@ set(sources-f90
 set(sources-fpp
   blasroutines.F90
   broydenmixer.F90
+  confinement.F90
   lapackroutines.F90
   simplemixer.F90)
 
diff --git a/slateratom/lib/confinement.F90 b/slateratom/lib/confinement.F90
new file mode 100644
index 00000000..777dcfaf
--- /dev/null
+++ b/slateratom/lib/confinement.F90
@@ -0,0 +1,606 @@
+#:include 'common.fypp'
+
+!> Module that builds up various supervectors.
+module confinement
+
+  use common_accuracy, only : dp
+  use utilities, only : fak
+  use core_overlap, only : v
+  use density, only : basis_times_basis_times_r2
+
+  implicit none
+  private
+
+  public :: TConf, TConfInp, confType
+  public :: TPowerConf, TPowerConf_init
+  public :: TWsConf, TWsConf_init
+
+
+  !> Generic wrapper for confinement potentials.
+  type, abstract :: TConf
+  contains
+
+    procedure(getPotOnGrid), deferred :: getPotOnGrid
+    procedure(getSupervec), deferred :: getSupervec
+
+  end type TConf
+
+
+  abstract interface
+
+    !> Tabulates the (shell-resolved) confinement potential on a given grid.
+    subroutine getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf)
+      import :: TConf, dp
+      implicit none
+
+      !> instance
+      class(TConf), intent(in) :: this
+
+      !> maximum angular momentum
+      integer, intent(in) :: max_l
+
+      !> number of numerical integration points
+      integer, intent(in) :: num_mesh_points
+
+      !> numerical integration abcissas
+      real(dp), intent(in) :: abcissa(:)
+
+      !> confinement potential on grid
+      real(dp), intent(out) :: vconf(:, 0:)
+
+    end subroutine getPotOnGrid
+
+
+    !> Calculates supervector matrix elements of the confining potential.
+    subroutine getSupervec(this, max_l, num_mesh_points, abcissa, weight, num_alpha, alpha,&
+        & poly_order, vconf, vconf_matrix)
+      import :: TConf, dp
+      implicit none
+
+      !> instance
+      class(TConf), intent(in) :: this
+
+      !> maximum angular momentum
+      integer, intent(in) :: max_l
+
+      !> number of numerical integration points
+      integer, intent(in) :: num_mesh_points
+
+      !> numerical integration abcissas
+      real(dp), intent(in) :: abcissa(:)
+
+      !> numerical integration weights
+      real(dp), intent(in) :: weight(:)
+
+      !> number of exponents in each shell
+      integer, intent(in) :: num_alpha(0:)
+
+      !> basis exponents
+      real(dp), intent(in) :: alpha(0:,:)
+
+      !> highest polynomial order + l in each shell
+      integer, intent(in) :: poly_order(0:)
+
+      !> Woods-Saxon potential on grid
+      real(dp), intent(in) :: vconf(:, 0:)
+
+      !> Woods-Saxon confinement supervector
+      real(dp), intent(out) :: vconf_matrix(0:,:,:)
+
+    end subroutine getSupervec
+
+  end interface
+
+
+  !> Input for the Power confinement.
+  type :: TPowerConfInp
+
+    !> confinement radii (power compression)
+    real(dp) :: r0(0:4)
+
+    !> power of confinement (power compression)
+    real(dp) :: power(0:4)
+
+  end type TPowerConfInp
+
+
+  !> Input for the Woods-Saxon confinement.
+  type :: TWsConfInp
+
+    !> onset radii (Woods-Saxon compression)
+    real(dp) :: rOnset(0:4)
+
+    !> cutoff of confinement (Woods-Saxon compression)
+    real(dp) :: rCut(0:4)
+
+    !> potential well height of confinement (Woods-Saxon compression)
+    real(dp) :: vMax(0:4)
+
+  end type TWsConfInp
+
+
+  !> Input for the Power confinement.
+  type :: TConfInp
+
+    !> Power confinement input structure
+    type(TPowerConfInp) :: power
+
+    !> Woods-Saxon confinement input structure
+    type(TWsConfInp) :: ws
+
+  end type TConfInp
+
+
+  !> Power confinement.
+  type, extends(TConf) :: TPowerConf
+    private
+
+    !> confinement radii (power compression)
+    real(dp) :: r0(0:4)
+
+    !> power of confinement (power compression)
+    real(dp) :: power(0:4)
+
+  contains
+
+    procedure :: getPotOnGrid => TPowerConf_getPotOnGrid
+    procedure :: getSupervec => TPowerConf_getSupervec
+
+  end type TPowerConf
+
+
+  !> Woods-Saxon confinement.
+  type, extends(TConf) :: TWsConf
+    private
+
+    !> onset radii (Woods-Saxon compression)
+    real(dp) :: rOnset(0:4)
+
+    !> cutoff of confinement (Woods-Saxon compression)
+    real(dp) :: rCut(0:4)
+
+    !> potential well height of confinement (Woods-Saxon compression)
+    real(dp) :: vMax(0:4)
+
+  contains
+
+    procedure :: getPotOnGrid => TWsConf_getPotOnGrid
+    procedure :: getSupervec => TWsConf_getSupervec
+
+  end type TWsConf
+
+
+  !> Enumerator for type of confinement potential.
+  type :: TConfEnum
+
+    !> no compression
+    integer :: none = 0
+
+    !> power compression
+    integer :: power = 1
+
+    !> Woods-Saxon compression
+    integer :: ws = 2
+
+  end type TConfEnum
+
+  !> Container for enumerated types of confinement potentials.
+  type(TConfEnum), parameter :: confType = TConfEnum()
+
+
+contains
+
+  !> Initializes a TPowerConf instance.
+  subroutine TPowerConf_init(this, input)
+
+    !> instance
+    type(TPowerConf), intent(out) :: this
+
+    !> input data
+    type(TPowerConfInp), intent(inout) :: input
+
+    this%r0(0:) = input%r0
+    this%power(0:) = input%power
+
+  end subroutine TPowerConf_init
+
+
+  !> Initializes a TWsConf instance.
+  subroutine TWsConf_init(this, input)
+
+    !> instance
+    type(TWsConf), intent(out) :: this
+
+    !> input data
+    type(TWsConfInp), intent(inout) :: input
+
+    this%rOnset(0:) = input%rOnset
+    this%rCut(0:) = input%rCut
+    this%vMax(0:) = input%vMax
+
+  end subroutine TWsConf_init
+
+
+  !> Tabulates the (shell-resolved) Power confinement potential on the grid.
+  subroutine TPowerConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf)
+
+    !> instance
+    class(TPowerConf), intent(in) :: this
+
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of numerical integration points
+    integer, intent(in) :: num_mesh_points
+
+    !> numerical integration abcissas
+    real(dp), intent(in) :: abcissa(:)
+
+    !> Woods-Saxon potential on grid
+    real(dp), intent(out) :: vconf(:, 0:)
+
+    !! auxiliary variables
+    integer :: ii, ll
+
+    @:ASSERT(size(vconf, dim=1) == num_mesh_points)
+    @:ASSERT(size(vconf, dim=2) == max_l + 1)
+
+    do ll = 0, max_l
+      do ii = 1, num_mesh_points
+        vconf(ii, ll) = getPowerPotential(abcissa(ii), this%r0(ll), this%power(ll))
+      end do
+    end do
+
+  end subroutine TPowerConf_getPotOnGrid
+
+
+  !> Tabulates the (shell-resolved) Woods-Saxon confinement potential on the grid.
+  subroutine TWsConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf)
+
+    !> instance
+    class(TWsConf), intent(in) :: this
+
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of numerical integration points
+    integer, intent(in) :: num_mesh_points
+
+    !> numerical integration abcissas
+    real(dp), intent(in) :: abcissa(:)
+
+    !> Woods-Saxon potential on grid
+    real(dp), intent(out) :: vconf(:, 0:)
+
+    !! auxiliary variables
+    integer :: ii, ll
+
+    @:ASSERT(size(vconf, dim=1) == num_mesh_points)
+    @:ASSERT(size(vconf, dim=2) == max_l + 1)
+
+    do ll = 0, max_l
+      do ii = 1, num_mesh_points
+        vconf(ii, ll) = getWSPotential(abcissa(ii), this%rOnset(ll), this%rCut(ll), this%vMax(ll))
+      end do
+    end do
+
+  end subroutine TWsConf_getPotOnGrid
+
+
+  !> Tabulates the (shell-resolved) Power confinement potential on the grid.
+  subroutine TPowerConf_getSupervec(this, max_l, num_mesh_points, abcissa, weight, num_alpha,&
+      & alpha, poly_order, vconf, vconf_matrix)
+
+    !> instance
+    class(TPowerConf), intent(in) :: this
+
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of numerical integration points
+    integer, intent(in) :: num_mesh_points
+
+    !> numerical integration abcissas
+    real(dp), intent(in) :: abcissa(:)
+
+    !> numerical integration weights
+    real(dp), intent(in) :: weight(:)
+
+    !> number of exponents in each shell
+    integer, intent(in) :: num_alpha(0:)
+
+    !> basis exponents
+    real(dp), intent(in) :: alpha(0:,:)
+
+    !> highest polynomial order + l in each shell
+    integer, intent(in) :: poly_order(0:)
+
+    !> Power potential on grid
+    real(dp), intent(in) :: vconf(:, 0:)
+
+    !> Power confinement supervector
+    real(dp), intent(out) :: vconf_matrix(0:,:,:)
+
+    @:ASSERT(size(vconf, dim=1) == num_mesh_points)
+    @:ASSERT(size(vconf, dim=2) == max_l + 1)
+
+    @:ASSERT(size(vconf_matrix, dim=1) == max_l + 1)
+
+    call build_dft_conf_matrix(max_l, num_alpha, poly_order, alpha, num_mesh_points, abcissa,&
+        & weight, vconf, vconf_matrix)
+
+  end subroutine TPowerConf_getSupervec
+
+
+  !> Tabulates the (shell-resolved) Woods-Saxon confinement potential on the grid.
+  subroutine TWsConf_getSupervec(this, max_l, num_mesh_points, abcissa, weight, num_alpha, alpha,&
+      & poly_order, vconf, vconf_matrix)
+
+    !> instance
+    class(TWsConf), intent(in) :: this
+
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of numerical integration points
+    integer, intent(in) :: num_mesh_points
+
+    !> numerical integration abcissas
+    real(dp), intent(in) :: abcissa(:)
+
+    !> numerical integration weights
+    real(dp), intent(in) :: weight(:)
+
+    !> number of exponents in each shell
+    integer, intent(in) :: num_alpha(0:)
+
+    !> basis exponents
+    real(dp), intent(in) :: alpha(0:,:)
+
+    !> highest polynomial order + l in each shell
+    integer, intent(in) :: poly_order(0:)
+
+    !> Woods-Saxon potential on grid
+    real(dp), intent(in) :: vconf(:, 0:)
+
+    !> Woods-Saxon confinement supervector
+    real(dp), intent(out) :: vconf_matrix(0:,:,:)
+
+    @:ASSERT(size(vconf, dim=1) == num_mesh_points)
+    @:ASSERT(size(vconf, dim=2) == max_l + 1)
+
+    @:ASSERT(size(vconf_matrix, dim=1) == max_l + 1)
+
+    call build_dft_conf_matrix(max_l, num_alpha, poly_order, alpha, num_mesh_points, abcissa,&
+        & weight, vconf, vconf_matrix)
+
+  end subroutine TWsConf_getSupervec
+
+
+  !> Builds DFT confinement matrix to be added to the Fock matrix by calculating the single matrix
+  !! elements and putting them together.
+  subroutine build_dft_conf_matrix(max_l, num_alpha, poly_order, alpha, num_mesh_points, abcissa,&
+      & weight, vconf, conf_matrix)
+
+    !> maximum angular momentum
+    integer, intent(in) :: max_l
+
+    !> number of exponents in each shell
+    integer, intent(in) :: num_alpha(0:)
+
+    !> highest polynomial order + l in each shell
+    integer, intent(in) :: poly_order(0:)
+
+    !> basis exponents
+    real(dp), intent(in) :: alpha(0:,:)
+
+    !> number of numerical integration points
+    integer, intent(in) :: num_mesh_points
+
+    !> numerical integration abcissas
+    real(dp), intent(in) :: abcissa(:)
+
+    !> numerical integration weights
+    real(dp), intent(in) :: weight(:)
+
+    !> confinement potential on grid
+    real(dp), intent(in) :: vconf(:, 0:)
+
+    !> DFT confinement matrix
+    real(dp), intent(out) :: conf_matrix(0:,:,:)
+
+    !> single matrix element of the confinement potential
+    real(dp) :: conf_matrixelement
+
+    !> auxiliary variables
+    integer :: ii, jj, kk, ll, mm, ss, tt, start
+
+    conf_matrix(:,:,:) = 0.0_dp
+    conf_matrixelement = 0.0_dp
+
+    do ii = 0, max_l
+      ss = 0
+      do jj = 1, num_alpha(ii)
+        do kk = 1, poly_order(ii)
+          ss = ss + 1
+
+          tt = ss - 1
+          do ll = jj, num_alpha(ii)
+
+            start = 1
+            if (ll == jj) start = kk
+
+            do mm = start, poly_order(ii)
+              tt = tt + 1
+
+              call dft_conf_matrixelement(num_mesh_points, weight, abcissa, vconf(:, ii),&
+                  & alpha(ii, jj), kk, alpha(ii, ll), mm, ii, conf_matrixelement)
+
+              conf_matrix(ii, ss, tt) = conf_matrixelement
+              conf_matrix(ii, tt, ss) = conf_matrixelement
+
+            end do
+          end do
+        end do
+      end do
+    end do
+
+  end subroutine build_dft_conf_matrix
+
+
+  !> Calculates a single matrix element of the exchange correlation potential.
+  pure subroutine dft_conf_matrixelement(num_mesh_points, weight, abcissa, vconf, alpha1, poly1,&
+      & alpha2, poly2, ll, conf_matrixelement)
+
+    !> number of numerical integration points
+    integer, intent(in) :: num_mesh_points
+
+    !> numerical integration weights
+    real(dp), intent(in) :: weight(:)
+
+    !> numerical integration abcissas
+    real(dp), intent(in) :: abcissa(:)
+
+    !> confinement potential on grid
+    real(dp), intent(in) :: vconf(:)
+
+    !> basis exponent of 1st basis
+    real(dp), intent(in) :: alpha1
+
+    !> highest polynomial order in 1st basis shell
+    integer, intent(in) :: poly1
+
+    !> basis exponent of 2nd basis
+    real(dp), intent(in) :: alpha2
+
+    !> highest polynomial order in 2nd basis shell
+    integer, intent(in) :: poly2
+
+    !> angular momentum
+    integer, intent(in) :: ll
+
+    !> single matrix element of the exchange correlation potential
+    real(dp), intent(out) :: conf_matrixelement
+
+    !> stores product of two basis functions and r^2
+    real(dp) :: basis
+
+    !> auxiliary variable
+    integer :: ii
+
+    conf_matrixelement = 0.0_dp
+
+    do ii = 1, num_mesh_points
+      basis = basis_times_basis_times_r2(alpha1, poly1, alpha2, poly2, ll, abcissa(ii))
+      conf_matrixelement = conf_matrixelement + weight(ii) * vconf(ii) * basis
+    end do
+
+  end subroutine dft_conf_matrixelement
+
+
+  !> Returns the Power potential for a given radial distance.
+  pure function getPowerPotential(rr, r0, power) result(pot)
+
+    !> radial distance
+    real(dp), intent(in) :: rr
+
+    !> confinement radius
+    real(dp), intent(in) :: r0
+
+    !> confinement power
+    real(dp), intent(in) :: power
+
+    !> Power potential at evaluation point
+    real(dp) :: pot
+
+    pot = (rr / r0)**power
+
+  end function getPowerPotential
+
+
+  !> Returns the Woods-Saxon potential for a given radial distance.
+  pure function getWSPotential(rr, rOnset, rCut, vMax) result(pot)
+
+    !> radial distance
+    real(dp), intent(in) :: rr
+
+    !> onset radius
+    real(dp), intent(in) :: rOnset
+
+    !> cutoff radius
+    real(dp), intent(in) :: rCut
+
+    !> potential well height
+    real(dp), intent(in) :: vMax
+
+    !> Woods-Saxon potential at evaluation point
+    real(dp) :: pot
+
+    ! case: rr <= rOnset
+    pot = 0.0_dp
+
+    if (rr > rOnset .and. rr < rCut) then
+      pot = vMax / (rCut - rOnset) * (rr - rOnset)
+    elseif (rr >= rCut) then
+      pot = vMax
+    end if
+
+  end function getWSPotential
+
+
+  ! !> Calculates analytic matrix elements of confining potential.
+  ! !! No checking for power, e.g. power==0 or power<0 etc. !
+  ! subroutine getConfPower_analytical(this, max_l, num_alpha, alpha, poly_order, vconf_matrix)
+
+  !   !> instance
+  !   class(TPowerConf), intent(in) :: this
+
+  !   !> maximum angular momentum
+  !   integer, intent(in) :: max_l
+
+  !   !> number of exponents in each shell
+  !   integer, intent(in) :: num_alpha(0:)
+
+  !   !> basis exponents
+  !   real(dp), intent(in) :: alpha(0:,:)
+
+  !   !> highest polynomial order + l in each shell
+  !   integer, intent(in) :: poly_order(0:)
+
+  !   !> confinement supervector
+  !   real(dp), intent(out) :: vconf_matrix(0:,:,:)
+
+  !   !! temporary storage
+  !   real(dp) :: alpha1
+
+  !   !! auxiliary variables
+  !   integer :: ii, jj, kk, ll, mm, nn, oo, nlp, nlq
+
+  !   vconf_matrix(:,:,:) = 0.0_dp
+
+  !   do ii = 0, max_l
+  !     if (this%power(ii) > 1.0e-06_dp) then
+  !       nn = 0
+  !       do jj = 1, num_alpha(ii)
+  !         do ll = 1, poly_order(ii)
+  !           nn = nn + 1
+  !           oo = 0
+  !           nlp = ll + ii
+  !           do kk = 1, num_alpha(ii)
+  !             alpha1 = 0.5_dp * (alpha(ii, jj) + alpha(ii, kk))
+  !             do mm = 1, poly_order(ii)
+  !               oo = oo + 1
+  !               nlq = mm + ii
+  !               vconf_matrix(ii, nn, oo) = 1.0_dp / sqrt(v(alpha(ii, jj), 2 * nlp)&
+  !                   & * v(alpha(ii, kk), 2 * nlq)) / (this%r0(ii) * 2.0_dp)**this%power(ii)&
+  !                   & * v(alpha1, nlp + nlq + this%power(ii))
+  !             end do
+  !           end do
+  !         end do
+  !       end do
+  !     end if
+  !   end do
+
+  ! end subroutine TConf_getConfPower_analytical
+
+end module confinement
diff --git a/slateratom/lib/confinement.f90 b/slateratom/lib/confinement.f90
deleted file mode 100644
index 9e6c55ee..00000000
--- a/slateratom/lib/confinement.f90
+++ /dev/null
@@ -1,118 +0,0 @@
-!> Module that builds up various supervectors.
-module confinement
-
-  use common_accuracy, only : dp
-  use utilities, only : fak
-  use core_overlap, only : v
-
-  implicit none
-  private
-
-  public :: TConf, confType
-
-
-  !> Confinement potential structure.
-  type :: TConf
-
-    !> type of confinement potential (0: none, 1: power, 2: Woods-Saxon)
-    !! at the moment different shells are compressed by the same type of potential
-    integer :: type = 0
-
-    !> confinement radii (power compression)
-    real(dp) :: r0(0:4)
-
-    !> power of confinement (power compression)
-    real(dp) :: power(0:4)
-
-    !> onset radii (Woods-Saxon compression)
-    real(dp) :: onset(0:4)
-
-    !> cutoff of confinement (Woods-Saxon compression)
-    real(dp) :: cutoff(0:4)
-
-    !> potential well height of confinement (Woods-Saxon compression)
-    real(dp) :: vmax(0:4)
-
-  contains
-
-    procedure :: getConfPower => TConf_getConfPower
-
-  end type TConf
-
-
-  !> Enumerator for type of confinement potential.
-  type :: TConfEnum
-
-    !> no compression
-    integer :: none = 0
-
-    !> power compression
-    integer :: power = 1
-
-    !> Woods-Saxon compression
-    integer :: ws = 2
-
-  end type TConfEnum
-
-  !> Container for enumerated types of confinement potentials.
-  type(TConfEnum), parameter :: confType = TConfEnum()
-
-
-contains
-
-  !> Calculates analytic matrix elements of confining potential.
-  !! No checking for power, e.g. power==0 or power<0 etc. !
-  subroutine TConf_getConfPower(this, max_l, num_alpha, alpha, poly_order, vconf)
-
-    !> instance
-    class(TConf), intent(in) :: this
-
-    !> maximum angular momentum
-    integer, intent(in) :: max_l
-
-    !> number of exponents in each shell
-    integer, intent(in) :: num_alpha(0:)
-
-    !> basis exponents
-    real(dp), intent(in) :: alpha(0:,:)
-
-    !> highest polynomial order + l in each shell
-    integer, intent(in) :: poly_order(0:)
-
-    !> confinement supervector
-    real(dp), intent(out) :: vconf(0:,:,:)
-
-    !! temporary storage
-    real(dp) :: alpha1
-
-    !! auxiliary variables
-    integer :: ii, jj, kk, ll, mm, nn, oo, nlp, nlq
-
-    vconf(:,:,:) = 0.0_dp
-
-    do ii = 0, max_l
-      if (this%power(ii) > 1.0e-06_dp) then
-        nn = 0
-        do jj = 1, num_alpha(ii)
-          do ll = 1, poly_order(ii)
-            nn = nn + 1
-            oo = 0
-            nlp = ll + ii
-            do kk = 1, num_alpha(ii)
-              alpha1 = 0.5_dp * (alpha(ii, jj) + alpha(ii, kk))
-              do mm = 1, poly_order(ii)
-                oo = oo + 1
-                nlq = mm + ii
-                vconf(ii, nn, oo) = 1.0_dp / sqrt(v(alpha(ii, jj), 2 * nlp)&
-                    & * v(alpha(ii, kk), 2 * nlq)) / (this%r0(ii) * 2.0_dp)**this%power(ii)&
-                    & * v(alpha1, nlp + nlq + this%power(ii))
-              end do
-            end do
-          end do
-        end do
-      end if
-    end do
-
-  end subroutine TConf_getConfPower
-
-end module confinement
diff --git a/slateratom/lib/globals.f90 b/slateratom/lib/globals.f90
index 1874eba4..92e50a35 100644
--- a/slateratom/lib/globals.f90
+++ b/slateratom/lib/globals.f90
@@ -5,13 +5,16 @@ module globals
   use mixer, only : TMixer, TMixer_init, TMixer_reset, mixerTypes
   use broydenmixer, only : TBroydenMixer, TBroydenMixer_init
   use simplemixer, only : TSimpleMixer, TSimpleMixer_init
-  use confinement, only : TConf
+  use confinement, only : TConfInp
 
   implicit none
 
 
-  !> confinement potential
-  type(TConf) :: conf
+  !> type of confinement potential
+  integer :: conf_type = 0
+
+  !> confinement potential input
+  type(TConfInp) :: confInp
 
   !> basis exponents
   real(dp) :: alpha(0:4, 10)
@@ -64,8 +67,11 @@ module globals
   !> kinetic supervector
   real(dp), allocatable :: tt(:,:,:)
 
+  !> confinement potential on grid
+  real(dp), allocatable :: vconf(:,:)
+
   !> confinement supervector
-  real(dp), allocatable :: vconf(:,:,:)
+  real(dp), allocatable :: vconf_matrix(:,:,:)
 
   !> coulomb supermatrix
   real(dp), allocatable :: jj(:,:,:,:,:,:)
@@ -217,7 +223,8 @@ subroutine allocate_globals()
 
     allocate(uu(0:max_l, problemsize, problemsize))
     allocate(tt(0:max_l, problemsize, problemsize))
-    allocate(vconf(0:max_l, problemsize, problemsize))
+    allocate(vconf(num_mesh_points, 0:max_l))
+    allocate(vconf_matrix(0:max_l, problemsize, problemsize))
     allocate(ff(2, 0:max_l, problemsize, problemsize))
     allocate(pot_old(2, 0:max_l, problemsize, problemsize))
     allocate(pot_new(2, 0:max_l, problemsize, problemsize))
diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.f90
index 3feedae9..16e9b994 100644
--- a/slateratom/lib/input.f90
+++ b/slateratom/lib/input.f90
@@ -4,7 +4,7 @@ module input
   use common_accuracy, only : dp
   use common_poisson, only : TBeckeGridParams
 
-  use confinement, only : TConf, confType
+  use confinement, only : TConfInp, confType
   use xcfunctionals, only : xcFunctional
 
   implicit none
@@ -17,9 +17,9 @@ module input
 
   !> Reads in all properties, except for occupation numbers.
   subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min_alpha,&
-      & max_alpha, num_alpha, tAutoAlphas, alpha, conf, num_occ, num_power, num_alphas, xcnr,&
-      & tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta,&
-      & grid_params)
+      & max_alpha, num_alpha, tAutoAlphas, alpha, conf_type, confInp, num_occ, num_power,&
+      & num_alphas, xcnr, tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega,&
+      & camAlpha, camBeta, grid_params)
 
     !> nuclear charge, i.e. atomic number
     integer, intent(out) :: nuc
@@ -54,8 +54,11 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     !> basis exponents
     real(dp), intent(out) :: alpha(0:4, 10)
 
-    !> confinement potential
-    type(TConf), intent(out) :: conf
+    !> type of confinement potential
+    integer, intent(out) :: conf_type
+
+    !> confinement potential input
+    type(TConfInp), intent(out) :: confInp
 
     !> maximal occupied shell
     integer, intent(out) :: num_occ
@@ -157,30 +160,30 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     last_conf_type = 0
     do ii = 0, max_l
       write(*, '(A,I3)') 'Enter confinement ID (0: none, 1: power, 2: WS) for l=', ii
-      read(*,*) conf%type
+      read(*,*) conf_type
       if (ii > 0) then
-        if (conf%type /= last_conf_type) then
+        if (conf_type /= last_conf_type) then
           error stop 'At the moment different shells are supposed to be compressed with the same&
               & type of potential'
         end if
       end if
-      last_conf_type = conf%type
+      last_conf_type = conf_type
     end do
 
-    select case (conf%type)
+    select case (conf_type)
     case(confType%none)
       continue
     case(confType%power)
       write(*, '(A)') 'Enter parameters r_0 and power'
       do ii = 0, max_l
         write(*, '(A,I3)') 'l=', ii
-        read(*,*) conf%r0(ii), conf%power(ii)
+        read(*,*) confInp%power%r0(ii), confInp%power%power(ii)
       end do
     case(confType%ws)
       write(*, '(A)') 'Enter parameters Ronset, Rcut and Vmax'
       do ii = 0, max_l
         write(*, '(A,I3)') 'l=', ii
-        read(*,*) conf%onset(ii), conf%cutoff(ii), conf%vmax(ii)
+        read(*,*) confInp%ws%rOnset(ii), confInp%ws%rCut(ii), confInp%ws%vMax(ii)
       end do
     case default
       error stop 'Invalid confinement potential.'
@@ -293,7 +296,8 @@ end subroutine read_input_2
 
   !> Echos gathered input to stdout.
   subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha,&
-      & conf, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const)
+      & conf_type, confInp, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points,&
+      & xalpha_const)
 
     !> nuclear charge, i.e. atomic number
     integer, intent(in) :: nuc
@@ -319,8 +323,11 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a
     !> basis exponents
     real(dp), intent(in) :: alpha(0:,:)
 
-    !> confinement potential
-    type(TConf), intent(in) :: conf
+    !> type of confinement potential
+    integer, intent(in) :: conf_type
+
+    !> confinement potential input
+    type(TConfInp), intent(in) :: confInp
 
     !> occupation numbers
     real(dp), intent(in) :: occ(:,0:,:)
@@ -415,15 +422,15 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a
 
     write(*, '(A)') ' '
     do ii = 0, max_l
-      select case (conf%type)
+      select case (conf_type)
       case(confType%none)
         write(*, '(A,I3,A)') 'l= ', ii, ' no confinement'
       case(confType%power)
-        write(*, '(A,I3,A,E15.7,A,E15.7)') 'l= ', ii, ', r0= ', conf%r0(ii),&
-            & ' power= ', conf%power(ii)
+        write(*, '(A,I3,A,E15.7,A,E15.7)') 'l= ', ii, ', r0= ', confInp%power%r0(ii),&
+            & ' power= ', confInp%power%power(ii)
       case(confType%ws)
-        write(*, '(A,I3,A,E15.7,A,E15.7,A,E15.7)') 'l= ', ii, ', Ronset= ', conf%onset(ii),&
-            & ' Rcutoff= ', conf%cutoff(ii), ' Vmax= ', conf%vmax(ii)
+        write(*, '(A,I3,A,E15.7,A,E15.7,A,E15.7)') 'l= ', ii, ', Ronset= ', confInp%ws%rOnset(ii),&
+            & ' Rcutoff= ', confInp%ws%rCut(ii), ' Vmax= ', confInp%ws%vMax(ii)
       end select
     end do
 
diff --git a/slateratom/prog/CMakeLists.txt b/slateratom/prog/CMakeLists.txt
index 5693fe3d..945ddd4a 100644
--- a/slateratom/prog/CMakeLists.txt
+++ b/slateratom/prog/CMakeLists.txt
@@ -1,9 +1,21 @@
+set(projectdir ${PROJECT_SOURCE_DIR})
+
+#
+# General options for all targets
+#
+set(fypp_flags ${FYPP_BUILD_FLAGS} ${FYPP_CONFIG_FLAGS})
+list(APPEND fypp_flags -I${projectdir}/common/include -DRELEASE="'${RELEASE}'")
+
 set(sources-f90
-  cmdargs.f90
-  main.f90)
+  cmdargs.f90)
+
+set(sources-fpp
+  main.F90)
 
-add_executable(slateratom ${sources-f90})
+skprogs_preprocess("${FYPP}" "${fypp_flags}" "F90" "f90" "${sources-fpp}" sources-f90-preproc)
+
+add_executable(slateratom ${sources-f90} ${sources-f90-preproc})
 
 target_link_libraries(slateratom skprogs-slateratom)
-  
+
 install(TARGETS slateratom EXPORT skprogs-targets DESTINATION ${CMAKE_INSTALL_BINDIR})
diff --git a/slateratom/prog/main.f90 b/slateratom/prog/main.F90
similarity index 78%
rename from slateratom/prog/main.f90
rename to slateratom/prog/main.F90
index 3bf52d23..9ac24388 100644
--- a/slateratom/prog/main.f90
+++ b/slateratom/prog/main.F90
@@ -1,3 +1,5 @@
+#:include 'common.fypp'
+
 program HFAtom
 
   use common_accuracy, only : dp
@@ -5,7 +7,7 @@ program HFAtom
   use integration, only : gauss_chebyshev_becke_mesh
   use input, only : read_input_1, read_input_2, echo_input
   use core_overlap, only : overlap, nuclear, kinetic
-  use confinement, only : confType
+  use confinement, only : TConf, confType, TPowerConf, TPowerConf_init, TWsConf, TWsConf_init
   use coulomb_hfex, only : coulomb, hfex, hfex_lr
   use densitymatrix, only : densmatrix
   use hamiltonian, only : build_hamiltonian
@@ -50,13 +52,17 @@ program HFAtom
   !! holds parameters, defining a Becke integration grid
   type(TBeckeGridParams) :: grid_params
 
+  !! general confinement potential
+  class(TConf), allocatable :: conf
+
   ! deactivate average potential calculation for now
   isAvgPotNeeded = .false.
 
   call parse_command_arguments()
   call read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min_alpha, max_alpha,&
-      & num_alpha, tAutoAlphas, alpha, conf, num_occ, num_power, num_alphas, xcnr, tPrintEigvecs,&
-      & tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta, grid_params)
+      & num_alpha, tAutoAlphas, alpha, conf_type, confInp, num_occ, num_power, num_alphas, xcnr,&
+      & tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta,&
+      & grid_params)
 
   problemsize = num_power * num_alphas
 
@@ -73,8 +79,8 @@ program HFAtom
   if (nuc > 36) num_mesh_points = 1250
   if (nuc > 54) num_mesh_points = 1500
 
-  call echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha, conf, occ,&
-      & num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const)
+  call echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha, conf_type,&
+      & confInp, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const)
 
   ! allocate global stuff and zero out
   call allocate_globals()
@@ -91,16 +97,21 @@ program HFAtom
   call nuclear(uu, max_l, num_alpha, alpha, poly_order)
   call kinetic(tt, max_l, num_alpha, alpha, poly_order)
 
-  select case (conf%type)
-  case(confType%none)
-    vconf(:,:,:) = 0.0_dp
+  select case (conf_type)
   case(confType%power)
-    call conf%getConfPower(max_l, num_alpha, alpha, poly_order, vconf)
+    @:CREATE_CLASS(conf, TPowerConf, TPowerConf_init, confInp%power)
   case(confType%ws)
-    error stop 'Woods-Saxon compression not yet supported.'
-    ! call conf%getConfWS()
+    @:CREATE_CLASS(conf, TWsConf, TWsConf_init, confInp%ws)
   end select
 
+  if (allocated(conf)) then
+    call conf%getPotOnGrid(max_l, num_mesh_points, abcissa, vconf)
+    call conf%getSupervec(max_l, num_mesh_points, abcissa, weight, num_alpha, alpha, poly_order,&
+        & vconf, vconf_matrix)
+  else
+    vconf_matrix(0:,:,:) = 0.0_dp
+  end if
+
   ! test for linear dependency
   call diagonalize_overlap(max_l, num_alpha, poly_order, ss)
 
@@ -134,7 +145,7 @@ program HFAtom
   pot_old(:,:,:,:) = 0.0_dp
 
   ! kinetic energy, nuclear-electron, and confinement matrix elements which are constant during SCF
-  call build_hamiltonian(pMixer, 0, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha,&
+  call build_hamiltonian(pMixer, 0, tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l, num_alpha,&
       & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old,&
       & pot_new, tZora, ff, camAlpha, camBeta)
 
@@ -158,19 +169,20 @@ program HFAtom
         & dz, xcnr, omega, camAlpha, camBeta, rho, drho, ddrho, vxc, exc, xalpha_const)
 
     ! build Fock matrix and get total energy during SCF
-    call build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha,&
-        & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old,&
-        & pot_new, tZora, ff, camAlpha, camBeta)
+    call build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l,&
+        & num_alpha, poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha,&
+        & pot_old, pot_new, tZora, ff, camAlpha, camBeta)
 
     if (tZora) then
-      call getTotalEnergyZora(tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha, poly_order,&
-          & problemsize, xcnr, num_mesh_points, weight, abcissa, rho, exc, vxc, eigval_scaled, occ,&
-          & camAlpha, camBeta, kinetic_energy, nuclear_energy, coulomb_energy, exchange_energy,&
-          & x_en_2, conf_energy, total_ene)
+      call getTotalEnergyZora(tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l, num_alpha,&
+          & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, rho, exc, vxc,&
+          & eigval_scaled, occ, camAlpha, camBeta, kinetic_energy, nuclear_energy, coulomb_energy,&
+          & exchange_energy, x_en_2, conf_energy, total_ene)
     else
-      call getTotalEnergy(tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha, poly_order,&
-          & xcnr, num_mesh_points, weight, abcissa, rho, exc, camAlpha, camBeta, kinetic_energy,&
-          & nuclear_energy, coulomb_energy, exchange_energy, x_en_2, conf_energy, total_ene)
+      call getTotalEnergy(tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l, num_alpha,&
+          & poly_order, xcnr, num_mesh_points, weight, abcissa, rho, exc, camAlpha, camBeta,&
+          & kinetic_energy, nuclear_energy, coulomb_energy, exchange_energy, x_en_2, conf_energy,&
+          & total_ene)
     end if
 
     call check_convergence(pot_old, pot_new, max_l, problemsize, scftol, iScf, change_max,&
@@ -215,10 +227,10 @@ program HFAtom
   end if
 
   if (tZora) then
-    call getTotalEnergyZora(tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha, poly_order,&
-        & problemsize, xcnr, num_mesh_points, weight, abcissa, rho, exc, vxc, eigval_scaled, occ,&
-        & camAlpha, camBeta, zora_ekin, nuclear_energy, coulomb_energy, exchange_energy, x_en_2,&
-        & conf_energy, total_ene)
+    call getTotalEnergyZora(tt, uu, nuc, vconf_matrix, jj, kk, kk_lr, pp, max_l, num_alpha,&
+        & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, rho, exc, vxc,&
+        & eigval_scaled, occ, camAlpha, camBeta, zora_ekin, nuclear_energy, coulomb_energy,&
+        & exchange_energy, x_en_2, conf_energy, total_ene)
   end if
 
   write(*, '(A,E20.12)') 'Potential Matrix Elements converged to ', change_max