Fix formatting and suggest minor code changes

sktools/src/sktools/calculators/slateratom.py

@@ -16,25 +16,13 @@
 
 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,
-    "tpss": 11,
-    "scan": 12,
-    "r2scan": 13,
-    "r4scan": 14,
-    "task": 15,
-    "task+cc": 16
-}
-
-SUPPORTED_MIXERS = {1: 'simple', 2: "broyden"}
+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"
@@ -79,7 +67,8 @@ class SlaterAtomSettings(sc.ClassDict):
         Maximal power for every angular momentum.
     """
 
-    def __init__(self, exponents, maxpowers, scftol, maxscfiter, mixer_id, mixing_parameter):
+    def __init__(self, exponents, maxpowers, scftol, maxscfiter, mixer_id,
+                 mixing_parameter):
         super().__init__()
         self.exponents = exponents
         self.maxpowers = maxpowers
@@ -102,7 +91,8 @@ def fromhsd(cls, root, query):
             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)
+        return cls(exponents, maxpowers, scftol, maxscfiter, mixer_id,
+                   mixing_parameter)
 
     def __eq__(self, other):
         if not isinstance(other, SlaterAtomSettings):
@@ -193,8 +183,9 @@ def __init__(self, settings, atomconfig, functional, compressions):
             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"
+        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):
@@ -357,7 +348,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} mixer, mixing factor".format(
-            self._settings.mixer, self._settings.mixing_parameter, self._COMMENT))
+            self._settings.mixer, self._settings.mixing_parameter,
+            self._COMMENT))
 
         # Occupations
         for ll, occperl in enumerate(self._atomconfig.occupations):
@@ -561,8 +553,8 @@ def get_density012(self):
         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]
+            grid = oc.RadialGrid(datarray[:, 0], datarray[:, 1])
+            density = datarray[:, 2:5]
             if datarray.shape[1] == 8:
                 density = np.column_stack([density, datarray[:, 7]])
 

sktools/src/sktools/xcfunctionals.py

@@ -188,6 +188,7 @@ def fromhsd(cls, root, query):
         myself.beta = 1.0
         return myself
 
+
 class XCTPSS(sc.ClassDict):
     """mGGA TPSS xc-functional."""
 
@@ -241,7 +242,8 @@ def fromhsd(cls, root, query):
         myself = cls()
         myself.type = "task"
         return myself
-
+
+
 class XCTASK_CC(sc.ClassDict):
     """mGGA TASK+CC xc-functional."""
 
@@ -252,6 +254,7 @@ def fromhsd(cls, root, query):
         myself.type = "task+cc"
         return myself
 
+
 class XCLocal(sc.ClassDict):
     '''Local xc-functional.'''
 
@@ -298,20 +301,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,
-    "tpss": XCTPSS,
-    "scan": XCSCAN,
-    "r2scan": XCR2SCAN,
-    "r4scan": XCR4SCAN,
-    "task": XCTASK,
-    "task+cc": XCTASK_CC
+    '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
 }

sktwocnt/lib/twocnt.f90

@@ -22,11 +22,11 @@ module twocnt
   use, intrinsic :: iso_c_binding, only : c_size_t
 
   use xc_f03_lib_m, only : xc_f03_func_t, xc_f03_func_init, xc_f03_func_end, xc_f03_lda_vxc,&
-  & xc_f03_gga_vxc, xc_f03_mgga_vxc, XC_LDA_X, XC_LDA_X_YUKAWA, XC_LDA_C_PW, XC_GGA_X_PBE,&
-  & XC_GGA_C_PBE, XC_GGA_X_B88, XC_GGA_C_LYP, XC_GGA_X_SFAT_PBE, XC_MGGA_X_TPSS, XC_MGGA_C_TPSS,&
-  & XC_MGGA_X_SCAN, XC_MGGA_C_SCAN,  XC_MGGA_X_R2SCAN, XC_MGGA_C_R2SCAN, XC_MGGA_X_R4SCAN,&
-  & XC_MGGA_X_TASK, XC_HYB_GGA_XC_B3LYP, XC_HYB_GGA_XC_CAMY_B3LYP, XC_MGGA_C_CC, &
-  & XC_UNPOLARIZED, xc_f03_func_set_ext_params
+      & xc_f03_gga_vxc, xc_f03_mgga_vxc, XC_LDA_X, XC_LDA_X_YUKAWA, XC_LDA_C_PW, XC_GGA_X_PBE,&
+      & XC_GGA_C_PBE, XC_GGA_X_B88, XC_GGA_C_LYP, XC_GGA_X_SFAT_PBE, XC_MGGA_X_TPSS,&
+      & XC_MGGA_C_TPSS, XC_MGGA_X_SCAN, XC_MGGA_C_SCAN,  XC_MGGA_X_R2SCAN, XC_MGGA_C_R2SCAN,&
+      & XC_MGGA_X_R4SCAN, XC_MGGA_X_TASK, XC_HYB_GGA_XC_B3LYP, XC_HYB_GGA_XC_CAMY_B3LYP,&
+      & XC_MGGA_C_CC,  XC_UNPOLARIZED, xc_f03_func_set_ext_params
 
   implicit none
   private
@@ -62,7 +62,7 @@ module twocnt
     !> atomic potential on grid
     type(TGridorb2) :: pot
 
-    !> atomic density and 1st/2nd derivative on grid
+    !> atomic density, 1st/2nd derivative and kinetic energy density on grid
     type(TGridorb2) :: rho, drho, ddrho, tau
 
   end type TAtomdata
@@ -505,7 +505,7 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
     !! radial grid-orbital portion and 1st/2nd derivative for all basis functions of atom 2
     real(dp), allocatable :: radval2(:,:), radval2p(:,:), radval2pp(:,:)
 
-    !! total potential and electron density of two atoms
+    !! total potential, electron density and kinetic energy density of two atoms
     real(dp), allocatable :: potval(:), densval(:), tauval(:), taupotval(:)
 
     !! atomic 1st density derivatives of atom 1
@@ -526,7 +526,7 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
     !! number of integration points
     integer :: nGrid
 
-    !> number of density grid points, compatible with libxc signature
+    !! number of density grid points, compatible with libxc signature
     integer(c_size_t) :: nGridLibxc
 
     !!  orbital indices/angular momenta on the two atoms and interaction type
@@ -540,8 +540,8 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
 
     !! libxc related objects
     real(dp), allocatable :: vxc(:), vx(:), vx_sr(:), vc(:)
-    real(dp), allocatable :: rhor(:), sigma(:), tau(:), vxcsigma(:), vxsigma(:), vxsigma_sr(:), vcsigma(:), vxtau(:), vctau(:)
-    real(dp), allocatable :: divvxc(:), divvx(:), divvc(:)
+    real(dp), allocatable :: rhor(:), sigma(:), tau(:), vxcsigma(:), vxsigma(:), vxsigma_sr(:)
+    real(dp), allocatable :: vcsigma(:), vxtau(:), vctau(:), divvxc(:), divvx(:), divvc(:)
 
     r1 => grid1(:, 1)
     theta1 => grid1(:, 2)
@@ -591,19 +591,18 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
         densval2p = atom2%drho%getValue(r2)
         ! care about correct 4pi normalization of density and compute sigma
         sigma = getLibxcSigma(densval1p, densval2p, dots)
-        if (xcFunctional%isMGGA(iXC)) then
-          allocate(tauval(nGrid))
-          tauval(:) = atom1%tau%getValue(r1) + atom2%tau%getValue(r2)
-          tau = getLibxcTau(tauval)
-
-          allocate(vxtau(nGrid))
-          vxtau(:) = 0.0_dp
-          if (iXC /= xcFunctional%MGGA_TASK) then
-             allocate(vctau(nGrid))
-             vctau(:) = 0.0_dp
-          end if
-       end if
       end if
+
+      if (xcFunctional%isMGGA(iXC)) then
+        allocate(tauval(nGrid))
+        tauval(:) = atom1%tau%getValue(r1) + atom2%tau%getValue(r2)
+        tau = getLibxcTau(tauval)
+        allocate(vxtau(nGrid), source=0.0_dp)
+        if (iXC /= xcFunctional%MGGA_TASK) then
+          allocate(vctau(nGrid), source=0.0_dp)
+        end if
+      end if
+
       if (tGlobalHybrid .or. tCam) then
         allocate(vxc(nGrid))
         vxc(:) = 0.0_dp
@@ -632,19 +631,20 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
         call getDivergence(nRad, nAng, densval1p, densval2p, r1, r2, theta1, theta2, vxsigma, divvx)
         call getDivergence(nRad, nAng, densval1p, densval2p, r1, r2, theta1, theta2, vcsigma, divvc)
         potval = vx + vc + divvx + divvc
-      ! 10: MGGA-TPSS, 11: MGGA-SCAN, 12: MGGA-r2SCAN, 13: MGGA-r4SCAN,
-      ! 14: MGGA-TASK, 15: MGGA-TASK+CC 
+      ! 10: MGGA-TPSS, 11: MGGA-SCAN, 12: MGGA-r2SCAN, 13: MGGA-r4SCAN 14: MGGA-TASK,
+      ! 15: MGGA-TASK+CC
       case(10:15)
         ! MGGA exchange
-        call xc_f03_mgga_vxc(xcfunc_x, nGridLibxc, rhor(1), sigma(1), lapl(1), tau(1),&
-        & vx(1), vxsigma(1), vlapl(1), vxtau(1))
+        call xc_f03_mgga_vxc(xcfunc_x, nGridLibxc, rhor(1), sigma(1), lapl(1), tau(1), vx(1),&
+            & vxsigma(1), vlapl(1), vxtau(1))
         call getDivergence(nRad, nAng, densval1p, densval2p, r1, r2, theta1, theta2, vxsigma, divvx)
         ! 11: MGGA-TPSS, 12: MGGA-SCAN, 13: MGGA-r2SCAN, 14: MGGA-TASK
         if (iXC /= xcFunctional%MGGA_TASK) then
           ! MGGA correlation
-          call xc_f03_mgga_vxc(xcfunc_c, nGridLibxc, rhor(1), sigma(1), lapl(1), tau(1),&
-          & vc(1), vcsigma(1), vlapl(1), vctau(1))
-          call getDivergence(nRad, nAng, densval1p, densval2p, r1, r2, theta1, theta2, vcsigma, divvc)
+          call xc_f03_mgga_vxc(xcfunc_c, nGridLibxc, rhor(1), sigma(1), lapl(1), tau(1), vc(1),&
+              & vcsigma(1), vlapl(1), vctau(1))
+          call getDivergence(nRad, nAng, densval1p, densval2p, r1, r2, theta1, theta2, vcsigma,&
+              & divvc)
           potval = vx + vc + divvx + divvc
           taupotval = vxtau + vctau
         else
@@ -714,7 +714,7 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
       ! calculate SK-quantities
       ! Hamiltonian
       integ1 = getHamiltonian(radval1(:, i1), radval2(:, i2), radval2p(:, i2), radval2pp(:, i2),&
-          & r2, l2, spherval1, spherval2, potval, taupotval, weights)
+          & r2, l2, spherval1, spherval2, potval, weights, pot_tau=taupotval)
       ! overlap integral: \sum_{r,\Omega} R_1(r) Y_1(\Omega) R_2(r) Y_2(\Omega) weight
       integ2 = getOverlap(radval1(:, i1), radval2(:, i2), spherval1, spherval2, weights)
       ! total density: \int (|\phi_1|^2 + |\phi_2|^2)
@@ -795,17 +795,18 @@ pure function getLibxcSigma(drho1, drho2, dots) result(sigma)
 
   end function getLibxcSigma
 
-  !> Calculates libXC renormalized density superposition of dimer.
-  pure function getLibxcTau(tau) result(tau_libxc)
 
-    !> superposition of atomic densities of atom 1 and atom 2
+  !> Calculates libXC renormalized kinetic energy density superposition of dimer.
+  pure function getLibxcTau(tau) result(rtau)
+
+    !> superposition of atomic kinetic energy densities of atom 1 and atom 2
     real(dp), intent(in) :: tau(:)
 
-    !> renormalized density
-    real(dp), allocatable :: tau_libxc(:)
+    !> renormalized kinetic energy density
+    real(dp), allocatable :: rtau(:)
 
-    ! renorm rho (incoming quantities are 4pi normed)
-    tau_libxc = tau * rec4pi
+    ! renorm tau (incoming quantities are 4pi normed)
+    rtau = tau * rec4pi
 
   end function getLibxcTau
 
@@ -969,8 +970,8 @@ end function getdensity
 
 
   !> Calculates Hamiltonian for a fixed orbital and interaction configuration.
-  pure function getHamiltonian(rad1, rad2, rad2p, rad2pp, r2, l2, spher1,&
-  & spher2, pot, pot_tau, weights) result(res)
+  pure function getHamiltonian(rad1, rad2, rad2p, rad2pp, r2, l2, spher1, spher2, pot, weights,&
+      & pot_tau) result(res)
 
     !> radial grid-orbital portion of atom 1 and atom 2
     real(dp), intent(in) :: rad1(:), rad2(:)
@@ -990,16 +991,16 @@ pure function getHamiltonian(rad1, rad2, rad2p, rad2pp, r2, l2, spher1,&
     !> effective potential on grid
     real(dp), intent(in) :: pot(:)
 
-    !> effective potential on grid
-    real(dp), intent(in), allocatable :: pot_tau(:)
-
     !> integration weights
     real(dp), intent(in) :: weights(:)
 
+    !> kinetic energy density on grid
+    real(dp), intent(in), optional :: pot_tau(:)
+
     !! resulting Hamiltonian matrix element
     real(dp) :: res
 
-    if (allocated(pot_tau)) then
+    if (present(pot_tau)) then
       res = sum((rad1 * spher1)&
           & * ((- 0.5_dp * rad2pp&
           & - rad2p / r2&

sktwocnt/lib/xcfunctionals.f90

@@ -109,26 +109,29 @@ pure function TXcFunctionalsEnum_isGGA(this, xcnr) result(isGGA)
 
   end function TXcFunctionalsEnum_isGGA
 
+
   pure function TXcFunctionalsEnum_isMGGA(this, xcnr) result(isMGGA)
 
-  !> Class instance
-  class(TXcFunctionalsEnum), intent(in) :: this
+    !> Class instance
+    class(TXcFunctionalsEnum), intent(in) :: this
 
-  !> identifier of exchange-correlation type
-  integer, intent(in) :: xcnr
+    !> identifier of exchange-correlation type
+    integer, intent(in) :: xcnr
 
-  !> True, if xc-functional index corresponds to a GGA functional
-  logical :: isMGGA
+    !> True, if xc-functional index corresponds to a GGA functional
+    logical :: isMGGA
 
-  isMGGA = .false.
+    isMGGA = .false.
 
-  if (xcnr == this%MGGA_TPSS .or. xcnr == this%MGGA_SCAN&
-  & .or. xcnr == this%MGGA_r2SCAN .or. xcnr == this%MGGA_r4SCAN&
-  & .or. xcnr == this%MGGA_TASK .or. xcnr == this%MGGA_TASK_CC)&
-  & isMGGA = .true.
+    if (xcnr == this%MGGA_TPSS .or. xcnr == this%MGGA_SCAN .or. xcnr == this%MGGA_r2SCAN&
+        & .or. xcnr == this%MGGA_r4SCAN .or. xcnr == this%MGGA_TASK&
+        & .or. xcnr == this%MGGA_TASK_CC) then
+      isMGGA = .true.
+    end if
 
   end function TXcFunctionalsEnum_isMGGA
 
+
   pure function TXcFunctionalsEnum_isLongRangeCorrected(this, xcnr) result(isLongRangeCorrected)
 
     !> Class instance