implement SAP guess into slateratom

.gitignore

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

sktools/src/sktools/calculators/slateratom.py

@@ -21,6 +21,9 @@
                          'lcy-bnl' : 6, 'pbe0' : 7, 'b3lyp' : 8,
                          'camy-b3lyp' : 9, 'camy-pbeh' : 10}
 
+GUESS_TYPES = {'Thomas-Fermi': 1,
+               'SAP': 2}
+
 INPUT_FILE = "slateratom.in"
 STDOUT_FILE = "output"
 DEFAULT_BINARY = "slateratom"
@@ -64,11 +67,12 @@ class SlaterAtomSettings(sc.ClassDict):
         Maximal power for every angular momentum.
     """
 
-    def __init__(self, exponents, maxpowers, scftol, maxscfiter):
+    def __init__(self, exponents, maxpowers, scfguess, scftol, maxscfiter):
         super().__init__()
         self.exponents = exponents
         self.maxpowers = maxpowers
         self.scftol = scftol
+        self.scfguess = scfguess
         self.maxscfiter = maxscfiter
 
     @classmethod
@@ -77,11 +81,13 @@ def fromhsd(cls, root, query):
         exponents = sc.get_shellvalues_list(node, query, conv.float1)
         node = query.getchild(root, "maxpowers")
         maxpowers = sc.get_shellvalues_list(node, query, conv.int0)
+        scfguess = query.getvalue(
+            root, "scfguess", defvalue="Thomas-Fermi")
         scftol = query.getvalue(
             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)
+        return cls(exponents, maxpowers, scfguess, scftol, maxscfiter)
 
     def __eq__(self, other):
         if not isinstance(other, SlaterAtomSettings):
@@ -100,6 +106,8 @@ def __eq__(self, other):
             for ii in range(len(myexps)):
                 if abs(myexps[ii] - otherexps[ii]) > sc.INPUT_FLOAT_TOLERANCE:
                     return False
+        if (self.scfguess != other.scfguess):
+            return False
         if (abs(self.scftol - other.scftol)
                 > sc.INPUT_FLOAT_TOLERANCE):
             return False
@@ -159,6 +167,10 @@ def __init__(self, settings, atomconfig, functional, compressions):
         if len(settings.maxpowers) != atomconfig.maxang + 1:
             msg = "Slateratom: Missing STO max. powers for some shells"
             raise sc.SkgenException(msg)
+
+        if self._settings.scfguess not in GUESS_TYPES.keys():
+            msg = "SlaterAtom: unsupported guess (must be either Thomas-Fermi or SAP)"
+            raise sc.SkgenException(msg)
 
         if self._settings.scftol <= 0.0:
             msg = "Slateratom: SCF tolerance must be >0.0 a.u."
@@ -365,6 +377,8 @@ def write(self, workdir):
 
         out.append("{:s} \t\t{:s} write eigenvectors".format(
             self._LOGICALSTRS[False], self._COMMENT))
+        out.append("{} \t\t\t{:s} SCF guess".format(GUESS_TYPES[self._settings.scfguess],
+                                                     self._COMMENT))
         out.append("{} {:g} \t\t\t{:s} broyden mixer, mixing factor".format(
             2, 0.1, self._COMMENT))
 

slateratom/lib/CMakeLists.txt

@@ -23,6 +23,7 @@ set(sources-f90
   mixer.f90
   numerical_differentiation.f90
   output.f90
+  sap.f90
   total_energy.f90
   utilities.f90
   zora_routines.f90)

slateratom/lib/dft.f90

@@ -13,7 +13,7 @@ module dft
   implicit none
   private
 
-  public :: dft_start_pot, density_grid, dft_exc_energy, dft_vxc_energy
+  public :: thomas_fermi_start_pot, density_grid, dft_exc_energy, dft_vxc_energy
   public :: dft_exc_matrixelement, xalpha
   public :: check_accuracy
 
@@ -25,7 +25,7 @@ module dft
   !!
   !! Generalized Thomas-Fermi atomic potential as published by R. Latter,
   !! Phys. Rev. 99, 510 (1955) eqn. 5/9 and implemented in Dirk Porezags scfatom.
-  pure subroutine dft_start_pot(abcissa, num_mesh_points, nuc, vxc)
+  pure subroutine thomas_fermi_start_pot(abcissa, num_mesh_points, nuc, vxc)
 
     !> numerical integration abcissas
     real(dp), intent(in) :: abcissa(:)
@@ -59,7 +59,7 @@ pure subroutine dft_start_pot(abcissa, num_mesh_points, nuc, vxc)
 
     end do
 
-  end subroutine dft_start_pot
+  end subroutine thomas_fermi_start_pot
 
 
   !> Calculate and store density and density derivatives on radial grid.

slateratom/lib/globals.f90

@@ -178,6 +178,9 @@ module globals
   !> identifier of mixer
   integer :: mixnr
 
+  !> SCF guess ID
+  integer :: scfGuess
+
   !> mixer instance
   type(TMixer), allocatable :: pMixer
 

slateratom/lib/hamiltonian.f90

@@ -17,7 +17,7 @@ module hamiltonian
 contains
 
   !> Main driver routine for Fock matrix build-up. Also calls mixer with potential matrix.
-  subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l,&
+  subroutine build_hamiltonian(pMixer, iScf, scfGuess, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l,&
       & num_alpha, poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha,&
       & pot_old, pot_new, tZora, ff, camAlpha, camBeta)
 
@@ -27,6 +27,9 @@ subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp
     !> current SCF step
     integer, intent(in) :: iScf
 
+    !> SCF guess ID
+    integer, intent(in) :: scfGuess
+
     !> kinetic supervector
     real(dp), intent(in) :: tt(0:,:,:)
 
@@ -176,8 +179,13 @@ subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp
     end if
 
     ! build mixer input
-    pot_new(1, :,:,:) = - real(nuc, dp) * uu + j_matrix - k_matrix(1, :,:,:)
-    pot_new(2, :,:,:) = - real(nuc, dp) * uu + j_matrix - k_matrix(2, :,:,:)
+    if (iScf == 0 .and. scfGuess == 2) then
+      pot_new(1, :,:,:) = -k_matrix(1, :,:,:)
+      pot_new(2, :,:,:) = -k_matrix(2, :,:,:)
+    else
+      pot_new(1, :,:,:) = -real(nuc, dp) * uu + j_matrix - k_matrix(1, :,:,:)
+      pot_new(2, :,:,:) = -real(nuc, dp) * uu + j_matrix - k_matrix(2, :,:,:)
+    end if
 
     ! mixer
     allocate(pot_diff, mold=pot_old)

slateratom/lib/input.F90

@@ -20,8 +20,8 @@ module input
   !> Reads in all properties, except for occupation numbers.
   subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min_alpha,&
       & max_alpha, num_alpha, tAutoAlphas, alpha, conf_type, confInp, num_occ, num_power,&
-      & num_alphas, xcnr, tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega,&
-      & camAlpha, camBeta, grid_params)
+      & num_alphas, xcnr, tPrintEigvecs, tZora, mixnr, mixing_factor, scfGuess, xalpha_const,&
+      & omega,camAlpha, camBeta, grid_params)
 
     !> nuclear charge, i.e. atomic number
     integer, intent(out) :: nuc
@@ -85,6 +85,9 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
 
     !> mixing factor
     real(dp), intent(out) :: mixing_factor
+
+    !> SCF guess
+    integer, intent(out) :: scfGuess
 
     !> exchange parameter for X-Alpha exchange
     real(dp), intent(out) :: xalpha_const
@@ -251,6 +254,10 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     write(*, '(A)') 'Print Eigenvectors ? .true./.false.'
     read(*,*) tPrintEigvecs
 
+    write(*, '(A)') 'Enter SCF guess:&
+        & 1: Thomas-Fermi, 2: superposition of atomic potentials (SAP)'
+    read(*,*) scfGuess
+
     write(*, '(A)') 'Enter mixer and mixing parameter <1:&
         & 1: Simple mixer, 2: Broyden mixer'
     read(*,*) mixnr, mixing_factor