Reimplement purification-based density matrix extrapolation

* The new code seems to work
* Reference:
  * The Journal of Chemical Physics 123, 044107 (2005)
  * The Journal of Chemical Physics 132, 124104 (2010)

Source/Fortran/GeometryOptimizationModule.F90

@@ -8,9 +8,9 @@ MODULE GeometryOptimizationModule
   USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, &
        & DestructMatrixMemoryPool
   USE PSMatrixAlgebraModule, ONLY : MatrixMultiply, MatrixNorm, &
-       & IncrementMatrix, MatrixTrace, ScaleMatrix
+       & IncrementMatrix, ScaleMatrix, DotMatrix
   USE PSMatrixModule, ONLY : Matrix_ps, DestructMatrix, ConstructEmptyMatrix, &
-       & FillMatrixIdentity, PrintMatrixInformation, CopyMatrix
+       & PrintMatrixInformation, CopyMatrix
   USE SolverParametersModule, ONLY : SolverParameters_t, PrintParameters, &
        & DestructSolverParameters
   USE SquareRootSolversModule, ONLY : SquareRoot, InverseSquareRoot
@@ -40,18 +40,13 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, nel, NewDensity,&
     !! Local Matrices
     TYPE(Matrix_ps) :: WorkingDensity
     TYPE(Matrix_ps) :: WorkingOverlap
-    TYPE(Matrix_ps) :: AddBranch, SubtractBranch
-    TYPE(Matrix_ps) :: TempMat1, TempMat2
-    TYPE(Matrix_ps) :: Identity
+    TYPE(Matrix_ps) :: TempMat
     !! Local Variables
-    REAL(NTREAL) :: subtract_trace
-    REAL(NTREAL) :: add_trace
     REAL(NTREAL) :: trace_value
     REAL(NTREAL) :: norm_value
     !! Temporary Variables
     TYPE(MatrixMemoryPool_p) :: pool1
     INTEGER :: outer_counter
-    INTEGER :: total_iterations
 
     !! Optional Parameters
     IF (PRESENT(solver_parameters_in)) THEN
@@ -72,8 +67,6 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, nel, NewDensity,&
     CALL ConstructEmptyMatrix(NewDensity, PreviousDensity)
     CALL ConstructEmptyMatrix(WorkingDensity, PreviousDensity)
     CALL ConstructEmptyMatrix(WorkingOverlap, PreviousDensity)
-    CALL ConstructEmptyMatrix(Identity, PreviousDensity)
-    CALL FillMatrixIdentity(Identity)
 
     !! Compute the working hamiltonian.
     CALL CopyMatrix(PreviousDensity, WorkingDensity)
@@ -85,83 +78,52 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, nel, NewDensity,&
             & solver_parameters%BalancePermutation, memorypool_in=pool1)
        CALL PermuteMatrix(WorkingOverlap, WorkingOverlap, &
             & solver_parameters%BalancePermutation, memorypool_in=pool1)
-       CALL PermuteMatrix(Identity, Identity, &
-            & solver_parameters%BalancePermutation, memorypool_in=pool1)
     END IF
 
-    !! Finish Setup
-    CALL CopyMatrix(WorkingDensity, NewDensity)
-    CALL CopyMatrix(WorkingDensity, AddBranch)
-    CALL CopyMatrix(WorkingDensity, SubtractBranch)
-
     !! Iterate
     IF (solver_parameters%be_verbose) THEN
        CALL WriteHeader("Iterations")
        CALL EnterSubLog
     END IF
     outer_counter = 1
-    norm_value = solver_parameters%converge_diff + 1.0_NTREAL
     DO outer_counter = 1,solver_parameters%max_iterations
-       !! Figure Out Sigma Value. After which, XnS is stored in TempMat
-       IF (outer_counter .GT. 1) THEN
-          CALL MatrixMultiply(AddBranch, WorkingOverlap, TempMat1, &
-               & threshold_in=solver_parameters%threshold, memory_pool_in=pool1)
-          CALL MatrixTrace(TempMat1, add_trace)
-          CALL MatrixMultiply(SubtractBranch, WorkingOverlap, TempMat2, &
-               & threshold_in=solver_parameters%threshold, memory_pool_in=pool1)
-          CALL MatrixTrace(TempMat2, subtract_trace)
-          IF (ABS(nel - add_trace) .GT. ABS(nel - subtract_trace)) THEN
-             !! Subtract Branch
-             trace_value = subtract_trace
-             CALL IncrementMatrix(AddBranch, NewDensity, -1.0_NTREAL)
-             norm_value = MatrixNorm(NewDensity)
-             CALL CopyMatrix(AddBranch, NewDensity)
-             CALL CopyMatrix(TempMat2, TempMat1)
-          ELSE
-             !! Add Branch
-             trace_value = add_trace
-             CALL IncrementMatrix(SubtractBranch, NewDensity, -1.0_NTREAL)
-             norm_value = MatrixNorm(NewDensity)
-             CALL CopyMatrix(SubtractBranch, NewDensity)
-          END IF
-       ELSE
-          CALL MatrixMultiply(NewDensity, WorkingOverlap, TempMat1, &
-               & threshold_in=solver_parameters%threshold, memory_pool_in=pool1)
+       !! Xn+1 = Xn S1 Xn
+       CALL MatrixMultiply(WorkingDensity, WorkingOverlap, TempMat, &
+            & threshold_in=solver_parameters%threshold, memory_pool_in=pool1)
+       CALL MatrixMultiply(TempMat, WorkingDensity, NewDensity, &
+            & threshold_in=solver_parameters%threshold, memory_pool_in=pool1)
+
+       !! Figure Out Sigma Value
+       CALL DotMatrix(WorkingDensity, WorkingOverlap, trace_value)
+
+       !! Xn+1 = 2 Xn - Xn S1 Xn
+       IF (nel * 0.5_NTREAL .GT. trace_value) THEN
+          CALL ScaleMatrix(NewDensity, -1.0_NTREAL)
+          CALL IncrementMatrix(WorkingDensity, NewDensity, 2.0_NTREAL)
        END IF
 
-       IF (solver_parameters%be_verbose .AND. outer_counter .GT. 1) THEN
-          CALL WriteListElement(key="Round", value=outer_counter-1)
+       !! Check Convergence
+       CALL IncrementMatrix(NewDensity, WorkingDensity, -1.0_NTREAL)
+       norm_value = MatrixNorm(WorkingDensity)
+
+       IF (solver_parameters%be_verbose) THEN
+          CALL WriteListElement(key="Round", value=outer_counter)
           CALL EnterSubLog
           CALL WriteElement(key="Convergence", value=norm_value)
           CALL WriteElement(key="Trace", value=trace_value)
-          CALL WriteElement(key="AddTrace", value=add_trace)
-          CALL WriteElement(key="SubtractTrace", value=subtract_trace)
           CALL ExitSubLog
        END IF
 
        IF (norm_value .LE. solver_parameters%converge_diff) THEN
           EXIT
        END IF
 
-       !! Compute (I - XnS)Xn
-       CALL IncrementMatrix(Identity, TempMat1, -1.0_NTREAL)
-       CALL ScaleMatrix(TempMat1, -1.0_NTREAL)
-       CALL MatrixMultiply(TempMat1, NewDensity, TempMat2, &
-            & threshold_in=solver_parameters%threshold, memory_pool_in=pool1)
-
-       !! Subtracted Version Xn - (I - XnS)Xn
-       CALL CopyMatrix(NewDensity, SubtractBranch)
-       CALL IncrementMatrix(TempMat2, SubtractBranch, -1.0_NTREAL)
-
-       !! Added Version Xn + (I - XnS)Xn
-       CALL CopyMatrix(NewDensity, AddBranch)
-       CALL IncrementMatrix(TempMat2, AddBranch)
-
+       !! Xn = Xn+1
+       CALL CopyMatrix(NewDensity, WorkingDensity)
     END DO
-    total_iterations = outer_counter-1
     IF (solver_parameters%be_verbose) THEN
        CALL ExitSubLog
-       CALL WriteElement(key="Total_Iterations", value=total_iterations)
+       CALL WriteElement(key="Total_Iterations", value=outer_counter)
        CALL PrintMatrixInformation(NewDensity)
     END IF
 
@@ -178,11 +140,7 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, nel, NewDensity,&
     !! Cleanup
     CALL DestructMatrix(WorkingDensity)
     CALL DestructMatrix(WorkingOverlap)
-    CALL DestructMatrix(Identity)
-    CALL DestructMatrix(TempMat1)
-    CALL DestructMatrix(TempMat2)
-    CALL DestructMatrix(AddBranch)
-    CALL DestructMatrix(SubtractBranch)
+    CALL DestructMatrix(TempMat)
     CALL DestructMatrixMemoryPool(pool1)
     CALL DestructSolverParameters(solver_parameters)