diff --git a/input/methods b/input/methods index c22322a..1e67d49 100644 --- a/input/methods +++ b/input/methods @@ -2,14 +2,14 @@ T F F F # MP2* MP3 MP2-F12 F F F -# CCD DCD CCSD CCSD(T) - F F F F -# drCCD rCCD lCCD pCCD - T T F F +# CCD pCCD DCD CCSD CCSD(T) + F F F F F +# drCCD rCCD crCCD lCCD + F F F F # CIS* CIS(D) CID CISD FCI F F F F F -# RPA* RPAx* ppRPA - T T F +# RPA* RPAx* crRPA ppRPA + F F F T # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW diff --git a/input/options b/input/options index a4093ba..f54eee2 100644 --- a/input/options +++ b/input/options @@ -11,8 +11,8 @@ # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 256 0.00001 T 5 T 0.00367493 F F F F F # ACFDT: AC Kx XBS - F F T + T T F # BSE: BSE dBSE dTDA evDyn - T T T F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/phenol.xyz b/mol/phenol.xyz new file mode 100644 index 0000000..f7871b3 --- /dev/null +++ b/mol/phenol.xyz @@ -0,0 +1,15 @@ +13 +Phenol; structure form HCP92 revised structure; l +C 4.555420 5.661760 4.489060 +C 4.584960 2.843420 4.498910 +C 3.378760 3.548270 4.498890 +C 3.359200 4.943960 4.500070 +C 5.758650 4.948670 4.502800 +C 5.789840 3.550210 4.500330 +H 4.622340 1.760920 4.495600 +H 3.631910 7.321310 4.523190 +H 2.474610 2.963630 4.498170 +H 2.384290 5.419790 4.496670 +H 6.695040 5.492670 4.498160 +H 6.727200 3.021210 4.497160 +O 4.537700 7.024010 4.500450 diff --git a/src/AOtoMO/AOtoMO_integral_transform.f90 b/src/AOtoMO/AOtoMO_integral_transform.f90 index dabeee8..fbb0d4b 100644 --- a/src/AOtoMO/AOtoMO_integral_transform.f90 +++ b/src/AOtoMO/AOtoMO_integral_transform.f90 @@ -79,7 +79,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI do nu=1,nBas ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l) enddo -! print*,i,k,j,l,ERI_MO_basis(i,j,k,l) +! write(11,'(I5,I5,I5,I5,F16.10)') i,j,k,l,ERI_MO_basis(i,j,k,l) enddo enddo enddo diff --git a/src/AOtoMO/AOtoMO_transform.f90 b/src/AOtoMO/AOtoMO_transform.f90 index 7919084..a4cd5a3 100644 --- a/src/AOtoMO/AOtoMO_transform.f90 +++ b/src/AOtoMO/AOtoMO_transform.f90 @@ -9,10 +9,19 @@ subroutine AOtoMO_transform(nBas,c,A) integer,intent(in) :: nBas double precision,intent(in) :: c(nBas,nBas) + integer :: i,j + ! Output variables double precision,intent(inout):: A(nBas,nBas) A = matmul(transpose(c),matmul(A,c)) +! do j=1,nBas +! do i=1,nBas +! write(10,'(I5,I5,F16.10)') i,j,A(i,j) +! enddo +! enddo + + end subroutine AOtoMO_transform diff --git a/src/CC/crCCD.f90 b/src/CC/crCCD.f90 index 36cd10d..eb2e0f7 100644 --- a/src/CC/crCCD.f90 +++ b/src/CC/crCCD.f90 @@ -38,21 +38,8 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF double precision,allocatable :: eV(:) double precision,allocatable :: delta_OOVV(:,:,:,:) - double precision,allocatable :: OOOO(:,:,:,:) double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OVOV(:,:,:,:) - double precision,allocatable :: OVVO(:,:,:,:) - double precision,allocatable :: VVVV(:,:,:,:) - - double precision,allocatable :: X1(:,:,:,:) - double precision,allocatable :: X2(:,:) - double precision,allocatable :: X3(:,:) - double precision,allocatable :: X4(:,:,:,:) - - double precision,allocatable :: u(:,:,:,:) - double precision,allocatable :: v(:,:,:,:) - double precision,allocatable :: r2r(:,:,:,:) - double precision,allocatable :: r2l(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) @@ -66,7 +53,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF write(*,*) write(*,*)'**************************************' - write(*,*)'| crossed-ring CCD calculation |' + write(*,*)'| Crossed-ring CCD calculation |' write(*,*)'**************************************' write(*,*) @@ -105,13 +92,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Create integral batches - allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV),OVVO(nO,nV,nV,nO)) + allocate(OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV)) - OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) - OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) - OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas) - OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) - VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas, 1:nO ,nO+1:nBas) deallocate(dbERI) @@ -129,9 +113,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Initialization - allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) - allocate(r2r(nO,nO,nV,nV),r2l(nO,nO,nV,nV)) - allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV)) + allocate(r2(nO,nO,nV,nV)) Conv = 1d0 nSCF = 0 @@ -159,23 +141,9 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Compute residual -! Form linear array + call form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2) - call form_u(nC,nO,nV,nR,OOOO,VVVV,OVOV,t2,u) - -! Form interemediate arrays - - call form_X(nC,nO,nV,nR,OOVV,t2,X1,X2,X3,X4) - -! Form quadratic array - - call form_v(nC,nO,nV,nR,X1,X2,X3,X4,t2,v) - - call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2r) - - call form_ladder_r(nC,nO,nV,nR,OOOO,OOVV,VVVV,t2,r2l) - - r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + u(:,:,:,:) + v(:,:,:,:) - r2r(:,:,:,:) - r2l(:,:,:,:) + r2(:,:,:,:) = OOVV(:,:,:,:) - delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) ! Check convergence @@ -183,7 +151,7 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Update amplitudes - t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + t2(:,:,:,:) = t2(:,:,:,:) + r2(:,:,:,:)/delta_OOVV(:,:,:,:) ! Compute correlation energy @@ -227,10 +195,10 @@ subroutine crCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)' crossed-ring CCD energy ' + write(*,*)' crossed-ring CCD energy ' write(*,*)'----------------------------------------------------' - write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD - write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD + write(*,'(1X,A30,1X,F15.10)')' E(crCCD) = ',ECCD + write(*,'(1X,A30,1X,F15.10)')' Ec(crCCD) = ',EcCCD write(*,*)'----------------------------------------------------' write(*,*) diff --git a/src/CC/drCCD.f90 b/src/CC/drCCD.f90 index 0d21165..063c4b0 100644 --- a/src/CC/drCCD.f90 +++ b/src/CC/drCCD.f90 @@ -39,7 +39,6 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:) - double precision,allocatable :: VOOV(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) @@ -85,11 +84,10 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Create integral batches - allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) + allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO)) OOVV(:,:,:,:) = sERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) OVVO(:,:,:,:) = sERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) - VOOV(:,:,:,:) = sERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas) deallocate(sERI) @@ -135,7 +133,7 @@ subroutine drCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF ! Compute residual - call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) + call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) diff --git a/src/CC/form_crossed_ring_r.f90 b/src/CC/form_crossed_ring_r.f90 new file mode 100644 index 0000000..0009e38 --- /dev/null +++ b/src/CC/form_crossed_ring_r.f90 @@ -0,0 +1,70 @@ +subroutine form_crossed_ring_r(nC,nO,nV,nR,OVOV,OOVV,t2,r2) + +! Form residuals for crossed-ring CCD + + implicit none + +! Input variables + + integer,intent(in) :: nC,nO,nV,nR + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: OVOV(nO,nV,nO,nV) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + double precision,allocatable :: y(:,:,:,:) + +! Output variables + + double precision,intent(out) :: r2(nO,nO,nV,nV) + + r2(:,:,:,:) = 0d0 + + allocate(y(nV,nO,nO,nV)) + + y(:,:,:,:) = 0d0 + + do i=nC+1,nO + do b=1,nV-nR + do l=nC+1,nO + do d=1,nV-nR + do k=nC+1,nO + do c=1,nV-nR + y(b,i,l,d) = y(b,i,l,d) + t2(i,k,c,b)*OOVV(k,l,c,d) + end do + end do + end do + end do + end do + end do + + do i=nC+1,nO + do j=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + + do k=nC+1,nO + do c=1,nV-nR + + r2(i,j,a,b) = r2(i,j,a,b) - OVOV(k,b,i,c)*t2(k,j,a,c) - OVOV(k,a,j,c)*t2(i,k,c,b) + + end do + end do + + do l=nC+1,nO + do d=1,nV-nR + + r2(i,j,a,b) = r2(i,j,a,b) - y(b,i,l,d)*t2(l,j,a,d) + + end do + end do + + end do + end do + end do + end do + +end subroutine form_crossed_ring_r diff --git a/src/CC/form_ring_r.f90 b/src/CC/form_ring_r.f90 index 2757650..8b75138 100644 --- a/src/CC/form_ring_r.f90 +++ b/src/CC/form_ring_r.f90 @@ -1,4 +1,4 @@ -subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) +subroutine form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) ! Form residuals for ring CCD @@ -9,7 +9,6 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) integer,intent(in) :: nC,nO,nV,nR double precision,intent(in) :: t2(nO,nO,nV,nV) double precision,intent(in) :: OVVO(nO,nV,nV,nO) - double precision,intent(in) :: VOOV(nV,nO,nO,nV) double precision,intent(in) :: OOVV(nO,nO,nV,nV) ! Local variables @@ -50,7 +49,7 @@ subroutine form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) do k=nC+1,nO do c=1,nV-nR - r2(i,j,a,b) = r2(i,j,a,b) + VOOV(a,k,i,c)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c) + r2(i,j,a,b) = r2(i,j,a,b) + OVVO(k,a,c,i)*t2(k,j,c,b) + OVVO(k,b,c,j)*t2(i,k,a,c) end do end do diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index 346f715..14dcb0c 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -40,7 +40,6 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, double precision,allocatable :: OOVV(:,:,:,:) double precision,allocatable :: OVVO(:,:,:,:) - double precision,allocatable :: VOOV(:,:,:,:) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) @@ -93,11 +92,10 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, ! Create integral batches - allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VOOV(nV,nO,nO,nV)) + allocate(OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO)) OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO ) - VOOV(:,:,:,:) = dbERI(nO+1:nBas, 1:nO , 1:nO ,nO+1:nBas) deallocate(dbERI) @@ -143,7 +141,7 @@ subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF, ! Compute residual - call form_ring_r(nC,nO,nV,nR,OVVO,VOOV,OOVV,t2,r2) + call form_ring_r(nC,nO,nV,nR,OVVO,OOVV,t2,r2) r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) diff --git a/src/LR/linear_response_B_pp.f90 b/src/LR/linear_response_B_pp.f90 index 3760373..46c151e 100644 --- a/src/LR/linear_response_B_pp.f90 +++ b/src/LR/linear_response_B_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp) +subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B_pp) ! Compute the B matrix of the pp channel @@ -10,7 +10,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B_pp integer,intent(in) :: ispin integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables diff --git a/src/LR/linear_response_C_pp_od.f90 b/src/LR/linear_response_C_pp_od.f90 new file mode 100644 index 0000000..82b22eb --- /dev/null +++ b/src/LR/linear_response_C_pp_od.f90 @@ -0,0 +1,91 @@ +subroutine linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,C_pp) + +! Compute the C matrix of the pp channel (without diagonal term) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,external :: Kronecker_delta + + integer :: a,b,c,d,ab,cd + +! Output variables + + double precision,intent(out) :: C_pp(nVV,nVV) + +! Build C matrix for the singlet manifold + + if(ispin == 1) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = lambda*(ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + end do + end do + end do + end do + + end if + +! Build C matrix for the triplet manifold, or alpha-alpha block, or in the spin-orbital basis + + if(ispin == 2 .or. ispin == 4) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = lambda*(ERI(a,b,c,d) - ERI(a,b,d,c)) + + end do + end do + end do + end do + + end if + +! Build the alpha-beta block of the C matrix + + if(ispin == 3) then + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = lambda*ERI(a,b,c,d) + + end do + end do + end do + end do + + end if + +end subroutine linear_response_C_pp_od diff --git a/src/LR/linear_response_D_pp_od.f90 b/src/LR/linear_response_D_pp_od.f90 new file mode 100644 index 0000000..d9c902d --- /dev/null +++ b/src/LR/linear_response_D_pp_od.f90 @@ -0,0 +1,91 @@ +subroutine linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,D_pp) + +! Compute the D matrix of the pp channel (without the diagonal term) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,external :: Kronecker_delta + + integer :: i,j,k,l,ij,kl + +! Output variables + + double precision,intent(out) :: D_pp(nOO,nOO) + +! Build the D matrix for the singlet manifold + + if(ispin == 1) then + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + + D_pp(ij,kl) = lambda*(ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + end do + end do + end do + end do + + end if + +! Build the D matrix for the triplet manifold, the alpha-alpha block, or in the spin-orbital basis + + if(ispin == 2 .or. ispin == 4) then + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + D_pp(ij,kl) = lambda*(ERI(i,j,k,l) - ERI(i,j,l,k)) + + end do + end do + end do + end do + + end if + +! Build the alpha-beta block of the D matrix + + if(ispin == 3) then + + ij = 0 + do i=nC+1,nO + do j=nC+1,nO + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=nC+1,nO + kl = kl + 1 + + D_pp(ij,kl) = lambda*ERI(i,j,k,l) + + end do + end do + end do + end do + + end if + +end subroutine linear_response_D_pp_od diff --git a/src/LR/linear_response_pp.f90 b/src/LR/linear_response_pp.f90 index 2d9873b..a1decea 100644 --- a/src/LR/linear_response_pp.f90 +++ b/src/LR/linear_response_pp.f90 @@ -75,7 +75,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om else - call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,B) + call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B) ! Diagonal blocks diff --git a/src/MBPT/G0T0.f90 b/src/MBPT/G0T0.f90 index ae946c0..a838c40 100644 --- a/src/MBPT/G0T0.f90 +++ b/src/MBPT/G0T0.f90 @@ -218,6 +218,11 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Free memory + + deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, & + Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t) + ! Compute the BSE correlation energy via the adiabatic connection if(doACFDT) then @@ -234,13 +239,13 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eG0T0,EcAC) + call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI_MO,eHF,eG0T0,EcAC) if(exchange_kernel) then EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) end if diff --git a/src/MBPT/dynamic_Tmatrix_A.f90 b/src/MBPT/dynamic_Tmatrix_A.f90 index 31e75e2..55d4c44 100644 --- a/src/MBPT/dynamic_Tmatrix_A.f90 +++ b/src/MBPT/dynamic_Tmatrix_A.f90 @@ -76,7 +76,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O 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,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) end do A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi @@ -89,7 +89,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O end do do kl=1,nOO - eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + 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 end do diff --git a/src/MBPT/evGT.f90 b/src/MBPT/evGT.f90 index 2309e94..12bf2a5 100644 --- a/src/MBPT/evGT.f90 +++ b/src/MBPT/evGT.f90 @@ -281,13 +281,13 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC) + call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI_MO,eGT,eGT,EcAC) if(exchange_kernel) then EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) end if diff --git a/src/MBPT/qsGT.f90 b/src/MBPT/qsGT.f90 index 1eeff46..1d8e32e 100644 --- a/src/MBPT/qsGT.f90 +++ b/src/MBPT/qsGT.f90 @@ -381,7 +381,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGT,eGT,EcAC) + call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + ERI_MO,eGT,eGT,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/MBPT/static_Tmatrix_TA.f90 index ac4ed6e..c3f2583 100644 --- a/src/MBPT/static_Tmatrix_TA.f90 +++ b/src/MBPT/static_Tmatrix_TA.f90 @@ -50,7 +50,7 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r do kl=1,nOO ! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) - chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) + chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) enddo TA(ia,jb) = TA(ia,jb) + 1d0*lambda*chi diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/MBPT/static_Tmatrix_TB.f90 index 3613980..d502c89 100644 --- a/src/MBPT/static_Tmatrix_TB.f90 +++ b/src/MBPT/static_Tmatrix_TB.f90 @@ -49,8 +49,8 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r enddo do kl=1,nOO -! 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)*Omega2(kl)/Omega2(kl)**2 + eta**2 +! 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)*Omega2(kl)/Omega2(kl)**2 + eta**2 enddo TB(ia,jb) = TB(ia,jb) + 1d0*lambda*chi diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 984550f..edaa4ab 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -9,10 +9,10 @@ program QuAcK logical :: dostab logical :: doKS logical :: doMP2,doMP3,doMP2F12 - logical :: doCCD,doDCD,doCCSD,doCCSDT - logical :: do_drCCD,do_rCCD,do_lCCD,do_crCCD,do_pCCD + logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical :: do_drCCD,do_rCCD,do_crCCD,do_lCCD logical :: doCIS,doCIS_D,doCID,doCISD,doFCI - logical :: doRPA,doRPAx,doppRPA + logical :: doRPA,doRPAx,docrRPA,doppRPA logical :: doADC logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW @@ -91,8 +91,6 @@ program QuAcK double precision :: start_CISD ,end_CISD ,t_CISD double precision :: start_FCI ,end_FCI ,t_FCI double precision :: start_RPA ,end_RPA ,t_RPA - double precision :: start_RPAx ,end_RPAx ,t_RPAx - double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA double precision :: start_ADC ,end_ADC ,t_ADC double precision :: start_GF2 ,end_GF2 ,t_GF2 double precision :: start_GF3 ,end_GF3 ,t_GF3 @@ -160,17 +158,17 @@ program QuAcK ! Which calculations do you want to do? - call read_methods(doRHF,doUHF,doKS,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - doRPA,doRPAx,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & - doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT, & + call read_methods(doRHF,doUHF,doKS,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + doRPA,doRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW, & + doufG0W0,doufGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read options for methods @@ -444,7 +442,7 @@ 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) end if end if @@ -645,6 +643,24 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform crossed-ring CCD calculation +!------------------------------------------------------------------------ + + if(do_crCCD) then + + call cpu_time(start_CCD) + call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & + ERI_MO,ENuc,ERHF,eHF) + + call cpu_time(end_CCD) + + t_CCD = end_CCD - start_CCD + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Perform ladder CCD calculation !------------------------------------------------------------------------ @@ -662,32 +678,11 @@ program QuAcK end if -!------------------------------------------------------------------------ -! Perform crossed-ring CCD calculation -!------------------------------------------------------------------------ - - do_crCCD = .false. - - if(do_crCCD) then - - call cpu_time(start_CCD) - call ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) - call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & - ERI_MO,ENuc,ERHF,eHF) - - call cpu_time(end_CCD) - - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds' - write(*,*) - - end if - !------------------------------------------------------------------------ ! Perform pair CCD calculation !------------------------------------------------------------------------ - if(do_pCCD) then + if(dopCCD) then call cpu_time(start_CCD) call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) @@ -787,7 +782,7 @@ program QuAcK if(doRPAx) then - call cpu_time(start_RPAx) + call cpu_time(start_RPA) if(unrestricted) then call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & @@ -798,10 +793,26 @@ program QuAcK call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) end if - call cpu_time(end_RPAx) + call cpu_time(end_RPA) - t_RPAx = end_RPAx - start_RPAx - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPAx,' seconds' + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute cr-RPA excitations +!------------------------------------------------------------------------ + + if(docrRPA) then + + call cpu_time(start_RPA) + call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + call cpu_time(end_RPA) + + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if @@ -812,12 +823,12 @@ program QuAcK if(doppRPA) then - call cpu_time(start_ppRPA) - call ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) - call cpu_time(end_ppRPA) + call cpu_time(start_RPA) + call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + call cpu_time(end_RPA) - t_ppRPA = end_ppRPA - start_ppRPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_ppRPA,' seconds' + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds' write(*,*) end if diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 3b77e7c..8ddd885 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,14 +1,14 @@ -subroutine read_methods(doRHF,doUHF,doKS,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - doRPA,doRPAx,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW, & - doufG0W0,doufGW, & - doG0T0,doevGT,doqsGT, & +subroutine read_methods(doRHF,doUHF,doKS,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + doRPA,doRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW, & + doufG0W0,doufGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read desired methods @@ -19,10 +19,10 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & logical,intent(out) :: doRHF,doUHF,doKS,doMOM logical,intent(out) :: doMP2,doMP3,doMP2F12 - logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT - logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD + logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI - logical,intent(out) :: doRPA,doRPAx,doppRPA + logical,intent(out) :: doRPA,doRPAx,docrRPA,doppRPA logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW logical,intent(out) :: doG0T0,doevGT,doqsGT @@ -48,14 +48,15 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doMP2F12 = .false. doCCD = .false. + dopCCD = .false. doDCD = .false. doCCSD = .false. doCCSDT = .false. do_drCCD = .false. do_rCCD = .false. + do_crCCD = .false. do_lCCD = .false. - do_pCCD = .false. doCIS = .false. doCIS_D = .false. @@ -65,6 +66,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doRPA = .false. doRPAx = .false. + docrRPA = .false. doppRPA = .false. doG0F2 = .false. @@ -105,19 +107,20 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read CC methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') doCCD = .true. - if(answer2 == 'T') doDCD = .true. - if(answer3 == 'T') doCCSD = .true. - if(answer4 == 'T') doCCSDT = .true. + read(1,*) answer1,answer2,answer3,answer4,answer5 + if(answer1 == 'T') doCCD = .true. + if(answer2 == 'T') dopCCD = .true. + if(answer3 == 'T') doDCD = .true. + if(answer4 == 'T') doCCSD = .true. + if(answer5 == 'T') doCCSDT = .true. ! Read weird CC methods read(1,*) read(1,*) answer1,answer2,answer3,answer4 if(answer1 == 'T') do_drCCD = .true. if(answer2 == 'T') do_rCCD = .true. - if(answer3 == 'T') do_lCCD = .true. - if(answer4 == 'T') do_pCCD = .true. + if(answer3 == 'T') do_crCCD = .true. + if(answer4 == 'T') do_lCCD = .true. ! Read excited state methods @@ -131,10 +134,11 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & if(doCIS_D) doCIS = .true. read(1,*) - read(1,*) answer1,answer2,answer3 + read(1,*) answer1,answer2,answer3,answer4 if(answer1 == 'T') doRPA = .true. if(answer2 == 'T') doRPAx = .true. - if(answer3 == 'T') doppRPA = .true. + if(answer3 == 'T') docrRPA = .true. + if(answer4 == 'T') doppRPA = .true. ! Read Green function methods diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index c41d3d7..3245e34 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -1,5 +1,4 @@ -subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t, & +subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & ERI,eT,eGT,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix @@ -27,32 +26,10 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple integer,intent(in) :: nR integer,intent(in) :: nS - integer,intent(in) :: nOOs - integer,intent(in) :: nOOt - integer,intent(in) :: nVVs - integer,intent(in) :: nVVt - double precision,intent(in) :: eT(nBas) double precision,intent(in) :: eGT(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega1s(nVVs) - double precision,intent(in) :: X1s(nVVs,nVVs) - double precision,intent(in) :: Y1s(nOOs,nVVs) - double precision,intent(in) :: Omega2s(nOOs) - double precision,intent(in) :: X2s(nVVs,nOOs) - double precision,intent(in) :: Y2s(nOOs,nOOs) - double precision,intent(in) :: rho1s(nBas,nBas,nVVs) - double precision,intent(in) :: rho2s(nBas,nBas,nOOs) - double precision,intent(in) :: Omega1t(nVVt) - double precision,intent(in) :: X1t(nVVt,nVVt) - double precision,intent(in) :: Y1t(nOOt,nVVt) - double precision,intent(in) :: Omega2t(nOOt) - double precision,intent(in) :: X2t(nVVt,nOOt) - double precision,intent(in) :: Y2t(nOOt,nOOt) - double precision,intent(in) :: rho1t(nBas,nBas,nVVt) - double precision,intent(in) :: rho2t(nBas,nBas,nOOt) - ! Local variables integer :: ispin @@ -62,6 +39,9 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision :: lambda double precision,allocatable :: Ec(:,:) + integer :: nOOs,nOOt + integer :: nVVs,nVVt + double precision :: EcRPA(nspin) double precision,allocatable :: TA(:,:) double precision,allocatable :: TB(:,:) @@ -69,14 +49,40 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) + ! Output variables double precision,intent(out) :: EcAC(nspin) +! Useful quantities + +! nOOs = nO*nO +! nVVs = nV*nV + + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 + + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 + ! Memory allocation - allocate(Ec(nAC,nspin)) + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(Ec(nAC,nspin)) ! Antisymmetrized kernel version @@ -115,31 +121,31 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple TA(:,:) = 0d0 TB(:,:) = 0d0 - if(doXBS) then +! if(doXBS) then - isp_T = 1 - iblock = 3 +! isp_T = 1 +! iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) +! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & +! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) - call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) - isp_T = 2 - iblock = 4 +! isp_T = 2 +! iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) +! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & +! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) - call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) - end if +! end if call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -180,31 +186,36 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) - if(doXBS) then + ! Initialize T matrix - isp_T = 1 - iblock = 3 + TA(:,:) = 0d0 + TB(:,:) = 0d0 + +! if(doXBS) then + +! isp_T = 1 +! iblock = 3 ! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI, & ! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(isp_T)) - call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) - isp_T = 2 - iblock = 4 +! isp_T = 2 +! iblock = 4 ! call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI, & ! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(isp_T)) - call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) +! call static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) +! if(.not.TDA) call static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) - end if +! end if call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/RPA/ACFDT_cr.f90 b/src/RPA/ACFDT_cr.f90 new file mode 100644 index 0000000..d82c77e --- /dev/null +++ b/src/RPA/ACFDT_cr.f90 @@ -0,0 +1,173 @@ +subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC) + +! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem +! for the crossed-ring contribution + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doXBS + logical,intent(in) :: exchange_kernel + logical,intent(in) :: dRPA + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: BSE + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: isp_W + integer :: iAC + double precision :: lambda + double precision,allocatable :: Ec(:,:) + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) + +! Output variables + + double precision,intent(out) :: EcAC(nspin) + +! Memory allocation + + allocate(Ec(nAC,nspin)) + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + +! Antisymmetrized kernel version + + if(exchange_kernel) then + + write(*,*) + write(*,*) '*** Exchange kernel version ***' + write(*,*) + + end if + + EcAC(:) = 0d0 + Ec(:,:) = 0d0 + +! Compute (singlet) RPA screening + + isp_W = 1 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + +! Singlet manifold + + if(singlet) then + + ispin = 1 + + write(*,*) '--------------' + write(*,*) 'Singlet states' + write(*,*) '--------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = -rAC(iAC) + + if(doXBS) then + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) +! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) + + end if + + call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & + rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & + ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin)) + + if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +! Triplet manifold + + if(triplet) then + + ispin = 2 + + write(*,*) '--------------' + write(*,*) 'Triplet states' + write(*,*) '--------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = -rAC(iAC) + + if(doXBS) then + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + + end if + + call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & + rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = -0.5d0*dot_product(wAC,Ec(:,ispin)) + + if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +end subroutine ACFDT_cr diff --git a/src/RPA/ACFDT_pp.f90 b/src/RPA/ACFDT_pp.f90 new file mode 100644 index 0000000..540c77c --- /dev/null +++ b/src/RPA/ACFDT_pp.f90 @@ -0,0 +1,145 @@ +subroutine ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + +! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for pp sector + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: singlet + logical,intent(in) :: triplet + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: iAC + double precision :: lambda + double precision,allocatable :: Ec(:,:) + + integer :: nOOs,nOOt + integer :: nVVs,nVVt + + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) + double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) + +! Output variables + + double precision,intent(out) :: EcAC(nspin) + +! Useful quantities + + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 + + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 + +! Memory allocation + + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & + rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & + rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + allocate(Ec(nAC,nspin)) + +! Antisymmetrized kernel version + + EcAC(:) = 0d0 + Ec(:,:) = 0d0 + +! Singlet manifold + + if(singlet) then + + ispin = 1 + + write(*,*) '--------------' + write(*,*) 'Singlet states' + write(*,*) '--------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = rAC(iAC) + + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,e,ERI,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcAC(ispin)) + + call ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOOs,nVVs,X1s,Y1s,X2s,Y2s,Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +! Triplet manifold + + if(triplet) then + + ispin = 2 + + write(*,*) '--------------' + write(*,*) 'Triplet states' + write(*,*) '--------------' + write(*,*) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' + write(*,*) '-----------------------------------------------------------------------------------' + + do iAC=1,nAC + + lambda = rAC(iAC) + + ! Initialize T matrix + + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,e,ERI,Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcAC(ispin)) + + call ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOOt,nVVt,X1t,Y1t,X2t,Y2t,Ec(iAC,ispin)) + + write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) + + end do + + EcAC(ispin) = 1.5d0*dot_product(wAC,Ec(:,ispin)) + + write(*,*) '-----------------------------------------------------------------------------------' + write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) + + end if + +end subroutine ACFDT_pp diff --git a/src/RPA/ACFDT_pp_correlation_energy.f90 b/src/RPA/ACFDT_pp_correlation_energy.f90 new file mode 100644 index 0000000..28991f8 --- /dev/null +++ b/src/RPA/ACFDT_pp_correlation_energy.f90 @@ -0,0 +1,51 @@ +subroutine ACFDT_pp_correlation_energy(ispin,nBas,nC,nO,nV,nR,nS,ERI,nOO,nVV,X1,Y1,X2,Y2,EcAC) + +! Compute the correlation energy via the adiabatic connection formula for the pp sector + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + double precision,allocatable :: B(:,:) + double precision,allocatable :: C(:,:) + double precision,allocatable :: D(:,:) + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: EcAC + +! Memory allocation + + allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO)) + +! Build pp matrices + + call linear_response_B_pp (ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,B) + call linear_response_C_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,C) + call linear_response_D_pp_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,D) + +! Compute Tr(K x P_lambda) + + EcAC = trace_matrix(nVV,matmul(transpose(X1),matmul(B,Y1)) + matmul(transpose(Y1),matmul(transpose(B),X1))) & + + trace_matrix(nVV,matmul(transpose(X1),matmul(C,X1)) + matmul(transpose(Y1),matmul(D,Y1))) & + + trace_matrix(nOO,matmul(transpose(X2),matmul(B,Y2)) + matmul(transpose(Y2),matmul(transpose(B),X2))) & + + trace_matrix(nOO,matmul(transpose(X2),matmul(C,X2)) + matmul(transpose(Y2),matmul(D,Y2))) & + - trace_matrix(nVV,C) - trace_matrix(nOO,D) + +end subroutine ACFDT_pp_correlation_energy + diff --git a/src/RPA/RPA.f90 b/src/RPA/RPA.f90 index 3d6cb63..0dce1a7 100644 --- a/src/RPA/RPA.f90 +++ b/src/RPA/RPA.f90 @@ -87,12 +87,12 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, endif - if(exchange_kernel) then +! if(exchange_kernel) then - EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(2) +! EcRPA(1) = 0.5d0*EcRPA(1) +! EcRPA(2) = 1.5d0*EcRPA(2) - end if +! end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/RPAx.f90 b/src/RPA/RPAx.f90 index 6b096e4..f47351d 100644 --- a/src/RPA/RPAx.f90 +++ b/src/RPA/RPAx.f90 @@ -13,8 +13,8 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet - double precision,intent(in) :: eta logical,intent(in) :: triplet + double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -89,12 +89,12 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR endif - if(exchange_kernel) then +! if(exchange_kernel) then EcRPAx(1) = 0.5d0*EcRPAx(1) EcRPAx(2) = 1.5d0*EcRPAx(2) - end if +! end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/ehRPA.f90 b/src/RPA/crRPA.f90 similarity index 66% rename from src/RPA/ehRPA.f90 rename to src/RPA/crRPA.f90 index 103f5c3..7671595 100644 --- a/src/RPA/ehRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -1,7 +1,7 @@ -subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & +subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI,dipole_int,eHF) -! Perform random phase approximation calculation with exchange (aka TDHF) +! Crossed-ring channel of the random phase approximation implicit none include 'parameters.h' @@ -42,7 +42,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n write(*,*) write(*,*)'***********************************************************' - write(*,*)'| Random phase approximation calculation: eh channel |' + write(*,*)'| Random phase approximation calculation: cr channel |' write(*,*)'***********************************************************' write(*,*) @@ -70,7 +70,7 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin)) + call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -83,48 +83,48 @@ subroutine ehRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('ehRPA@HF ',ispin,nS,Omega(:,ispin)) + call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif - if(exchange_kernel) then +! if(exchange_kernel) then EcRPAx(1) = 0.5d0*EcRPAx(1) EcRPAx(2) = 1.5d0*EcRPAx(2) - end if +! end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (singlet) =',EcRPAx(1) - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy (triplet) =',EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA correlation energy =',EcRPAx(1) + EcRPAx(2) - write(*,'(2X,A50,F20.10)') 'Tr@ehRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (singlet) =',EcRPAx(1) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy (triplet) =',EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA correlation energy =',EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@crRPA total energy =',ENuc + ERHF + EcRPAx(1) + EcRPAx(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) ! Compute the correlation energy via the adiabatic connection -! if(doACFDT) then + if(doACFDT) then -! write(*,*) '-------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of ehRPA correlation energy' -! write(*,*) '-------------------------------------------------------' -! write(*,*) + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of crRPA correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) -! call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) + call ACFDT_cr(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1) -! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2) -! write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@crRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if -end subroutine ehRPA +end subroutine crRPA diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index a40e297..3e0b18a 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) +subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Perform pp-RPA calculation @@ -8,6 +8,7 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Input variables logical,intent(in) :: TDA + logical,intent(in) :: doACFDT logical,intent(in) :: singlet logical,intent(in) :: triplet integer,intent(in) :: nBas @@ -23,16 +24,18 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Local variables integer :: ispin - integer :: nOO - integer :: nVV - double precision,allocatable :: Omega1(:,:) - double precision,allocatable :: X1(:,:,:) - double precision,allocatable :: Y1(:,:,:) - double precision,allocatable :: Omega2(:,:) - double precision,allocatable :: X2(:,:,:) - double precision,allocatable :: Y2(:,:,:) + integer :: nS + integer :: nOOs,nOOt + integer :: nVVs,nVVt + double precision,allocatable :: Omega1s(:),Omega1t(:) + double precision,allocatable :: X1s(:,:),X1t(:,:) + double precision,allocatable :: Y1s(:,:),Y1t(:,:) + double precision,allocatable :: Omega2s(:),Omega2t(:) + double precision,allocatable :: X2s(:,:),X2t(:,:) + double precision,allocatable :: Y2s(:,:),Y2t(:,:) double precision :: Ec_ppRPA(nspin) + double precision :: EcAC(nspin) ! Hello world @@ -45,6 +48,25 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Initialization Ec_ppRPA(:) = 0d0 + EcAC(:) = 0d0 + +! Useful quantities + + nS = nO*nV + + nOOs = nO*(nO+1)/2 + nVVs = nV*(nV+1)/2 + + nOOt = nO*(nO-1)/2 + nVVt = nV*(nV-1)/2 + + ! Memory allocation + + allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & + Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) + + allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & + Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) ! Singlet manifold @@ -52,25 +74,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ispin = 1 - ! Useful quantities + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, & + Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin)) - nOO = nO*(nO+1)/2 - nVV = nV*(nV+1)/2 - - ! Memory allocation - - allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & - Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - - call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & - Ec_ppRPA(ispin)) - - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) - - deallocate(Omega1,X1,Y1,Omega2,X2,Y2) + call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s) + call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s) endif @@ -80,26 +88,11 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ispin = 2 - ! Useful quantities + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, & + Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin)) - nOO = nO*(nO-1)/2 - nVV = nV*(nV-1)/2 - - ! Memory allocation - - allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & - Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - - - call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & - Ec_ppRPA(ispin)) - - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) - - deallocate(Omega1,X1,Y1,Omega2,X2,Y2) + call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t) + call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t) endif @@ -112,4 +105,27 @@ subroutine ppRPA(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Compute the correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '---------------------------------------------------------' + write(*,*) 'Adiabatic connection version of pp-RPA correlation energy' + write(*,*) '---------------------------------------------------------' + write(*,*) + + call ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end subroutine ppRPA