Add support for mGGA, SAP guess, orbital rotation gradient, DIIS mixer

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

common/lib/accuracy.F90

@@ -3,7 +3,7 @@
 !! Not all routines use the string length specifications to set their character string lengths.
 module common_accuracy
 
-  use, intrinsic :: iso_fortran_env, only : real64
+  use, intrinsic :: iso_fortran_env, only : real64, real128
 
   implicit none
   private

common/lib/gridgenerator.F90

@@ -20,7 +20,7 @@ subroutine gengrid1_1(radQuad, rm, coordtrans, grid, weights)
 
     interface
       !> General interface for utilized coordinate transformation.
-      subroutine coordtrans(oldc, rm, newc, jacobi)
+      pure subroutine coordtrans(oldc, rm, newc, jacobi)
         use common_accuracy, only : dp
 
         !> old coordinate

common/lib/partition.F90

@@ -32,7 +32,7 @@ pure function partitionFunc(r1, r2, dist, partparams) result(res)
       real(dp), intent(in) :: partparams(:)
 
       !! resulting value of the partition function, between [0,1]
-      real(dp) :: res
+      real(dp):: res
 
     end function partitionFunc
 

common/lib/poisson.F90

@@ -11,11 +11,10 @@ module common_poisson
 
   use common_accuracy, only : dp
   use common_constants, only : pi, pi_hlf
-
+  use common_coordtrans, only : coordtrans_radial_becke1, coordtrans_radial_becke2
   use common_quadratures, only : TQuadrature, TQuadrature2D, gauss_chebyshev_quadrature
   use common_quadratures, only : lebedev_laikov_quadrature
   use common_gridgenerator, only : gengrid1_1, gengrid1_3, gengrid2_3
-  use common_coordtrans, only : coordtrans_radial_becke1, coordtrans_radial_becke2
   use common_partition, only : partition_becke_homo
   use common_anglib, only : initGaunt
   use common_finitedifferences, only : makeHelmholzFDMatrix7P, makePoissonFDMatrix7P

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, 'ywb97m': 16}
 
 INPUT_FILE = "sktwocnt.in"
 STDOUT_FILE = "output"
@@ -133,7 +135,7 @@ def _store_atomdata(self, workdir, atomdata, iatom):
                                                 iatom)
         xcn = self._functional.type
         if xcn in ('lcy-bnl', 'lcy-pbe', 'pbe0', 'b3lyp', 'camy-b3lyp',
-                   'camy-pbeh'):
+                   'camy-pbeh', 'ywb97m'):
             atomfiles.dens_wavefuncs = self._store_dens_wavefuncs(
                 workdir, atomdata.dens_wavefuncs, iatom)
         atomfiles.occshells = atomdata.occshells
@@ -227,7 +229,7 @@ def _write_twocnt_gridinfo(self, fp):
             fp.write("{:f}\n".format(self._functional.alpha))
             fp.write("{:s}\n".format(becke))
         # CAM functionals
-        elif self._functional.type in ('camy-b3lyp', 'camy-pbeh'):
+        elif self._functional.type in ('camy-b3lyp', 'camy-pbeh', 'ywb97m'):
             becke = '2000 194 11 1.0'
             fp.write("{:f} {:f} {:f}\n".format(self._functional.omega,
                                                self._functional.alpha,
@@ -243,7 +245,7 @@ def _write_twocnt_integration_parameters(self, fp):
 
     def _write_twocnt_atom_block(self, fp, atomfiles):
         if self._functional.type in ('lcy-bnl', 'lcy-pbe', 'pbe0', 'b3lyp',
-                                     'camy-b3lyp', 'camy-pbeh'):
+                                     'camy-b3lyp', 'camy-pbeh', 'ywb97m'):
             fp.write("{:d} {:d}\n".format(len(atomfiles.wavefuncs),
                                           len(atomfiles.dens_wavefuncs)))
         else:
@@ -253,7 +255,7 @@ def _write_twocnt_atom_block(self, fp, atomfiles):
             fp.write("'{}' {:d}\n".format(wavefuncfile, ll))
 
         if self._functional.type in ('lcy-bnl', 'lcy-pbe', 'pbe0', 'b3lyp',
-                                     'camy-b3lyp', 'camy-pbeh'):
+                                     'camy-b3lyp', 'camy-pbeh', 'ywb97m'):
             occdict = {}
             for xx in atomfiles.occshells:
                 occdict[xx[0]] = xx[1]

sktools/src/sktools/calculators/slateratom.py

@@ -19,7 +19,11 @@
 
 SUPPORTED_FUNCTIONALS = {'lda' : 2, 'pbe' : 3, 'blyp' : 4, 'lcy-pbe' : 5,
                          'lcy-bnl' : 6, 'pbe0' : 7, 'b3lyp' : 8,
-                         'camy-b3lyp' : 9, 'camy-pbeh' : 10}
+                         'camy-b3lyp' : 9, 'camy-pbeh' : 10, "tpss": 11,
+                         'scan': 12, 'r2scan': 13, 'r4scan': 14, 'task': 15,
+                         'task+cc': 16, 'ywb97m': 17}
+
+SUPPORTED_MIXERS = {1: 'simple', 2: 'broyden', 3: 'diis'}
 
 INPUT_FILE = "slateratom.in"
 STDOUT_FILE = "output"
@@ -64,12 +68,15 @@ 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):
@@ -81,7 +88,12 @@ 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):
@@ -168,16 +180,25 @@ 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.0 < self._settings.mixing_parameter <= 1.0:
+            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]
 
-            if xcfkey in ('lcy-pbe', 'lcy-bnl', 'camy-b3lyp', 'camy-pbeh'):
+            if xcfkey in ('lcy-pbe', 'lcy-bnl', 'camy-b3lyp', 'camy-pbeh', 'ywb97m'):
                 self._omega = functional.omega
             else:
                 self._omega = None
 
-            if xcfkey in ('camy-b3lyp', 'camy-pbeh'):
+            if xcfkey in ('camy-b3lyp', 'camy-pbeh', 'ywb97m'):
                 self._alpha = functional.alpha
                 self._beta = functional.beta
             elif xcfkey == 'pbe0':
@@ -294,7 +315,7 @@ def write(self, workdir):
                 "2000 194 11 1.0 \t{:s} Becke integrator settings"
                 .format(self._COMMENT)]
         # CAM functionals
-        elif xctype in ('camy-b3lyp', 'camy-pbeh'):
+        elif xctype in ('camy-b3lyp', 'camy-pbeh', 'ywb97m'):
             out += [
                 "{:g} {:g} {:g} \t{:s} ".format(
                     self._omega, self._alpha, self._beta, self._COMMENT) + \
@@ -365,8 +386,9 @@ def write(self, workdir):
 
         out.append("{:s} \t\t{:s} write eigenvectors".format(
             self._LOGICALSTRS[False], self._COMMENT))
-        out.append("{} {:g} \t\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):
@@ -566,19 +588,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

@@ -97,6 +97,46 @@ def fromhsd(cls, root, query):
         return myself
 
 
+class XCYWB97M(sc.ClassDict):
+    '''Range-separated YwB97M xc-functional.
+
+    Attributes
+    ----------
+    omega (float): range-separation parameter
+    alpha (float): fraction of the global exact HF exchange
+    beta (float): determines (alpha + beta) fraction of long-range HF exchange
+    '''
+
+    @classmethod
+    def fromhsd(cls, root, query):
+        '''Creates instance from a HSD-node and with given query object.'''
+        omega, child = query.getvalue(root, 'omega', conv.float0,
+                                      returnchild=True)
+        if omega <= 0.0:
+            raise hsd.HSDInvalidTagValueException(
+                msg='Invalid rs-parameter {:f}'.format(omega),
+                node=child)
+
+        alpha, child = query.getvalue(root, 'alpha', conv.float0,
+                                      returnchild=True)
+
+        beta, child = query.getvalue(root, 'beta', conv.float0,
+                                     returnchild=True)
+
+        if not alpha + beta > 0.0:
+            raise hsd.HSDInvalidTagValueException(
+                msg='Invalid CAM-parameter combination alpha={:f}, beta={:f}!\n'
+                .format(alpha, beta) +
+                'Should satisfy alpha + beta > 0.0', node=child)
+
+        myself = cls()
+        myself.type = 'ywb97m'
+        myself.omega = omega
+        myself.alpha = alpha
+        myself.beta = beta
+        return myself
+
+
 class XCCAMYPBEH(sc.ClassDict):
     '''Range-separated CAMY-PBEh xc-functional.
 
@@ -189,6 +229,72 @@ def fromhsd(cls, root, query):
         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.'''
 
@@ -244,5 +350,12 @@ def fromhsd(cls, root, query):
     'pbe0': XCPBE0,
     'b3lyp': XCB3LYP,
     'camy-b3lyp': XCCAMYB3LYP,
-    'camy-pbeh': XCCAMYPBEH
+    'camy-pbeh': XCCAMYPBEH,
+    'tpss': XCTPSS,
+    'scan': XCSCAN,
+    'r2scan': XCR2SCAN,
+    'r4scan': XCR4SCAN,
+    'task': XCTASK,
+    'task+cc': XCTASK_CC,
+    'ywb97m': XCYWB97M
 }