diff --git a/src/CI/CIInitial.f90 b/src/CI/CIInitial.f90 index 9b5c18e..ceb41f6 100644 --- a/src/CI/CIInitial.f90 +++ b/src/CI/CIInitial.f90 @@ -413,6 +413,7 @@ function CIInitial_calculateEnergyOne( n, thisA, thisB ) result (auxCIenergy) end function CIInitial_calculateEnergyOne + function CIInitial_calculateEnergyTwo( n, thisA, thisB ) result (auxCIenergy) implicit none integer(8) :: thisA(:), thisB(:) @@ -521,7 +522,7 @@ function CIInitial_calculateEnergyTwo( n, thisA, thisB ) result (auxCIenergy) auxCIenergy = auxCIenergy + & CIcore_instance%fourCenterIntegrals(i,j)%values(auxIndex, 1) - + end if end do end do diff --git a/src/CI/CIJadamilu.f90 b/src/CI/CIJadamilu.f90 index 9e055b7..aa285e4 100644 --- a/src/CI/CIJadamilu.f90 +++ b/src/CI/CIJadamilu.f90 @@ -33,7 +33,7 @@ subroutine CIJadamilu_buildCouplingMatrix() integer(8), allocatable :: indexConfB(:) integer(1), allocatable :: couplingOrder(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies coupling = 0 !! allocate arrays @@ -79,7 +79,7 @@ subroutine CIJadamilu_buildCouplingOrderList() integer :: ssize integer, allocatable :: cilevel(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies ssize = 1 do i = 1, numberOfSpecies @@ -428,7 +428,7 @@ subroutine matvec2 ( nx, v, w, iter) if ( abs(v(i) ) >= tol) nonzero = nonzero + 1 end do - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies allocate ( indexConf ( numberOfSpecies, nproc ) ) allocate ( auxindexConf ( numberOfSpecies, nproc ) ) @@ -654,7 +654,7 @@ function CIJadamilu_buildMatrixRecursion2(nproc, s, indexConf, auxindexConf, cc, integer :: auxcilevel(:,:) integer, allocatable :: counter(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies allocate (counter(numberOfSpecies)) counter = 0 @@ -709,6 +709,7 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) integer :: numberOfSpecies, s integer, allocatable :: stringAinB(:) integer(4) :: coupling + integer(4) :: coupling2 integer(4) :: ssize,auxcoupling(3) !! 0,1,2 integer(8) :: indexConfA(:) integer(8), allocatable :: indexConfB(:) @@ -723,7 +724,7 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) !!$ CIcore_instance%timeA(1) = omp_get_wtime() - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies do i = 1, numberOfSpecies @@ -746,23 +747,27 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) do ci = 1, size(CIcore_instance%numberOfStrings(i)%values, dim = 1) do b = 1 + ssize , CIcore_instance%numberOfStrings(i)%values(ci) + ssize - !b = ssize + bb - do p = CIcore_instance%numberOfCoreOrbitals%values(i)+1, & - CIcore_instance%numberOfOccupiedOrbitals%values(i) - !do p = 1, & + !do p = CIcore_instance%numberOfCoreOrbitals%values(i)+1, & ! CIcore_instance%numberOfOccupiedOrbitals%values(i) + ! stringAinB(p) = CIcore_instance%orbitals(i)%values( & + ! CIcore_instance%strings(i)%values(p,a),b) + !end do - stringAinB(p) = CIcore_instance%orbitals(i)%values( & - CIcore_instance%strings(i)%values(p,a),b) + !coupling = CIcore_instance%numberOfOccupiedOrbitals%values(i) - sum ( stringAinB ) - CIcore_instance%numberOfCoreOrbitals%values(i) + + !coupling = 0 + !!$omp simd + !do p = CIcore_instance%numberOfCoreOrbitals%values(i)+1, CIcore_instance%numberOfOrbitals%values(i) + ! coupling = coupling + CIcore_instance%orbitals(i)%values(p,a) * CIcore_instance%orbitals(i)%values(p,b) + !end do + !coupling = CIcore_instance%numberOfOccupiedOrbitals%values(i) - coupling - CIcore_instance%numberOfCoreOrbitals%values(i) - !stringBinA(p) = CIcore_instance%orbitals(i)%values( & - ! CIcore_instance%strings(i)%values(p,b),a) - end do - coupling = CIcore_instance%numberOfOccupiedOrbitals%values(i) - sum ( stringAinB ) - & - CIcore_instance%numberOfCoreOrbitals%values(i) + coupling = CIcore_instance%numberOfOccupiedOrbitals%values(i) - sum(CIcore_instance%orbitals(i)%values(:,a) * CIcore_instance%orbitals(i)%values(:,b)) ! - CIcore_instance%numberOfCoreOrbitals%values(i) - ! coupling = CIcore_instance%numberOfOccupiedOrbitals%values(i) - sum ( stringAinB ) + !!$omp simd + !coupling = sum(abs(CIcore_instance%orbitals(i)%values(:,a) - CIcore_instance%orbitals(i)%values(:,b))) + !coupling = coupling / 2 if ( coupling <= 2 ) then @@ -823,6 +828,7 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) do ci = 1, size(CIcore_instance%numberOfStrings(i)%values, dim = 1) !! 1 is always zero cilevel(i) = ci - 1 + if ( CIcore_instance%nCouplingOneTwo(i,nn)%values( 2,ci ) == 0 ) cycle auxos = CIJadamilu_buildRowRecursionFirstOne( i, indexConfA, indexConfB, nn, cilevel ) end do @@ -832,7 +838,6 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) !!$ CIcore_instance%timeB(3) = omp_get_wtime() !!$ CIcore_instance%timeA(4) = omp_get_wtime() - !$omp atomic w(c) = w(c) + vc*CIcore_instance%diagonalHamiltonianMatrix%values(c) !$omp end atomic @@ -849,7 +854,7 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) do ci = 1, size(CIcore_instance%numberOfStrings(i)%values, dim = 1) !! 1 is always zero cilevel(i) = ci - 1 - + if ( CIcore_instance%nCouplingOneTwo(i,nn)%values( 2,ci ) == 0 ) cycle do u = 1, CIcore_instance%sizeciorderlist if ( sum(abs(cilevel - & CIcore_instance%ciorderlist( CIcore_instance%auxciorderlist(u), :))) == 0 ) then @@ -877,7 +882,7 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) do ci = 1, size(CIcore_instance%numberOfStrings(i)%values, dim = 1) !! 1 is always zero cilevel(i) = ci - 1 - + if ( CIcore_instance%nCouplingOneTwo(i,nn)%values( 3,ci ) == 0 ) cycle do u = 1, CIcore_instance%sizeCiOrderList if ( sum(abs(cilevel - & CIcore_instance%ciOrderList( CIcore_instance%auxciOrderList(u), :))) == 0 ) then @@ -911,7 +916,9 @@ subroutine CIJadamilu_buildRow( nn, indexConfA, c, w, vc, cilevelA) do ci = 1, size(CIcore_instance%numberOfStrings(i)%values, dim = 1) !! 1 is always zero cilevel(i) = ci - 1 + if ( CIcore_instance%nCouplingOneTwo(i,nn)%values( 2,ci ) == 0 ) cycle do cj = 1, size(CIcore_instance%numberOfStrings(j)%values, dim = 1) !! 1 is always zero + if ( CIcore_instance%nCouplingOneTwo(j,nn)%values( 2,cj ) == 0 ) cycle cilevel(j) = cj - 1 do u = 1, CIcore_instance%sizeCiOrderList if ( sum(abs(cilevel - & @@ -984,7 +991,7 @@ recursive function CIJadamilu_buildRowRecursionSecondOne( ii, indexConfB, w, vc real(8) :: CIenergy integer :: cilevel(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies ci = cilevel(ii) + 1 ssize = CIcore_instance%nCouplingSize(ii,nn)%values( 2,ci ) @@ -1030,7 +1037,7 @@ function CIJadamilu_buildRowRecursionSecondTwoCal( ii, indexConfA, indexConfB, w real(8) :: CIenergy integer :: cilevel(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies ci = cilevel(ii) + 1 ssize = CIcore_instance%nCouplingSize(ii,nn)%values( 3,ci ) @@ -1075,7 +1082,7 @@ function CIJadamilu_buildRowRecursionSecondTwoGet( ii, indexConfA, indexConfB, w real(8) :: CIenergy integer :: cilevel(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies ci = cilevel(ii) + 1 ssize = CIcore_instance%nCouplingSize(ii,nn)%values( 3,ci ) @@ -1107,20 +1114,22 @@ end function CIJadamilu_buildRowRecursionSecondTwoGet function CIJadamilu_buildRowRecursionSecondTwoDiff( ii, jj, indexConfB, w, vc, dd, nn, cilevel, u ) result (os) implicit none + integer, intent(in) :: ii, nn, u, jj + integer, intent(in) :: cilevel(:) + integer(8), intent(out) :: dd(:) + real(8), intent(in) :: vc + integer(8), intent(inout) :: indexConfB(:) + real(8), intent(inout) :: w(:) integer(8) :: ai,aj,d, aai, aaj - integer :: ii, nn, ci, u, k, jj, cj - integer :: ssizei, ssizej + integer :: ci, k, cj + integer(8) :: ssizei, ssizej + integer(8) :: dd_i_shift, dd_j_shift integer :: bi, bj, factor, factori integer :: auxIndex1, auxIndex2, auxIndex integer :: os,numberOfSpecies - integer(8) :: indexConfB(:) - integer(8) :: dd(:) - real(8) :: vc - real(8) :: w(:) real(8) :: CIenergy - integer :: cilevel(:) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies ci = cilevel(ii) + 1 cj = cilevel(jj) + 1 ssizei = CIcore_instance%nCouplingSize(ii,nn)%values( 2,ci ) @@ -1131,11 +1140,16 @@ function CIJadamilu_buildRowRecursionSecondTwoDiff( ii, jj, indexConfB, w, vc, d CIcore_instance%ciOrderSize1(u,k) )* CIcore_instance%ciOrderSize2(u,k) end do + dd_i_shift = - CIcore_instance%numberOfStrings2(ii)%values(ci) + & + CIcore_instance%ciOrderSize1(u,ii) + + dd_j_shift = - CIcore_instance%numberOfStrings2(jj)%values(cj) + & + CIcore_instance%ciOrderSize1(u,jj) + do aai = 1, CIcore_instance%nCouplingOneTwo(ii,nn)%values( 2,ci ) ai = ssizei + aai indexConfB(ii) = CIcore_instance%couplingMatrix(ii,nn)%values(ai, 2) - dd(ii) = (indexConfB(ii) - CIcore_instance%numberOfStrings2(ii)%values(ci) + & - CIcore_instance%ciOrderSize1(u,ii) )* CIcore_instance%ciOrderSize2(u,ii) + dd(ii) = (indexConfB(ii) + dd_i_shift )* CIcore_instance%ciOrderSize2(u,ii) bi = indexConfB(ii) factori = CIcore_instance%couplingMatrixFactorOne(ii,nn)%values(bi) @@ -1146,8 +1160,7 @@ function CIJadamilu_buildRowRecursionSecondTwoDiff( ii, jj, indexConfB, w, vc, d aj = ssizej + aaj indexConfB(jj) = CIcore_instance%couplingMatrix(jj,nn)%values(aj, 2) - dd(jj) = (indexConfB(jj) - CIcore_instance%numberOfStrings2(jj)%values(cj) + & - CIcore_instance%ciOrderSize1(u,jj) )* CIcore_instance%ciOrderSize2(u,jj) + dd(jj) = (indexConfB(jj) + dd_j_shift )* CIcore_instance%ciOrderSize2(u,jj) d = sum(dd) !CIenergy = vc*CIcore_calculateEnergyTwoDiff ( ii, jj, indexConfB, nn ) diff --git a/src/CI/CIOrder.f90 b/src/CI/CIOrder.f90 index 398bf1f..8485034 100644 --- a/src/CI/CIOrder.f90 +++ b/src/CI/CIOrder.f90 @@ -57,6 +57,14 @@ subroutine CIOrder_settingCILevel() CIcore_instance%maxCILevel = sum(CIcore_instance%CILevel) + case ( "SCI" ) !! same as FCI + + do i=1, numberOfSpecies + CIcore_instance%CILevel(i) = CIcore_instance%numberOfOccupiedOrbitals%values(i) + end do + + CIcore_instance%maxCILevel = sum(CIcore_instance%CILevel) + case ( "CIS" ) do i=1, numberOfSpecies diff --git a/src/CI/CISCI.f90 b/src/CI/CISCI.f90 new file mode 100644 index 0000000..dd1ce79 --- /dev/null +++ b/src/CI/CISCI.f90 @@ -0,0 +1,1027 @@ +module CISCI_ + use Exception_ + use Matrix_ + use Vector_ + use MolecularSystem_ + use Configuration_ + use MolecularSystem_ + use String_ + use IndexMap_ + use InputCI_ + use CIcore_ + use CIJadamilu_ + use CIInitial_ + use omp_lib + implicit none + + type, public :: CISCI + type (Vector8) :: amplitudeCore + type (Vector8) :: amplitudeCore2 + type (Vector8) :: coefficientCore + type (Matrix) :: coefficientTarget + type (Vector8) :: auxcoefficientTarget + type (Vector8) :: diagonalTarget + type (Vector8) :: eigenValues + integer :: coreSpaceSize + integer :: targetSpaceSize + real(8) :: PT2energy + end type CISCI + + type(CISCI) :: CISCI_instance + +contains + + subroutine CISCI_show() + + write (6,*) "" + write (6,"(T2,A62)") " SELECTED CONFIGURATION INTERACTION (SCI): " + write (6,"(T2,A62)") " Adaptive Sampling CI (ASCI) " + write (6,"(T2,A62)") " Deterministic Algorithm " + write (6,"(T2,A62)") " Based on 10.1063/1.4955109 " + + write(6,*) "" + write(6,*) " Diagonalizer for target space hamiltonian : ", trim(String_getUppercase((CONTROL_instance%CI_DIAGONALIZATION_METHOD))) + write(6,*) "=============================================================" + write(6,*) "M. BOLLHÖFER AND Y. NOTAY, JADAMILU:" + write(6,*) "a software code for computing selected eigenvalues of " + write(6,*) "large sparse symmetric matrices, " + write(6,*) "Computer Physics Communications, vol. 177, pp. 951-964, 2007." + write(6,*) "=============================================================" + + write (6,*) "" + write (6,"(T2,A,F14.5,A3 )") "Estimated memory needed: ", float(CIcore_instance%numberOfConfigurations*3*8)/(1024**3) , " GB" + write (6,*) "" + + end subroutine CISCI_show + + !! Allocating arrays + subroutine CISCI_constructor() + implicit none + + type(Configuration) :: auxConfigurationA, auxConfigurationB + integer :: a,b,c,aa,bb,i + real(8) :: CIenergy + integer :: nproc + + CISCI_instance%coreSpaceSize = CONTROL_instance%CI_SCI_CORE_SPACE + CISCI_instance%targetSpaceSize = CONTROL_instance%CI_SCI_TARGET_SPACE + + !! copy and destroying diagonal vector... because of the intial matrix subroutine + call Vector_constructor8 ( CIcore_instance%diagonalHamiltonianMatrix, & + CIcore_instance%numberOfConfigurations, 0.0_8 ) + CIcore_instance%diagonalHamiltonianMatrix%values = CIcore_instance%diagonalHamiltonianMatrix2%values + call Vector_destructor8 ( CIcore_instance%diagonalHamiltonianMatrix2 ) + + + !!This will carry the index changes after sorting configurations + call Vector_constructorInteger8 ( CIcore_instance%auxIndexCIMatrix, & + CIcore_instance%numberOfConfigurations, 0_8 ) !hmm + + do a = 1, CIcore_instance%numberOfConfigurations + CIcore_instance%auxIndexCIMatrix%values(a)= a + end do + + !! auxiliary array to get the index vector to build a configuration. get the configurations for the hamiltonian matrix in the core and target space + call CISCI_getInitialIndexes( CIcore_instance%coreConfigurations, CIcore_instance%coreConfigurationsLevel, CISCI_instance%coreSpaceSize ) + call CISCI_getInitialIndexes( CIcore_instance%targetConfigurations, CIcore_instance%targetConfigurationsLevel, CISCI_instance%targetSpaceSize ) + + !! arrays for CISCI + call Vector_constructor8 ( CISCI_instance%amplitudeCore, CIcore_instance%numberOfConfigurations, 0.0_8) + call Vector_constructor8 ( CISCI_instance%coefficientCore, int(CISCI_instance%coreSpaceSize,8), 0.0_8) + call Matrix_constructor ( CISCI_instance%coefficientTarget, int(CISCI_instance%targetSpaceSize,8), & + int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + call Vector_constructor8 ( CISCI_instance%auxcoefficientTarget, int(CISCI_instance%targetSpaceSize,8), 0.0_8) !! meh... just to avoid changing everything from matrix to vector format + call Vector_constructor8 ( CISCI_instance%diagonalTarget, int(CISCI_instance%targetSpaceSize,8), 0.0_8) !! Jadamilu requires to store diagonal vector (in target space) + call Vector_constructor8 ( CISCI_instance%eigenValues, 15_8, 0.0_8) !! store the eigenvalues per macro iterations + + end subroutine CISCI_constructor + + !! main part + subroutine CISCI_run() + implicit none + integer(8) :: i, j, ii, jj + integer :: k ! macro SCI iteration + real(8) :: timeA(15), timeB(15) + real(8) :: timeAA, timeBB + type(Vector8) :: eigenValuesTarget + real(8) :: currentEnergy + + currentEnergy = HartreeFock_instance%totalEnergy + call Vector_constructor8 ( eigenValuesTarget, int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + + !! HF determinant coefficient + CISCI_instance%coefficientCore%values(1) = 1.0_8 + CISCI_instance%coefficientTarget%values(1,1) = 1.0_8 + CIcore_instance%eigenVectors%values(1,1) = 1.0_8 + CISCI_instance%eigenValues%values(1) = HartreeFock_instance%totalEnergy + + write (6,*) "" + write (6,"(T2,A29 )") "Starting SCI macro iterations " + write (6,*) "" + do k = 2, 15 + +!$ timeA(k) = omp_get_wtime() + + !! calculating the amplitudes in core space. This is the pertubation guess of CI eigenvector + CISCI_instance%amplitudeCore%values = 0.0_8 + call CISCI_core_amplitudes ( CISCI_instance%amplitudeCore%values, CIcore_instance%numberOfConfigurations, & + CISCI_instance%coefficientCore%values, CISCI_instance%coreSpaceSize, currentEnergy ) + + !! setting the HF again... because the HF amplitude diverges + if ( k == 2 ) CISCI_instance%amplitudeCore%values(1) = 1.0_8 + + !! resetting the original index changes + do i = 1, CIcore_instance%numberOfConfigurations + CIcore_instance%auxIndexCIMatrix%values(i)= i + end do + + !! getting the target absolute largest coefficients + call Vector_sortElementsAbsolute8( CISCI_instance%amplitudeCore, & + CIcore_instance%auxIndexCIMatrix, int( CISCI_instance%targetSpaceSize ,8) ) + + !! recover the configurations for the hamiltonian matrix in the target space + call CISCI_getInitialIndexes( CIcore_instance%targetConfigurations, CIcore_instance%targetConfigurationsLevel, CISCI_instance%targetSpaceSize ) + + !! storing only the largest diagonal elements (for jadamilu) + do i = 1, CISCI_instance%targetSpaceSize + ii = CIcore_instance%auxIndexCIMatrix%values(i) + !! storing only the largest diagonal elements (for jadamilu) + CISCI_instance%diagonalTarget%values(i) = CIcore_instance%diagonalHamiltonianMatrix%values(ii) + enddo + + !! using the amplitued as the initial coeff guess, after that, use the previous diganolized eigenvectors in target space + if ( k == 2 ) then + do i = 1, CISCI_instance%targetSpaceSize + ii = CIcore_instance%auxIndexCIMatrix%values(i) + CISCI_instance%coefficientTarget%values(i,1) = CISCI_instance%amplitudeCore%values(ii) + enddo + else + do i = 1, CISCI_instance%targetSpaceSize + ii = CIcore_instance%auxIndexCIMatrix%values(i) + CISCI_instance%coefficientTarget%values(i,1) = CIcore_instance%eigenVectors%values(ii,1) + enddo + end if + + !! eigenvalue guess + eigenValuesTarget%values(1) = currentEnergy + + !! diagonalize in target space + call CISCI_jadamiluInterface( int(CISCI_instance%targetSpaceSize,8), & + 1_8, & + eigenValuesTarget, & + CISCI_instance%coefficientTarget, timeAA, timeBB ) + + !! saving the eigenvectors coeff to an aux vector. Only ground state + CISCI_instance%auxcoefficientTarget%values(:) = CISCI_instance%coefficientTarget%values(:,1) + + !! updating the full eigenvector with the solution in the target space + CIcore_instance%eigenVectors%values(:,1) = 0.0_8 + do i = 1, CISCI_instance%targetSpaceSize + ii = CIcore_instance%auxIndexCIMatrix%values(i) + CIcore_instance%eigenVectors%values(ii,1) = CISCI_instance%coefficientTarget%values(i,1) + end do + + !! convergence criteria + CISCI_instance%eigenValues%values(k) = eigenValuesTarget%values(1) + if ( abs( CISCI_instance%eigenValues%values(k) - currentEnergy ) < 1.0E-5 ) then + write (6,"(T2,A10,I4,A8,F25.12)") "SCI Iter: ", k , " Energy: ", CISCI_instance%eigenValues%values(k) + !$ timeB(k) = omp_get_wtime() + exit + end if + + !! getting the core absolute largest coefficients + call Vector_sortElementsAbsolute8( CISCI_instance%auxcoefficientTarget, & + CIcore_instance%auxIndexCIMatrix, int( CISCI_instance%coreSpaceSize ,8) ) + + !! recover the configurations for the hamiltonian matrix in the core space + call CISCI_getInitialIndexes( CIcore_instance%coreConfigurations, CIcore_instance%coreConfigurationsLevel, CISCI_instance%coreSpaceSize ) + !! call CISCI_getInitialIndexes( CIcore_instance%fullConfigurations, CIcore_instance%fullConfigurationsLevel , int(CIcore_instance%numberOfConfigurations,4) ) + + !! storing only the largest coefficients, and rearraing the next eigenvector guess + do i = 1, CISCI_instance%coreSpaceSize + CISCI_instance%coefficientCore%values(i) = CISCI_instance%auxcoefficientTarget%values(i) + enddo + + !! restart amplitudes for next run + CISCI_instance%amplitudeCore%values = 0.0_8 + write (6,"(T2,A10,I4,A8,F25.12)") "SCI Iter: ", k , " Energy: ", CISCI_instance%eigenValues%values(k) + +!$ timeB(k) = omp_get_wtime() + !! updating new reference + currentEnergy = eigenValuesTarget%values(1) + + enddo !k + + !! summary of the macro iteration + write (6,*) "" + write (6,"(T2,A95 )") " Selected CI (SCI) summary " + write (6,"(T2,A95 )") "Iter Ground-State Energy Correlation Energy Energy Diff. Time(s) " + do k = 2, 15 + write (6,"(T2,I2, F25.12, F25.12, F25.12, F16.4 )") k-1, CISCI_instance%eigenValues%values(k), & + CISCI_instance%eigenValues%values(k) - HartreeFock_instance%totalEnergy, & + CISCI_instance%eigenValues%values(k) - CISCI_instance%eigenValues%values(k-1), & + timeB(k) - timeA(k) + if ( abs( CISCI_instance%eigenValues%values(k) - CISCI_instance%eigenValues%values(k-1) ) < 1.0E-5 ) then + CIcore_instance%eigenvalues%values(1) = CISCI_instance%eigenValues%values(k) + exit + endif + enddo !k + + !! calculating PT2 correction. A pertuberd estimation of configurations not include in the target space + call CISCI_PT2 ( CIcore_instance%eigenVectors%values(:,1), CISCI_instance%amplitudeCore%values, & + CISCI_instance%targetSpaceSize, CIcore_instance%numberOfConfigurations, & + CISCI_instance%eigenValues%values(k), CISCI_instance%PT2energy ) + + write (6,*) "" + write (6,"(T2,A,F25.12)") "CI-PT2 energy correction :", CISCI_instance%PT2energy + + end subroutine CISCI_run + + + !> + !! @brief Muestra informacion del objeto + !! + !! @param this + !< + !! Map the indexes of initial CI matrix to the complete matrix. + subroutine CISCI_getInitialIndexes( auxConfigurationMatrix, auxConfigurationLevel, auxMatrixSize ) + implicit none + + type(imatrix) :: auxConfigurationMatrix + type(imatrix) :: auxConfigurationLevel + integer :: auxMatrixSize + integer(8) :: a,b,c + integer :: u,v + integer :: ci + integer :: i, j, ii, jj + integer :: s, numberOfSpecies, auxnumberOfSpecies + real(8) :: timeA, timeB + integer(1) :: coupling + integer(8) :: numberOfConfigurations + real(8) :: CIenergy + integer(8), allocatable :: indexConf(:) + integer, allocatable :: cilevel(:) + +!$ timeA = omp_get_wtime() + + numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + + s = 0 + c = 0 + + call Matrix_constructorInteger ( auxConfigurationMatrix, int( numberOfSpecies,8), & + int(auxMatrixSize,8), 0 ) + call Matrix_constructorInteger ( auxConfigurationLevel, int( numberOfSpecies,8), & + int(auxMatrixSize,8), 0 ) + + !! call recursion + allocate ( cilevel ( numberOfSpecies ) ) + allocate ( indexConf ( numberOfSpecies ) ) + + s = 0 + c = 0 + indexConf = 0 + cilevel = 0 + + do ci = 1, CIcore_instance%sizeCiOrderList + cilevel(:) = CIcore_instance%ciOrderList( CIcore_instance%auxciOrderList(ci), :) + s = 0 + auxnumberOfSpecies = CISCI_getIndexesRecursion( auxConfigurationMatrix, auxConfigurationLevel, auxMatrixSize, s, numberOfSpecies, indexConf, c, cilevel ) + end do + + deallocate ( indexConf ) + deallocate ( cilevel ) + +!$ timeB = omp_get_wtime() +!$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for getting sorted indexes : ", timeB - timeA ," (s)" + + end subroutine CISCI_getInitialIndexes + + +recursive function CISCI_getIndexesRecursion(auxConfigurationMatrix, auxConfigurationLevel, auxMatrixSize, s, numberOfSpecies, indexConf, c, cilevel) result (os) + implicit none + + type(imatrix) :: auxConfigurationMatrix + type(imatrix) :: auxConfigurationLevel + integer :: auxMatrixSize + integer(8) :: a,b,c + integer :: u,v + integer :: i, j, ii, jj + integer :: s, ss, numberOfSpecies + integer :: os,is + integer :: size1, size2 + integer(8) :: indexConf(:) + integer(1) :: coupling + integer :: ssize + integer :: cilevel(:) + + is = s + 1 + if ( is < numberOfSpecies ) then + i = cilevel(is) + 1 + ssize = CIcore_instance%numberOfStrings2(is)%values(i) + + do a = 1, CIcore_instance%numberOfStrings(is)%values(i) + indexConf(is) = ssize + a + os = CISCI_getIndexesRecursion( auxConfigurationMatrix, auxConfigurationLevel, auxMatrixSize, is, numberOfSpecies, indexConf, c, cilevel) + + end do + else + os = is + i = cilevel(is) + 1 + ssize = CIcore_instance%numberOfStrings2(is)%values(i) + + do a = 1, CIcore_instance%numberOfStrings(is)%values(i) + c = c + 1 + indexConf(is) = ssize + a + do u = 1, auxMatrixSize + if ( c == CIcore_instance%auxIndexCIMatrix%values(u) ) then + do ss = 1, numberOfSpecies + auxConfigurationMatrix%values(ss,u) = indexConf(ss) + auxConfigurationLevel%values(ss,u) = cilevel(ss) !? check... + end do + exit + end if + end do + end do + end if + + end function CISCI_getIndexesRecursion + + + subroutine CISCI_core_amplitudes ( amplitudeCore, numberOfConfigurations, coefficientCore, SCICoreSpaceSize, oldEnergy ) + + implicit none + + integer(4) SCICoreSpaceSize + integer(8) numberOfConfigurations + real(8) amplitudeCore ( numberOfConfigurations ) + real(8) coefficientCore ( SCICoreSpaceSize ) + real(8) :: CIEnergy + integer(8) :: nonzero + integer(8) :: i, j, ia, ib, ii, jj, iii, jjj + integer(4) :: nproc, n, nn + real(8) :: wi + real(8) :: timeA, timeB + real(8) :: tol + integer(4) :: iter, size1, size2 + integer :: ci + integer :: auxSize + integer(8) :: a,b,c, aa + integer :: s, numberOfSpecies + integer(8), allocatable :: indexConfA(:) !! ncore, species + integer, allocatable :: cilevel(:) + real(8) :: diagEnergy + real(8) :: oldEnergy + real(8) :: shift + +!$ timeA = omp_get_wtime() + call omp_set_num_threads(omp_get_max_threads()) + nproc = omp_get_max_threads() + shift = 1E-6 !! to avoid divergence +! shift = 0.0_8 + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies + + allocate ( indexConfA ( numberOfSpecies ) ) + allocate ( cilevel ( numberOfSpecies ) ) + + cilevel = 0 + indexConfA = 0 + + !$omp parallel & + !$omp& private( n, a, aa, cilevel, indexConfA, diagEnergy ),& + !$omp& shared( amplitudeCore, coefficientCore, oldEnergy ) + !$omp do schedule (static) + do a = 1, SCICoreSpaceSize + aa = CIcore_instance%auxIndexCIMatrix%values(a) + cilevel = CIcore_instance%coreConfigurationsLevel%values(:,a) + indexConfA = CIcore_instance%coreConfigurations%values(:,a) + n = OMP_GET_THREAD_NUM() + 1 + + !! using jadamilu subroutine to calculate all configurations coupled to configuration a + call CIJadamilu_buildRow( n, indexConfA, aa, amplitudeCore, coefficientCore(a), cilevel ) + + !! removing diagonal term + amplitudeCore(aa) = amplitudeCore(aa) - CIcore_instance%diagonalHamiltonianMatrix%values(aa) * coefficientCore(a) + + end do + !$omp end do nowait + !$omp end parallel + + do b = 1, CIcore_instance%numberOfConfigurations + amplitudeCore(b) = amplitudeCore(b) / ( CIcore_instance%diagonalHamiltonianMatrix%values(b) - oldEnergy + shift ) + end do + + CIcore_instance%pindexConf = 0 +!$ timeB = omp_get_wtime() + deallocate ( cilevel ) + deallocate ( indexConfA ) +!$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for calculating SCI amplitudes : ", timeB - timeA ," (s)" + + end subroutine CISCI_core_amplitudes + + subroutine CISCI_jadamiluInterface(n, maxeig, eigenValues, eigenVectors, timeA, timeB) + implicit none + external DPJDREVCOM + integer(8) :: maxnev + real(8) :: CIenergy + integer(8) :: nproc + type(Vector8), intent(inout) :: eigenValues + type(Matrix), intent(inout) :: eigenVectors + +! N: size of the problem +! MAXEIG: max. number of wanteg eig (NEIG<=MAXEIG) +! MAXSP: max. value of MADSPACE + integer(8) :: n, maxeig, MAXSP + integer(8) :: LX + real(8), allocatable :: EIGS(:), RES(:), X(:)!, D(:) +! arguments to pass to the routines + integer(8) :: NEIG, MADSPACE, ISEARCH, NINIT + integer(8) :: JA(1), IA(1) + integer(8) :: ICNTL(5) + integer(8) :: ITER, IPRINT, INFO + real(8) :: SIGMA, TOL, GAP, MEM, DROPTOL, SHIFT + integer(8) :: NDX1, NDX2, NDX3 + integer(8) :: IJOB! some local variables + integer(8) :: auxSize + integer(4) :: size1,size2 + 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 + + LX = N*(3*MAXSP+MAXEIG+1)+4*MAXSP*MAXSP + + if ( allocated ( eigs ) ) deallocate ( eigs ) + allocate ( eigs ( maxeig ) ) + eigs = 0.0_8 + if ( allocated ( res ) ) deallocate ( res ) + allocate ( res ( maxeig ) ) + res = 0.0_8 + if ( allocated ( x ) ) deallocate ( x ) + allocate ( x ( lx ) ) + x = 0.0_8 + +! set input variables + IPRINT = 0 ! standard report on standard output + ISEARCH = 1 ! we want the smallest eigenvalues + NEIG = maxeig ! number of wanted eigenvalues + !NINIT = 0 ! no initial approximate eigenvectors + NINIT = NEIG ! initial approximate eigenvectors + MADSPACE = maxsp ! desired size of the search space + ITER = 1000*NEIG ! maximum number of iteration steps + TOL = CONTROL_instance%CI_CONVERGENCE !1.0d-4 ! tolerance for the eigenvector residual + TOL = 1e-3 !1.0d-4 ! tolerance for the eigenvector residual, for ASCI this can be higher + + NDX1 = 0 + NDX2 = 0 + MEM = 0 + +! additional parameters set to default + ICNTL(1)=0 + ICNTL(2)=0 + ICNTL(3)=0 + ICNTL(4)=0 + ICNTL(5)=1 + + IJOB=0 + + JA(1) = -1 + IA(1) = -1 + + ! set initial eigenpairs + do j = 1, n + X(j) = eigenVectors%values(j,1) + end do + + do i = 1, CONTROL_instance%NUMBER_OF_CI_STATES + EIGS(i) = eigenValues%values(i) + end do + + DROPTOL = 1E-4 + + SIGMA = EIGS(1) + gap = 0 + SHIFT = 0!EIGS(1) + + do i = 1, CONTROL_instance%NUMBER_OF_CI_STATES + write(6,"(T2,A5,I4,2X,A10,F20.10,2X,A11,F20.10)") "State", i, "Eigenvalue", eigs( i ), "Eigenvector", x((i-1)*n + i) + end do + + iiter = 0 + +10 CALL DPJDREVCOM( N, CISCI_instance%diagonalTarget%values , JA, IA, EIGS, RES, X, LX, NEIG, & + SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL, & + SHIFT, DROPTOL, MEM, ICNTL, & + IJOB, NDX1, NDX2, IPRINT, INFO, GAP) + if (CONTROL_instance%CI_JACOBI ) then + fullMatrix = .false. + else + fullMatrix = .true. + end if + +!! your private matrix-vector multiplication + iiter = iiter +1 + IF (IJOB.EQ.1) THEN + call CISCI_matvec ( N, X(NDX1), X(NDX2), iiter) + GOTO 10 + END IF + + !! saving the eigenvalues + eigenValues%values = EIGS + + !! saving the eigenvectors + k = 0 + do j = 1, maxeig + do i = 1, N + k = k + 1 + eigenVectors%values(i,j) = X(k) + end do + end do + +! release internal memory and discard preconditioner + CALL PJDCLEANUP + if ( allocated ( x ) ) deallocate ( x ) + +!$ timeB = omp_get_wtime() + + end subroutine CISCI_jadamiluInterface + + subroutine CISCI_matvec ( nx, v, w, iter) + + !******************************************************************************* + !! AV computes w <- A * V where A is a discretized Laplacian. + ! Parameters: + ! Input, integer NX, the length of the vectors. + ! Input, real V(NX), the vector to be operated on by A. + ! Output, real W(NX), the result of A*V. + ! + implicit none + + integer(8) nx + real(8) v(nx) + real(8) w(nx) + integer(4) :: iter + integer(8) :: a,b,aa,bb + integer(8) :: nonzero, nonzerow + real(8) :: tol + integer :: uu,vv + integer :: i, ii, jj, n + integer :: numberOfSpecies + real(8) :: timeA, timeB + real(8) :: CIenergy + real(8) :: CIenergy2 + integer(1) :: coupling + integer(1), allocatable :: orbitalsA(:,:), orbitalsB(:,:), couplingS(:) + integer :: initialCIMatrixSize + integer :: nproc + integer(8), allocatable :: indexConfA(:) + integer(8), allocatable :: indexConfB(:) + + numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + call omp_set_num_threads(omp_get_max_threads()) + nproc = omp_get_max_threads() + + allocate ( couplingS ( numberOfSpecies ) ) + allocate ( indexConfA ( numberOfSpecies ) ) + allocate ( indexConfB ( numberOfSpecies ) ) + allocate (orbitalsA ( maxval(CIcore_instance%numberOfOrbitals%values(:)), numberOfSpecies)) + allocate (orbitalsB ( maxval(CIcore_instance%numberOfOrbitals%values(:)), numberOfSpecies)) + + + nonzero = 0 + nonzerow = 0 + w = 0.0_8 + tol = CONTROL_instance%CI_MATVEC_TOLERANCE + do a = 1 , nx + if ( abs(v(a) ) >= tol) nonzero = nonzero + 1 + end do + +!$ timeA= omp_get_wtime() + + !$omp parallel & + !$omp& private( a, aa, b, bb, uu, vv, coupling, couplingS, CIenergy, i, ii, jj, orbitalsA, orbitalsB),& + !$omp& firstprivate( indexConfA, indexConfB ),& + !$omp& shared( v, nx, numberOfSpecies, tol ) reduction (+:w) + !$omp do schedule (dynamic) + do a = 1, nx + + indexConfA = CIcore_instance%targetConfigurations%values(:,a) + + orbitalsA = 0 + do i = 1, numberOfSpecies + do uu = 1, CIcore_instance%numberOfOccupiedOrbitals%values(i) + orbitalsA( CIcore_instance%strings(i)%values(uu,indexConfA(i) ), i ) = 1 + end do + end do + + do b = a, nx + + couplingS = 0 + indexConfB(:) = CIcore_instance%targetConfigurations%values(:,b) + + orbitalsB = 0 + do i = 1, numberOfSpecies + do vv = 1, CIcore_instance%numberOfOccupiedOrbitals%values(i) + orbitalsB( CIcore_instance%strings(i)%values(vv,indexConfB(i) ), i ) = 1 + end do + end do + + do i = 1, numberOfSpecies + couplingS(i) = couplingS(i) + & + CIcore_instance%numberOfOccupiedOrbitals%values(i) - sum ( orbitalsA(:,i) * orbitalsB(:,i) ) + 1 + end do + + coupling = product(couplingS) + + select case (coupling) + + case(1) + CIenergy = CISCI_instance%diagonalTarget%values(a) + w(a) = w(a) + CIEnergy*v(a) + + case(2) + do i = 1, numberOfSpecies + if ( couplingS(i) == 2 ) ii = i + end do + + CIenergy = CISCI_calculateEnergyOne ( ii, indexConfA, indexConfB ) + w(b) = w(b) + CIenergy * v(a) + w(a) = w(a) + CIenergy * v(b) + + case(3) + do i = 1, numberOfSpecies + if ( couplingS(i) == 3 ) ii = i + end do + + CIenergy = CISCI_calculateEnergyTwoSame ( ii, indexConfA, indexConfB ) + w(b) = w(b) + CIenergy * v(a) + w(a) = w(a) + CIenergy * v(b) + + case(4) + do i = 1, numberOfSpecies + if ( couplingS(i) == 2 ) then + ii = i + exit + end if + end do + do i = ii+1, numberOfSpecies + if ( couplingS(i) == 2 ) jj = i + end do + CIenergy = CISCI_calculateEnergyTwoDiff ( ii, jj, indexConfA, indexConfB ) + w(b) = w(b) + CIenergy * v(a) + w(a) = w(a) + CIenergy * v(b) + + end select + + end do !b + end do !a + !$omp end do nowait + !$omp end parallel + +!$ timeB = omp_get_wtime() + !! to check how dense is the w vector + do a = 1 , nx + if ( abs(w(a) ) >= tol) nonzerow = nonzerow + 1 + end do + !stop + deallocate (orbitalsA ) + deallocate (orbitalsB ) + deallocate ( indexConfB ) + deallocate ( indexConfA ) + deallocate ( couplingS ) + +!$ write(*,"(A,I2,A,E10.3,A2,I12,I12)") " ", iter, " ", timeB -timeA ," ", nonzero, nonzerow + return + + end subroutine CISCI_matvec + + function CISCI_calculateEnergyOne( ii, thisA, thisB ) result (auxCIenergy) + implicit none + integer(8) :: thisA(:), thisB(:) + integer(8) :: a, b + integer :: i,j,s,n, nn,ii + integer :: l,k,z,kk,ll + integer :: factor, factor2, auxOcc, AA, BB + logical(1) :: equalA, equalB + integer :: auxnumberOfOtherSpecieSpatialOrbitals + integer(8) :: auxIndex1, auxIndex11, auxIndex2, auxIndex + integer :: diffOrb(2), otherdiffOrb(2) !! to avoid confusions + real(8) :: auxCIenergy + + auxCIenergy = 0.0_8 + factor = 1 + + !! copy a + a = thisA(ii) + b = thisB(ii) + + diffOrb = 0 + + do kk = 1, CIcore_instance%occupationNumber(ii) + if ( CIcore_instance%orbitals(ii)%values( & + CIcore_instance%strings(ii)%values(kk,a),b) == 0 ) then + diffOrb(1) = CIcore_instance%strings(ii)%values(kk,a) + AA = kk + exit + end if + end do + + do kk = 1, CIcore_instance%occupationNumber(ii) + if ( CIcore_instance%orbitals(ii)%values( & + CIcore_instance%strings(ii)%values(kk,b),a) == 0 ) then + diffOrb(2) = CIcore_instance%strings(ii)%values(kk,b) + BB = kk + exit + end if + end do + + factor = (-1)**(AA-BB) + + !One particle terms + + auxCIenergy= auxCIenergy + CIcore_instance%twoCenterIntegrals(ii)%values( diffOrb(1), diffOrb(2) ) + + !! save the different orbitals + + auxIndex1= CIcore_instance%twoIndexArray(ii)%values( diffOrb(1), diffOrb(2)) + + do ll=1, CIcore_instance%occupationNumber( ii ) !! the same orbitals pair are excluded by the exchange + + l = CIcore_instance%strings(ii)%values(ll,b) !! or a + + auxIndex2 = CIcore_instance%twoIndexArray(ii)%values( l,l) + auxIndex = CIcore_instance%fourIndexArray(ii)%values( auxIndex1, auxIndex2 ) + + auxCIenergy = auxCIenergy + CIcore_instance%fourCenterIntegrals(ii,ii)%values(auxIndex, 1) + + auxIndex = CIcore_instance%fourIndexArray(ii)%values( & + CIcore_instance%twoIndexArray(ii)%values(diffOrb(1),l), & + CIcore_instance%twoIndexArray(ii)%values(l,diffOrb(2)) ) + + auxCIenergy = auxCIenergy + & + MolecularSystem_instance%species(ii)%kappa*CIcore_instance%fourCenterIntegrals(ii,ii)%values(auxIndex, 1) + + end do + + !end if + do j=1, ii - 1 !! avoid ii, same species + + b = thisB(j) + + auxnumberOfOtherSpecieSpatialOrbitals = CIcore_instance%numberOfSpatialOrbitals2%values(j) + auxIndex11 = auxnumberOfOtherSpecieSpatialOrbitals * (auxIndex1 - 1 ) + + do ll=1, CIcore_instance%occupationNumber( j ) + + l = CIcore_instance%strings(j)%values(ll,b) + + auxIndex = auxIndex11 + CIcore_instance%twoIndexArray(j)%values( l,l) + + auxCIenergy = auxCIenergy + & + CIcore_instance%fourCenterIntegrals(ii,j)%values(auxIndex, 1) + + end do + + end do + + do j= ii + 1, MolecularSystem_instance%numberOfQuantumSpecies!! avoid ii, same species + + b = thisB(j) + + auxnumberOfOtherSpecieSpatialOrbitals = CIcore_instance%numberOfSpatialOrbitals2%values(j) + + auxIndex11 = auxnumberOfOtherSpecieSpatialOrbitals * (auxIndex1 - 1 ) + + do ll=1, CIcore_instance%occupationNumber( j ) + + l = CIcore_instance%strings(j)%values(ll,b) + + auxIndex = auxIndex11 + CIcore_instance%twoIndexArray(j)%values( l,l) + + auxCIenergy = auxCIenergy + & + CIcore_instance%fourCenterIntegrals(ii,j)%values(auxIndex, 1) + end do + + end do + + auxCIenergy= auxCIenergy * factor + + end function CISCI_calculateEnergyOne + + function CISCI_calculateEnergyTwoSame( ii, thisA, thisB ) result (auxCIenergy) + implicit none + integer(8) :: a, b + integer :: ii + integer :: kk,z + integer :: factor, AA(2), BB(2) + integer(8) :: thisA(:), thisB(:) + integer(8) :: auxIndex + integer :: diffOrbA(2), diffOrbB(2) !! to avoid confusions + real(8) :: auxCIenergy + + a = thisA(ii) + b = thisB(ii) + !diffOrbA = 0 + !diffOrbB = 0 + z = 0 + auxCIenergy = 0.0_8 + + do kk = 1, CIcore_instance%occupationNumber(ii) + if ( CIcore_instance%orbitals(ii)%values( & + CIcore_instance%strings(ii)%values(kk,a),b) == 0 ) then + z = z + 1 + diffOrbA(z) = CIcore_instance%strings(ii)%values(kk,a) + AA(z) = kk + if ( z == 2 ) exit + end if + end do + + z = 0 + do kk = 1, CIcore_instance%occupationNumber(ii) + if ( CIcore_instance%orbitals(ii)%values( & + CIcore_instance%strings(ii)%values(kk,b),a) == 0 ) then + z = z + 1 + diffOrbB(z) = CIcore_instance%strings(ii)%values(kk,b) + BB(z) = kk + if ( z == 2 ) exit + end if + end do + + factor = (-1)**(AA(1)-BB(1) + AA(2) - BB(2) ) + auxIndex = CIcore_instance%fourIndexArray(ii)%values( & + CIcore_instance%twoIndexArray(ii)%values(& + diffOrbA(1),diffOrbB(1)),& + CIcore_instance%twoIndexArray(ii)%values(& + diffOrbA(2),diffOrbB(2)) ) + + auxCIenergy = CIcore_instance%fourCenterIntegrals(ii,ii)%values(auxIndex, 1) + + auxIndex = CIcore_instance%fourIndexArray(ii)%values( & + CIcore_instance%twoIndexArray(ii)%values(& + diffOrbA(1),diffOrbB(2)),& + CIcore_instance%twoIndexArray(ii)%values(& + diffOrbA(2),diffOrbB(1)) ) + auxCIenergy = auxCIenergy + & + MolecularSystem_instance%species(ii)%kappa*CIcore_instance%fourCenterIntegrals(ii,ii)%values(auxIndex, 1) + + auxCIenergy= auxCIenergy * factor + + end function CISCI_calculateEnergyTwoSame + + function CISCI_calculateEnergyTwoDiff( ii, jj, thisA, thisB ) result (auxCIenergy) + implicit none + integer(8) :: a, b + integer :: ii, jj + integer :: kk,z + integer :: factori, factorj, AA, BB + integer(8) :: thisA(:), thisB(:) + integer(8) :: auxIndex, auxIndex1, auxIndex2 + integer :: diffOrb(2) + real(8) :: auxCIenergy + + a = thisA(ii) + b = thisB(ii) + + diffOrb = 0 + + do kk = 1, CIcore_instance%occupationNumber(ii) + if ( CIcore_instance%orbitals(ii)%values( & + CIcore_instance%strings(ii)%values(kk,a),b) == 0 ) then + diffOrb(1) = CIcore_instance%strings(ii)%values(kk,a) + AA = kk + exit + end if + end do + + do kk = 1, CIcore_instance%occupationNumber(ii) + if ( CIcore_instance%orbitals(ii)%values( & + CIcore_instance%strings(ii)%values(kk,b),a) == 0 ) then + diffOrb(2) = CIcore_instance%strings(ii)%values(kk,b) + BB = kk + exit + end if + end do + + factori = (-1)**(AA-BB) + auxIndex1= CIcore_instance%twoIndexArray(ii)%values( diffOrb(1), diffOrb(2)) + auxIndex1 = CIcore_instance%numberOfSpatialOrbitals2%values(jj) * (auxIndex1 - 1 ) + + a = thisA(jj) + b = thisB(jj) + + diffOrb = 0 + + do kk = 1, CIcore_instance%occupationNumber(jj) + if ( CIcore_instance%orbitals(jj)%values( & + CIcore_instance%strings(jj)%values(kk,a),b) == 0 ) then + diffOrb(1) = CIcore_instance%strings(jj)%values(kk,a) + AA = kk + exit + end if + end do + + do kk = 1, CIcore_instance%occupationNumber(jj) + if ( CIcore_instance%orbitals(jj)%values( & + CIcore_instance%strings(jj)%values(kk,b),a) == 0 ) then + diffOrb(2) = CIcore_instance%strings(jj)%values(kk,b) + BB = kk + exit + end if + end do + + factorj = (-1)**(AA-BB) + + auxIndex2= CIcore_instance%twoIndexArray(jj)%values( diffOrb(1), diffOrb(2)) + auxIndex = auxIndex1 + auxIndex2 + + auxCIenergy = factori * factorj *CIcore_instance%fourCenterIntegrals(ii,jj)%values(auxIndex, 1) + + end function CISCI_calculateEnergyTwoDiff + + subroutine CISCI_PT2 ( coefficients, auxenergyCorrection, SCITargetSpaceSize, numberOfConfigurations, refEnergy, energyCorrection ) + + !******************************************************************************* + !! AV computes w <- A * V where A is a discretized Laplacian. + ! Parameters: + ! Input, integer NX, the length of the vectors. + ! Input, real V(NX), the vector to be operated on by A. + ! Output, real W(NX), the result of A*V. + ! + implicit none + integer :: SCITargetSpaceSize + integer(8) numberOfConfigurations + real(8) :: coefficients ( numberOfConfigurations ) + real(8) :: auxenergyCorrection ( numberOfConfigurations ) + real(8) :: refEnergy + real(8) :: energyCorrection + real(8) :: CIEnergy + integer(8) :: nonzero + integer(8) :: i, j, ia, ib, ii, jj, iii, jjj + integer(4) :: nproc, n, nn + real(8) :: wi + real(8) :: timeA, timeB + real(8) :: tol + integer(4) :: iter, size1, size2 + integer :: ci + integer :: auxSize + integer(8) :: a,b,c, aa, bb + integer :: s, numberOfSpecies + integer(8), allocatable :: indexConfA(:) !! ncore, species + integer, allocatable :: cilevel(:) + real(8) :: diagEnergy + real(8) :: oldEnergy + real(8) :: shift + + call omp_set_num_threads(omp_get_max_threads()) + nproc = omp_get_max_threads() +!$ timeA = omp_get_wtime() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies + + allocate ( indexConfA ( numberOfSpecies ) ) + allocate ( cilevel ( numberOfSpecies ) ) + + cilevel = 0 + indexConfA = 0 + energyCorrection = 0.0_8 + auxenergyCorrection = 0.0_8 + + !$omp parallel & + !$omp& private( n, a, aa, cilevel, indexConfA ),& + !$omp& shared( auxenergyCorrection, coefficients, refEnergy ) + !$omp do schedule (static) + do a = 1, SCITargetSpaceSize + aa = CIcore_instance%auxIndexCIMatrix%values(a) + cilevel = CIcore_instance%targetConfigurationsLevel%values(:,a) + indexConfA = CIcore_instance%targetConfigurations%values(:,a) + n = OMP_GET_THREAD_NUM() + 1 + + !! using jadamilu subroutine to calculate all configurations coupled to configuration a + call CIJadamilu_buildRow( n, indexConfA, aa, auxenergyCorrection, coefficients(aa), cilevel ) + + !! removing the contributions from configurations within the target space + do b = 1, SCITargetSpaceSize + bb = CIcore_instance%auxIndexCIMatrix%values(b) + auxenergyCorrection(bb) = 0.0_8 + enddo + + end do + !$omp end do nowait + !$omp end parallel + + do b = 1, numberOfConfigurations + energyCorrection = energyCorrection + auxenergyCorrection(b) **2 / ( refEnergy - CIcore_instance%diagonalHamiltonianMatrix%values(b) ) + enddo + + + CIcore_instance%pindexConf = 0 +!$ timeB = omp_get_wtime() + deallocate ( cilevel ) + deallocate ( indexConfA ) +!$ write(*,"(A,E10.3)") "Time for CI-PT2 correction: ", timeB -timeA + + end subroutine CISCI_PT2 + +end module CISCI_ diff --git a/src/CI/CIcore.f90 b/src/CI/CIcore.f90 index aacc126..5adb3a3 100644 --- a/src/CI/CIcore.f90 +++ b/src/CI/CIcore.f90 @@ -21,6 +21,7 @@ module CIcore_ type(vector8) :: initialEigenValues integer(8) :: numberOfConfigurations integer :: nproc + integer :: numberOfQuantumSpecies type(ivector) :: numberOfCoreOrbitals type(ivector) :: numberOfOccupiedOrbitals type(ivector) :: numberOfOrbitals @@ -33,8 +34,8 @@ module CIcore_ type(matrix), allocatable :: twoCenterIntegrals(:) type(imatrix8), allocatable :: twoIndexArray(:) type(imatrix8), allocatable :: fourIndexArray(:) - type(imatrix), allocatable :: strings(:) !! species, conf, occupations - type(imatrix1), allocatable :: orbitals(:) !! species, conf, occupations + type(imatrix), allocatable :: strings(:) !! species, conf, occupations. index for occupied orbitals, e.g. 1 2 5 6 + type(imatrix1), allocatable :: orbitals(:) !! species, conf, occupations. array with 1 for occupied and 0 unoccupied orb, e.g. 1 1 0 0 1 1 integer, allocatable :: sumstrings(:) !! species type(ivector), allocatable :: auxstring(:,:) !! species, occupations type(ivector8), allocatable :: numberOfStrings(:) !! species, excitation level, number of strings @@ -64,6 +65,13 @@ module CIcore_ integer :: ncouplingOrderTwoDiff type(imatrix) :: auxConfigurations !! species, configurations for initial hamiltonian + type(imatrix) :: coreConfigurations !! species, configurations for core SCI space + type(imatrix) :: targetConfigurations !! species, configurations for target SCI space + type(imatrix) :: fullConfigurations !! species, configurations for target SCI space + type(imatrix) :: coreConfigurationsLevel !! species, configurations for CI level of core SCI space + type(imatrix) :: targetConfigurationsLevel !! species, configurations for CI level target SCI space + type(imatrix) :: fullConfigurationsLevel !! species, configurations for CI level target SCI space + type(configuration), allocatable :: configurations(:) integer(2), allocatable :: auxconfs(:,:,:) ! nconf, species, occupation type (Vector8) :: diagonalHamiltonianMatrix @@ -74,7 +82,8 @@ module CIcore_ integer, allocatable :: recursionVector1(:) integer, allocatable :: recursionVector2(:) integer, allocatable :: CILevel(:) - integer, allocatable :: pindexConf(:,:) + integer, allocatable :: pindexConf(:,:) !! save previous configuration to avoid unneccesary calculations + integer :: maxCILevel type (Matrix) :: initialHamiltonianMatrix type (Matrix) :: initialHamiltonianMatrix2 @@ -137,7 +146,8 @@ subroutine CIcore_constructor(level) call Vector_getFromFile(unit=wfnUnit, binary=.true., value=HartreeFock_instance%puntualInteractionEnergy, & arguments=["PUNTUALINTERACTIONENERGY"]) - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + CIcore_instance%numberOfQuantumSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies CIcore_instance%numberOfSpecies = numberOfSpecies @@ -287,7 +297,6 @@ subroutine CIcore_constructor(level) !!Even occupation number = beta end do - call Configuration_globalConstructor() close(wfnUnit) @@ -338,7 +347,7 @@ function CIcore_getIndex ( indexConf ) result ( output ) integer(8) :: output, ssize integer :: i,j, numberOfSpecies - numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + numberOfSpecies = CIcore_instance%numberOfQuantumSpecies output = 0 !! simplify!! do i = 1, numberOfSpecies diff --git a/src/CI/CImod.f90 b/src/CI/CImod.f90 index 5116436..aabcf03 100644 --- a/src/CI/CImod.f90 +++ b/src/CI/CImod.f90 @@ -44,6 +44,7 @@ module CImod_ use CIDiag_ use CIFullMatrix_ use CIInitial_ + use CISCI_ use CIJadamilu_ use CIOrder_ use CIStrings_ @@ -92,13 +93,14 @@ module CImod_ subroutine CImod_run() implicit none integer :: i,j,m, numberOfSpecies - integer :: a + integer :: a, ms real(8) :: timeA, timeB real(8), allocatable :: eigenValues(:) real(8) :: ecorr ! select case ( trim(CIcore_instance%level) ) numberOfSpecies = MolecularSystem_getNumberOfQuantumSpecies() + ms = CONTROL_instance%CI_MADSPACE write (*,*) "" write (*,*) "===============================================" @@ -178,263 +180,323 @@ subroutine CImod_run() call Vector_constructor8 ( CIcore_instance%eigenvalues, & int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8 ) - select case (trim(String_getUppercase(CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - ! case ("ARPACK") + if ( CONTROL_instance%CONFIGURATION_INTERACTION_LEVEL /= "SCI" ) then - ! write (*,*) "This method was removed" + select case (trim(String_getUppercase(CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - case ("JADAMILU") + ! case ("ARPACK") - write (*,*) "Building Strings..." - call CIStrings_buildStrings() + ! write (*,*) "This method was removed" - write (*,*) "Building CI level table..." - call CIOrder_buildCIOrderList() + case ("JADAMILU") - call CIJadamilu_buildCouplingMatrix() - call CIJadamilu_buildCouplingOrderList() + write (*,*) "Building Strings..." + call CIStrings_buildStrings() - write (*,*) "Building diagonal..." - call CIDiag_buildDiagonal() + write (*,*) "Building CI level table..." + call CIOrder_buildCIOrderList() - write (*,*) "Building initial hamiltonian..." - call CIInitial_buildInitialCIMatrix2() + call CIJadamilu_buildCouplingMatrix() + call CIJadamilu_buildCouplingOrderList() - call Matrix_constructor (CIcore_instance%eigenVectors, & - int(CIcore_instance%numberOfConfigurations,8), & - int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + write (*,*) "Building diagonal..." + call CIDiag_buildDiagonal() - if ( CONTROL_instance%CI_LOAD_EIGENVECTOR ) then - call CImod_loadEigenVector (CIcore_instance%eigenvalues, & - CIcore_instance%eigenVectors) - end if + write (*,*) "Building initial hamiltonian..." + call CIInitial_buildInitialCIMatrix2() - write(*,*) "" - write(*,*) "Diagonalizing hamiltonian..." - write(*,*) " Using : ", trim(String_getUppercase((CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - write(*,*) "=============================================================" - write(*,*) "M. BOLLHÖFER AND Y. NOTAY, JADAMILU:" - write(*,*) " a software code for computing selected eigenvalues of " - write(*,*) " large sparse symmetric matrices, " - write(*,*) "Computer Physics Communications, vol. 177, pp. 951-964, 2007." - write(*,*) "=============================================================" + 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 + if ( CONTROL_instance%CI_LOAD_EIGENVECTOR ) then + call CImod_loadEigenVector (CIcore_instance%eigenvalues, & + CIcore_instance%eigenVectors) + end if - call Vector_constructor ( CIcore_instance%groundStateEnergies, 30, 0.0_8) - call Vector_constructor ( CIcore_instance%DDCISDTiming, 30, 0.0_8) + write(*,*) "" + write(*,*) "Diagonalizing hamiltonian..." + write(*,*) " Using : ", trim(String_getUppercase((CONTROL_instance%CI_DIAGONALIZATION_METHOD))) + write(*,*) "=============================================================" + write(*,*) "M. BOLLHÖFER AND Y. NOTAY, JADAMILU:" + write(*,*) " a software code for computing selected eigenvalues of " + write(*,*) " large sparse symmetric matrices, " + write(*,*) "Computer Physics Communications, vol. 177, pp. 951-964, 2007." + write(*,*) "=============================================================" 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,"(T2,A,F14.5,A3 )") "Estimated memory needed: ", & + float(CIcore_instance%numberOfConfigurations*( 2 + (3*ms + CONTROL_instance%NUMBER_OF_CI_STATES + 1) + 4*ms*ms)*8)/(1024**3) , " GB" write (6,*) "" - ecorr = 0.0_8 + !! diagonal correction. See 10.1016/j.chemphys.2007.07.001 + if ( CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT == "CISD") then - do i = 2, 31 + 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 + + 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 + !! 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) + 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 + !! 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 + 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 + 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. + !! 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 + if ( abs( CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i) ) <= 1e-6) exit - end do + 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 + 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 + if ( CONTROL_instance%CI_SAVE_EIGENVECTOR ) then + call CImod_saveEigenVector () + end if - else !! standard CI, no diagonal correction + 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 ) + 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 () + if ( CONTROL_instance%CI_SAVE_EIGENVECTOR ) then + call CImod_saveEigenVector () + end if + end if - end if + case ("DSYEVX") + + write (*,*) "Building Strings..." + call CIStrings_buildStrings() - case ("DSYEVX") + write (*,*) "Building CI level table..." + call CIOrder_buildCIOrderList() - write (*,*) "Building Strings..." - call CIStrings_buildStrings() + write (*,*) "Building diagonal..." + call CIDiag_buildDiagonal() - write (*,*) "Building CI level table..." - call CIOrder_buildCIOrderList() + call Matrix_constructor (CIcore_instance%eigenVectors, & + int(CIcore_instance%numberOfConfigurations,8), & + int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) - write (*,*) "Building diagonal..." - call CIDiag_buildDiagonal() + write(*,*) "" + write(*,*) "Diagonalizing hamiltonian..." + write(*,*) " Using : ", trim(String_getUppercase((CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - call Matrix_constructor (CIcore_instance%eigenVectors, & - int(CIcore_instance%numberOfConfigurations,8), & - int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + write (6,*) "" + write (6,"(T2,A,F14.5,A3 )") "Estimated memory needed: ", & + float((CIcore_instance%numberOfConfigurations**2 + 2 )*8)/(1024**3) , " GB" + write (6,*) "" - write(*,*) "" - write(*,*) "Diagonalizing hamiltonian..." - write(*,*) " Using : ", trim(String_getUppercase((CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - !! diagonal correction. See 10.1016/j.chemphys.2007.07.001 - if ( CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT == "CISD") then + !! 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%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) " + 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 + ecorr = 0.0_8 - do i = 2, 31 + do i = 2, 31 - call CIFullMatrix_buildHamiltonianMatrix( timeA, timeB) + 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 + do a = 2, CIcore_instance%numberOfConfigurations + CIcore_instance%hamiltonianMatrix%values(a,a) = CIcore_instance%hamiltonianMatrix%values(a,a) + ecorr + end do - call Matrix_eigen_select (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & - int(1), int(CONTROL_instance%NUMBER_OF_CI_STATES), & - eigenVectors = CIcore_instance%eigenVectors, & - flags = int(SYMMETRIC,4)) + call Matrix_eigen_select (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & + int(1), int(CONTROL_instance%NUMBER_OF_CI_STATES), & + eigenVectors = CIcore_instance%eigenVectors, & + flags = int(SYMMETRIC,4)) - ecorr = CIcore_instance%eigenvalues%values(1) - HartreeFock_instance%totalEnergy - CIcore_instance%groundStateEnergies%values(i) = CIcore_instance%eigenvalues%values(1) + 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 + 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 + if ( abs( CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i) ) <= 1e-6) exit - end do + end do - else !! standard CI, 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)" + call CIFullMatrix_buildHamiltonianMatrix(timeA, timeB) +!$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for building Hamiltonian Matrix : ", timeB - timeA ," (s)" - call Matrix_eigen_select (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & - int(1), int(CONTROL_instance%NUMBER_OF_CI_STATES), & - eigenVectors = CIcore_instance%eigenVectors, & - flags = int(SYMMETRIC,4)) - - end if - - !! deallocate transformed integrals - deallocate(CIcore_instance%twoCenterIntegrals) - deallocate(CIcore_instance%fourCenterIntegrals) + call Matrix_eigen_select (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & + int(1), int(CONTROL_instance%NUMBER_OF_CI_STATES), & + eigenVectors = CIcore_instance%eigenVectors, & + flags = int(SYMMETRIC,4)) - case ("DSYEVR") + end if - write (*,*) "Building Strings..." - call CIStrings_buildStrings() + !! deallocate transformed integrals + deallocate(CIcore_instance%twoCenterIntegrals) + deallocate(CIcore_instance%fourCenterIntegrals) - write (*,*) "Building CI level table..." - call CIOrder_buildCIOrderList() + case ("DSYEVR") - write (*,*) "Building diagonal..." - call CIDiag_buildDiagonal() + write (*,*) "Building Strings..." + call CIStrings_buildStrings() - call Matrix_constructor (CIcore_instance%eigenVectors, & - int(CIcore_instance%numberOfConfigurations,8), & - int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + write (*,*) "Building CI level table..." + call CIOrder_buildCIOrderList() - !! diagonal correction. See 10.1016/j.chemphys.2007.07.001 - if ( CONTROL_instance%CI_DIAGONAL_DRESSED_SHIFT == "CISD") then + write (*,*) "Building diagonal..." + call CIDiag_buildDiagonal() - call Vector_constructor ( CIcore_instance%groundStateEnergies, 30, 0.0_8) + call Matrix_constructor (CIcore_instance%eigenVectors, & + int(CIcore_instance%numberOfConfigurations,8), & + int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + write(*,*) "" + write(*,*) "Diagonalizing hamiltonian..." + write(*,*) " Using : ", trim(String_getUppercase((CONTROL_instance%CI_DIAGONALIZATION_METHOD))) 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,"(T2,A,F14.5,A3 )") "Estimated memory needed: ", & + float((CIcore_instance%numberOfConfigurations**2 + 2 )*8)/(1024**3) , " GB" write (6,*) "" - write (6,"(T2,A95 )") "Iter Ground-State Energy Correlation Energy Energy Diff. Time(s) " - ecorr = 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) " - do i = 2, 31 + ecorr = 0.0_8 + + do i = 2, 31 - call CIFullMatrix_buildHamiltonianMatrix( timeA, timeB) + 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 + 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 Matrix_eigen_dsyevr (CIcore_instance%hamiltonianMatrix, CIcore_instance%eigenvalues, & + 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) - - 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 + end if - if ( abs( CIcore_instance%groundStateEnergies%values(i-1) - CIcore_instance%groundStateEnergies%values(i) ) <= 1e-6) exit + !! deallocate transformed integrals + deallocate(CIcore_instance%twoCenterIntegrals) + deallocate(CIcore_instance%fourCenterIntegrals) - end do + case default - else !! standard CI, no diagonal correction + call CImod_exception( ERROR, "CImod run", "Diagonalization method not implemented") - 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 + end select + + + elseif ( CONTROL_instance%CONFIGURATION_INTERACTION_LEVEL == "SCI" ) then - !! deallocate transformed integrals - deallocate(CIcore_instance%twoCenterIntegrals) - deallocate(CIcore_instance%fourCenterIntegrals) + select case (trim(String_getUppercase(CONTROL_instance%CI_DIAGONALIZATION_METHOD))) - case default + case ("JADAMILU") - call CImod_exception( ERROR, "CImod run", "Diagonalization method not implemented") + write (*,*) "Building Strings..." + call CIStrings_buildStrings() + write (*,*) "Building CI level table..." + call CIOrder_buildCIOrderList() - end select + call CIJadamilu_buildCouplingMatrix() + call CIJadamilu_buildCouplingOrderList() + + write (*,*) "Building diagonal..." + call CIDiag_buildDiagonal() + + call CISCI_show() + + write (*,*) "Allocating arrays for SCI ..." + call CISCI_constructor() + + call Matrix_constructor (CIcore_instance%eigenVectors, & + int(CIcore_instance%numberOfConfigurations,8), & + int(CONTROL_instance%NUMBER_OF_CI_STATES,8), 0.0_8) + + call CISCI_run() + + case default + + call CImod_exception( ERROR, "CImod run", "Diagonalization method not implemented for SCI") + + end select + + end if write(*,*) "" write(*,*) "-----------------------------------------------" @@ -790,6 +852,15 @@ subroutine CImod_show() write (6,"(T8,A19, F25.12)") "\delta E(Q) = ", davidsonCorrection write (6,"(T8,A19, F25.12)") "E(CISDTQ) ESTIMATE ", HartreeFock_instance%totalEnergy +& CIcorrection + davidsonCorrection + + else if ( CIcore_instance%level == "SCI" ) then + + write(*,"(A)") "" + write (6,"(T2,A34)") "EPSTEIN-NESBET PT2 CORRECTION:" + write(*,"(A)") "" + write (6,"(T8,A19, F25.12)") "E_PT2 :", CISCI_instance%PT2energy + write (6,"(T8,A19, F25.12)") "E_SCI + E_PT2 :", CIcore_instance%eigenvalues%values(1) + CISCI_instance%PT2energy + else write(*,"(A)") "" @@ -798,8 +869,6 @@ subroutine CImod_show() end if - else - end if end subroutine CImod_show diff --git a/src/core/CONTROL.f90 b/src/core/CONTROL.f90 index a09c0d8..b85c6dd 100644 --- a/src/core/CONTROL.f90 +++ b/src/core/CONTROL.f90 @@ -205,6 +205,8 @@ module CONTROL_ logical :: CI_BUILD_FULL_MATRIX integer :: CI_MADSPACE logical :: CI_NATURAL_ORBITALS + integer :: CI_SCI_CORE_SPACE + integer :: CI_SCI_TARGET_SPACE !!*************************************************************************** !! Non-orthogonal CI @@ -545,6 +547,8 @@ module CONTROL_ logical :: LowdinParameters_CIBuildFullMatrix integer :: LowdinParameters_CIMadSpace logical :: LowdinParameters_CINaturalOrbitals + integer :: LowdinParameters_CISCICoreSpace + integer :: LowdinParameters_CISCITargetSpace !!*************************************************************************** !! Non-orthogonal CI @@ -884,6 +888,10 @@ module CONTROL_ LowdinParameters_CINaturalOrbitals, & LowdinParameters_CIPrintEigenVectorsFormat, & LowdinParameters_CIPrintThreshold, & + LowdinParameters_CISCICoreSpace, & + LowdinParameters_CISCITargetSpace, & + + !!*************************************************************************** !! Non-orthogonal CI @@ -1585,6 +1593,8 @@ subroutine CONTROL_start() CONTROL_instance%CI_NATURAL_ORBITALS=.FALSE. CONTROL_instance%CI_PRINT_EIGENVECTORS_FORMAT = "OCCUPIED" CONTROL_instance%CI_PRINT_THRESHOLD = 1E-1 + CONTROL_instance%CI_SCI_CORE_SPACE = 100 + CONTROL_instance%CI_SCI_TARGET_SPACE = 10000 !!*************************************************************************** !! Non-orthogonal CI @@ -1974,6 +1984,10 @@ subroutine CONTROL_load(unit) CONTROL_instance%CI_NATURAL_ORBITALS= LowdinParameters_CINaturalOrbitals CONTROL_instance%CI_PRINT_EIGENVECTORS_FORMAT = LowdinParameters_CIPrintEigenVectorsFormat CONTROL_instance%CI_PRINT_THRESHOLD = LowdinParameters_CIPrintThreshold + CONTROL_instance%CI_SCI_CORE_SPACE = LowdinParameters_CISCICoreSpace + CONTROL_instance%CI_SCI_TARGET_SPACE = LowdinParameters_CISCITargetSpace + + !!*************************************************************************** !! Non-orthogonal CI @@ -2336,9 +2350,10 @@ subroutine CONTROL_save( unit, lastStep, firstStep ) LowdinParameters_CIBuildFullMatrix = CONTROL_instance%CI_BUILD_FULL_MATRIX LowdinParameters_CIMadSpace = CONTROL_instance%CI_MADSPACE LowdinParameters_CINaturalOrbitals = CONTROL_instance%CI_NATURAL_ORBITALS - LowdinParameters_CIPrintEigenVectorsFormat = CONTROL_instance%CI_PRINT_EIGENVECTORS_FORMAT LowdinParameters_CIPrintThreshold = CONTROL_instance%CI_PRINT_THRESHOLD + LowdinParameters_CISCICoreSpace = CONTROL_instance%CI_SCI_CORE_SPACE + LowdinParameters_CISCITargetSpace = CONTROL_instance%CI_SCI_TARGET_SPACE !!*************************************************************************** !! Non-orthogonal CI @@ -2673,6 +2688,8 @@ subroutine CONTROL_copy(this, otherThis) otherThis%CI_NATURAL_ORBITALS = this%CI_NATURAL_ORBITALS otherThis%CI_PRINT_EIGENVECTORS_FORMAT = this%CI_PRINT_EIGENVECTORS_FORMAT otherThis%CI_PRINT_THRESHOLD = this%CI_PRINT_THRESHOLD + otherThis%CI_SCI_CORE_SPACE = this%CI_SCI_CORE_SPACE + otherThis%CI_SCI_TARGET_SPACE = this%CI_SCI_TARGET_SPACE !!*************************************************************************** !! Non-orthogonal CI diff --git a/src/core/Vector.f90 b/src/core/Vector.f90 index b199c8a..e74a7f9 100644 --- a/src/core/Vector.f90 +++ b/src/core/Vector.f90 @@ -91,6 +91,8 @@ module Vector_ Vector_reverseSortElements, & Vector_reverseSortElements8, & Vector_reverseSortElements8Int, & + Vector_reverseSortElementsAbsolute8, & + Vector_sortElementsAbsolute8, & Vector_swapElements, & Vector_getSize, & Vector_getElement, & @@ -1512,6 +1514,114 @@ subroutine Vector_reverseSortElements8(this,indexVector,m) end subroutine Vector_reverseSortElements8 + subroutine Vector_reverseSortElementsAbsolute8(this,indexVector,m) + type(Vector8) :: this + type(IVector8), optional :: indexVector + integer(8), optional :: m + integer(8) i,j,n + + n = Vector_getSize8(this) + if ( .not. present (indexVector) ) then + do i=1,n + do j=i+1,n + if ( abs(this%values(j)).lt. abs(this%values(i)) ) then + call Vector_swapElements8( this, i, j ) + end if + end do + end do + else + + if ( .not. present (m) ) then + + do i=1,n + indexVector%values(i) = i + end do + + do i=1,n + do j=i+1,n + if ( abs(this%values(j)).lt. abs(this%values(i)) ) then + call Vector_swapElements8( this, i, j ) + call Vector_swapIntegerElements8( indexVector, i, j ) + end if + end do + end do + else + + do i=1,n + indexVector%values(i) = i + end do + + do i=1,m + do j=i+1,n + if ( abs(this%values(j)).lt.abs(this%values(i))) then + call Vector_swapElements8( this, i, j ) + call Vector_swapIntegerElements8( indexVector, i, j ) + end if + end do + end do + end if + end if + + end subroutine Vector_reverseSortElementsAbsolute8 + + subroutine Vector_sortElementsAbsolute8(this,indexVector,m) + type(Vector8) :: this + type(IVector8), optional :: indexVector + integer(8), optional :: m + integer(8) i,j,n + real(8) :: timeA, timeB + +!$ timeA = omp_get_wtime() + + n = Vector_getSize8(this) + if ( .not. present (indexVector) ) then + do i=1,n + do j=i+1,n + if ( abs(this%values(j)).gt. abs(this%values(i)) ) then + call Vector_swapElements8( this, i, j ) + end if + end do + end do + else + + if ( .not. present (m) ) then + + do i=1,n + indexVector%values(i) = i + end do + + do i=1,n + do j=i+1,n + if ( abs(this%values(j)).gt. abs(this%values(i)) ) then + call Vector_swapElements8( this, i, j ) + call Vector_swapIntegerElements8( indexVector, i, j ) + end if + end do + end do + else + + !do i=1,n + ! indexVector%values(i) = i + !end do + + do i=1,m + do j=i+1,n + if ( abs(this%values(j)).gt.abs(this%values(i))) then + call Vector_swapElements8( this, i, j ) + call Vector_swapIntegerElements8( indexVector, i, j ) + end if + end do + end do + end if + end if + +!$ timeB = omp_get_wtime() +!$ write(*,"(A,E10.3,A4)") "** TOTAL Elapsed Time for sorting the vector : ", timeB - timeA ," (s)" + + end subroutine Vector_sortElementsAbsolute8 + + + subroutine Vector_reverseSortElements8Int(this,indexVector,m) type(IVector8) :: this type(IVector8), optional :: indexVector diff --git a/test/PsH.SCI.lowdin b/test/PsH.SCI.lowdin new file mode 100644 index 0000000..77ee9bc --- /dev/null +++ b/test/PsH.SCI.lowdin @@ -0,0 +1,35 @@ +GEOMETRY + e-(H) SHARON-E-6S2P 0.00 0.00 0.00 addParticles=1 + H dirac 0.00 0.00 0.00 + e+ SHARON-E+6S2P 0.00 0.00 0.00 +END GEOMETRY + +TASKS + method = "UHF" + configurationInteractionLevel ="SCI" + !configurationInteractionLevel ="CISD" +END TASKS + +CONTROL + readCoefficients=F + numberOfCIstates=1 + CINaturalOrbitals=T + CIStatesToPrint = 1 + CIdiagonalizationMethod = "JADAMILU" + !CIPrintEigenVectorsFormat = "NONE" + !CIPrintEigenVectorsFormat = "OCCUPIED" + !CIPrintEigenVectorsFormat = "ORBITALS" + !CIPrintThreshold = 5e-2 + buildTwoParticlesMatrixForOneParticle=T + CISCICoreSpace = 100 + CISCITargetSpace = 500 + CIMatVecTolerance = 1e-10 + HFprintEigenVectors = "ALL" +END CONTROL + +INPUT_CI + species="E-ALPHA" core=0 active=0 + species="E-BETA" core=0 active=0 + species="POSITRON" core=0 active=0 +END INPUT_CI + diff --git a/test/PsH.SCI.py b/test/PsH.SCI.py new file mode 100644 index 0000000..a98410a --- /dev/null +++ b/test/PsH.SCI.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python +from __future__ import print_function +import os +import sys +from colorstring import * + +if len(sys.argv)==2: + lowdinbin = sys.argv[1] +else: + lowdinbin = "lowdin2" + +testName = sys.argv[0][:-3] +inputName = testName + ".lowdin" +outputName = testName + ".out" + +# Reference values and tolerance + +refValues = { +"HF energy" : [-0.666783062050,1E-8], +"E_SCI+PT2" : [-0.743335783194,1E-6], +} + +testValues = dict(refValues) #copy +for value in testValues: #reset + testValues[value] = 0 #reset + +# Run calculation + +status = os.system(lowdinbin + " -i " + inputName) + +if status: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output = open(outputName, "r") +outputRead = output.readlines() + +# Values +for i in range(0,len(outputRead)): + line = outputRead[i] + if "TOTAL ENERGY =" in line: + testValues["HF energy"] = float(line.split()[3]) + if "E_SCI + E_PT2 :" in line: + testValues["E_SCI+PT2"] = float(line.split()[4]) + + +passTest = True + +for value in refValues: + diffValue = abs(refValues[value][0] - testValues[value]) + if ( diffValue <= refValues[value][1] ): + passTest = passTest * True + else : + passTest = passTest * False + print("%s %.8f %.8f %.2e" % ( value, refValues[value][0], testValues[value], diffValue)) + +if passTest : + print(testName + str_green(" ... OK")) +else: + print(testName + str_red(" ... NOT OK")) + sys.exit(1) + +output.close() diff --git a/test/runtest.sh b/test/runtest.sh index aaedd1e..f91fdf9 100644 --- a/test/runtest.sh +++ b/test/runtest.sh @@ -9,7 +9,7 @@ fi echo $lowdinbin for testfile in `ls *.py`; do - echo $testfile + #echo $testfile python3 $testfile $lowdinbin status=$((status + $?)) done