diff --git a/Source/Fortran/GeometryOptimizationModule.F90 b/Source/Fortran/GeometryOptimizationModule.F90 index 31df82d9..a040def8 100644 --- a/Source/Fortran/GeometryOptimizationModule.F90 +++ b/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,57 +78,39 @@ 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 @@ -143,25 +118,12 @@ SUBROUTINE PurificationExtrapolate(PreviousDensity, Overlap, nel, NewDensity,& 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)