add new SCF diagnostic quantities for slateratom

.gitignore

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

slateratom/lib/CMakeLists.txt

@@ -31,6 +31,7 @@ set(sources-fpp
   blasroutines.F90
   broydenmixer.F90
   confinement.F90
+  diismixer.F90
   input.F90
   lapackroutines.F90
   simplemixer.F90

slateratom/lib/broydenmixer.F90

@@ -48,11 +48,11 @@ module broydenmixer
     !> Weights for prev. iterations
     real(dp), allocatable :: ww(:)
 
-    !> Charge difference in last iteration
-    real(dp), allocatable :: qDiffLast(:)
+    !> Input increment in last iteration
+    real(dp), allocatable :: incrementLast(:)
 
-    !> Input charge in last iteration
-    real(dp), allocatable :: qInpLast(:)
+    !> Input quantity that is being mixed, last iteration
+    real(dp), allocatable :: inputLast(:)
 
     !> Storage for the "a" matrix
     real(dp), allocatable :: aa(:,:)
@@ -70,7 +70,7 @@ module broydenmixer
 
   !> Creates a Broyden mixer instance.
   !! The weight associated with an iteration is calculated as weigthFac/ww where ww is the Euclidean
-  !! norm of the charge difference vector. If the calculated weigth is outside of the
+  !! norm of the quantity difference vector. If the calculated weigth is outside of the
   !! [minWeight, maxWeight] region it is replaced with the appropriate boundary value.
   subroutine TBroydenMixer_init(this, mIter, mixParam, omega0, minWeight, maxWeight, weightFac)
 
@@ -107,8 +107,8 @@ subroutine TBroydenMixer_init(this, mIter, mixParam, omega0, minWeight, maxWeigh
     this%maxWeight = maxWeight
     this%weightFac = weightFac
     allocate(this%ww(mIter-1))
-    allocate(this%qInpLast(this%nElem))
-    allocate(this%qDiffLast(this%nElem))
+    allocate(this%inputLast(this%nElem))
+    allocate(this%incrementLast(this%nElem))
     allocate(this%aa(mIter-1, mIter-1))
     allocate(this%dF(this%nElem, mIter - 1))
     allocate(this%uu(this%nElem, mIter - 1))
@@ -129,10 +129,10 @@ subroutine TBroydenMixer_reset(this, nElem)
 
     if (nElem /= this%nElem) then
       this%nElem = nElem
-      deallocate(this%qInpLast)
-      deallocate(this%qDiffLast)
-      allocate(this%qInpLast(this%nElem))
-      allocate(this%qDiffLast(this%nElem))
+      deallocate(this%inputLast)
+      deallocate(this%incrementLast)
+      allocate(this%inputLast(this%nElem))
+      allocate(this%incrementLast(this%nElem))
       deallocate(this%dF)
       allocate(this%dF(this%nElem, this%mIter - 1))
       deallocate(this%uu)
@@ -153,24 +153,24 @@ end subroutine TBroydenMixer_reset
   !!   The restriction arises from the assumption that the dot-products of density matrices are
   !!   real-valued (imaginary parts add up to zero due to the hermitian property) and the linear
   !!   system of equations remains real-valued.
-  subroutine TBroydenMixer_mix(this, qInpResult, qDiff)
+  subroutine TBroydenMixer_mix(this, inputResult, increment)
 
     !> The Broyden mixer
     type(TBroydenMixer), intent(inout) :: this
 
     !> Input charges on entry, mixed charges on exit
-    real(dp), intent(inout) :: qInpResult(:)
+    real(dp), intent(inout) :: inputResult(:)
 
     !> Charge difference between output and input charges
-    real(dp), intent(in) :: qDiff(:)
+    real(dp), intent(in) :: increment(:)
 
     this%iIter = this%iIter + 1
     if (this%iIter > this%mIter) then
       error stop "Broyden mixer: Maximal nr. of steps exceeded"
     end if
 
-    call modifiedBroydenMixing(qInpResult, this%qInpLast, this%qDiffLast, this%aa,&
-        & this%ww, this%iIter, qDiff, this%alpha, this%omega0, this%minWeight, this%maxWeight,&
+    call modifiedBroydenMixing(inputResult, this%inputLast, this%incrementLast, this%aa,&
+        & this%ww, this%iIter, increment, this%alpha, this%omega0, this%minWeight, this%maxWeight,&
         & this%weightFac, this%nElem, this%dF, this%uu)
 
   end subroutine TBroydenMixer_mix

slateratom/lib/diagonalizations.f90

@@ -14,7 +14,7 @@ module diagonalizations
 
   !> Diagonalizes overlap matrix to check for linear dependency of basis set.
   !! Implicitely LAPACK's dsyev is called.
-  subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
+  subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss, invsqrt_ss)
 
     !> maximum angular momentum
     integer, intent(in) :: max_l
@@ -28,14 +28,17 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
     !> overlap supervector
     real(dp), intent(in) :: ss(0:, :,:)
 
+    !> inv. sqrt. of the overlap supervector
+    real(dp), intent(inout) :: invsqrt_ss(0:, :,:)
+
     !! overlap matrices
     real(dp), allocatable :: overlap(:,:)
 
     !! eigenvalues of overlap matrices
     real(dp), allocatable :: eigenvalues(:)
 
     !> auxiliary variables
-    integer :: ll, diagsize
+    integer :: ll, diagsize, ii
 
     do ll = 0, max_l
 
@@ -46,7 +49,7 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
 
       overlap = ss(ll, :,:)
 
-      call heev(overlap, eigenvalues, 'U', 'N')
+      call heev(overlap, eigenvalues, 'U', 'V')
 
       write(*, '(A,I3,A,E16.8)') 'Smallest eigenvalue of overlap for l= ', ll, ' : ', eigenvalues(1)
 
@@ -57,6 +60,12 @@ subroutine diagonalize_overlap(max_l, num_alpha, poly_order, ss)
         stop
       end if
 
+      ! compute S^(-1/2)
+      do ii = 1, diagsize
+        invsqrt_ss(ll, ii, ii) = eigenvalues(ii)**(-1.0_dp / 2)
+      end do
+      invsqrt_ss(ll,:,:) = matmul(overlap, matmul(invsqrt_ss(ll,:,:), transpose(overlap)))
+
       deallocate(overlap, eigenvalues)
 
     end do

slateratom/lib/diismixer.F90

@@ -0,0 +1,219 @@
+#:include 'common.fypp'
+
+!> Contains the DIIS mixer
+!! The DIIS mixing is done by building a weighted combination over the previous input charges to
+!! minimise the residue of the error.
+!! Only a specified number of previous charge vectors are considered.
+!! The modification based on from Kovalenko et al. (J. Comput. Chem., 20: 928-936 1999) and Patrick
+!! Briddon to add a contribution from the gradient vector as well is also used.
+!! In order to use the mixer you have to create and reset it.
+!!
+!! Code is adapted from DFTB+.
+module diismixer
+  use common_accuracy, only : dp
+  use lapackroutines, only : gesv
+  implicit none
+
+  private
+
+  public :: TDiisMixer, TDiisMixer_init, TDiisMixer_mix, TDiisMixer_reset
+
+  !> Contains the necessary data for an DIIS mixer.
+  type TDiisMixer
+    private
+
+    !> Initial mixing parameter
+    real(dp) :: initMixParam
+
+    !> Max. nr. of stored prev. vectors
+    integer :: mPrevVector
+
+    !> Nr. of stored previous vectors
+    integer :: iPrevVector
+
+    !> Nr. of elements in the vectors
+    integer :: nElem
+
+    !> Index for the storage
+    integer :: indx
+
+    !> Stored previous input quantities
+    real(dp), allocatable :: prevInput(:,:)
+
+    !> Stored differences of previous input quantities
+    real(dp), allocatable :: prevIncrement(:,:)
+
+    !> Stored prev. error vectors
+    real(dp), allocatable :: prevErrorVec(:,:)
+
+    !> True if DIIS used from iteration 2 as well as mixing
+    logical :: tFromStart
+
+    !> Alpha factor to add in new information
+    real(dp) :: alpha
+
+    contains
+      procedure :: reset => TDiisMixer_reset
+      procedure :: mix1D => TDiisMixer_mix
+  end type TDiisMixer
+
+contains
+
+  !> Initializes a DIIS mixer instance.
+  subroutine TDiisMixer_init(this, iGenerations, initMixParam, tFromStart)
+
+    !> Pointer to an initialized DIIS mixer on exit
+    type(TDiisMixer), intent(out) :: this
+
+    !> Number of generations to consider (including current)
+    integer, intent(in) :: iGenerations
+
+    !> Damping parameter for the first mixing steps
+    real(dp), intent(in) :: initMixParam
+
+    !> Use DIIS from step 2 onwards?
+    logical, intent(in) :: tFromStart
+
+    @:ASSERT(iGenerations >= 2)
+
+    this%nElem = 0
+    this%mPrevVector = iGenerations
+
+    allocate(this%prevInput(this%nElem, this%mPrevVector))
+    allocate(this%prevErrorVec(this%nElem, this%mPrevVector))
+    allocate(this%prevIncrement(this%nElem, this%mPrevVector))
+
+    this%initMixParam = initMixParam
+    this%tFromStart = tFromStart
+
+  end subroutine TDiisMixer_init
+
+
+  !> Makes the mixer ready for a new SCC cycle.
+  subroutine TDiisMixer_reset(this, nElem)
+
+    !> DIIS mixer instance
+    class(TDiisMixer), intent(inout) :: this
+
+    !> Nr. of elements in the vectors to mix
+    integer, intent(in) :: nElem
+
+    @:ASSERT(nElem > 0)
+
+    if (nElem /= this%nElem) then
+      this%nElem = nElem
+      deallocate(this%prevInput)
+      deallocate(this%prevErrorVec)
+      deallocate(this%prevIncrement)
+      allocate(this%prevInput(this%nElem, this%mPrevVector))
+      allocate(this%prevErrorVec(this%nElem, this%mPrevVector))
+      allocate(this%prevIncrement(this%nElem, this%mPrevVector))
+    end if
+    this%iPrevVector = 0
+    this%indx = 0
+
+  end subroutine TDiisMixer_reset
+
+
+  !> Mixes quantities according to the DIIS method.
+  subroutine TDiisMixer_mix(this, inputResult, inputIncrement, errorVector)
+
+    !> Pointer to the diis mixer
+    class(TDiisMixer), intent(inout) :: this
+
+    !> Input quantity on entry, mixed quantity on exit.
+    real(dp), intent(inout) :: inputResult(:)
+
+    !> Increment vector between input and output quantities
+    real(dp), intent(in) :: inputIncrement(:)
+
+    !> Error metric vector
+    real(dp), intent(in) :: errorVector(:)
+
+
+    real(dp), allocatable :: aa(:,:), bb(:,:)
+    integer :: ii, jj
+
+    @:ASSERT(size(inputResult) == this%nElem)
+    @:ASSERT(size(errorVector) == this%nElem)
+
+    if (this%iPrevVector < this%mPrevVector) then
+      this%iPrevVector = this%iPrevVector + 1
+    end if
+
+    call storeVectors(this%prevInput, this%prevIncrement, this%prevErrorVec, this%indx,&
+        & inputResult, inputIncrement, errorVector, this%mPrevVector)
+    if (this%tFromStart .or. this%iPrevVector == this%mPrevVector) then
+
+      allocate(aa(this%iPrevVector + 1, this%iPrevVector + 1))
+      allocate(bb(this%iPrevVector + 1, 1))
+
+      aa(:,:) = 0.0_dp
+      bb(:,:) = 0.0_dp
+
+      ! (due to the hermitian property of our density matrices, the dot-product below is real)
+      do ii = 1, this%iPrevVector
+        do jj = 1, this%iPrevVector
+          aa(ii, jj) = dot_product(this%prevErrorVec(:, ii), this%prevErrorVec(:, jj))
+        end do
+      end do
+      aa(this%iPrevVector + 1, 1:this%iPrevVector) = -1.0_dp
+      aa(1:this%iPrevVector, this%iPrevVector + 1) = -1.0_dp
+
+      bb(this%iPrevVector + 1, 1) = -1.0_dp
+
+      ! Solve DIIS system of linear equations
+      call gesv(aa, bb)
+
+      inputResult(:) = 0.0_dp
+      do ii = 1, this%iPrevVector
+        inputResult(:) = inputResult + bb(ii, 1) * (this%prevInput(:, ii) + this%prevIncrement(:, ii))
+      end do
+
+    end if
+
+    if (this%iPrevVector < this%mPrevVector) then
+      ! First few iterations return simple mixed vector
+      inputResult(:) = inputResult + this%initMixParam * inputIncrement(:)
+    end if
+
+  end subroutine TDiisMixer_mix
+
+
+  !> Stores a vector pair in a limited storage.
+  !! If the stack is full, oldest vector pair is overwritten.
+  subroutine storeVectors(prevInp, prevIncrement, prevErrorVector, indx, input, increment,&
+        & errorVector, mPrevVector)
+
+    !> Contains previous vectors of the first type
+    real(dp), intent(inout) :: prevInp(:,:)
+
+    !> Contains previous vectors of the second type
+    real(dp), intent(inout) :: prevErrorVector(:,:)
+
+    !> Contains previous differences of vectors of first type
+    real(dp), intent(inout) :: prevIncrement(:,:)
+
+    !> Indexing of data
+    integer, intent(inout) :: indx
+
+    !> New first vector
+    real(dp), intent(in) :: input(:)
+
+    !> New increment of first vector
+    real(dp), intent(in) :: increment(:)
+
+    !> New second vector
+    real(dp), intent(in) :: errorVector(:)
+
+    !> Size of the stacks.
+    integer, intent(in) :: mPrevVector
+
+    indx = mod(indx, mPrevVector) + 1
+    prevInp(:, indx) = input
+    prevIncrement(:, indx) = increment
+    prevErrorVector(:, indx) = errorVector
+
+  end subroutine storeVectors
+
+end module diismixer