From 647f27fb620daa90f81b6342cc59f8bdb8739982 Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 12 Sep 2024 22:01:06 +0200 Subject: [PATCH] ok with RG0T0pp now --- src/GT/GG0T0pp.f90 | 72 ------------------------ src/GT/RGTpp_excitation_density.f90 | 2 +- src/LR/ppLR.f90 | 3 +- src/LR/ppLR_B.f90 | 2 +- src/LR/ppLR_C.f90 | 2 +- src/LR/ppLR_D.f90 | 2 +- src/RPA/crACFDT.f90 | 7 +++ src/RPA/crRRPA.f90 | 5 -- src/RPA/phACFDT.f90 | 18 ++++++ src/RPA/phRRPA.f90 | 5 -- src/RPA/phRRPAx.f90 | 5 -- src/RPA/ppACFDT.f90 | 7 +++ src/RPA/ppRRPA.f90 | 5 -- src/RPA/sort_ppRPA.f90 | 85 ----------------------------- 14 files changed, 38 insertions(+), 182 deletions(-) diff --git a/src/GT/GG0T0pp.f90 b/src/GT/GG0T0pp.f90 index f7ced74..c0012d4 100644 --- a/src/GT/GG0T0pp.f90 +++ b/src/GT/GG0T0pp.f90 @@ -152,8 +152,6 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d end if -! call GGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOO,nVV,eHF,eGT,Om1,rho1,Om2,rho2) - !---------------------------------------------- ! Dump results !---------------------------------------------- @@ -172,76 +170,6 @@ subroutine GG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d call print_GG0T0pp(nOrb,nO,eHF,ENuc,EGHF,Sig,Z,eGT,EcGM,EcRPA) -! Perform BSE calculation - -! if(dophBSE) then - -! call GGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV, & -! Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,dipole_int,eHF,eGT,EcBSE) - -! if(exchange_kernel) then - -! EcBSE = 0.5d0*EcBSE - -! end if - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@phBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! Compute the BSE correlation energy via the adiabatic connection - -! if(doACFDT) then - -! write(*,*) '--------------------------------------------------------' -! write(*,*) 'Adiabatic connection version of phBSE correlation energy' -! write(*,*) '--------------------------------------------------------' -! write(*,*) - -! if(doXBS) then - -! write(*,*) '*** scaled screening version (XBS) ***' -! write(*,*) - -! end if - -! call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,eta,nOrb,nC,nO,nV,nR,nS, & -! nOO,nVV,Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,eHF,eGT,EcBSE) - -! if(exchange_kernel) then - -! EcBSE = 0.5d0*EcBSE - -! end if - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - -! end if - -! if(doppBSE) then - -! call GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV, & -! Om1,X1,Y1,Om2,X2,Y2,rho1,rho2,ERI,dipole_int,eHF,eGT,EcBSE) - -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF correlation energy = ',EcBSE,' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0T0pp@GHF total energy = ',ENuc + EGHF + EcBSE,' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) - -! end if - ! Testing zone if(dotest) then diff --git a/src/GT/RGTpp_excitation_density.f90 b/src/GT/RGTpp_excitation_density.f90 index 821db82..e894d02 100644 --- a/src/GT/RGTpp_excitation_density.f90 +++ b/src/GT/RGTpp_excitation_density.f90 @@ -130,7 +130,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho ! Triplet manifold !---------------------------------------------- - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then dim_1 = (nBas - nO) * (nBas - nO - 1) / 2 dim_2 = nO * (nO - 1) / 2 diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 63f5a97..6d62908 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -111,6 +111,7 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) end if ! Compute the RPA correlation energy + EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp)) EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp) EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) @@ -119,7 +120,7 @@ subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' endif - deallocate(M, Z, Om) + deallocate(M,Z,Om) end subroutine diff --git a/src/LR/ppLR_B.f90 b/src/LR/ppLR_B.f90 index bfb3734..30a093d 100644 --- a/src/LR/ppLR_B.f90 +++ b/src/LR/ppLR_B.f90 @@ -51,7 +51,7 @@ subroutine ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) ! Build the B matrix for the triplet or alpha-alpha manifold - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then ab = 0 do a=nO+1,nOrb-nR diff --git a/src/LR/ppLR_C.f90 b/src/LR/ppLR_C.f90 index 764a0ed..fbc3900 100644 --- a/src/LR/ppLR_C.f90 +++ b/src/LR/ppLR_C.f90 @@ -108,7 +108,7 @@ subroutine ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) ! Build C matrix for the triplet or alpha-alpha manifold - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then !$OMP PARALLEL & !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nOrb,nR) & !$OMP PRIVATE(c,d,a,b,ab,cd) & diff --git a/src/LR/ppLR_D.f90 b/src/LR/ppLR_D.f90 index ee759a3..3744f5c 100644 --- a/src/LR/ppLR_D.f90 +++ b/src/LR/ppLR_D.f90 @@ -58,7 +58,7 @@ subroutine ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) ! Build the D matrix for the triplet or alpha-alpha manifold - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then ij = 0 do i=nC+1,nO diff --git a/src/RPA/crACFDT.f90 b/src/RPA/crACFDT.f90 index 9e82467..296d9f5 100644 --- a/src/RPA/crACFDT.f90 +++ b/src/RPA/crACFDT.f90 @@ -48,6 +48,13 @@ subroutine crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS, allocate(Ec(nAC,nspin)) allocate(Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) +! Hello world + + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of crRPA correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) + ! Antisymmetrized kernel version if(exchange_kernel) then diff --git a/src/RPA/crRRPA.f90 b/src/RPA/crRRPA.f90 index 01febce..e2c64df 100644 --- a/src/RPA/crRRPA.f90 +++ b/src/RPA/crRRPA.f90 @@ -113,11 +113,6 @@ subroutine crRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, if(doACFDT) then - write(*,*) '-------------------------------------------------------' - write(*,*) 'Adiabatic connection version of crRPA correlation energy' - write(*,*) '-------------------------------------------------------' - write(*,*) - call crACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,eHF,EcRPA) write(*,*) diff --git a/src/RPA/phACFDT.f90 b/src/RPA/phACFDT.f90 index 2cb024c..c6e9e2d 100644 --- a/src/RPA/phACFDT.f90 +++ b/src/RPA/phACFDT.f90 @@ -42,6 +42,24 @@ subroutine phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS, double precision,intent(out) :: EcAC(nspin) +! Hello world + + if(dRPA) then + + write(*,*) '--------------------------------------------------------' + write(*,*) 'Adiabatic connection version of phRPA correlation energy' + write(*,*) '--------------------------------------------------------' + write(*,*) + + else + + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of crRPA correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) + + end if + ! Memory allocation allocate(Ec(nAC,nspin)) diff --git a/src/RPA/phRRPA.f90 b/src/RPA/phRRPA.f90 index 27a12ab..7a13d8f 100644 --- a/src/RPA/phRRPA.f90 +++ b/src/RPA/phRRPA.f90 @@ -117,11 +117,6 @@ subroutine phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO, if(doACFDT) then - write(*,*) '--------------------------------------------------------' - write(*,*) 'Adiabatic connection version of phRPA correlation energy' - write(*,*) '--------------------------------------------------------' - write(*,*) - call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,eHF,EcRPA) write(*,*) diff --git a/src/RPA/phRRPAx.f90 b/src/RPA/phRRPAx.f90 index 2f27eb5..491ea12 100644 --- a/src/RPA/phRRPAx.f90 +++ b/src/RPA/phRRPAx.f90 @@ -118,11 +118,6 @@ subroutine phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO if(doACFDT) then - write(*,*) '---------------------------------------------------------' - write(*,*) 'Adiabatic connection version of phRPAx correlation energy' - write(*,*) '---------------------------------------------------------' - write(*,*) - call phACFDT(exchange_kernel,dRPA,TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcRPA) write(*,*) diff --git a/src/RPA/ppACFDT.f90 b/src/RPA/ppACFDT.f90 index ddb606e..0bd32d6 100644 --- a/src/RPA/ppACFDT.f90 +++ b/src/RPA/ppACFDT.f90 @@ -47,6 +47,13 @@ subroutine ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,e,EcAC) double precision,intent(out) :: EcAC(nspin) +! Hello world + + write(*,*) '--------------------------------------------------------' + write(*,*) 'Adiabatic connection version of ppRPA correlation energy' + write(*,*) '--------------------------------------------------------' + write(*,*) + ! Memory allocation allocate(Ec(nAC,nspin)) diff --git a/src/RPA/ppRRPA.f90 b/src/RPA/ppRRPA.f90 index 2fe9658..24e17aa 100644 --- a/src/RPA/ppRRPA.f90 +++ b/src/RPA/ppRRPA.f90 @@ -132,11 +132,6 @@ subroutine ppRRPA(dotest,TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF, if(doACFDT) then - write(*,*) '--------------------------------------------------------' - write(*,*) 'Adiabatic connection version of ppRPA correlation energy' - write(*,*) '--------------------------------------------------------' - write(*,*) - call ppACFDT(TDA,singlet,triplet,nBas,nC,nO,nV,nR,ERI,eHF,EcRPA) write(*,*) diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 5c6714a..24c6f59 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -113,91 +113,6 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) call set_order(Z2,order2,nOO+nVV,nOO) end if - - -! Orthogonalize eigenvectors - -! deg1 = 1 -! ab_start = 1 - -! do ab=2,nVV - -! if(abs(Om1(ab-1) - Om1(ab)) < 1d-6) then - -! deg1 = deg1 + 1 - -! if(ab == nVV) then - -! ab_end = ab - -! allocate(S1(deg1,deg1),O1(deg1,deg1)) - -! S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end))) -! call orthogonalization_matrix(1,deg1,S1,O1) -! Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1) - -! deallocate(S1,O1) - -! end if - -! else - -! ab_end = ab - 1 - -! allocate(S1(deg1,deg1),O1(deg1,deg1)) - -! S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end))) -! call orthogonalization_matrix(1,deg1,S1,O1) -! Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1) - -! deallocate(S1,O1) - -! deg1 = 1 -! ab_start = ab - -! end if -! end do - -! deg2 = 1 -! ij_start = 1 - -! do ij=2,nOO - -! if(abs(Om2(ij-1) - Om2(ij)) < 1d-6) then - -! deg2 = deg2 + 1 - -! if(ij == nOO) then - -! ij_end = ij -! -! allocate(S2(deg2,deg2),O2(deg2,deg2)) -! -! S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end))) -! call orthogonalization_matrix(1,deg2,S2,O2) -! Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2) -! -! deallocate(S2,O2) - -! end if - -! else - -! ij_end = ij - 1 - -! allocate(S2(deg2,deg2),O2(deg2,deg2)) - -! S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end))) -! call orthogonalization_matrix(1,deg2,S2,O2) -! Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2) - -! deallocate(S2,O2) - -! deg2 = 1 -! ij_start = ij - -! end if -! end do allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO)) allocate(tmp1(nOO+nVV,nVV),tmp2(nOO+nVV,nOO))