From 5f459c10b9bcf3d3be66cd2e7a7223c89c713b8a Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Thu, 4 Jul 2024 15:13:33 +0200
Subject: [PATCH 1/8] Hand over WS compression to slateratom

---
 sktools/src/sktools/calculators/slateratom.py |  75 ++++++++---
 sktools/src/sktools/common.py                 |   2 +
 sktools/src/sktools/compressions.py           |  73 +++++++++++
 slateratom/lib/CMakeLists.txt                 |   1 +
 slateratom/lib/confinement.f90                | 118 ++++++++++++++++++
 slateratom/lib/core_overlap.f90               |  61 +--------
 slateratom/lib/globals.f90                    |   7 +-
 slateratom/lib/input.f90                      |  78 ++++++++----
 slateratom/prog/main.f90                      |  23 ++--
 9 files changed, 324 insertions(+), 114 deletions(-)
 create mode 100644 slateratom/lib/confinement.f90

diff --git a/sktools/src/sktools/calculators/slateratom.py b/sktools/src/sktools/calculators/slateratom.py
index 5659ddb6..f429d284 100644
--- a/sktools/src/sktools/calculators/slateratom.py
+++ b/sktools/src/sktools/calculators/slateratom.py
@@ -10,8 +10,9 @@
 import sktools.hsd.converter as conv
 import sktools.common as sc
 from sktools.taggedfile import TaggedFile
-import sktools.compressions
+import sktools.compressions as skcomp
 import sktools.radial_grid as oc
+import sktools.xcfunctionals as xc
 
 
 LOGGER = logging.getLogger('slateratom')
@@ -134,7 +135,7 @@ class SlateratomInput:
     functional : str
         DFT functional type ('lda', 'pbe', 'blyp', 'lcpbe', 'lcbnl')
     compressions : list
-        List of PowerCompression objects. Either empty (no compression applied)
+        List of Compression objects. Either empty (no compression applied)
         or has a compression object for every angular momentum of the atom.
     settings : SlaterAtom
         Further detailed settings of the program.
@@ -193,16 +194,26 @@ def __init__(self, settings, atomconfig, functional, compressions):
         if compressions is None:
             compressions = []
         for comp in compressions:
-            if not isinstance(comp, sktools.compressions.PowerCompression):
-                msg = "Invalid compressiont type {} for slateratom".format(
+            if not isinstance(comp, (skcomp.PowerCompression,
+                                     skcomp.WoodsSaxonCompression)):
+                msg = "Invalid compression type {} for slateratom".format(
                     comp.__class__.__name__)
                 raise sc.SkgenException(msg)
         maxang = atomconfig.maxang
         ncompr = len(compressions)
-        if ncompr and ncompr != maxang + 1:
+        if ncompr > 0 and ncompr != maxang + 1:
             msg = "Invalid number of compressions" \
                 "(expected {:d}, got {:d})".format(maxang + 1, ncompr)
             raise sc.SkgenException(msg)
+
+        # block different compression types of different shells for now
+        if ncompr > 0:
+            compids = [comp.compid for comp in compressions]
+            if not len(set(compids)) == 1:
+                msg = "At the moment, shells may only be compressed by the " \
+                    + "same type of potential."
+                raise sc.SkgenException(msg)
+
         self._compressions = compressions
         myrelativistics = sc.RELATIVISTICS_NONE, sc.RELATIVISTICS_ZORA
         if atomconfig.relativistics not in myrelativistics:
@@ -213,7 +224,7 @@ def __init__(self, settings, atomconfig, functional, compressions):
     def isXCFunctionalSupported(self, functional):
         '''Checks if the given xc-functional is supported by the calculator,
            in particular: checks if AVAILABLE_FUNCTIONALS intersect with
-           sktools.xcfunctionals.XCFUNCTIONALS
+           xc.XCFUNCTIONALS
 
         Args:
             functional: xc-functional, defined in xcfunctionals.py
@@ -224,8 +235,8 @@ def isXCFunctionalSupported(self, functional):
 
         tmp = []
         for xx in SUPPORTED_FUNCTIONALS:
-            if xx in sktools.xcfunctionals.XCFUNCTIONALS:
-                tmp.append(sktools.xcfunctionals.XCFUNCTIONALS[xx])
+            if xx in xc.XCFUNCTIONALS:
+                tmp.append(xc.XCFUNCTIONALS[xx])
 
         return bool(functional.__class__ in tmp)
 
@@ -297,14 +308,42 @@ def write(self, workdir):
 
         # Compressions
         if len(self._compressions) == 0:
-            out += ["{:g} {:d} \t\t{:s} Compr. radius and power ({:s})".format(
-                1e30, 0, self._COMMENT, sc.ANGMOM_TO_SHELL[ll])
+            # no compression for all shells
+            out += ["{:d} \t\t\t{:s} Compression ID".format(
+                skcomp.SUPPORTED_COMPRESSIONS['nocompression'],
+                self._COMMENT) + " ({:s})".format(sc.ANGMOM_TO_SHELL[ll])
                     for ll in range(maxang + 1)]
         else:
-            out += ["{:g} {:g} \t\t{:s} Compr. radius and power ({:s})".format(
-                compr.radius, compr.power, self._COMMENT,
-                sc.ANGMOM_TO_SHELL[ll])
-                    for ll, compr in enumerate(self._compressions)]
+            # define the type of compression for each shell
+            for ll, compr in enumerate(self._compressions):
+                if compr.compid == \
+                   skcomp.SUPPORTED_COMPRESSIONS['powercompression']:
+                    out += ["{:d} \t\t{:s} Compression ID".format(
+                        skcomp.SUPPORTED_COMPRESSIONS['powercompression'],
+                        self._COMMENT) \
+                            + " ({:s})".format(sc.ANGMOM_TO_SHELL[ll])]
+                elif compr.compid == \
+                     skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression']:
+                    out += ["{:d} \t\t{:s} Compression ID".format(
+                        skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression'],
+                        self._COMMENT) \
+                            + " ({:s})".format(sc.ANGMOM_TO_SHELL[ll])]
+                else:
+                    msg = 'Invalid compression type.'
+                    raise sc.SkgenException(msg)
+            # provide the compression parametrization for each shell
+            for ll, compr in enumerate(self._compressions):
+                if compr.compid == \
+                   skcomp.SUPPORTED_COMPRESSIONS['powercompression']:
+                    out += ["{:g} {:g} \t\t{:s} Compr. radius and power ({:s})"
+                            .format(compr.radius, compr.power, self._COMMENT,
+                                    sc.ANGMOM_TO_SHELL[ll])]
+                elif compr.compid == \
+                     skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression']:
+                    out += ["{:g} {:g} {:g} \t\t{:s}".format(
+                        compr.onset, compr.cutoff, compr.vmax, self._COMMENT) \
+                            + " Compr. onset, cutoff and vmax ({:s})"
+                            .format(sc.ANGMOM_TO_SHELL[ll])]
 
         out += ["{:d} \t\t\t{:s} nr. of occupied shells ({:s})".format(
             len(occ), self._COMMENT, sc.ANGMOM_TO_SHELL[ll])
@@ -313,9 +352,9 @@ def write(self, workdir):
         # STO powers and exponents
         exponents = self._settings.exponents
         maxpowers = self._settings.maxpowers
-        out += ["{:d} {:d} \t\t\t{:s} nr. of exponents, max. power ({:s})".format(
-            len(exponents[ll]), maxpowers[ll], self._COMMENT,
-            sc.ANGMOM_TO_SHELL[ll])
+        out += ["{:d} {:d} \t\t\t{:s} nr. of exponents, max. power ({:s})"
+                .format(len(exponents[ll]), maxpowers[ll], self._COMMENT,
+                        sc.ANGMOM_TO_SHELL[ll])
                 for ll in range(maxang + 1)]
         out.append("{:s} \t\t{:s} automatic exponent generation".format(
             self._LOGICALSTRS[False], self._COMMENT))
@@ -326,7 +365,7 @@ def write(self, workdir):
 
         out.append("{:s} \t\t{:s} write eigenvectors".format(
             self._LOGICALSTRS[False], self._COMMENT))
-        out.append("{} {:g} \t\t{:s} broyden mixer, mixing factor".format(
+        out.append("{} {:g} \t\t\t{:s} broyden mixer, mixing factor".format(
             2, 0.1, self._COMMENT))
 
         # Occupations
diff --git a/sktools/src/sktools/common.py b/sktools/src/sktools/common.py
index 7c97d78d..8e3aa313 100644
--- a/sktools/src/sktools/common.py
+++ b/sktools/src/sktools/common.py
@@ -528,6 +528,8 @@ def is_shelf_file_matching(shelf_file, mydict):
         db = shelve.open(shelf_file, 'r')
     except dbm.error:
         return False
+    if not type(db) is type(mydict):
+        return False
     match = True
     for key in mydict:
         match = key in db and db[key] == mydict[key]
diff --git a/sktools/src/sktools/compressions.py b/sktools/src/sktools/compressions.py
index 62a1eb05..ec54c48b 100644
--- a/sktools/src/sktools/compressions.py
+++ b/sktools/src/sktools/compressions.py
@@ -39,6 +39,8 @@ def fromhsd(cls, root, query):
                 node=child)
 
         myself = cls()
+        myself.compid = SUPPORTED_COMPRESSIONS[
+            myself.__class__.__name__.lower()]
         myself.power = power
         myself.radius = radius
 
@@ -62,9 +64,80 @@ def __eq__(self, other):
         return power_ok and radius_ok
 
 
+class WoodsSaxonCompression(sc.ClassDict):
+    '''Compression by the Woods Saxon potential.
+
+    Attributes
+    ----------
+    onset : float
+        Onset radius of the compression
+    cutoff : float
+        Cutoff radius of the compression
+    vmax : float
+        Potential well depth/height
+    '''
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        '''Creates instance from a HSD-node and with given query object.'''
+
+        onset, child = query.getvalue(root, 'onset', conv.float0,
+                                      returnchild=True)
+        if onset < 0.0:
+            msg = 'Invalid onset radius {:f}'.format(onset)
+            raise hsd.HSDInvalidTagValueException(msg=msg, node=child)
+
+        cutoff, child = query.getvalue(root, 'cutoff', conv.float0,
+                                       returnchild=True)
+        if cutoff <= onset:
+            msg = 'Invalid cutoff radius {:f}'.format(cutoff)
+            raise hsd.HSDInvalidTagValueException(msg=msg, node=child)
+
+        vmax, child = query.getvalue(
+            root, 'vmax', conv.float0, defvalue=100.0, returnchild=True)
+        if vmax <= 0.0:
+            msg = 'Invalid potential well height {:f}'.format(vmax)
+            raise hsd.HSDInvalidTagValueException(msg=msg, node=child)
+
+        myself = cls()
+        myself.compid = SUPPORTED_COMPRESSIONS[
+            myself.__class__.__name__.lower()]
+        myself.onset = onset
+        myself.cutoff = cutoff
+        myself.vmax = vmax
+
+        return myself
+
+
+    def tohsd(self, root, query, parentname=None):
+        ''''''
+
+        if parentname is None:
+            mynode = root
+        else:
+            mynode = query.setchild(root, 'WoodsSaxonCompression')
+
+        query.setchildvalue(mynode, 'onset', conv.float0, self.onset)
+        query.setchildvalue(mynode, 'cutoff', conv.float0, self.cutoff)
+        query.setchildvalue(mynode, 'vmax', conv.float0, self.vmax)
+
+
+    def __eq__(self, other):
+        onset_ok = abs(self.onset - other.onset) < 1e-8
+        cutoff_ok = abs(self.cutoff - other.cutoff) < 1e-8
+        vmax_ok = abs(self.vmax - other.vmax) < 1e-8
+        return onset_ok and cutoff_ok and vmax_ok
+
+
 # Registered compressions with corresponding hsd name as key
 COMPRESSIONS = {
     'powercompression': PowerCompression,
+    'woodssaxoncompression': WoodsSaxonCompression,
+}
+SUPPORTED_COMPRESSIONS = {
+    'nocompression': 0,
+    'powercompression': 1,
+    'woodssaxoncompression': 2,
 }
 
 
diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt
index 78399778..9dea55e3 100644
--- a/slateratom/lib/CMakeLists.txt
+++ b/slateratom/lib/CMakeLists.txt
@@ -9,6 +9,7 @@ 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
diff --git a/slateratom/lib/confinement.f90 b/slateratom/lib/confinement.f90
new file mode 100644
index 00000000..9e6c55ee
--- /dev/null
+++ b/slateratom/lib/confinement.f90
@@ -0,0 +1,118 @@
+!> 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/core_overlap.f90 b/slateratom/lib/core_overlap.f90
index 258386c6..91e9d523 100644
--- a/slateratom/lib/core_overlap.f90
+++ b/slateratom/lib/core_overlap.f90
@@ -7,7 +7,7 @@ module core_overlap
   implicit none
   private
 
-  public :: overlap, kinetic, nuclear, moments, v, confinement
+  public :: overlap, kinetic, nuclear, moments, v
 
 
   interface v
@@ -192,65 +192,6 @@ pure subroutine kinetic(tt, max_l, num_alpha, alpha, poly_order)
   end subroutine kinetic
 
 
-  !> Calculates analytic matrix elements of confining potential.
-  !! No checking for power, e.g. power==0 or power<0 etc. !
-  pure subroutine confinement(vconf, max_l, num_alpha, alpha, poly_order, conf_r0, conf_power)
-
-    !> confinement supervector
-    real(dp), intent(out) :: vconf(0:,:,:)
-
-    !> 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 radii
-    real(dp), intent(in) :: conf_r0(0:)
-
-    !> power of confinement
-    real(dp), intent(in) :: conf_power(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 (conf_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)) / (conf_r0(ii) * 2.0_dp)**conf_power(ii)&
-                    & * v(alpha1, nlp + nlq + conf_power(ii))
-              end do
-            end do
-          end do
-        end do
-      end if
-    end do
-
-  end subroutine confinement
-
-
   !> Calculates arbitrary moments of electron distribution, e.g. expectation values of <r>, <r^2>
   !! etc.; this is implemented analytically for arbitrary powers.
   pure subroutine moments(moment, max_l, num_alpha, alpha, poly_order, cof, power)
diff --git a/slateratom/lib/globals.f90 b/slateratom/lib/globals.f90
index 6687874d..1874eba4 100644
--- a/slateratom/lib/globals.f90
+++ b/slateratom/lib/globals.f90
@@ -5,14 +5,13 @@ 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
 
   implicit none
 
-  !> confinement radii
-  real(dp) :: conf_r0(0:4)
 
-  !> power of confinement
-  real(dp) :: conf_power(0:4)
+  !> confinement potential
+  type(TConf) :: conf
 
   !> basis exponents
   real(dp) :: alpha(0:4, 10)
diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.f90
index 51104403..3feedae9 100644
--- a/slateratom/lib/input.f90
+++ b/slateratom/lib/input.f90
@@ -4,6 +4,7 @@ module input
   use common_accuracy, only : dp
   use common_poisson, only : TBeckeGridParams
 
+  use confinement, only : TConf, confType
   use xcfunctionals, only : xcFunctional
 
   implicit none
@@ -16,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_r0, conf_power, 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, 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
@@ -53,11 +54,8 @@ 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 radii
-    real(dp), intent(out) :: conf_r0(0:4)
-
-    !> power of confinement
-    real(dp), intent(out) :: conf_power(0:4)
+    !> confinement potential
+    type(TConf), intent(out) :: conf
 
     !> maximal occupied shell
     integer, intent(out) :: num_occ
@@ -98,6 +96,9 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     !> holds parameters, defining a Becke integration grid
     type(TBeckeGridParams), intent(out) :: grid_params
 
+    !! confinement type of last query
+    integer :: last_conf_type
+
     !! auxiliary variables
     integer :: ii, jj
 
@@ -153,12 +154,38 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
       stop
     end if
 
-    write(*, '(A)') 'Enter Confinement: r_0 and integer power, power=0.0 -> off'
+    last_conf_type = 0
     do ii = 0, max_l
-      write(*, '(A,I3)') 'l=', ii
-      read(*,*) conf_r0(ii), conf_power(ii)
+      write(*, '(A,I3)') 'Enter confinement ID (0: none, 1: power, 2: WS) for l=', ii
+      read(*,*) conf%type
+      if (ii > 0) 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
     end do
 
+    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)
+      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)
+      end do
+    case default
+      error stop 'Invalid confinement potential.'
+    end select
+
     write(*, '(A)') 'Enter number of occupied shells for each angular momentum.'
     do ii = 0, max_l
       write(*, '(A,I3)') 'l=', ii
@@ -188,7 +215,8 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
       end do
     else
       do ii = 0, max_l
-        write(*, '(A,I3,A,I3,A)') 'Enter ', num_alpha(ii), 'exponents for l=', ii, ', one per line.'
+        write(*, '(A,I3,A,I3,A)') 'Enter ', num_alpha(ii),&
+            & ' exponents for l=', ii, ', one per line.'
         do jj = 1, num_alpha(ii)
           read(*,*) alpha(ii, jj)
         end do
@@ -265,8 +293,7 @@ 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_r0, conf_power, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points,&
-      & xalpha_const)
+      & conf, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const)
 
     !> nuclear charge, i.e. atomic number
     integer, intent(in) :: nuc
@@ -292,11 +319,8 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a
     !> basis exponents
     real(dp), intent(in) :: alpha(0:,:)
 
-    !> confinement radii
-    real(dp), intent(in) :: conf_r0(0:)
-
-    !> power of confinement
-    real(dp), intent(in) :: conf_power(0:)
+    !> confinement potential
+    type(TConf), intent(in) :: conf
 
     !> occupation numbers
     real(dp), intent(in) :: occ(:,0:,:)
@@ -391,12 +415,16 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a
 
     write(*, '(A)') ' '
     do ii = 0, max_l
-      if (conf_power(ii) > 1.0e-06_dp) then
-        write(*, '(A,I3,A,E15.7,A,E15.7)') 'l= ', ii, ', r0= ', conf_r0(ii), ' power= ',&
-            & conf_power(ii)
-      else
-        write(*, '(A,I3,A)') 'l= ', ii, ' no confinement '
-      end if
+      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)
+      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)
+      end select
     end do
 
     write(*, '(A)') ' '
diff --git a/slateratom/prog/main.f90 b/slateratom/prog/main.f90
index a4143149..3bf52d23 100644
--- a/slateratom/prog/main.f90
+++ b/slateratom/prog/main.f90
@@ -4,7 +4,8 @@ program HFAtom
   use common_message, only : error
   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, confinement
+  use core_overlap, only : overlap, nuclear, kinetic
+  use confinement, only : confType
   use coulomb_hfex, only : coulomb, hfex, hfex_lr
   use densitymatrix, only : densmatrix
   use hamiltonian, only : build_hamiltonian
@@ -54,9 +55,8 @@ program HFAtom
 
   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_r0, conf_power, num_occ, num_power, num_alphas, xcnr,&
-      & tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta,&
-      & grid_params)
+      & num_alpha, tAutoAlphas, alpha, conf, 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 +73,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_r0,&
-      & conf_power, 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, occ,&
+      & num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const)
 
   ! allocate global stuff and zero out
   call allocate_globals()
@@ -90,7 +90,16 @@ program HFAtom
   call overlap(ss, max_l, num_alpha, alpha, poly_order)
   call nuclear(uu, max_l, num_alpha, alpha, poly_order)
   call kinetic(tt, max_l, num_alpha, alpha, poly_order)
-  call confinement(vconf, max_l, num_alpha, alpha, poly_order, conf_r0, conf_power)
+
+  select case (conf%type)
+  case(confType%none)
+    vconf(:,:,:) = 0.0_dp
+  case(confType%power)
+    call conf%getConfPower(max_l, num_alpha, alpha, poly_order, vconf)
+  case(confType%ws)
+    error stop 'Woods-Saxon compression not yet supported.'
+    ! call conf%getConfWS()
+  end select
 
   ! test for linear dependency
   call diagonalize_overlap(max_l, num_alpha, poly_order, ss)

From 424c38fa2d9d5da7c46db9f381236174622a52b1 Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Tue, 16 Jul 2024 10:55:41 +0200
Subject: [PATCH 2/8] Implement Woods-Saxon confinement matrix elements

---
 slateratom/lib/CMakeLists.txt          |   2 +-
 slateratom/lib/confinement.F90         | 606 +++++++++++++++++++++++++
 slateratom/lib/confinement.f90         | 118 -----
 slateratom/lib/globals.f90             |  17 +-
 slateratom/lib/input.f90               |  47 +-
 slateratom/prog/CMakeLists.txt         |  20 +-
 slateratom/prog/{main.f90 => main.F90} |  64 +--
 7 files changed, 700 insertions(+), 174 deletions(-)
 create mode 100644 slateratom/lib/confinement.F90
 delete mode 100644 slateratom/lib/confinement.f90
 rename slateratom/prog/{main.f90 => main.F90} (78%)

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

From f10a9b4aa3dde01117a029e8e7ab29b68d56a822 Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Tue, 16 Jul 2024 17:47:43 +0200
Subject: [PATCH 3/8] Rework WS potential in terms of W, a and r0

---
 sktools/src/sktools/calculators/slateratom.py |  4 +-
 sktools/src/sktools/compressions.py           | 54 +++++++------
 slateratom/lib/CMakeLists.txt                 |  2 +-
 slateratom/lib/confinement.F90                | 78 +++++++++----------
 slateratom/lib/{input.f90 => input.F90}       | 13 +++-
 5 files changed, 74 insertions(+), 77 deletions(-)
 rename slateratom/lib/{input.f90 => input.F90} (97%)

diff --git a/sktools/src/sktools/calculators/slateratom.py b/sktools/src/sktools/calculators/slateratom.py
index f429d284..4fd25de3 100644
--- a/sktools/src/sktools/calculators/slateratom.py
+++ b/sktools/src/sktools/calculators/slateratom.py
@@ -341,8 +341,8 @@ def write(self, workdir):
                 elif compr.compid == \
                      skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression']:
                     out += ["{:g} {:g} {:g} \t\t{:s}".format(
-                        compr.onset, compr.cutoff, compr.vmax, self._COMMENT) \
-                            + " Compr. onset, cutoff and vmax ({:s})"
+                        compr.ww, compr.aa, compr.r0, self._COMMENT) \
+                            + " Compr. height, slope, half-height radius ({:s})"
                             .format(sc.ANGMOM_TO_SHELL[ll])]
 
         out += ["{:d} \t\t\t{:s} nr. of occupied shells ({:s})".format(
diff --git a/sktools/src/sktools/compressions.py b/sktools/src/sktools/compressions.py
index ec54c48b..c2c69b4d 100644
--- a/sktools/src/sktools/compressions.py
+++ b/sktools/src/sktools/compressions.py
@@ -69,42 +69,40 @@ class WoodsSaxonCompression(sc.ClassDict):
 
     Attributes
     ----------
-    onset : float
-        Onset radius of the compression
-    cutoff : float
-        Cutoff radius of the compression
-    vmax : float
-        Potential well depth/height
+    ww : float
+        Height (W) of the compression
+    aa : float
+        Slope (a) of the compression
+    r0 : float
+        Half-height radius of the compression
     '''
 
     @classmethod
     def fromhsd(cls, root, query):
         '''Creates instance from a HSD-node and with given query object.'''
 
-        onset, child = query.getvalue(root, 'onset', conv.float0,
-                                      returnchild=True)
-        if onset < 0.0:
-            msg = 'Invalid onset radius {:f}'.format(onset)
+        ww, child = query.getvalue(
+            root, 'w', conv.float0, defvalue=100.0, returnchild=True)
+        if ww < 0.0:
+            msg = 'Invalid potential height (W) {:f}'.format(ww)
             raise hsd.HSDInvalidTagValueException(msg=msg, node=child)
 
-        cutoff, child = query.getvalue(root, 'cutoff', conv.float0,
-                                       returnchild=True)
-        if cutoff <= onset:
-            msg = 'Invalid cutoff radius {:f}'.format(cutoff)
+        aa, child = query.getvalue(root, 'a', conv.float0, returnchild=True)
+        if aa <= 0.0:
+            msg = 'Invalid potential slope (a) {:f}'.format(aa)
             raise hsd.HSDInvalidTagValueException(msg=msg, node=child)
 
-        vmax, child = query.getvalue(
-            root, 'vmax', conv.float0, defvalue=100.0, returnchild=True)
-        if vmax <= 0.0:
-            msg = 'Invalid potential well height {:f}'.format(vmax)
+        r0, child = query.getvalue(root, 'r0', conv.float0, returnchild=True)
+        if r0 < 0.0:
+            msg = 'Invalid potential half-height radius {:f}'.format(r0)
             raise hsd.HSDInvalidTagValueException(msg=msg, node=child)
 
         myself = cls()
         myself.compid = SUPPORTED_COMPRESSIONS[
             myself.__class__.__name__.lower()]
-        myself.onset = onset
-        myself.cutoff = cutoff
-        myself.vmax = vmax
+        myself.ww = ww
+        myself.aa = aa
+        myself.r0 = r0
 
         return myself
 
@@ -117,16 +115,16 @@ def tohsd(self, root, query, parentname=None):
         else:
             mynode = query.setchild(root, 'WoodsSaxonCompression')
 
-        query.setchildvalue(mynode, 'onset', conv.float0, self.onset)
-        query.setchildvalue(mynode, 'cutoff', conv.float0, self.cutoff)
-        query.setchildvalue(mynode, 'vmax', conv.float0, self.vmax)
+        query.setchildvalue(mynode, 'w', conv.float0, self.ww)
+        query.setchildvalue(mynode, 'a', conv.float0, self.aa)
+        query.setchildvalue(mynode, 'r0', conv.float0, self.r0)
 
 
     def __eq__(self, other):
-        onset_ok = abs(self.onset - other.onset) < 1e-8
-        cutoff_ok = abs(self.cutoff - other.cutoff) < 1e-8
-        vmax_ok = abs(self.vmax - other.vmax) < 1e-8
-        return onset_ok and cutoff_ok and vmax_ok
+        ww_ok = abs(self.ww - other.ww) < 1e-8
+        aa_ok = abs(self.aa - other.aa) < 1e-8
+        r0_ok = abs(self.r0 - other.r0) < 1e-8
+        return ww_ok and aa_ok and r0_ok
 
 
 # Registered compressions with corresponding hsd name as key
diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt
index c23d9ad2..1c47c778 100644
--- a/slateratom/lib/CMakeLists.txt
+++ b/slateratom/lib/CMakeLists.txt
@@ -18,7 +18,6 @@ set(sources-f90
   diagonalizations.f90
   globals.f90
   hamiltonian.f90
-  input.f90
   integration.f90
   lapack.f90
   mixer.f90
@@ -33,6 +32,7 @@ set(sources-fpp
   blasroutines.F90
   broydenmixer.F90
   confinement.F90
+  input.F90
   lapackroutines.F90
   simplemixer.F90)
 
diff --git a/slateratom/lib/confinement.F90 b/slateratom/lib/confinement.F90
index 777dcfaf..8ea748d4 100644
--- a/slateratom/lib/confinement.F90
+++ b/slateratom/lib/confinement.F90
@@ -107,14 +107,14 @@ end subroutine getSupervec
   !> Input for the Woods-Saxon confinement.
   type :: TWsConfInp
 
-    !> onset radii (Woods-Saxon compression)
-    real(dp) :: rOnset(0:4)
+    !> potential heights (W)
+    real(dp) :: ww(0:4)
 
-    !> cutoff of confinement (Woods-Saxon compression)
-    real(dp) :: rCut(0:4)
+    !> potential slopes (a)
+    real(dp) :: aa(0:4)
 
-    !> potential well height of confinement (Woods-Saxon compression)
-    real(dp) :: vMax(0:4)
+    !> half-height radii (r0)
+    real(dp) :: r0(0:4)
 
   end type TWsConfInp
 
@@ -135,10 +135,10 @@ end subroutine getSupervec
   type, extends(TConf) :: TPowerConf
     private
 
-    !> confinement radii (power compression)
+    !> confinement radii
     real(dp) :: r0(0:4)
 
-    !> power of confinement (power compression)
+    !> power of confinement
     real(dp) :: power(0:4)
 
   contains
@@ -153,14 +153,14 @@ end subroutine getSupervec
   type, extends(TConf) :: TWsConf
     private
 
-    !> onset radii (Woods-Saxon compression)
-    real(dp) :: rOnset(0:4)
+    !> potential heights (W)
+    real(dp) :: ww(0:4)
 
-    !> cutoff of confinement (Woods-Saxon compression)
-    real(dp) :: rCut(0:4)
+    !> potential slopes (a)
+    real(dp) :: aa(0:4)
 
-    !> potential well height of confinement (Woods-Saxon compression)
-    real(dp) :: vMax(0:4)
+    !> half-height radii (r0)
+    real(dp) :: r0(0:4)
 
   contains
 
@@ -214,9 +214,9 @@ subroutine TWsConf_init(this, input)
     !> input data
     type(TWsConfInp), intent(inout) :: input
 
-    this%rOnset(0:) = input%rOnset
-    this%rCut(0:) = input%rCut
-    this%vMax(0:) = input%vMax
+    this%ww(0:) = input%ww
+    this%aa(0:) = input%aa
+    this%r0(0:) = input%r0
 
   end subroutine TWsConf_init
 
@@ -239,16 +239,14 @@ subroutine TPowerConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf)
     !> Woods-Saxon potential on grid
     real(dp), intent(out) :: vconf(:, 0:)
 
-    !! auxiliary variables
-    integer :: ii, ll
+    !! angular momentum
+    integer :: 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
+      vconf(:, ll) = getPowerPotential(abcissa, this%r0(ll), this%power(ll))
     end do
 
   end subroutine TPowerConf_getPotOnGrid
@@ -272,16 +270,14 @@ subroutine TWsConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf)
     !> Woods-Saxon potential on grid
     real(dp), intent(out) :: vconf(:, 0:)
 
-    !! auxiliary variables
-    integer :: ii, ll
+    !! angular momentum
+    integer :: 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
+      vconf(:, ll) = getWSPotential(abcissa, this%ww(ll), this%aa(ll), this%r0(ll))
     end do
 
   end subroutine TWsConf_getPotOnGrid
@@ -499,7 +495,7 @@ end subroutine dft_conf_matrixelement
 
 
   !> Returns the Power potential for a given radial distance.
-  pure function getPowerPotential(rr, r0, power) result(pot)
+  elemental impure function getPowerPotential(rr, r0, power) result(pot)
 
     !> radial distance
     real(dp), intent(in) :: rr
@@ -513,37 +509,35 @@ pure function getPowerPotential(rr, r0, power) result(pot)
     !> Power potential at evaluation point
     real(dp) :: pot
 
+    @:ASSERT(rr >= 0.0_dp)
+
     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)
+  !! see J. Chem. Theory Comput. 12, 1, 53-64 (2016) eqn. 4.
+  elemental impure function getWSPotential(rr, ww, aa, r0) result(pot)
 
     !> radial distance
     real(dp), intent(in) :: rr
 
-    !> onset radius
-    real(dp), intent(in) :: rOnset
+    !> potential height (W)
+    real(dp), intent(in) :: ww
 
-    !> cutoff radius
-    real(dp), intent(in) :: rCut
+    !> potential slope (a)
+    real(dp), intent(in) :: aa
 
-    !> potential well height
-    real(dp), intent(in) :: vMax
+    !> half-height radius (r0)
+    real(dp), intent(in) :: r0
 
     !> Woods-Saxon potential at evaluation point
     real(dp) :: pot
 
-    ! case: rr <= rOnset
-    pot = 0.0_dp
+    @:ASSERT(rr >= 0.0_dp)
 
-    if (rr > rOnset .and. rr < rCut) then
-      pot = vMax / (rCut - rOnset) * (rr - rOnset)
-    elseif (rr >= rCut) then
-      pot = vMax
-    end if
+    pot = ww / (1.0 + exp(-aa * (rr - r0)))
 
   end function getWSPotential
 
diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.F90
similarity index 97%
rename from slateratom/lib/input.f90
rename to slateratom/lib/input.F90
index 16e9b994..138e3b94 100644
--- a/slateratom/lib/input.f90
+++ b/slateratom/lib/input.F90
@@ -1,3 +1,5 @@
+#:include 'common.fypp'
+
 !> Module that reads input values from stdin.
 module input
 
@@ -180,10 +182,13 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
         read(*,*) confInp%power%r0(ii), confInp%power%power(ii)
       end do
     case(confType%ws)
-      write(*, '(A)') 'Enter parameters Ronset, Rcut and Vmax'
+      write(*, '(A)') 'Enter parameters compr. height, slope and half-height radius'
       do ii = 0, max_l
         write(*, '(A,I3)') 'l=', ii
-        read(*,*) confInp%ws%rOnset(ii), confInp%ws%rCut(ii), confInp%ws%vMax(ii)
+        read(*,*) confInp%ws%ww(ii), confInp%ws%aa(ii), confInp%ws%r0(ii)
+        @:ASSERT(confInp%ws%ww(ii) >= 0.0)
+        @:ASSERT(confInp%ws%aa(ii) > 0.0)
+        @:ASSERT(confInp%ws%r0(ii) >= 0.0)
       end do
     case default
       error stop 'Invalid confinement potential.'
@@ -429,8 +434,8 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a
         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= ', confInp%ws%rOnset(ii),&
-            & ' Rcutoff= ', confInp%ws%rCut(ii), ' Vmax= ', confInp%ws%vMax(ii)
+        write(*, '(A,I3,A,E15.7,A,E15.7,A,E15.7)') 'l= ', ii, ', W= ', confInp%ws%ww(ii),&
+            & ' a= ', confInp%ws%aa(ii), ' r0= ', confInp%ws%r0(ii)
       end select
     end do
 

From 60c892b519c1ab1bffb3d3f5daad9872c40926d9 Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Tue, 16 Jul 2024 17:58:35 +0200
Subject: [PATCH 4/8] Add example to mio/skdef.hsd

---
 examples/mio/skdef.hsd | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/examples/mio/skdef.hsd b/examples/mio/skdef.hsd
index b0ce861e..6e03fc0e 100644
--- a/examples/mio/skdef.hsd
+++ b/examples/mio/skdef.hsd
@@ -85,6 +85,12 @@ AtomParameters {
     DftbAtom {
       ShellResolved = No
       DensityCompression = PowerCompression { Power = 2; Radius = 2.5 }
+      # Alternatively: Woods-Saxon compression potential
+      # (see J. Chem. Theory Comput. 12, 1, 53-64 (2016) eqn. 4.)
+      # With default W = 100.0:
+      # DensityCompression = WoodsSaxonCompression { a = 2.0; r0 = 6.0 }
+      # or with manual W:
+      # DensityCompression = WoodsSaxonCompression { W = 300.0; a = 2.0; r0 = 6.0 }
       WaveCompressions = SingleAtomCompressions {
 	S = PowerCompression { Power = 2; Radius = 3.0 }
       }

From 2f14d28032f21936ff54cb3627dbbca436f4ad62 Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Wed, 17 Jul 2024 11:06:21 +0200
Subject: [PATCH 5/8] Add regression test covering the WS compression

---
 .../Non-Relativistic/{ => Power}/_C-C.skf     |  0
 .../Non-Relativistic/{ => Power}/_C-O.skf     |  0
 .../Non-Relativistic/{ => Power}/_O-C.skf     |  0
 .../Non-Relativistic/{ => Power}/_O-O.skf     |  0
 .../Non-Relativistic/{ => Power}/config       |  0
 .../Non-Relativistic/{ => Power}/skdef.hsd    |  0
 .../LDA-PW91/Non-Relativistic/WS/_C-C.skf     | 92 +++++++++++++++++++
 .../LDA-PW91/Non-Relativistic/WS/config       |  1 +
 .../LDA-PW91/Non-Relativistic/WS/skdef.hsd    | 74 +++++++++++++++
 test/prog/sktable/tests                       |  3 +-
 10 files changed, 169 insertions(+), 1 deletion(-)
 rename test/prog/sktable/LDA-PW91/Non-Relativistic/{ => Power}/_C-C.skf (100%)
 rename test/prog/sktable/LDA-PW91/Non-Relativistic/{ => Power}/_C-O.skf (100%)
 rename test/prog/sktable/LDA-PW91/Non-Relativistic/{ => Power}/_O-C.skf (100%)
 rename test/prog/sktable/LDA-PW91/Non-Relativistic/{ => Power}/_O-O.skf (100%)
 rename test/prog/sktable/LDA-PW91/Non-Relativistic/{ => Power}/config (100%)
 rename test/prog/sktable/LDA-PW91/Non-Relativistic/{ => Power}/skdef.hsd (100%)
 create mode 100644 test/prog/sktable/LDA-PW91/Non-Relativistic/WS/_C-C.skf
 create mode 100644 test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config
 create mode 100644 test/prog/sktable/LDA-PW91/Non-Relativistic/WS/skdef.hsd

diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/_C-C.skf b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_C-C.skf
similarity index 100%
rename from test/prog/sktable/LDA-PW91/Non-Relativistic/_C-C.skf
rename to test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_C-C.skf
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/_C-O.skf b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_C-O.skf
similarity index 100%
rename from test/prog/sktable/LDA-PW91/Non-Relativistic/_C-O.skf
rename to test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_C-O.skf
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/_O-C.skf b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_O-C.skf
similarity index 100%
rename from test/prog/sktable/LDA-PW91/Non-Relativistic/_O-C.skf
rename to test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_O-C.skf
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/_O-O.skf b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_O-O.skf
similarity index 100%
rename from test/prog/sktable/LDA-PW91/Non-Relativistic/_O-O.skf
rename to test/prog/sktable/LDA-PW91/Non-Relativistic/Power/_O-O.skf
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/config b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/config
similarity index 100%
rename from test/prog/sktable/LDA-PW91/Non-Relativistic/config
rename to test/prog/sktable/LDA-PW91/Non-Relativistic/Power/config
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/skdef.hsd b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/skdef.hsd
similarity index 100%
rename from test/prog/sktable/LDA-PW91/Non-Relativistic/skdef.hsd
rename to test/prog/sktable/LDA-PW91/Non-Relativistic/Power/skdef.hsd
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/_C-C.skf b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/_C-C.skf
new file mode 100644
index 00000000..6b60e8cd
--- /dev/null
+++ b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/_C-C.skf
@@ -0,0 +1,92 @@
+0.020000 69
+   0.000000000000E+00  -1.991451862433E-01  -5.008031450092E-01  -4.388315590227E-02   0.000000000000E+00   3.643020862183E-01   3.643020862183E-01   0.000000000000E+00   2.000000000000E+00   2.000000000000E+00
+   1.201000000000E+01 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -6.381760667102E-01  -1.457225383447E+00   0.000000000000E+00  -2.781529807256E-01  -8.318367775611E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   5.572139974757E-01   8.404174953383E-01   0.000000000000E+00  -1.905803575596E-01   8.094163478854E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -5.837314388851E-01  -1.429923183678E+00   0.000000000000E+00  -2.486977954925E-01  -8.057290724036E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   5.329323332977E-01   8.308907927640E-01   0.000000000000E+00  -1.993756789262E-01   8.016932616479E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -5.304822244232E-01  -1.402644467587E+00   0.000000000000E+00  -2.189403014961E-01  -7.824076380814E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   5.084716804364E-01   8.211978085221E-01   0.000000000000E+00  -2.081949854654E-01   7.940460008169E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -4.784955599301E-01  -1.375424366584E+00   0.000000000000E+00  -1.890194583538E-01  -7.616582659754E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.838707508287E-01   8.113488238503E-01   0.000000000000E+00  -2.170212085455E-01   7.864756723071E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -4.278291849715E-01  -1.348295714043E+00   0.000000000000E+00  -1.590631634127E-01  -7.432748695947E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.591674096031E-01   8.013540366786E-01   0.000000000000E+00  -2.258374760168E-01   7.789825131351E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -3.785319908731E-01  -1.321289112944E+00   0.000000000000E+00  -1.291886172414E-01  -7.270600733583E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.343986302243E-01   7.912235486829E-01   0.000000000000E+00  -2.346271774869E-01   7.715659951136E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -3.306445689180E-01  -1.294433007484E+00   0.000000000000E+00  -9.950272526593E-02  -7.128256222100E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.096004547687E-01   7.809673533215E-01   0.000000000000E+00  -2.433740223936E-01   7.642249219751E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.841997537569E-01  -1.267753757839E+00   0.000000000000E+00  -7.010252715559E-02  -7.003926401807E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.848079588999E-01   7.705953248032E-01   0.000000000000E+00  -2.520620912822E-01   7.569575191556E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.392231584231E-01  -1.241275717315E+00   0.000000000000E+00  -4.107564657073E-02  -6.895917620483E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.600552211412E-01   7.601172079331E-01   0.000000000000E+00  -2.606758807177E-01   7.497615165227E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.957336978732E-01  -1.215021311211E+00   0.000000000000E+00  -1.250075478206E-02  -6.802631588775E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.353752960671E-01   7.495426087830E-01   0.000000000000E+00  -2.692003422690E-01   7.426342243672E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.537440984636E-01  -1.189011116767E+00   0.000000000000E+00   1.555195751402E-02  -6.722564753587E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.108001910672E-01   7.388809861370E-01   0.000000000000E+00  -2.776209160059E-01   7.355726030042E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.132613912126E-01  -1.163263943646E+00   0.000000000000E+00   4.302030512984E-02  -6.654306944109E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.863608463640E-01   7.281416436623E-01   0.000000000000E+00  -2.859235589456E-01   7.285733263414E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -7.428738709134E-02  -1.137796914454E+00   0.000000000000E+00   6.984963435550E-02  -6.596539424110E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.620871180001E-01   7.173337227586E-01   0.000000000000E+00  -2.940947688739E-01   7.216328397808E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -3.681913294126E-02  -1.112625544847E+00   0.000000000000E+00   9.599235365450E-02  -6.548032465922E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.380077635382E-01   7.064661960417E-01   0.000000000000E+00  -3.021216039532E-01   7.147474128151E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -8.493469272953E-04  -1.087763822859E+00   0.000000000000E+00   1.214074722414E-01  -6.507642545890E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.141504302464E-01   6.955478614191E-01   0.000000000000E+00  -3.099916985090E-01   7.079131866784E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.363316728233E-02  -1.063224287081E+00   0.000000000000E+00   1.460601492153E-01  -6.474309247389E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.905416455727E-01   6.845873367190E-01   0.000000000000E+00  -3.176932753679E-01   7.011262173972E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   6.664312816847E-02  -1.039018103420E+00   0.000000000000E+00   1.699212554858E-01  -6.447051945662E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.672068097348E-01   6.735930548359E-01   0.000000000000E+00  -3.252151550938E-01   6.943825145756E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   9.819840617581E-02  -1.015155140164E+00   0.000000000000E+00   1.929669503337E-01  -6.424966338368E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.441701902785E-01   6.625732593598E-01   0.000000000000E+00  -3.325467624471E-01   6.876780762349E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.283196758107E-01  -9.916440411512E-01   0.000000000000E+00   2.151782740909E-01  -6.407220876626E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.214549184798E-01   6.515360006571E-01   0.000000000000E+00  -3.396781303658E-01   6.810089200070E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.570300899489E-01  -9.684922968449E-01   0.000000000000E+00   2.365407580943E-01  -6.393053143369E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   9.908298748523E-02   6.404891323766E-01   0.000000000000E+00  -3.465999017406E-01   6.743711109690E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.843549760981E-01  -9.457063131805E-01   0.000000000000E+00   2.570440527715E-01  -6.381766218752E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   7.707525210303E-02   6.294403083537E-01   0.000000000000E+00  -3.533033292340E-01   6.677607863830E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.103215541896E-01  -9.232914780520E-01   0.000000000000E+00   2.766815744410E-01  -6.372725066125E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   5.545143017344E-02   6.183969798905E-01   0.000000000000E+00  -3.597802733660E-01   6.611741775903E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.349586753171E-01  -9.012522253536E-01   0.000000000000E+00   2.954501711562E-01  -6.365352966526E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.423010545896E-02   6.073663933909E-01   0.000000000000E+00  -3.660231990690E-01   6.546076292918E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.582965806941E-01  -8.795920965074E-01   0.000000000000E+00   3.133498076952E-01  -6.359128024740E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   1.342873200677E-02   5.963555883316E-01   0.000000000000E+00  -3.720251708888E-01   6.480576164255E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   2.803666799779E-01  -8.583137994394E-01   0.000000000000E+00   3.303832695908E-01  -6.353579765553E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -6.936360055855E-03   5.853713955533E-01   0.000000000000E+00  -3.777798469921E-01   6.415207588408E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.012013479961E-01  -8.374192649818E-01   0.000000000000E+00   3.465558859085E-01  -6.348285834957E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.684995732541E-02   5.744204358554E-01   0.000000000000E+00  -3.832814721179E-01   6.349938339491E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.208337388333E-01  -8.169097007063E-01   0.000000000000E+00   3.618752703227E-01  -6.342868817604E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -4.629795455255E-02   5.635091188830E-01   0.000000000000E+00  -3.885248695958E-01   6.284737875179E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.392976161719E-01  -7.967856422082E-01   0.000000000000E+00   3.763510799036E-01  -6.336993178767E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -6.526734472286E-02   5.526436422931E-01   0.000000000000E+00  -3.935054325365E-01   6.219577427610E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.566271987543E-01  -7.770470018782E-01   0.000000000000E+00   3.899947909200E-01  -6.330362336431E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -8.374620780993E-02   5.418299911891E-01   0.000000000000E+00  -3.982191142866E-01   6.154430078642E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.728570198272E-01  -7.576931152134E-01   0.000000000000E+00   4.028194908779E-01  -6.322715866867E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.017236982220E-01   5.310739378152E-01   0.000000000000E+00  -4.026624182258E-01   6.089270820754E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   3.880217994516E-01  -7.387227847312E-01   0.000000000000E+00   4.148396859594E-01  -6.313826845098E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.191900309643E-01   5.203810415005E-01   0.000000000000E+00  -4.068323869760E-01   6.024076604755E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.021563286043E-01  -7.201343215593E-01   0.000000000000E+00   4.260711229924E-01  -6.303499320091E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.361364665392E-01   5.097566488469E-01   0.000000000000E+00  -4.107265910792E-01   5.958826375365E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.152953640604E-01  -7.019255847814E-01   0.000000000000E+00   4.365306250702E-01  -6.291565923171E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.525552946120E-01   4.992058941510E-01   0.000000000000E+00  -4.143431171952E-01   5.893501095672E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.274735331198E-01  -6.840940186252E-01   0.000000000000E+00   4.462359399491E-01  -6.277885607157E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.684398164677E-01   4.887337000562E-01   0.000000000000E+00  -4.176805558613E-01   5.828083761339E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.387252473247E-01  -6.666366875815E-01   0.000000000000E+00   4.552056003733E-01  -6.262341512891E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.837843262935E-01   4.783447784258E-01   0.000000000000E+00  -4.207379888527E-01   5.762559405403E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.490846244013E-01  -6.495503095445E-01   0.000000000000E+00   4.634587955139E-01  -6.244838959280E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -1.985840913198E-01   4.680436314328E-01   0.000000000000E+00  -4.235149761748E-01   5.696915094414E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.585854177437E-01  -6.328312870641E-01   0.000000000000E+00   4.710152527495E-01  -6.225303552540E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.128353308613E-01   4.578345528602E-01   0.000000000000E+00  -4.260115427168E-01   5.631139916629E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.672609528388E-01  -6.164757367998E-01   0.000000000000E+00   4.778951290669E-01  -6.203679410077E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.265351942999E-01   4.477216296049E-01   0.000000000000E+00  -4.282281645922E-01   5.565224962892E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.751440701048E-01  -6.004795172627E-01   0.000000000000E+00   4.841189114093E-01  -6.179927494290E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.396817380574E-01   4.377087433806E-01   0.000000000000E+00  -4.301657551907E-01   5.499163300810E-01
+   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   4.822670736819E-01  -5.848382549295E-01   0.000000000000E+00   4.897073253495E-01  -6.154024051534E-01   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00   0.000000000000E+00  -2.522739016091E-01   4.277995726138E-01   0.000000000000E+00  -4.318256509620E-01   5.432949942779E-01
+
+
+Spline
+12 0.0553585
+112.9353346817185 2.801373701455403 -0.1119994835253462
+0.035 0.0375    0.204206 -35.71077211012958 2016.504000000031 24177.93762071238
+0.0375 0.04    0.12791 -25.17491577974109 2197.838532155373 -120889.6881035729
+0.04 0.0425    0.07682029999999999 -16.45240477090621 1291.165871378576 -57585.58520643491
+0.0425 0.045    0.0428593 -11.07630513663398 859.2739823303137 16659.22892930921
+0.045 0.04533    0.0207993 -6.467574682557872 984.2181993001326 -2167173.572075024
+0.04533 0.045334    0.0186943 -6.526006277016704 -1161.283637054166 353213222.4907721
+0.045334 0.046259    0.0186682 -6.518342311433599 3077.275032831984 -1324559.571220061
+0.046259 0.047184    0.0142234 -4.225362350069925 -598.3777773036936 561811.1110751317
+0.047184 0.0493131    0.0102476 -3.890262342340788 960.6480559297889 -100763.5210502349
+0.0493131 0.0503195    0.00534702 -1.169934109375229 317.0412179256228 -143026.9144497911
+0.0503195 0.0513259    0.00434492 -0.9663840979460291 -114.7856421811885 10348.58893883691
+0.0513259 0.0553585    0.00326664 -1.165980214261954 -83.5411824570522 -5782.515169399558 27636944.82683195 -3877959552.095367
+
+This SPLINE is just a DUMMY-SPLINE!!!!!!!!!!!!!!!
+
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config
new file mode 100644
index 00000000..fef984ce
--- /dev/null
+++ b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config
@@ -0,0 +1 @@
+C C
\ No newline at end of file
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/skdef.hsd b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/skdef.hsd
new file mode 100644
index 00000000..7505e52e
--- /dev/null
+++ b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/skdef.hsd
@@ -0,0 +1,74 @@
+SkdefVersion = 1
+
+Globals {
+  Superposition = density
+  XCFunctional = lda {}
+}
+
+AtomParameters {
+
+  C {
+    AtomConfig {
+      AtomicNumber = 6
+      Mass = 12.01
+      Occupations {
+	1S = 1.0 1.0
+	2S = 1.0 1.0
+	2P = 2.0 0.0
+      }
+      ValenceShells = 2s 2p
+      Relativistics = None
+    }
+    DftbAtom {
+      ShellResolved = No
+      DensityCompression = WoodsSaxonCompression { W = 10.0; a = 2.0; r0 = 7.0 }
+      WaveCompressions = SingleAtomCompressions {
+	S = WoodsSaxonCompression { W = 10.0; a = 2.0; r0 = 2.7 }
+	P = WoodsSaxonCompression { W = 10.0; a = 2.0; r0 = 2.7 }
+      }
+    }
+  }
+
+}
+
+
+OnecenterParameters {
+
+  $StandardDeltaFilling {
+    DeltaFilling = 0.01
+  }
+
+  C {
+    $StandardDeltaFilling
+    Calculator = SlaterAtom {
+    MaxScfIterations = 120
+    ScfTolerance = 1.0e-10
+      Exponents {
+	S = 0.5 1.14 2.62 6.0
+	P = 0.5 1.14 2.62 6.0
+      }
+      MaxPowers {
+	S = 3
+	P = 3
+      }
+    }
+  }
+
+}
+
+TwoCenterParameters {
+
+  $EqGrid = EquidistantGrid {
+      GridStart = 0.6
+      GridSeparation = 0.02
+      Tolerance = 5e-5
+      MaxDistance = -1.4
+  }
+
+  $SkTwocnt_300_150 = Sktwocnt {
+    IntegrationPoints = 300 150
+  }
+
+  C-C { Grid = $EqGrid; Calculator = $SkTwocnt_300_150 }
+
+}
diff --git a/test/prog/sktable/tests b/test/prog/sktable/tests
index 4ef363ba..af0e0790 100644
--- a/test/prog/sktable/tests
+++ b/test/prog/sktable/tests
@@ -3,7 +3,8 @@
 
 #:include 'common.fypp'
 
-LDA-PW91/Non-Relativistic
+LDA-PW91/Non-Relativistic/Power
+LDA-PW91/Non-Relativistic/WS
 GGA-PBE96/Non-Relativistic
 HYB-PBE0/Non-Relativistic
 HYB-B3LYP/Non-Relativistic

From 46a698453a39aaf807f98228621f25da6a22646c Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Thu, 18 Jul 2024 19:06:56 +0200
Subject: [PATCH 6/8] Apply the feedback of the code review

---
 slateratom/lib/confinement.F90 | 130 ++++++---------------------------
 slateratom/lib/input.F90       |   2 +
 2 files changed, 24 insertions(+), 108 deletions(-)

diff --git a/slateratom/lib/confinement.F90 b/slateratom/lib/confinement.F90
index 8ea748d4..9d8afcc6 100644
--- a/slateratom/lib/confinement.F90
+++ b/slateratom/lib/confinement.F90
@@ -20,8 +20,9 @@ module confinement
   type, abstract :: TConf
   contains
 
+    procedure, nopass :: getSupervec => TConf_getSupervec
+
     procedure(getPotOnGrid), deferred :: getPotOnGrid
-    procedure(getSupervec), deferred :: getSupervec
 
   end type TConf
 
@@ -123,10 +124,10 @@ end subroutine getSupervec
   type :: TConfInp
 
     !> Power confinement input structure
-    type(TPowerConfInp) :: power
+    type(TPowerConfInp), allocatable :: power
 
     !> Woods-Saxon confinement input structure
-    type(TWsConfInp) :: ws
+    type(TWsConfInp), allocatable :: ws
 
   end type TConfInp
 
@@ -144,7 +145,6 @@ end subroutine getSupervec
   contains
 
     procedure :: getPotOnGrid => TPowerConf_getPotOnGrid
-    procedure :: getSupervec => TPowerConf_getSupervec
 
   end type TPowerConf
 
@@ -165,7 +165,6 @@ end subroutine getSupervec
   contains
 
     procedure :: getPotOnGrid => TWsConf_getPotOnGrid
-    procedure :: getSupervec => TWsConf_getSupervec
 
   end type TWsConf
 
@@ -283,12 +282,10 @@ subroutine TWsConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf)
   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
+  !> Constructs the DFT confinement supervector to be added to the Fock matrix, by calculating the
+  !! single matrix elements and putting them together.
+  subroutine TConf_getSupervec(max_l, num_mesh_points, abcissa, weight, num_alpha, alpha,&
+      & poly_order, vconf, vconf_matrix)
 
     !> maximum angular momentum
     integer, intent(in) :: max_l
@@ -311,108 +308,25 @@ subroutine TPowerConf_getSupervec(this, max_l, num_mesh_points, abcissa, weight,
     !> highest polynomial order + l in each shell
     integer, intent(in) :: poly_order(0:)
 
-    !> Power potential on grid
+    !> confinement potential on grid
     real(dp), intent(in) :: vconf(:, 0:)
 
-    !> Power confinement supervector
+    !> 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:)
+    !> single matrix element of the confinement potential
+    real(dp) :: vconf_matrixelement
 
-    !> Woods-Saxon confinement supervector
-    real(dp), intent(out) :: vconf_matrix(0:,:,:)
+    !> auxiliary variables
+    integer :: ii, jj, kk, ll, mm, ss, tt, start
 
     @: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
+    vconf_matrix(:,:,:) = 0.0_dp
+    vconf_matrixelement = 0.0_dp
 
     do ii = 0, max_l
       ss = 0
@@ -430,10 +344,10 @@ subroutine build_dft_conf_matrix(max_l, num_alpha, poly_order, alpha, num_mesh_p
               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)
+                  & alpha(ii, jj), kk, alpha(ii, ll), mm, ii, vconf_matrixelement)
 
-              conf_matrix(ii, ss, tt) = conf_matrixelement
-              conf_matrix(ii, tt, ss) = conf_matrixelement
+              vconf_matrix(ii, ss, tt) = vconf_matrixelement
+              vconf_matrix(ii, tt, ss) = vconf_matrixelement
 
             end do
           end do
@@ -441,10 +355,10 @@ subroutine build_dft_conf_matrix(max_l, num_alpha, poly_order, alpha, num_mesh_p
       end do
     end do
 
-  end subroutine build_dft_conf_matrix
+  end subroutine TConf_getSupervec
 
 
-  !> Calculates a single matrix element of the exchange correlation potential.
+  !> Calculates a single matrix element of the confinement potential.
   pure subroutine dft_conf_matrixelement(num_mesh_points, weight, abcissa, vconf, alpha1, poly1,&
       & alpha2, poly2, ll, conf_matrixelement)
 
@@ -475,7 +389,7 @@ pure subroutine dft_conf_matrixelement(num_mesh_points, weight, abcissa, vconf,
     !> angular momentum
     integer, intent(in) :: ll
 
-    !> single matrix element of the exchange correlation potential
+    !> single matrix element of the confinement potential
     real(dp), intent(out) :: conf_matrixelement
 
     !> stores product of two basis functions and r^2
diff --git a/slateratom/lib/input.F90 b/slateratom/lib/input.F90
index 138e3b94..c3f01454 100644
--- a/slateratom/lib/input.F90
+++ b/slateratom/lib/input.F90
@@ -176,12 +176,14 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     case(confType%none)
       continue
     case(confType%power)
+      allocate(confInp%power)
       write(*, '(A)') 'Enter parameters r_0 and power'
       do ii = 0, max_l
         write(*, '(A,I3)') 'l=', ii
         read(*,*) confInp%power%r0(ii), confInp%power%power(ii)
       end do
     case(confType%ws)
+      allocate(confInp%ws)
       write(*, '(A)') 'Enter parameters compr. height, slope and half-height radius'
       do ii = 0, max_l
         write(*, '(A,I3)') 'l=', ii

From 2e866b7cdc05b3fe3040d04bd08b373767ecdab9 Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Tue, 20 Aug 2024 20:53:50 +0200
Subject: [PATCH 7/8] Fix shelf search bug

---
 sktools/src/sktools/common.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/sktools/src/sktools/common.py b/sktools/src/sktools/common.py
index 8e3aa313..15610dae 100644
--- a/sktools/src/sktools/common.py
+++ b/sktools/src/sktools/common.py
@@ -528,11 +528,12 @@ def is_shelf_file_matching(shelf_file, mydict):
         db = shelve.open(shelf_file, 'r')
     except dbm.error:
         return False
-    if not type(db) is type(mydict):
-        return False
     match = True
     for key in mydict:
-        match = key in db and db[key] == mydict[key]
+        try:
+            match = key in db and db[key] == mydict[key]
+        except KeyError:
+            match = False
         if not match:
             return False
     return True

From 36ab95c9dd560bf8e773d74d205253176fc9554e Mon Sep 17 00:00:00 2001
From: Tammo van der Heide <vanderhe@uni-bremen.de>
Date: Sun, 25 Aug 2024 21:38:22 +0200
Subject: [PATCH 8/8] Address the feedback of the code review

Co-authored-by: Ben Hourahine <benjamin.hourahine@strath.ac.uk>
---
 slateratom/lib/confinement.F90                | 58 +------------------
 .../CAMY-B3LYP/Non-Relativistic/config        |  2 +-
 .../sktable/CAMY-PBEh/Non-Relativistic/config |  2 +-
 .../sktable/GGA-PBE96/Non-Relativistic/config |  2 +-
 .../sktable/HYB-B3LYP/Non-Relativistic/config |  2 +-
 .../sktable/HYB-PBE0/Non-Relativistic/config  |  2 +-
 .../LDA-PW91/Non-Relativistic/Power/config    |  2 +-
 .../LDA-PW91/Non-Relativistic/WS/config       |  2 +-
 8 files changed, 8 insertions(+), 64 deletions(-)

diff --git a/slateratom/lib/confinement.F90 b/slateratom/lib/confinement.F90
index 9d8afcc6..7fd523ed 100644
--- a/slateratom/lib/confinement.F90
+++ b/slateratom/lib/confinement.F90
@@ -1,6 +1,6 @@
 #:include 'common.fypp'
 
-!> Module that builds up various supervectors.
+!> Module that builds up various confining potential supervectors.
 module confinement
 
   use common_accuracy, only : dp
@@ -455,60 +455,4 @@ elemental impure function getWSPotential(rr, ww, aa, r0) result(pot)
 
   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/test/prog/sktable/CAMY-B3LYP/Non-Relativistic/config b/test/prog/sktable/CAMY-B3LYP/Non-Relativistic/config
index bfa360b6..03fd460a 100644
--- a/test/prog/sktable/CAMY-B3LYP/Non-Relativistic/config
+++ b/test/prog/sktable/CAMY-B3LYP/Non-Relativistic/config
@@ -1 +1 @@
-C,N C,N
\ No newline at end of file
+C,N C,N
diff --git a/test/prog/sktable/CAMY-PBEh/Non-Relativistic/config b/test/prog/sktable/CAMY-PBEh/Non-Relativistic/config
index 57bc8b62..897a44bc 100644
--- a/test/prog/sktable/CAMY-PBEh/Non-Relativistic/config
+++ b/test/prog/sktable/CAMY-PBEh/Non-Relativistic/config
@@ -1 +1 @@
-O,N O,N
\ No newline at end of file
+O,N O,N
diff --git a/test/prog/sktable/GGA-PBE96/Non-Relativistic/config b/test/prog/sktable/GGA-PBE96/Non-Relativistic/config
index e5eac9cd..fe945713 100644
--- a/test/prog/sktable/GGA-PBE96/Non-Relativistic/config
+++ b/test/prog/sktable/GGA-PBE96/Non-Relativistic/config
@@ -1 +1 @@
-N,O N,O
\ No newline at end of file
+N,O N,O
diff --git a/test/prog/sktable/HYB-B3LYP/Non-Relativistic/config b/test/prog/sktable/HYB-B3LYP/Non-Relativistic/config
index bfa360b6..03fd460a 100644
--- a/test/prog/sktable/HYB-B3LYP/Non-Relativistic/config
+++ b/test/prog/sktable/HYB-B3LYP/Non-Relativistic/config
@@ -1 +1 @@
-C,N C,N
\ No newline at end of file
+C,N C,N
diff --git a/test/prog/sktable/HYB-PBE0/Non-Relativistic/config b/test/prog/sktable/HYB-PBE0/Non-Relativistic/config
index bfa360b6..03fd460a 100644
--- a/test/prog/sktable/HYB-PBE0/Non-Relativistic/config
+++ b/test/prog/sktable/HYB-PBE0/Non-Relativistic/config
@@ -1 +1 @@
-C,N C,N
\ No newline at end of file
+C,N C,N
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/config b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/config
index 2b6041a2..d380badc 100644
--- a/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/config
+++ b/test/prog/sktable/LDA-PW91/Non-Relativistic/Power/config
@@ -1 +1 @@
-C,O C,O
\ No newline at end of file
+C,O C,O
diff --git a/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config
index fef984ce..3aab270d 100644
--- a/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config
+++ b/test/prog/sktable/LDA-PW91/Non-Relativistic/WS/config
@@ -1 +1 @@
-C C
\ No newline at end of file
+C C