Implement mGGA, SAP guess, orbital gradient

.gitignore

@@ -11,3 +11,4 @@ __pycache__
 _gitmsg.saved.txt
 *.egg-info
 dist
+.vscode/

AUTHORS.rst

@@ -15,4 +15,6 @@ contributed to this package :
 
 * Christof Köhler (University of Bremen)
 
+* Alexander Maryewski (Karlsruhe Institute of Technology, Germany)
+
 * Thomas Niehaus (University of Lyon, France)

sktools/src/sktools/calculators/sktwocnt.py

@@ -11,7 +11,9 @@
 
 SUPPORTED_FUNCTIONALS = {'lda' : 1, 'pbe' : 2, 'blyp' : 3, 'lcy-pbe' : 4,
                          'lcy-bnl' : 5, 'pbe0' : 6, 'b3lyp' : 7,
-                         'camy-b3lyp' : 8, 'camy-pbeh' : 9}
+                         'camy-b3lyp' : 8, 'camy-pbeh' : 9, 'tpss': 10,
+                         'scan': 11, 'r2scan': 12, 'r4scan': 13, 'task': 14,
+                         'task+cc': 15}
 
 INPUT_FILE = "sktwocnt.in"
 STDOUT_FILE = "output"

sktools/src/sktools/calculators/slateratom.py

@@ -16,9 +16,25 @@
 
 LOGGER = logging.getLogger('slateratom')
 
-SUPPORTED_FUNCTIONALS = {'lda' : 2, 'pbe' : 3, 'blyp' : 4, 'lcy-pbe' : 5,
-                         'lcy-bnl' : 6, 'pbe0' : 7, 'b3lyp' : 8,
-                         'camy-b3lyp' : 9, 'camy-pbeh' : 10}
+SUPPORTED_FUNCTIONALS = {
+    "lda": 2,
+    "pbe": 3,
+    "blyp": 4,
+    "lcy-pbe": 5,
+    "lcy-bnl": 6,
+    "pbe0": 7,
+    "b3lyp": 8,
+    "camy-b3lyp": 9,
+    "camy-pbeh": 10,
+    "tpss": 11,
+    "scan": 12,
+    "r2scan": 13,
+    "r4scan": 14,
+    "task": 15,
+    "task+cc": 16
+}
+
+SUPPORTED_MIXERS = {1: 'simple', 2: "broyden"}
 
 INPUT_FILE = "slateratom.in"
 STDOUT_FILE = "output"
@@ -63,12 +79,14 @@ class SlaterAtomSettings(sc.ClassDict):
         Maximal power for every angular momentum.
     """
 
-    def __init__(self, exponents, maxpowers, scftol, maxscfiter):
+    def __init__(self, exponents, maxpowers, scftol, maxscfiter, mixer_id, mixing_parameter):
         super().__init__()
         self.exponents = exponents
         self.maxpowers = maxpowers
         self.scftol = scftol
         self.maxscfiter = maxscfiter
+        self.mixer = mixer_id
+        self.mixing_parameter = mixing_parameter
 
     @classmethod
     def fromhsd(cls, root, query):
@@ -80,7 +98,11 @@ def fromhsd(cls, root, query):
             root, "scftolerance", converter=conv.float0, defvalue=1.0e-10)
         maxscfiter = query.getvalue(
             root, "maxscfiterations", converter=conv.int0, defvalue=120)
-        return cls(exponents, maxpowers, scftol, maxscfiter)
+        mixer_id = query.getvalue(
+            root, "mixer", converter=conv.int0, defvalue=2)
+        mixing_parameter = query.getvalue(
+            root, "mixingparameter", converter=conv.float0, defvalue=0.1)
+        return cls(exponents, maxpowers, scftol, maxscfiter, mixer_id, mixing_parameter)
 
     def __eq__(self, other):
         if not isinstance(other, SlaterAtomSettings):
@@ -167,6 +189,14 @@ def __init__(self, settings, atomconfig, functional, compressions):
             msg = "Slateratom: Maximum number of SCF iterations must be >=1"
             raise sc.SkgenException(msg)
 
+        if self._settings.mixer not in SUPPORTED_MIXERS.keys():
+            msg = f"Slateratom: mixer {self._settings.mixer} not found"
+            raise sc.SkgenException(msg)
+
+        if not (0. <= self._settings.mixing_parameter <= 1.):
+            msg = "Slateratom: mixing parameter must lie within the [0, 1] interval"
+            raise sc.SkgenException(msg)
+
         if self.isXCFunctionalSupported(functional):
             xcfkey = functional.type
             self._functional = SUPPORTED_FUNCTIONALS[xcfkey]
@@ -326,8 +356,8 @@ 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(
-            2, 0.1, self._COMMENT))
+        out.append("{} {:g} \t\t{:s} mixer, mixing factor".format(
+            self._settings.mixer, self._settings.mixing_parameter, self._COMMENT))
 
         # Occupations
         for ll, occperl in enumerate(self._atomconfig.occupations):
@@ -527,19 +557,15 @@ def get_density012(self):
         density : GridData
            Grid data with the density and its first and second derivatives.
         """
-        fp = open(os.path.join(self._workdir, "dens.dat"), "r")
-        fp.readline()
-        fp.readline()
-        fp.readline()
-        fp.readline()
-        fp.readline()
-        ngrid = int(fp.readline())
-        # noinspection PyNoneFunctionAssignment,PyTypeChecker
-        dens = np.fromfile(fp, dtype=float, count=ngrid * 7, sep=" ")
-        fp.close()
-        dens.shape = (ngrid, 7)
-        grid = oc.RadialGrid(dens[:,0], dens[:,1])
-        density = dens[:,2:5]
+
+        with open(os.path.join(self._workdir, "dens.dat"), 'r') as handle:
+            lines_ = [x.split() for x in handle.readlines()]
+            datarray = np.asarray(lines_[6:], dtype=float)
+            grid = oc.RadialGrid(datarray[:,0], datarray[:,1])
+            density = datarray[:,2:5]
+            if datarray.shape[1] == 8:
+                density = np.column_stack([density, datarray[:, 7]])
+
         return oc.GridData(grid, density)
 
     def get_wavefunction012(self, ss, nn, ll):

sktools/src/sktools/xcfunctionals.py

@@ -188,6 +188,69 @@ def fromhsd(cls, root, query):
         myself.beta = 1.0
         return myself
 
+class XCTPSS(sc.ClassDict):
+    """mGGA TPSS xc-functional."""
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        """Creates instance from a HSD-node and with given query object."""
+        myself = cls()
+        myself.type = "tpss"
+        return myself
+
+
+class XCSCAN(sc.ClassDict):
+    """mGGA SCAN xc-functional."""
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        """Creates instance from a HSD-node and with given query object."""
+        myself = cls()
+        myself.type = "scan"
+        return myself
+
+
+class XCR2SCAN(sc.ClassDict):
+    """mGGA r2SCAN xc-functional."""
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        """Creates instance from a HSD-node and with given query object."""
+        myself = cls()
+        myself.type = "r2scan"
+        return myself
+
+
+class XCR4SCAN(sc.ClassDict):
+    """mGGA r4SCAN xc-functional."""
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        """Creates instance from a HSD-node and with given query object."""
+        myself = cls()
+        myself.type = "r4scan"
+        return myself
+
+
+class XCTASK(sc.ClassDict):
+    """mGGA TASK xc-functional."""
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        """Creates instance from a HSD-node and with given query object."""
+        myself = cls()
+        myself.type = "task"
+        return myself
+
+class XCTASK_CC(sc.ClassDict):
+    """mGGA TASK+CC xc-functional."""
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        """Creates instance from a HSD-node and with given query object."""
+        myself = cls()
+        myself.type = "task+cc"
+        return myself
 
 class XCLocal(sc.ClassDict):
     '''Local xc-functional.'''
@@ -235,14 +298,20 @@ def fromhsd(cls, root, query):
 
 # Registered xc-functionals with corresponding HSD name as key:
 XCFUNCTIONALS = {
-    'lcy-bnl': XCLCYBNL,
-    'lcy-pbe': XCLCYPBE,
-    'local': XCLocal,
-    'pbe': XCPBE,
-    'lda': XCLDA,
-    'blyp': XCBLYP,
-    'pbe0': XCPBE0,
-    'b3lyp': XCB3LYP,
-    'camy-b3lyp': XCCAMYB3LYP,
-    'camy-pbeh': XCCAMYPBEH
+    "lcy-bnl": XCLCYBNL,
+    "lcy-pbe": XCLCYPBE,
+    "local": XCLocal,
+    "pbe": XCPBE,
+    "lda": XCLDA,
+    "blyp": XCBLYP,
+    "pbe0": XCPBE0,
+    "b3lyp": XCB3LYP,
+    "camy-b3lyp": XCCAMYB3LYP,
+    "camy-pbeh": XCCAMYPBEH,
+    "tpss": XCTPSS,
+    "scan": XCSCAN,
+    "r2scan": XCR2SCAN,
+    "r4scan": XCR4SCAN,
+    "task": XCTASK,
+    "task+cc": XCTASK_CC
 }