From a13ab94a52cb075910afa83f841d4c5f7308e404 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 22 Mar 2020 22:22:47 +0100 Subject: [PATCH] FC in pCCD --- examples/molecule.F2 | 2 +- src/QuAcK/G0T0.f90 | 2 -- src/QuAcK/QuAcK.f90 | 8 +++--- src/QuAcK/pCCD.f90 | 64 +++++++++++++++++++++++++++----------------- 4 files changed, 46 insertions(+), 30 deletions(-) diff --git a/examples/molecule.F2 b/examples/molecule.F2 index 9c99be9..12bba91 100644 --- a/examples/molecule.F2 +++ b/examples/molecule.F2 @@ -2,4 +2,4 @@ 2 9 9 0 0 # Znuc x y z F 0. 0. 0. - F 0. 0. 2 + F 0. 0. 3.3 diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index bbeed71..8c43f4b 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -194,8 +194,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Compute the BSE correlation energy via the adiabatic connection - ! Compute the BSE correlation energy via the adiabatic connection if(doACFDT) then diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ad8acea..fd7eacf 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -447,7 +447,7 @@ program QuAcK if(do_pCCD) then call cpu_time(start_CCD) - call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nO(1),nV(1),ERI_MO_basis,ENuc,ERHF,eHF) + call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -513,7 +513,8 @@ program QuAcK if(doppRPA) then call cpu_time(start_ppRPA) - call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF) + call ppRPA(singlet_manifold,triplet_manifold, & + nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_ppRPA) t_ppRPA = end_ppRPA - start_ppRPA @@ -529,7 +530,8 @@ program QuAcK if(doADC) then call cpu_time(start_ADC) - call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF,nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO_basis) + call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & + nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO_basis) call cpu_time(end_ADC) t_ADC = end_ADC - start_ADC diff --git a/src/QuAcK/pCCD.f90 b/src/QuAcK/pCCD.f90 index ecd7f11..7c8d5bf 100644 --- a/src/QuAcK/pCCD.f90 +++ b/src/QuAcK/pCCD.f90 @@ -1,4 +1,4 @@ -subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) +subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) ! pair CCD module @@ -10,7 +10,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas,nO,nV + integer,intent(in) :: nBas,nC,nO,nV,nR double precision,intent(in) :: ENuc,ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -31,8 +31,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) double precision,allocatable :: OVVO(:,:) double precision,allocatable :: VVVV(:,:) - double precision,allocatable :: X(:,:) - double precision,allocatable :: Y(:,:) + double precision,allocatable :: y(:,:) double precision,allocatable :: r(:,:) double precision,allocatable :: t(:,:) @@ -51,8 +50,8 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) allocate(delta_OOVV(nO,nV)) - do i=1,nO - do a=1,nV + do i=nC+1,nO + do a=1,nV-nR delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i)) enddo enddo @@ -61,22 +60,22 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV)) - do i=1,nO - do j=1,nO + do i=nC+1,nO + do j=nC+1,nO OOOO(i,j) = ERI(i,i,j,j) end do end do - do i=1,nO - do a=1,nV + do i=nC+1,nO + do a=1,nV-nR OOVV(i,a) = ERI(i,i,nO+a,nO+a) OVOV(i,a) = ERI(i,nO+a,i,nO+a) OVVO(i,a) = ERI(i,nO+a,nO+a,i) end do end do - do a=1,nV - do b=1,nV + do a=1,nV-nR + do b=1,nV-nR VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b) end do end do @@ -89,7 +88,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) ! Initialization - allocate(r(nO,nV),X(nV,nV),Y(nO,nO)) + allocate(r(nO,nV),y(nO,nO)) Conv = 1d0 nSCF = 0 @@ -113,22 +112,34 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) ! Form intermediate array - X(:,:) = matmul(transpose(OOVV(:,:)),t(:,:)) - Y(:,:) = matmul(t(:,:),transpose(OOVV(:,:))) + y(:,:) = 0d0 + do i=nC+1,nO + do j=nC+1,nO + do b=1,nV-nR + y(i,j) = y(i,j) + OOVV(j,b)*t(i,b) + end do + end do + end do ! Compute residual - do i=1,nO - do a=1,nV - r(i,a) = - 2d0*(X(a,a) + Y(i,i))*t(i,a) + do i=nC+1,nO + do a=1,nV-nR + + r(i,a) = OOVV(i,a) + delta_OOVV(i,a)*t(i,a) & + - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a) + + do j=nC+1,nO + r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a) + end do + + do b=1,nV-nR + r(i,a) = r(i,a) - 2d0*OOVV(i,b)*t(i,b)*t(i,a) + VVVV(a,b)*t(i,b) + end do + end do end do - r(:,:) = r(:,:) + OOVV(:,:) + delta_OOVV(:,:)*t(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t(:,:))*t(:,:) & - + matmul(t(:,:),transpose(VVVV(:,:))) & - + matmul(transpose(OOOO(:,:)),t(:,:)) & - + matmul(Y(:,:),t) ! Check convergence @@ -140,7 +151,12 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nO,nV,ERI,ENuc,ERHF,eHF) ! Compute correlation energy - EcCCD = trace_matrix(nO,matmul(t(:,:),transpose(OOVV(:,:)))) + EcCCD = 0d0 + do i=nC+1,nO + do a=1,nV-nR + EcCCD = EcCCD + t(i,a)*OOVV(i,a) + end do + end do ! Dump results