diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index 1fcf56d..3e40966 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -117,7 +117,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! Memory allocation allocate(old_eh_Om(nS),old_ee_Om(nVV),old_hh_Om(nOO)) - allocate(eh_rho(nOrb,nOrb,nS+nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO)) + allocate(eh_rho(nOrb,nOrb,nS),ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO)) allocate(old_eh_Phi(nOrb,nOrb,nOrb,nOrb),old_pp_Phi(nOrb,nOrb,nOrb,nOrb)) ! Initialization @@ -259,9 +259,9 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if - ! Bpp(:,:) = Bpp(:,:) + pp_Gam_B(:,:) - ! Cpp(:,:) = Cpp(:,:) + pp_Gam_C(:,:) - ! Dpp(:,:) = Dpp(:,:) + pp_Gam_D(:,:) + Bpp(:,:) = Bpp(:,:) + pp_Gam_B(:,:) + Cpp(:,:) = Cpp(:,:) + pp_Gam_C(:,:) + Dpp(:,:) = Dpp(:,:) + pp_Gam_D(:,:) call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,ee_Om,X1,Y1,hh_Om,X2,Y2,EcRPA) call wall_time(end_t) @@ -328,7 +328,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, deallocate(eh_rho,ee_rho,hh_rho) ! TODO Once we will compute the blocks of kernel starting from the 4-tensors we can move the freeing up ! Memory allocation - allocate(eh_rho(nOrb,nOrb,nS+nS)) + allocate(eh_rho(nOrb,nOrb,nS)) allocate(ee_rho(nOrb,nOrb,nVV),hh_rho(nOrb,nOrb,nOO)) ! Build singlet eh integrals @@ -384,7 +384,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for pp reducible kernel =',t,' seconds' write(*,*) - ! alpha = 0.01d0 + ! alpha = 0.1d0 ! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:) ! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:) diff --git a/src/Parquet/G_eh_Phi.f90 b/src/Parquet/G_eh_Phi.f90 index 5ec86d3..d81c6fc 100644 --- a/src/Parquet/G_eh_Phi.f90 +++ b/src/Parquet/G_eh_Phi.f90 @@ -1,12 +1,12 @@ subroutine G_eh_Phi(nOrb,nC,nR,nS,eh_Om,eh_rho,eh_Phi) -! Compute irreducible vertex in the triplet pp channel +! Compute irreducible vertex in the eh channel implicit none ! Input variables integer,intent(in) :: nOrb,nC,nR,nS double precision,intent(in) :: eh_Om(nS) - double precision,intent(in) :: eh_rho(nOrb,nOrb,nS+nS) + double precision,intent(in) :: eh_rho(nOrb,nOrb,nS) ! Local variables integer :: p,q,r,s @@ -23,11 +23,11 @@ subroutine G_eh_Phi(nOrb,nC,nR,nS,eh_Om,eh_rho,eh_Phi) do q = nC+1, nOrb-nR do p = nC+1, nOrb-nR - ! do n=1,nS - ! eh_Phi(p,q,r,s) = eh_Phi(p,q,r,s) & - ! - eh_rho(r,p,n)*eh_rho(q,s,n)/eh_Om(n) & - ! - eh_rho(p,r,nS+n)*eh_rho(s,q,nS+n)/eh_Om(n) - ! end do + do n=1,nS + eh_Phi(p,q,r,s) = eh_Phi(p,q,r,s) & + - eh_rho(r,p,n)*eh_rho(q,s,n)/eh_Om(n) & + - eh_rho(p,r,n)*eh_rho(s,q,n)/eh_Om(n) + end do enddo enddo diff --git a/src/Parquet/G_screened_integrals.f90 b/src/Parquet/G_screened_integrals.f90 index f6c52bf..9358173 100644 --- a/src/Parquet/G_screened_integrals.f90 +++ b/src/Parquet/G_screened_integrals.f90 @@ -15,7 +15,7 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho double precision :: X,Y ! Output variables - double precision,intent(out) :: rho(nOrb,nOrb,nS+nS) + double precision,intent(out) :: rho(nOrb,nOrb,nS) rho(:,:,:) = 0d0 ! !$OMP PARALLEL & @@ -25,26 +25,28 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho ! !$OMP DO do q=nC+1,nOrb-nR do p=nC+1,nOrb-nR - - jb = 0 - do j=nC+1,nO - do b=nO+1,nOrb-nR - jb = jb + 1 - do ia=1,nS + do ia=1,nS + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) - rho(p,q,ia) = (ERI(q,j,p,b) - ERI(q,j,b,p)) * X & - + (- eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X - - rho(p,q,nS+ia) = (ERI(q,b,p,j) - ERI(q,b,j,p)) * Y & - + (- eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y + rho(p,q,ia) = rho(p,q,ia) & + + (ERI(q,j,p,b) - ERI(q,j,b,p) & + - 0d0*eh_Phi(q,j,b,p) + 0d0*pp_Phi(q,j,p,b)) * X & + + (ERI(q,b,p,j) - ERI(q,b,j,p) & + - 0d0*eh_Phi(q,b,j,p) + 0d0*pp_Phi(q,b,p,j)) * Y end do end do + end do + end do end do ! !$OMP END DO @@ -76,7 +78,6 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 integer :: a,b,c,d integer :: p,q integer :: ab,cd,ij,kl - double precision,external :: Kronecker_delta ! Output variables @@ -106,8 +107,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do c=nO+1,nOrb-nR do d=c+1,nOrb-nR cd = cd + 1 - rho1(p,q,ab) = ( ERI(p,q,c,d) - ERI(p,q,d,c) ) * X1(cd,ab) & - + ( eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X1(cd,ab) + rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & + + 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X1(cd,ab) end do end do @@ -115,8 +116,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho1(p,q,ab) = ( ERI(p,q,k,l) - ERI(p,q,l,k) ) * Y1(kl,ab) & - + ( eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y1(kl,ab) + rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & + + 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab) end do end do @@ -131,8 +132,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do c=nO+1,nOrb-nR do d=c+1,nOrb-nR cd = cd + 1 - rho2(p,q,ij) = ( ERI(p,q,c,d) - ERI(p,q,d,c) ) * X2(cd,ij) & - + ( eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X2(cd,ij) + rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & + + 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X2(cd,ij) end do end do @@ -140,8 +141,8 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho2(p,q,ij) = ( ERI(p,q,k,l) - ERI(p,q,l,k) ) * Y2(kl,ij) & - + ( eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y2(kl,ij) + rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & + + 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij) end do end do end do diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index dbe8e42..cb0e02f 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -117,7 +117,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, allocate(old_eh_sing_Om(nS),old_eh_trip_Om(nS)) allocate(old_ee_sing_Om(nVVs),old_hh_sing_Om(nOOs)) allocate(old_ee_trip_Om(nVVt),old_hh_trip_Om(nOOt)) - allocate(eh_sing_rho(nOrb,nOrb,nS+nS),eh_trip_rho(nOrb,nOrb,nS+nS)) + allocate(eh_sing_rho(nOrb,nOrb,nS),eh_trip_rho(nOrb,nOrb,nS)) allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs)) allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt)) allocate(old_eh_sing_Phi(nOrb,nOrb,nOrb,nOrb),old_eh_trip_Phi(nOrb,nOrb,nOrb,nOrb)) @@ -330,9 +330,9 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if - ! Bpp(:,:) = Bpp(:,:) + pp_sing_Gam_B(:,:) - ! Cpp(:,:) = Cpp(:,:) + pp_sing_Gam_C(:,:) - ! Dpp(:,:) = Dpp(:,:) + pp_sing_Gam_D(:,:) + Bpp(:,:) = Bpp(:,:) + pp_sing_Gam_B(:,:) + Cpp(:,:) = Cpp(:,:) + pp_sing_Gam_C(:,:) + Dpp(:,:) = Dpp(:,:) + pp_sing_Gam_D(:,:) call ppRLR(TDA,nOOs,nVVs,Bpp,Cpp,Dpp,ee_sing_Om,X1s,Y1s,hh_sing_Om,X2s,Y2s,EcRPA) call wall_time(end_t) @@ -388,9 +388,9 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, end if - ! Bpp(:,:) = Bpp(:,:) + pp_trip_Gam_B(:,:) - ! Cpp(:,:) = Cpp(:,:) + pp_trip_Gam_C(:,:) - ! Dpp(:,:) = Dpp(:,:) + pp_trip_Gam_D(:,:) + Bpp(:,:) = Bpp(:,:) + pp_trip_Gam_B(:,:) + Cpp(:,:) = Cpp(:,:) + pp_trip_Gam_C(:,:) + Dpp(:,:) = Dpp(:,:) + pp_trip_Gam_D(:,:) call ppRLR(TDA,nOOt,nVVt,Bpp,Cpp,Dpp,ee_trip_Om,X1t,Y1t,hh_trip_Om,X2t,Y2t,EcRPA) @@ -439,7 +439,7 @@ subroutine RParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, deallocate(eh_sing_rho,eh_trip_rho,ee_sing_rho,ee_trip_rho,hh_sing_rho,hh_trip_rho) ! TODO Once we will compute the blocks of kernel starting from the 4-tensors we can move the freeing up ! Memory allocation - allocate(eh_sing_rho(nOrb,nOrb,nS+nS),eh_trip_rho(nOrb,nOrb,nS+nS)) + allocate(eh_sing_rho(nOrb,nOrb,nS),eh_trip_rho(nOrb,nOrb,nS)) allocate(ee_sing_rho(nOrb,nOrb,nVVs),hh_sing_rho(nOrb,nOrb,nOOs)) allocate(ee_trip_rho(nOrb,nOrb,nVVt),hh_trip_rho(nOrb,nOrb,nOOt)) diff --git a/src/Parquet/R_eh_singlet_Phi.f90 b/src/Parquet/R_eh_singlet_Phi.f90 index 5971e26..91bb865 100644 --- a/src/Parquet/R_eh_singlet_Phi.f90 +++ b/src/Parquet/R_eh_singlet_Phi.f90 @@ -7,7 +7,7 @@ subroutine R_eh_singlet_Phi(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_sing_Phi) ! Input variables integer,intent(in) :: nOrb,nC,nR,nS double precision,intent(in) :: eh_sing_Om(nS) - double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS+nS) + double precision,intent(in) :: eh_sing_rho(nOrb,nOrb,nS) ! Local variables integer :: p,q,r,s @@ -24,11 +24,11 @@ subroutine R_eh_singlet_Phi(nOrb,nC,nR,nS,eh_sing_Om,eh_sing_rho,eh_sing_Phi) do q = nC+1, nOrb-nR do p = nC+1, nOrb-nR - ! do n=1,nS - ! eh_sing_Phi(p,q,r,s) = eh_sing_Phi(p,q,r,s) & - ! - eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) & - ! - eh_sing_rho(p,r,n+nS)*eh_sing_rho(s,q,n+nS)/eh_sing_Om(n) - ! end do + do n=1,nS + eh_sing_Phi(p,q,r,s) = eh_sing_Phi(p,q,r,s) & + - eh_sing_rho(r,p,n)*eh_sing_rho(q,s,n)/eh_sing_Om(n) & + - eh_sing_rho(p,r,n)*eh_sing_rho(s,q,n)/eh_sing_Om(n) + end do enddo enddo diff --git a/src/Parquet/R_eh_triplet_Phi.f90 b/src/Parquet/R_eh_triplet_Phi.f90 index a6ec43f..e6d8587 100644 --- a/src/Parquet/R_eh_triplet_Phi.f90 +++ b/src/Parquet/R_eh_triplet_Phi.f90 @@ -7,7 +7,7 @@ subroutine R_eh_triplet_Phi(nOrb,nC,nR,nS,eh_trip_Om,eh_trip_rho,eh_trip_Phi) ! Input variables integer,intent(in) :: nOrb,nC,nR,nS double precision,intent(in) :: eh_trip_Om(nS) - double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS+nS) + double precision,intent(in) :: eh_trip_rho(nOrb,nOrb,nS) ! Local variables integer :: p,q,r,s @@ -24,11 +24,11 @@ subroutine R_eh_triplet_Phi(nOrb,nC,nR,nS,eh_trip_Om,eh_trip_rho,eh_trip_Phi) do q = nC+1, nOrb-nR do p = nC+1, nOrb-nR - ! do n=1,nS - ! eh_trip_Phi(p,q,r,s) = eh_trip_Phi(p,q,r,s) & - ! - eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) & - ! - eh_trip_rho(p,r,n+nS)*eh_trip_rho(s,q,n+nS)/eh_trip_Om(n) - ! end do + do n=1,nS + eh_trip_Phi(p,q,r,s) = eh_trip_Phi(p,q,r,s) & + - eh_trip_rho(r,p,n)*eh_trip_rho(q,s,n)/eh_trip_Om(n) & + - eh_trip_rho(p,r,n)*eh_trip_rho(s,q,n)/eh_trip_Om(n) + end do enddo enddo diff --git a/src/Parquet/R_pp_singlet_Gam.f90 b/src/Parquet/R_pp_singlet_Gam.f90 index cd9987c..4bf6f88 100644 --- a/src/Parquet/R_pp_singlet_Gam.f90 +++ b/src/Parquet/R_pp_singlet_Gam.f90 @@ -35,7 +35,7 @@ subroutine R_pp_singlet_Gamma_D(nOrb,nC,nO,nOOs,eh_sing_Phi,eh_trip_Phi,pp_sing_ kl = kl +1 pp_sing_Gam_D(ij,kl) = 0.5d0*eh_sing_Phi(i,j,k,l) - 1.5d0*eh_trip_Phi(i,j,k,l) & - - 1.5d0*eh_sing_Phi(i,j,l,k) + 0.5d0*eh_trip_Phi(i,j,l,k) + + 0.5d0*eh_sing_Phi(i,j,l,k) - 1.5d0*eh_trip_Phi(i,j,l,k) pp_sing_Gam_D(ij,kl) = pp_sing_Gam_D(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) @@ -86,7 +86,7 @@ subroutine R_pp_singlet_Gamma_C(nOrb,nO,nR,nVVs,eh_sing_Phi,eh_trip_Phi,pp_sing_ pp_sing_Gam_C(ab,cd) = 0.5d0*eh_sing_Phi(a,b,c,d) - 1.5d0*eh_trip_Phi(a,b,c,d) & - - 1.5d0*eh_sing_Phi(a,b,d,c) + 0.5d0*eh_trip_Phi(a,b,d,c) + + 0.5d0*eh_sing_Phi(a,b,d,c) - 1.5d0*eh_trip_Phi(a,b,d,c) pp_sing_Gam_C(ab,cd) = pp_sing_Gam_C(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) @@ -136,7 +136,7 @@ subroutine R_pp_singlet_Gamma_B(nOrb,nC,nO,nR,nOOs,nVVs,eh_sing_Phi,eh_trip_Phi, ij = ij +1 pp_sing_Gam_B(ab,ij) = 0.5d0*eh_sing_Phi(a,b,i,j) - 1.5d0*eh_trip_Phi(a,b,i,j) & - - 1.5d0*eh_sing_Phi(a,b,j,i) + 0.5d0*eh_trip_Phi(a,b,j,i) + + 0.5d0*eh_sing_Phi(a,b,j,i) - 1.5d0*eh_trip_Phi(a,b,j,i) pp_sing_Gam_B(ab,ij) = pp_sing_Gam_B(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) diff --git a/src/Parquet/R_screened_integrals.f90 b/src/Parquet/R_screened_integrals.f90 index 4e2ad64..0b619d3 100644 --- a/src/Parquet/R_screened_integrals.f90 +++ b/src/Parquet/R_screened_integrals.f90 @@ -17,7 +17,7 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr double precision :: X,Y ! Output variables - double precision,intent(out) :: rho(nOrb,nOrb,nS+nS) + double precision,intent(out) :: rho(nOrb,nOrb,nS) rho(:,:,:) = 0d0 ! !$OMP PARALLEL & @@ -27,24 +27,23 @@ subroutine R_eh_singlet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr ! !$OMP DO do q=nC+1,nOrb-nR do p=nC+1,nOrb-nR - - jb = 0 - do j=nC+1,nO - do b=nO+1,nOrb-nR - jb = jb + 1 - do ia=1,nS + do ia=1,nS + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) - rho(p,q,ia) = (2d0*ERI(q,j,p,b) - ERI(q,j,b,p) & - - 0.5d0*eh_sing_Phi(q,j,b,p) - 1.5d0*eh_trip_Phi(q,j,b,p) & - + 0.5d0*pp_sing_Phi(q,j,p,b) + 1.5d0*pp_trip_Phi(q,j,p,b)) * X - - rho(p,q,nS+ia) = (2d0*ERI(q,b,p,j) - ERI(q,b,j,p) & - - 0.5d0*eh_sing_Phi(q,b,j,p) - 1.5d0*eh_trip_Phi(q,b,j,p) & - + 0.5d0*pp_sing_Phi(q,b,p,j) + 1.5d0*pp_trip_Phi(q,b,p,j)) * Y + rho(p,q,ia) = rho(p,q,ia) + (2d0*ERI(q,j,p,b) - ERI(q,j,b,p) & + - 0d0*0.5d0*eh_sing_Phi(q,j,b,p) - 0d0*1.5d0*eh_trip_Phi(q,j,b,p) & + + 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*1.5d0*pp_trip_Phi(q,j,p,b)) * X & + + (2d0*ERI(q,b,p,j) - ERI(q,b,j,p) & + - 0d0*0.5d0*eh_sing_Phi(q,b,j,p) - 0d0*1.5d0*eh_trip_Phi(q,b,j,p) & + + 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*1.5d0*pp_trip_Phi(q,b,p,j)) * Y end do @@ -77,7 +76,7 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr double precision :: X,Y ! Output variables - double precision,intent(out) :: rho(nOrb,nOrb,nS+nS) + double precision,intent(out) :: rho(nOrb,nOrb,nS) rho(:,:,:) = 0d0 ! !$OMP PARALLEL & @@ -87,28 +86,27 @@ subroutine R_eh_triplet_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_sing_Phi,eh_tr ! !$OMP DO do q=nC+1,nOrb-nR do p=nC+1,nOrb-nR - - jb = 0 - do j=nC+1,nO - do b=nO+1,nOrb-nR - jb = jb + 1 - do ia=1,nS + do ia=1,nS + + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 X = 0.5d0*(XpY(ia,jb) + XmY(ia,jb)) Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) - rho(p,q,ia) = (- ERI(q,j,b,p) & - - 0.5d0*eh_sing_Phi(q,j,b,p) + 0.5d0*eh_trip_Phi(q,j,b,p) & - - 0.5d0*pp_sing_Phi(q,j,p,b) + 0.5d0*pp_trip_Phi(q,j,p,b)) * X - - rho(p,q,nS+ia) = (- ERI(q,b,j,p) & - - 0.5d0*eh_sing_Phi(q,b,j,p) + 0.5d0*eh_trip_Phi(q,b,j,p) & - - 0.5d0*pp_sing_Phi(q,b,p,j) + 0.5d0*pp_trip_Phi(q,b,p,j)) * Y + rho(p,q,ia) = rho(p,q,ia) + (- ERI(q,j,b,p) & + - 0d0*0.5d0*eh_sing_Phi(q,j,b,p) + 0d0*0.5d0*eh_trip_Phi(q,j,b,p) & + - 0d0*0.5d0*pp_sing_Phi(q,j,p,b) + 0d0*0.5d0*pp_trip_Phi(q,j,p,b)) * X & + + (- ERI(q,b,j,p) & + - 0d0*0.5d0*eh_sing_Phi(q,b,j,p) + 0d0*0.5d0*eh_trip_Phi(q,b,j,p) & + - 0d0*0.5d0*pp_sing_Phi(q,b,p,j) + 0d0*0.5d0*pp_trip_Phi(q,b,p,j)) * Y end do - end do + end do end do @@ -175,9 +173,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, do d = c, nOrb-nR cd = cd + 1 - rho1(p,q,ab) = (ERI(p,q,c,d) + ERI(p,q,d,c) & - + 0.5d0*eh_sing_Phi(p,q,c,d) - 1.5d0*eh_trip_Phi(p,q,c,d) & - - 1.5d0*eh_sing_Phi(p,q,d,c) + 0.5d0*eh_trip_Phi(p,q,d,c))& + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) + ERI(p,q,d,c) & + + 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) & + + 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))& *X1(cd,ab)/sqrt(1d0 + Kronecker_delta(c,d)) end do ! d @@ -189,9 +187,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, kl = kl + 1 - rho1(p,q,ab) = (ERI(p,q,k,l) + ERI(p,q,l,k) & - + 0.5d0*eh_sing_Phi(p,q,k,l) - 1.5d0*eh_trip_Phi(p,q,k,l) & - - 1.5d0*eh_sing_Phi(p,q,l,k) + 0.5d0*eh_trip_Phi(p,q,l,k))& + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) + ERI(p,q,l,k) & + + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) & + + 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))& *Y1(kl,ab)/sqrt(1d0 + Kronecker_delta(k,l)) end do ! l end do ! k @@ -208,9 +206,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, do d = c, nOrb-nR cd = cd + 1 - rho2(p,q,ij) = (ERI(p,q,c,d) + ERI(p,q,d,c) & - + 0.5d0*eh_sing_Phi(p,q,c,d) - 1.5d0*eh_trip_Phi(p,q,c,d) & - - 1.5d0*eh_sing_Phi(p,q,d,c) + 0.5d0*eh_trip_Phi(p,q,d,c))& + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) + ERI(p,q,d,c) & + + 0d0*0.5d0*eh_sing_Phi(p,q,c,d) - 0d0*1.5d0*eh_trip_Phi(p,q,c,d) & + + 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*1.5d0*eh_trip_Phi(p,q,d,c))& *X2(cd,ij)/sqrt(1d0 + Kronecker_delta(c,d)) end do ! d end do ! c @@ -220,9 +218,9 @@ subroutine R_pp_singlet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, do l = k, nO kl = kl + 1 - rho2(p,q,ij) = (ERI(p,q,k,l) + ERI(p,q,l,k) & - + 0.5d0*eh_sing_Phi(p,q,k,l) - 1.5d0*eh_trip_Phi(p,q,k,l) & - - 1.5d0*eh_sing_Phi(p,q,l,k) + 0.5d0*eh_trip_Phi(p,q,l,k))& + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) + ERI(p,q,l,k) & + + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) - 0d0*1.5d0*eh_trip_Phi(p,q,k,l) & + + 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*1.5d0*eh_trip_Phi(p,q,l,k))& *Y2(kl,ij)/sqrt(1d0 + Kronecker_delta(k,l)) end do ! l end do ! k @@ -292,9 +290,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, do d = c+1, nOrb-nR cd = cd + 1 - rho1(p,q,ab) = (ERI(p,q,c,d) - ERI(p,q,d,c) & - + 0.5d0*eh_sing_Phi(p,q,c,d) + 0.5d0*eh_trip_Phi(p,q,c,d) & - - 0.5d0*eh_sing_Phi(p,q,d,c) - 0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab) + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c) & + + 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) & + - 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X1(cd,ab) end do ! d end do ! c @@ -305,9 +303,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, kl = kl + 1 - rho1(p,q,ab) = (ERI(p,q,k,l) - ERI(p,q,l,k) & - + 0.5d0*eh_sing_Phi(p,q,k,l) + 0.5d0*eh_trip_Phi(p,q,k,l) & - - 0.5d0*eh_sing_Phi(p,q,l,k) - 0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab) + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k) & + + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) & + - 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y1(kl,ab) end do ! l end do ! k end do ! b @@ -323,9 +321,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, do d = c+1, nOrb-nR cd = cd + 1 - rho2(p,q,ij) = (ERI(p,q,c,d) - ERI(p,q,d,c) & - + 0.5d0*eh_sing_Phi(p,q,c,d) + 0.5d0*eh_trip_Phi(p,q,c,d) & - - 0.5d0*eh_sing_Phi(p,q,d,c) - 0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij) + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c) & + + 0d0*0.5d0*eh_sing_Phi(p,q,c,d) + 0d0*0.5d0*eh_trip_Phi(p,q,c,d) & + - 0d0*0.5d0*eh_sing_Phi(p,q,d,c) - 0d0*0.5d0*eh_trip_Phi(p,q,d,c) )*X2(cd,ij) end do ! d end do ! c @@ -334,9 +332,9 @@ subroutine R_pp_triplet_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_sing_Phi, do l = k+1, nO kl = kl + 1 - rho2(p,q,ij) = (ERI(p,q,k,l) - ERI(p,q,l,k) & - + 0.5d0*eh_sing_Phi(p,q,k,l) + 0.5d0*eh_trip_Phi(p,q,k,l) & - - 0.5d0*eh_sing_Phi(p,q,l,k) - 0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij) + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k) & + + 0d0*0.5d0*eh_sing_Phi(p,q,k,l) + 0d0*0.5d0*eh_trip_Phi(p,q,k,l) & + - 0d0*0.5d0*eh_sing_Phi(p,q,l,k) - 0d0*0.5d0*eh_trip_Phi(p,q,l,k) )*Y2(kl,ij) end do ! l end do ! k