Rework WS potential in terms of W, a and r0

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(

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

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)
 

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 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 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
 

slateratom/lib/input.f90 → 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