PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:42 +01:00
Code Releases Activity

ok with RG0T0pp now

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-12 22:01:06 +02:00
parent 8ec7d3d170
commit 647f27fb62
14 changed files with 38 additions and 182 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
src/GT/GG0T0pp.f90
View File

@ -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

2
src/GT/RGTpp_excitation_density.f90
View File

@ -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

1
src/LR/ppLR.f90
View File

@ -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)

2
src/LR/ppLR_B.f90
View File

@ -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

2
src/LR/ppLR_C.f90
View File

@ -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) &

2
src/LR/ppLR_D.f90
View File

@ -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

7
src/RPA/crACFDT.f90
View File

@ -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

5
src/RPA/crRRPA.f90
View File

@ -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(*,*)

18
src/RPA/phACFDT.f90
View File

@ -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))

5
src/RPA/phRRPA.f90
View File

@ -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(*,*)

5
src/RPA/phRRPAx.f90
View File

@ -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(*,*)

7
src/RPA/ppACFDT.f90
View File

@ -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))

5
src/RPA/ppRRPA.f90
View File

@ -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(*,*)

85
src/RPA/sort_ppRPA.f90
View File

@ -114,91 +114,6 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
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))