Hand over WS compression to slateratom

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

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]

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

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

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