From 1a27f5a1a00528f879614240a737b9f2278edd09 Mon Sep 17 00:00:00 2001 From: jacharrym2 Date: Fri, 16 Aug 2024 17:36:16 +0200 Subject: [PATCH] Extended the DDCISD method to Jadamilu diagonalizer. --- src/CI/CIJadamilu.f90 | 6 +- src/CI/CIOrder.f90 | 6 ++ src/CI/CIcore.f90 | 1 + src/CI/CImod.f90 | 143 ++++++++++++++++++++++++++++++------- test/H-e+H-.DD-CISD.lowdin | 4 +- test/H-e+H-.DD-CISD.py | 2 +- 6 files changed, 131 insertions(+), 31 deletions(-) diff --git a/src/CI/CIJadamilu.f90 b/src/CI/CIJadamilu.f90 index f4b587c..9e055b7 100644 --- a/src/CI/CIJadamilu.f90 +++ b/src/CI/CIJadamilu.f90 @@ -221,7 +221,7 @@ recursive function CIJadamilu_buildCouplingOrderRecursion( s, numberOfSpecies, end function CIJadamilu_buildCouplingOrderRecursion - subroutine CIJadamilu_jadamiluInterface(n, maxeig, eigenValues, eigenVectors) + subroutine CIJadamilu_jadamiluInterface(n, maxeig, eigenValues, eigenVectors, timeA, timeB) implicit none external DPJDREVCOM integer(8) :: maxnev @@ -249,7 +249,9 @@ subroutine CIJadamilu_jadamiluInterface(n, maxeig, eigenValues, eigenVectors) integer(8) :: I,J,K,ii,jj,jjj integer(4) :: iiter logical :: fullMatrix + real(8) :: timeA, timeB +!$ timeA = omp_get_wtime() maxsp = CONTROL_instance%CI_MADSPACE !!if ( CONTROL_instance%CI_JACOBI ) then @@ -372,6 +374,8 @@ subroutine CIJadamilu_jadamiluInterface(n, maxeig, eigenValues, eigenVectors) CALL PJDCLEANUP if ( allocated ( x ) ) deallocate ( x ) +!$ timeB = omp_get_wtime() + end subroutine CIJadamilu_jadamiluInterface subroutine matvec2 ( nx, v, w, iter) diff --git a/src/CI/CIOrder.f90 b/src/CI/CIOrder.f90 index 5db139d..398bf1f 100644 --- a/src/CI/CIOrder.f90 +++ b/src/CI/CIOrder.f90 @@ -127,6 +127,12 @@ subroutine CIOrder_settingCILevel() end select + if ( CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT == "CISD" .and. trim(CIcore_instance%level) /= "CISD" ) then + + call CIOrder_exception( ERROR, "Configuration interactor constructor", "DDCISD shift are only valid for CISD level!") + + end if + end subroutine CIOrder_settingCILevel diff --git a/src/CI/CIcore.f90 b/src/CI/CIcore.f90 index 842580d..aacc126 100644 --- a/src/CI/CIcore.f90 +++ b/src/CI/CIcore.f90 @@ -27,6 +27,7 @@ module CIcore_ type(vector) :: numberOfSpatialOrbitals2 type(vector8) :: eigenvalues type(vector) :: groundStateEnergies + type(vector) :: DDCISDTiming type(vector) :: lambda !!Number of particles per orbital, module only works for 1 or 2 particles per orbital type(matrix), allocatable :: fourCenterIntegrals(:,:) type(matrix), allocatable :: twoCenterIntegrals(:) diff --git a/src/CI/CImod.f90 b/src/CI/CImod.f90 index 98990b5..5116436 100644 --- a/src/CI/CImod.f90 +++ b/src/CI/CImod.f90 @@ -200,7 +200,6 @@ subroutine CImod_run() write (*,*) "Building initial hamiltonian..." call CIInitial_buildInitialCIMatrix2() - !!call CIFullMatrix_buildHamiltonianMatrix() This should be modified to build the CI matrix in memory call Matrix_constructor (CIcore_instance%eigenVectors, & int(CIcore_instance%numberOfConfigurations,8), & @@ -221,15 +220,76 @@ subroutine CImod_run() write(*,*) "Computer Physics Communications, vol. 177, pp. 951-964, 2007." write(*,*) "=============================================================" + !! diagonal correction. See 10.1016/j.chemphys.2007.07.001 + if ( CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT == "CISD") then + + call Vector_constructor ( CIcore_instance%groundStateEnergies, 30, 0.0_8) + call Vector_constructor ( CIcore_instance%DDCISDTiming, 30, 0.0_8) + + write (6,*) "" + write (6,"(T2,A50, A12)") " ITERATIVE DIAGONAL DRESSED CISD SHIFT: " , CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT + write (6,"(T2,A62)") " ( Size-extensive correction) " + write (6,"(T2,A62)") " Based on 10.1016/j.chemphys.2007.07.001 and 10.1063/5.0182498" + write (6,*) "" + + ecorr = 0.0_8 - call CIJadamilu_jadamiluInterface(CIcore_instance%numberOfConfigurations, & - int(CONTROL_instance%NUMBER_OF_CI_STATES,8), & - CIcore_instance%eigenvalues, & - CIcore_instance%eigenVectors ) + do i = 2, 31 + + !! add the diagonal shift + do a = 2, CIcore_instance%numberOfConfigurations + CIcore_instance%diagonalHamiltonianMatrix%values(a) = CIcore_instance%diagonalHamiltonianMatrix%values(a) + ecorr + end do + + call CIJadamilu_jadamiluInterface(CIcore_instance%numberOfConfigurations, & + int(CONTROL_instance%NUMBER_OF_CI_STATES,8), & + CIcore_instance%eigenvalues, & + CIcore_instance%eigenVectors, timeA, timeB) + + !! restore the original diagonal + do a = 2, CIcore_instance%numberOfConfigurations + CIcore_instance%diagonalHamiltonianMatrix%values(a) = CIcore_instance%diagonalHamiltonianMatrix%values(a) - ecorr + end do + + ecorr = CIcore_instance%eigenvalues%values(1) - HartreeFock_instance%totalEnergy + CIcore_instance%groundStateEnergies%values(i) = CIcore_instance%eigenvalues%values(1) + CIcore_instance%DDCISDTiming%values(i) = timeB - timeA + + write (6,"(T2,I2, F25.12, F25.12, F25.12, F16.4 )") i-1, CIcore_instance%groundStateEnergies%values(i), ecorr, (CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i)) , timeB - timeA + + !! Restart ci matrix diagonalization from previous eigenvectors + CONTROL_instance%CI_LOAD_EIGENVECTOR = .True. + + if ( abs( CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i) ) <= 1e-6) exit + + end do + + + write (6,*) "" + write (6,"(T2,A42 )") " ITERATIVE DIAGONAL DRESSED CONVERGENCE " + write (6,"(T2,A95 )") "Iter Ground-State Energy Correlation Energy Energy Diff. Time(s) " + do i = 2, 31 + write (6,"(T2,I2, F25.12, F25.12, F25.12, F16.4 )") i-1, CIcore_instance%groundStateEnergies%values(i), ecorr, (CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i)) , CIcore_instance%DDCISDTiming%values(i) + if ( abs( CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i) ) <= 1e-6) exit + end do + + if ( CONTROL_instance%CI_SAVE_EIGENVECTOR ) then + call CImod_saveEigenVector () + end if + + else !! standard CI, no diagonal correction + + call CIJadamilu_jadamiluInterface(CIcore_instance%numberOfConfigurations, & + int(CONTROL_instance%NUMBER_OF_CI_STATES,8), & + CIcore_instance%eigenvalues, & + CIcore_instance%eigenVectors, timeA, timeB ) + + if ( CONTROL_instance%CI_SAVE_EIGENVECTOR ) then + call CImod_saveEigenVector () + end if - if ( CONTROL_instance%CI_SAVE_EIGENVECTOR ) then - call CImod_saveEigenVector () end if + case ("DSYEVX") write (*,*) "Building Strings..." @@ -241,9 +301,6 @@ subroutine CImod_run() write (*,*) "Building diagonal..." call CIDiag_buildDiagonal() - !write (*,*) "Building Hamiltonian..." - !call CIFullMatrix_buildHamiltonianMatrix() - call Matrix_constructor (CIcore_instance%eigenVectors, & int(CIcore_instance%numberOfConfigurations,8), & int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) @@ -289,7 +346,7 @@ subroutine CImod_run() end do - else !! no diagonal correction + else !! standard CI, no diagonal correction call CIFullMatrix_buildHamiltonianMatrix(timeA, timeB) !$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for building Hamiltonian Matrix : ", timeB - timeA ," (s)" @@ -301,10 +358,6 @@ subroutine CImod_run() end if -! call Matrix_eigen_select (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & -! 1, CONTROL_instance%NUMBER_OF_CI_STATES, & -! flags = SYMMETRIC, dm = CIcore_instance%numberOfConfigurations ) - !! deallocate transformed integrals deallocate(CIcore_instance%twoCenterIntegrals) deallocate(CIcore_instance%fourCenterIntegrals) @@ -320,26 +373,62 @@ subroutine CImod_run() write (*,*) "Building diagonal..." call CIDiag_buildDiagonal() - write (*,*) "Building Hamiltonian..." - call CIFullMatrix_buildHamiltonianMatrix( timeA, timeB) - call Matrix_constructor (CIcore_instance%eigenVectors, & int(CIcore_instance%numberOfConfigurations,8), & int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + !! diagonal correction. See 10.1016/j.chemphys.2007.07.001 + if ( CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT == "CISD") then + + call Vector_constructor ( CIcore_instance%groundStateEnergies, 30, 0.0_8) + + write (6,*) "" + write (6,"(T2,A50, A12)") " ITERATIVE DIAGONAL DRESSED CISD SHIFT: " , CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT + write (6,"(T2,A62)") " ( Size-extensive correction) " + write (6,"(T2,A62)") " Based on 10.1016/j.chemphys.2007.07.001 and 10.1063/5.0182498" + write (6,*) "" + write (6,"(T2,A95 )") "Iter Ground-State Energy Correlation Energy Energy Diff. Time(s) " + + ecorr = 0.0_8 + + do i = 2, 31 + + call CIFullMatrix_buildHamiltonianMatrix( timeA, timeB) + + do a = 2, CIcore_instance%numberOfConfigurations + CIcore_instance%hamiltonianMatrix%values(a,a) = CIcore_instance%hamiltonianMatrix%values(a,a) + ecorr + end do + + call Matrix_eigen_dsyevr (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & + 1, CONTROL_instance%NUMBER_OF_CI_STATES, & + eigenVectors = CIcore_instance%eigenVectors, & + flags = SYMMETRIC) + + ecorr = CIcore_instance%eigenvalues%values(1) - HartreeFock_instance%totalEnergy + CIcore_instance%groundStateEnergies%values(i) = CIcore_instance%eigenvalues%values(1) + + write (6,"(T2,I2, F25.12, F25.12, F25.12, F16.4 )") i-1, CIcore_instance%groundStateEnergies%values(i), ecorr, (CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i)) , timeB - timeA + + if ( abs( CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i) ) <= 1e-6) exit + + end do + + else !! standard CI, no diagonal correction + + call CIFullMatrix_buildHamiltonianMatrix(timeA, timeB) +!$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for building Hamiltonian Matrix : ", timeB - timeA ," (s)" + + call Matrix_eigen_dsyevr (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & + 1, CONTROL_instance%NUMBER_OF_CI_STATES, & + eigenVectors = CIcore_instance%eigenVectors, & + flags = SYMMETRIC) + + end if + !! deallocate transformed integrals deallocate(CIcore_instance%twoCenterIntegrals) deallocate(CIcore_instance%fourCenterIntegrals) - call Matrix_eigen_dsyevr (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & - 1, CONTROL_instance%NUMBER_OF_CI_STATES, & - eigenVectors = CIcore_instance%eigenVectors, & - flags = SYMMETRIC) - -! call Matrix_eigen_dsyevr (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & -! 1, CONTROL_instance%NUMBER_OF_CI_STATES, & -! flags = SYMMETRIC, dm = CIcore_instance%numberOfConfigurations ) - case default call CImod_exception( ERROR, "CImod run", "Diagonalization method not implemented") diff --git a/test/H-e+H-.DD-CISD.lowdin b/test/H-e+H-.DD-CISD.lowdin index ab79cd6..f2e1052 100644 --- a/test/H-e+H-.DD-CISD.lowdin +++ b/test/H-e+H-.DD-CISD.lowdin @@ -16,8 +16,8 @@ END TASKS CONTROL readCoefficients=F numberOfCIstates=1 - CIdiagonalizationMethod = "DSYEVX" - !CIdiagonalizationMethod = "JADAMILU" + !CIdiagonalizationMethod = "DSYEVX" + CIdiagonalizationMethod = "JADAMILU" CIdiagonalDressedShift = "CISD" END CONTROL diff --git a/test/H-e+H-.DD-CISD.py b/test/H-e+H-.DD-CISD.py index 8bc5a65..743fba8 100644 --- a/test/H-e+H-.DD-CISD.py +++ b/test/H-e+H-.DD-CISD.py @@ -17,7 +17,7 @@ refValues = { "HF energy" : [-1.165428966723,1E-8], -"CISD energy" : [-1.284366244580,1E-8], +"CISD energy" : [-1.284366244580,1E-7], } testValues = dict(refValues) #copy