From f6270a0ba5dccf09897efd607c17b1d98529205e Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 4 Jan 2022 11:39:33 +0100 Subject: [PATCH] CID and CISD --- input/methods | 4 +- input/options | 2 +- mol/h2.xyz | 2 +- src/CI/CID.f90 | 218 ++++++++++++++++++ src/CI/CISD.f90 | 342 +++++++++++++++++++---------- src/GT/dynamic_Tmatrix_A.f90 | 12 +- src/GT/static_Tmatrix_A.f90 | 4 +- src/GT/static_Tmatrix_B.f90 | 4 +- src/HF/RHF.f90 | 6 +- src/QuAcK/QuAcK.f90 | 14 +- src/utils/spatial_to_spin_fock.f90 | 29 +++ 11 files changed, 497 insertions(+), 140 deletions(-) create mode 100644 src/CI/CID.f90 create mode 100644 src/utils/spatial_to_spin_fock.f90 diff --git a/input/methods b/input/methods index 259d3f4..9f505d6 100644 --- a/input/methods +++ b/input/methods @@ -7,7 +7,7 @@ # drCCD rCCD crCCD lCCD F F F F # CIS* CIS(D) CID CISD FCI - F F F F F + F F T T F # RPA* RPAx* crRPA ppRPA F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 @@ -15,7 +15,7 @@ # G0W0* evGW* qsGW* ufG0W0 ufGW F F F F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index e8b50b8..d6417cb 100644 --- a/input/options +++ b/input/options @@ -15,6 +15,6 @@ # ACFDT: AC Kx XBS F F F # BSE: BSE dBSE dTDA evDyn - T F T F + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index bb00204..a4e936a 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 1.0 +H 0. 0. 2.000000 diff --git a/src/CI/CID.f90 b/src/CI/CID.f90 new file mode 100644 index 0000000..ca6a78d --- /dev/null +++ b/src/CI/CID.f90 @@ -0,0 +1,218 @@ +subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0) + +! Perform configuration interaction with doubles + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold + integer,intent(in) :: nBasin + integer,intent(in) :: nCin + integer,intent(in) :: nOin + integer,intent(in) :: nVin + integer,intent(in) :: nRin + double precision,intent(in) :: Fin(nBasin,nBasin) + double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin) + double precision,intent(in) :: E0 + +! Local variables + + integer :: nBas + integer :: nC + integer :: nO + integer :: nV + integer :: nR + + double precision,allocatable :: F(:,:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: ERI(:,:,:,:) + + logical :: dump_trans = .false. + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,kc,iajb,kcld + integer :: ishift,jshift + integer :: ispin + integer :: nD + integer :: nH + integer :: maxH + double precision,external :: Kronecker_delta + double precision,allocatable :: H(:,:) + double precision,allocatable :: ECID(:) + + double precision :: tmp + +! Hello world + + write(*,*) + write(*,*)'******************************************************' + write(*,*)'| Configuration Interaction with Singles and Doubles |' + write(*,*)'******************************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas = 2*nBasin + nC = 2*nCin + nO = 2*nOin + nV = 2*nVin + nR = 2*nRin + + allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas)) + + call spatial_to_spin_fock(nBasin,Fin,nBas,F) + call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI) + +! Antysymmetrize ERIs + + allocate(ERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,sERI,ERI) + + deallocate(sERI) + +! Compute CID matrix + + nD = (nO - nC)*(nO - nC - 1)/2*(nV - nR)*(nV - nR - 1)/2 + nH = 1 + nD + + write(*,*) 'nD = ',nD + write(*,*) 'nH = ',nH + write(*,*) + + maxH = min(nH,21) + + ! Memory allocation + + allocate(H(nH,nH),ECID(nH)) + + ! 00 block + + ishift = 0 + jshift = 0 + + H(ishift+1,jshift+1) = E0 + + print*,'00 block done...' + + ! 0D blocks + + ishift = 0 + jshift = 1 + + iajb = 0 + do i=nC+1,nO + do a=1,nV-nR + do j=i+1,nO + do b=a+1,nV-nR + + iajb = iajb + 1 + tmp = ERI(i,j,nO+a,nO+b) + + H(ishift+1,jshift+iajb) = tmp + H(jshift+iajb,ishift+1) = tmp + + end do + end do + end do + end do + + print*,'0D blocks done...' + + ! DD block + + ishift = 1 + jshift = 1 + + iajb = 0 + do i=nC+1,nO + do a=1,nV-nR + do j=i+1,nO + do b=a+1,nV-nR + + iajb = iajb + 1 + + kcld = 0 + do k=nC+1,nO + do c=1,nV-nR + do l=k+1,nO + do d=c+1,nV-nR + + kcld = kcld + 1 + tmp = & + E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + - F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & + - F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & + + F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + - F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & + + F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + + F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + - F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & + + F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + - F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + - F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + - F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + + F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + + F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + - F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + - ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) & + + ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + - ERI(nO+a,l,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & + - ERI(nO+a,k,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & + + ERI(nO+a,k,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + - ERI(nO+a,l,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & + + ERI(nO+a,l,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + + ERI(nO+a,k,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + - ERI(nO+a,k,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & + - ERI(nO+b,l,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,k) & + + ERI(nO+b,l,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,k) & + + ERI(nO+b,k,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,l) & + - ERI(nO+b,k,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,l) & + + ERI(nO+b,l,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,k) & + - ERI(nO+b,l,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,k) & + - ERI(nO+b,k,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,l) & + + ERI(nO+b,k,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,l) & + - ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + + ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) + + H(ishift+iajb,jshift+kcld) = tmp + + end do + end do + end do + end do + + end do + end do + end do + end do + + print*,'DD block done...' + + write(*,*) + write(*,*) 'Diagonalizing CID matrix...' + write(*,*) + + call diagonalize_matrix(nH,H,ECID) + + print*,'CID energies (au)' + call matout(maxH,1,ECID) + write(*,*) + + print*,'CID excitation energies (eV)' + call matout(maxH-1,1,(ECID(2:maxH)-ECID(1))*HaToeV) + write(*,*) + + if(dump_trans) then + print*,'Singlet CID transition vectors' + call matout(nH,nH,H) + write(*,*) + endif + +end subroutine CID diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90 index 611d940..1037d11 100644 --- a/src/CI/CISD.f90 +++ b/src/CI/CISD.f90 @@ -1,4 +1,4 @@ -subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF) +subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0) ! Perform configuration interaction with singles and doubles @@ -9,23 +9,42 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF) logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + integer,intent(in) :: nBasin + integer,intent(in) :: nCin + integer,intent(in) :: nOin + integer,intent(in) :: nVin + integer,intent(in) :: nRin + double precision,intent(in) :: Fin(nBasin,nBasin) + double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin) + double precision,intent(in) :: E0 ! Local variables + integer :: nBas + integer :: nC + integer :: nO + integer :: nV + integer :: nR + + double precision,allocatable :: F(:,:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: ERI(:,:,:,:) + logical :: dump_trans = .false. integer :: i,j,k,l integer :: a,b,c,d - integer :: ia,jb,iajb,kcld + integer :: ia,kc,iajb,kcld integer :: ishift,jshift integer :: ispin integer :: nS integer :: nD - integer :: nSD + integer :: nH + integer :: maxH double precision,external :: Kronecker_delta - double precision,allocatable :: H(:,:),Omega(:) + double precision,allocatable :: H(:,:) + double precision,allocatable :: ECISD(:) + + double precision :: tmp ! Hello world @@ -35,170 +54,257 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI,eHF) write(*,*)'******************************************************' write(*,*) -! Compute CIS matrix +! Spatial to spin orbitals - if(singlet_manifold) then + nBas = 2*nBasin + nC = 2*nCin + nO = 2*nOin + nV = 2*nVin + nR = 2*nRin - ispin = 1 + allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas)) - ! Dimensions + call spatial_to_spin_fock(nBasin,Fin,nBas,F) + call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI) - nS = (nO - nC)*(nV - nR) - nD = (nO - nC)*(nO - nC + 1)/2*(nV - nR)*(nV - nR + 1)/2 - nSD = 1 + nS + nD +! Antysymmetrize ERIs - print*,'nS = ',nS - print*,'nD = ',nD - print*,'nSD = ',nSD + allocate(ERI(nBas,nBas,nBas,nBas)) - ! Memory allocation + call antisymmetrize_ERI(2,nBas,sERI,ERI) - allocate(H(nSD,nSD),Omega(nSD)) - - ! 0D block + deallocate(sERI) - ishift = 0 - jshift = 1 + nS +! Compute CISD matrix - iajb = 0 + nS = (nO - nC)*(nV - nR) + nD = (nO - nC)*(nO - nC - 1)/2*(nV - nR)*(nV - nR - 1)/2 + nH = 1 + nS + nD - do i=nC+1,nO - do a=1,nV-nR - do j=i,nO - do b=a,nV-nR + write(*,*) 'nS = ',nS + write(*,*) 'nD = ',nD + write(*,*) 'nH = ',nH + write(*,*) - iajb = iajb + 1 - H(ishift+1,jshift+iajb) = ERI(i,j,nO+a,nO+b) + maxH = min(nH,21) + + ! Memory allocation + + allocate(H(nH,nH),ECISD(nH)) + + ! 00 block + + ishift = 0 + jshift = 0 + + H(ishift+1,jshift+1) = E0 + + print*,'00 block done...' + + ! 0S blocks + + ishift = 0 + jshift = 1 + + ia = 0 + do i=nC+1,nO + do a=1,nV-nR + + ia = ia + 1 + tmp = F(i,nO+a) + H(ishift+1,jshift+ia) = tmp + H(jshift+ia,ishift+1) = tmp + + end do + end do + + print*,'0S blocks done...' + + ! 0D blocks + + ishift = 0 + jshift = 1 + nS + + iajb = 0 + do i=nC+1,nO + do a=1,nV-nR + do j=i+1,nO + do b=a+1,nV-nR + + iajb = iajb + 1 + tmp = ERI(i,j,nO+a,nO+b) + + H(ishift+1,jshift+iajb) = tmp + H(jshift+iajb,ishift+1) = tmp - end do end do end do end do - - ! SS block + end do + + print*,'0D blocks done...' - ishift = 1 - jshift = 1 + ! SS block - ia = 0 - jb = 0 + ishift = 1 + jshift = 1 - do i=nC+1,nO - do a=1,nV-nR + ia = 0 + do i=nC+1,nO + do a=1,nV-nR - ia = ia + 1 + ia = ia + 1 + kc = 0 + do k=nC+1,nO + do c=1,nV-nR - do j=nC+1,nO - do b=1,nV-nR + kc = kc + 1 + tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + - F(i,k)*Kronecker_delta(a,c) & + + F(nO+a,nO+c)*Kronecker_delta(i,k) & + - ERI(nO+a,k,nO+c,i) - jb = jb + 1 + H(ishift+ia,jshift+kc) = tmp - H(ishift+ia,jshift+jb) & - = Kronecker_delta(i,j)*Kronecker_delta(a,b)*(eHF(nO+a) - eHF(i)) & - + ERI(nO+a,j,i,nO+b) - ERI(nO+a,j,nO+b,i) - - end do end do - end do + end do - - ! SD block + end do - ishift = 1 - jshift = 1 + nS + print*,'SS block done...' - ia = 0 - kcld = 0 + ! SD blocks - do i=nC+1,nO - do a=1,nV-nR + ishift = 1 + jshift = 1 + nS - ia = ia + 1 + ia = 0 + do i=nC+1,nO + do a=1,nV-nR - do k=nC+1,nO - do c=1,nV-nR - do l=k,nO - do d=c,nV-nR + ia = ia + 1 + kcld = 0 - kcld = kcld + 1 + do k=nC+1,nO + do c=1,nV-nR + do l=k+1,nO + do d=c+1,nV-nR - H(ishift+ia,jshift+kcld) & - = Kronecker_delta(i,k)*(ERI(nO+a,l,nO+c,nO+d) - ERI(nO+a,l,nO+d,nO+c)) & - - Kronecker_delta(i,l)*(ERI(nO+a,k,nO+c,nO+d) - ERI(nO+a,k,nO+d,nO+c)) & - - Kronecker_delta(a,c)*(ERI(k,l,i,nO+d) - ERI(k,l,nO+d,i)) & - + Kronecker_delta(a,d)*(ERI(k,l,i,nO+c) - ERI(k,l,nO+c,i)) + kcld = kcld + 1 + tmp = - F(l,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k) & + + F(l,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k) & + - F(k,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l) & + + F(k,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l) & + - ERI(k,l,nO+d,i)*Kronecker_delta(a,c) & + + ERI(k,l,nO+c,i)*Kronecker_delta(a,d) & + - ERI(nO+a,l,nO+c,nO+d)*Kronecker_delta(i,k) & + + ERI(nO+a,k,nO+c,nO+d)*Kronecker_delta(i,l) + + H(ishift+ia,jshift+kcld) = tmp + H(jshift+kcld,ishift+ia) = tmp - end do end do end do end do - end do + end do + end do - ! DD block + print*,'SD blocks done...' - ishift = 1 + nS - jshift = 1 + nS + ! DD block - iajb = 0 - kcld = 0 + ishift = 1 + nS + jshift = 1 + nS - do i=nC+1,nO - do a=1,nV-nR - do j=i,nO - do b=a,nV-nR + iajb = 0 + do i=nC+1,nO + do a=1,nV-nR + do j=i+1,nO + do b=a+1,nV-nR - iajb = iajb + 1 + iajb = iajb + 1 - do k=nC+1,nO - do c=1,nV-nR - do l=k,nO - do d=c,nV-nR - - kcld = kcld + 1 - -! H(ishift+iajb,jshift+kcld) & -! = Kronecker_delta(i,k)*(ERI(a,l,c,d) - ERI(a,l,d,c)) & -! - Kronecker_delta(i,l)*(ERI(a,k,c,d) - ERI(a,k,d,c)) & -! - Kronecker_delta(a,c)*(ERI(k,l,i,d) - ERI(k,l,d,i)) & -! + Kronecker_delta(a,d)*(ERI(k,l,i,c) - ERI(k,l,c,i)) - - end do + kcld = 0 + do k=nC+1,nO + do c=1,nV-nR + do l=k+1,nO + do d=c+1,nV-nR + + kcld = kcld + 1 + tmp = & + E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + - F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & + - F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & + + F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + - F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & + + F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + + F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + - F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & + + F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + - F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + - F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + - F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + + F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + + F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + - F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + - ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) & + + ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + - ERI(nO+a,l,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & + - ERI(nO+a,k,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & + + ERI(nO+a,k,nO+c,j)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + - ERI(nO+a,l,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & + + ERI(nO+a,l,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + + ERI(nO+a,k,nO+d,i)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + - ERI(nO+a,k,nO+c,i)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & + - ERI(nO+b,l,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,k) & + + ERI(nO+b,l,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,k) & + + ERI(nO+b,k,nO+d,j)*Kronecker_delta(a,c)*Kronecker_delta(i,l) & + - ERI(nO+b,k,nO+c,j)*Kronecker_delta(a,d)*Kronecker_delta(i,l) & + + ERI(nO+b,l,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,k) & + - ERI(nO+b,l,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,k) & + - ERI(nO+b,k,nO+d,i)*Kronecker_delta(a,c)*Kronecker_delta(j,l) & + + ERI(nO+b,k,nO+c,i)*Kronecker_delta(a,d)*Kronecker_delta(j,l) & + - ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + + ERI(nO+a,nO+b,nO+c,nO+d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) + + H(ishift+iajb,jshift+kcld) = tmp + end do end do end do - end do + end do end do end do + end do - call diagonalize_matrix(nSD,H,Omega) - call print_excitation('CISD ',ispin,nS,Omega) - - if(dump_trans) then - print*,'Singlet CISD transition vectors' - call matout(nSD,nSD,H) - write(*,*) - endif + print*,'DD block done...' + write(*,*) + write(*,*) 'Diagonalizing CISD matrix...' + write(*,*) + + call diagonalize_matrix(nH,H,ECISD) + + print*,'CISD energies (au)' + call matout(maxH,1,ECISD) + write(*,*) + + print*,'CISD excitation energies (eV)' + call matout(maxH-1,1,(ECISD(2:maxH)-ECISD(1))*HaToeV) + write(*,*) + + if(dump_trans) then + print*,'Singlet CISD transition vectors' + call matout(nH,nH,H) + write(*,*) endif -! if(triplet_manifold) then - -! ispin = 2 -! -! call diagonalize_matrix(nSD,H,Omega) -! call print_excitation('CISD ',ispin,nSD,Omega) - -! if(dump_trans) then -! print*,'Triplet CIS transition vectors' -! call matout(nSD,nSD,H) -! write(*,*) -! endif - -! endif - end subroutine CISD diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90 index 4cef4c2..9adef18 100644 --- a/src/GT/dynamic_Tmatrix_A.f90 +++ b/src/GT/dynamic_Tmatrix_A.f90 @@ -59,12 +59,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O do cd=1,nVV eps = + Omega1(cd) - chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) + chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2) end do do kl=1,nOO eps = - Omega2(kl) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2) end do A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi @@ -73,12 +73,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O do cd=1,nVV eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) - chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) + chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2) end do do kl=1,nOO eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2) end do A_dyn(ia,jb) = A_dyn(ia,jb) - 1d0*lambda*chi @@ -87,12 +87,12 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O do cd=1,nVV eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) - chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOO eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + 1d0*lambda*chi diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90 index f6a4817..0710ca2 100644 --- a/src/GT/static_Tmatrix_A.f90 +++ b/src/GT/static_Tmatrix_A.f90 @@ -47,13 +47,13 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh do cd=1,nVV eps = + Omega1(cd) ! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) - chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) + chi = chi + rho1(i,a,cd)*rho1(j,b,cd)*eps/(eps**2 + eta**2) enddo do kl=1,nOO eps = - Omega2(kl) ! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + chi = chi + rho2(i,a,kl)*rho2(j,b,kl)*eps/(eps**2 + eta**2) enddo TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi diff --git a/src/GT/static_Tmatrix_B.f90 b/src/GT/static_Tmatrix_B.f90 index 4103ecc..53561be 100644 --- a/src/GT/static_Tmatrix_B.f90 +++ b/src/GT/static_Tmatrix_B.f90 @@ -47,13 +47,13 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rh do cd=1,nVV eps = + Omega1(cd) ! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2 - chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) + chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) enddo do kl=1,nOO eps = - Omega2(kl) ! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 - chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2) + chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) enddo TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 05499a5..a4bf5fe 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,4 +1,4 @@ -subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,e,c,P,Vx) +subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx) ! Perform restricted Hartree-Fock calculation @@ -45,7 +45,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) double precision,allocatable :: cp(:,:) - double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) double precision,allocatable :: ON(:) @@ -56,6 +55,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T double precision,intent(out) :: c(nBas,nBas) double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: Vx(nBas) + double precision,intent(out) :: F(nBas,nBas) ! Hello world @@ -72,7 +72,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Memory allocation allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & - cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas),ON(nBas), & + cp(nBas,nBas),Fp(nBas,nBas),ON(nBas), & error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Guess coefficients and eigenvalues diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ee41859..48f2189 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -64,6 +64,8 @@ program QuAcK double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: dipole_int_aa(:,:,:) double precision,allocatable :: dipole_int_bb(:,:,:) + double precision,allocatable :: F_AO(:,:) + double precision,allocatable :: F_MO(:,:) double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz @@ -247,7 +249,7 @@ program QuAcK allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), & S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), & - dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin)) + dipole_int_AO(nBas,nBas,ncart),dipole_int_MO(nBas,nBas,ncart),Vxc(nBas,nspin),F_AO(nBas,nBas)) ! Read integrals @@ -291,7 +293,7 @@ program QuAcK call cpu_time(start_HF) call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) + nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) call cpu_time(end_HF) t_HF = end_HF - start_HF @@ -433,6 +435,7 @@ program QuAcK ! Memory allocation allocate(ERI_MO(nBas,nBas,nBas,nBas)) + allocate(F_MO(nBas,nBas)) ! Read and transform dipole-related integrals @@ -448,7 +451,8 @@ program QuAcK ket1 = 1 ket2 = 1 call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO) -! call AOtoMO_transform(nBas,cHF,T+V) + F_MO(:,:) = F_AO(:,:) + call AOtoMO_transform(nBas,cHF,F_MO) end if end if @@ -734,7 +738,7 @@ program QuAcK if(doCID) then call cpu_time(start_CID) -! call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF) + call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF) call cpu_time(end_CID) t_CID = end_CID - start_CID @@ -750,7 +754,7 @@ program QuAcK if(doCISD) then call cpu_time(start_CISD) -! call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF) + call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF) call cpu_time(end_CISD) t_CISD = end_CISD - start_CISD diff --git a/src/utils/spatial_to_spin_fock.f90 b/src/utils/spatial_to_spin_fock.f90 new file mode 100644 index 0000000..614c693 --- /dev/null +++ b/src/utils/spatial_to_spin_fock.f90 @@ -0,0 +1,29 @@ +subroutine spatial_to_spin_fock(nBas,F,nBas2,sF) + +! Convert Fock matrix from spatial to spin orbitals + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nBas2 + double precision,intent(in) :: F(nBas,nBas) + +! Local variables + + integer :: p,q + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: sF(nBas2,nBas2) + + do p=1,nBas2 + do q=1,nBas2 + + sF(p,q) = Kronecker_delta(mod(p,2),mod(q,2))*F((p+1)/2,(q+1)/2) + + enddo + enddo + +end subroutine spatial_to_spin_fock