Rework WS potential in terms of W, a and r0

Add example to mio/skdef.hsd
Add regression test covering the WS compression
Fix shelf search bug

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

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

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)