Skip to content

Commit 684c3f4

Browse files
author
lburth
committed
RPA@MOM-ref in real branch
1 parent 63145f7 commit 684c3f4

File tree

8 files changed

+386
-10
lines changed

8 files changed

+386
-10
lines changed

‎src/LR/CVS_phRLR.f90‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
subroutine CVS_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
2+
3+
! Compute linear response for complex unrestricted formalism using the c-product
4+
5+
implicit none
6+
include 'parameters.h'
7+
8+
! Input variables
9+
10+
logical,intent(in) :: TDA
11+
integer,intent(in) :: nS
12+
double precision,intent(in) :: Aph(nS,nS)
13+
double precision,intent(in) :: Bph(nS,nS)
14+
15+
! Local variables
16+
17+
double precision,external :: trace_matrix
18+
complex*16,allocatable :: RPA_matrix(:,:)
19+
complex*16,allocatable :: OmOmminus(:)
20+
21+
! Output variables
22+
23+
double precision,intent(out) :: EcRPA
24+
double precision,intent(out) :: Om(nS)
25+
double precision,intent(out) :: XpY(nS,nS)
26+
double precision,intent(out) :: XmY(nS,nS)
27+
28+
29+
30+
! Tamm-Dancoff approximation
31+
32+
if(TDA) then
33+
34+
XpY(:,:) = Aph(:,:)
35+
call diagonalize_matrix(nS,XpY,Om)
36+
XpY(:,:) = transpose(XpY(:,:))
37+
XmY(:,:) = XpY(:,:)
38+
39+
else
40+
41+
allocate(RPA_matrix(2*nS,2*nS),OmOmminus(2*nS))
42+
RPA_matrix(1:nS,1:nS) = cmplx(Aph(:,:),0d0,kind=8)
43+
RPA_matrix(1:nS,nS+1:2*nS) = cmplx(Bph(:,:),0d0,kind=8)
44+
RPA_matrix(nS+1:2*nS,1:nS) = cmplx(-Bph(:,:),0d0,kind=8)
45+
RPA_matrix(nS+1:2*nS,nS+1:2*nS) = cmplx(-Aph(:,:),0d0,kind=8)
46+
call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus)
47+
call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix)
48+
call complex_normalize_RPA(nS,RPA_matrix)
49+
50+
Om(:) = real(OmOmminus(1:nS))
51+
if(maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS))) > 1e-12) then
52+
call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.')
53+
write(*,*) "Maximal difference :", maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS)))
54+
end if
55+
if(minval(real(OmOmminus(:))) < 0d0) &
56+
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
57+
if(maxval(abs(aimag(OmOmminus(:))))>1d-8)&
58+
call print_warning('You may have instabilities in linear response: complex excitation eigenvalue !')
59+
XpY(:,:) = transpose(real(RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,1:nS)))
60+
XmY(:,:) = transpose(real(RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,1:nS)))
61+
62+
deallocate(RPA_matrix,OmOmminus)
63+
64+
end if
65+
66+
! Compute the RPA correlation energy
67+
68+
EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,Aph))
69+
70+
end subroutine

‎src/LR/CVS_phULR.f90‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
subroutine CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
2+
3+
! Compute linear response for complex unrestricted formalism using the c-product
4+
5+
implicit none
6+
include 'parameters.h'
7+
8+
! Input variables
9+
10+
logical,intent(in) :: TDA
11+
integer,intent(in) :: nSa
12+
integer,intent(in) :: nSb
13+
integer,intent(in) :: nSt
14+
double precision,intent(in) :: Aph(nSt,nSt)
15+
double precision,intent(in) :: Bph(nSt,nSt)
16+
17+
! Local variables
18+
19+
double precision,external :: trace_matrix
20+
complex*16,allocatable :: RPA_matrix(:,:)
21+
complex*16,allocatable :: OmOmminus(:)
22+
23+
! Output variables
24+
25+
double precision,intent(out) :: EcRPA
26+
double precision,intent(out) :: Om(nSt)
27+
double precision,intent(out) :: XpY(nSt,nSt)
28+
double precision,intent(out) :: XmY(nSt,nSt)
29+
30+
! Tamm-Dancoff approximation
31+
32+
if(TDA) then
33+
34+
XpY(:,:) = Aph(:,:)
35+
call diagonalize_matrix(nSt,XpY,Om)
36+
XpY(:,:) = transpose(XpY(:,:))
37+
XmY(:,:) = XpY(:,:)
38+
39+
else
40+
41+
allocate(RPA_matrix(2*nSt,2*nSt),OmOmminus(2*nSt))
42+
RPA_matrix(1:nSt,1:nSt) = cmplx(Aph(:,:),0d0,kind=8)
43+
RPA_matrix(1:nSt,nSt+1:2*nSt) = cmplx(Bph(:,:),0d0,kind=8)
44+
RPA_matrix(nSt+1:2*nSt,1:nSt) = cmplx(-Bph(:,:),0d0,kind=8)
45+
RPA_matrix(nSt+1:2*nSt,nSt+1:2*nSt) = cmplx(-Aph(:,:),0d0,kind=8)
46+
call complex_diagonalize_matrix_without_sort(2*nSt,RPA_matrix,OmOmminus)
47+
call complex_sort_eigenvalues_RPA(2*nSt,OmOmminus,RPA_matrix)
48+
call complex_normalize_RPA(nSt,RPA_matrix)
49+
50+
Om(:) = real(OmOmminus(1:nSt))
51+
if(maxval(abs(OmOmminus(1:nSt)+OmOmminus(nSt+1:2*nSt))) > 1e-12) then
52+
call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.')
53+
write(*,*) "Maximal difference :", maxval(abs(OmOmminus(1:nSt)+OmOmminus(nSt+1:2*nSt)))
54+
end if
55+
if(minval(real(OmOmminus(:))) < 0d0) &
56+
call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
57+
if(maxval(abs(aimag(OmOmminus(:))))>1d-8)&
58+
call print_warning('You may have instabilities in linear response: complex excitation eigenvalue !')
59+
XpY(:,:) = transpose(real(RPA_matrix(1:nSt,1:nSt) + RPA_matrix(nSt+1:2*nSt,1:nSt)))
60+
XmY(:,:) = transpose(real(RPA_matrix(1:nSt,1:nSt) - RPA_matrix(nSt+1:2*nSt,1:nSt)))
61+
62+
deallocate(RPA_matrix,OmOmminus)
63+
64+
end if
65+
66+
! Compute the RPA correlation energy
67+
68+
EcRPA = 0.5d0*(sum(Om) - trace_matrix(nSt,Aph))
69+
70+
end subroutine

‎src/RPA/CVS_phRRPA.f90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ subroutine CVS_phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC
119119
call CVS_phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,eHF,ERI,Aph)
120120
if(.not.TDA) call CVS_phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,ERI,Bph)
121121

122-
call phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
122+
call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
123123
call print_excitation_energies('phRPA@RHF','singlet',nSCVS,Om)
124124
call CVS_phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
125125

@@ -134,7 +134,7 @@ subroutine CVS_phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC
134134
call CVS_phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,eHF,ERI,Aph)
135135
if(.not.TDA) call CVS_phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,ERI,Bph)
136136

137-
call phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
137+
call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
138138
call print_excitation_energies('phRPA@RHF','triplet',nSCVS,Om)
139139
call CVS_phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
140140

‎src/RPA/CVS_phRRPAx.f90‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
subroutine CVS_phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,nCVS,FC,occupations,ENuc,ERHF,ERI,dipole_int,eHF)
2+
3+
! Perform a direct random phase approximation calculation
4+
5+
implicit none
6+
include 'parameters.h'
7+
include 'quadrature.h'
8+
9+
! Input variables
10+
11+
logical,intent(in) :: dotest
12+
13+
logical,intent(in) :: TDA
14+
logical,intent(in) :: doACFDT
15+
logical,intent(in) :: exchange_kernel
16+
logical,intent(in) :: singlet
17+
logical,intent(in) :: triplet
18+
integer,intent(in) :: nBas
19+
integer,intent(in) :: nC
20+
integer,intent(in) :: nO
21+
integer,intent(in) :: nV
22+
integer,intent(in) :: nR
23+
integer,intent(in) :: nS
24+
integer,intent(in) :: nCVS
25+
integer,intent(in) :: FC
26+
integer,intent(in) :: occupations(nO)
27+
double precision,intent(in) :: ENuc
28+
double precision,intent(in) :: ERHF
29+
double precision,intent(in) :: eHF(nBas)
30+
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
31+
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
32+
33+
! Local variables
34+
35+
integer :: ia,i
36+
integer :: ispin
37+
integer :: nSCVS, nFC
38+
logical :: dRPA,found
39+
double precision :: t1, t2
40+
double precision :: lambda
41+
double precision,allocatable :: Aph(:,:)
42+
double precision,allocatable :: Bph(:,:)
43+
double precision,allocatable :: Om(:)
44+
double precision,allocatable :: XpY(:,:)
45+
double precision,allocatable :: XmY(:,:)
46+
integer,allocatable :: virtuals(:)
47+
integer,allocatable :: occupations_fc(:)
48+
49+
double precision :: EcRPA(nspin)
50+
51+
! Hello world
52+
53+
write(*,*)
54+
write(*,*)'*********************************'
55+
write(*,*)'* Restricted ph-RPA Calculation *'
56+
write(*,*)'*********************************'
57+
write(*,*)
58+
59+
! CVS
60+
61+
print *, "No exications to the first", nCVS, "orbital(s) are considered."
62+
if(nC/=0) then
63+
print *, "Do not use PyDuck frozen core with CVS !"
64+
stop
65+
end if
66+
67+
! Frozen Core
68+
69+
nFC = MERGE(1,0,FC/=0)
70+
allocate(occupations_fc(nO-nFC))
71+
! remove FC from occupations
72+
do ispin=1,nspin
73+
occupations_fc(1:nO-nFC) = occupations(1:nO - nFC)
74+
found = .false.
75+
do i=1,nO-1
76+
if(.not. found) then
77+
if(occupations(i)==FC) then
78+
found = .true.
79+
occupations_fc(i) = occupations(i+1)
80+
else
81+
occupations_fc(i) = occupations(i)
82+
endif
83+
else
84+
occupations_fc(i) = occupations(i+1)
85+
endif
86+
enddo
87+
enddo
88+
print *, "Not Frozen orbitals:"
89+
print *,occupations_fc(1:nO-nFC)
90+
91+
! TDA
92+
93+
if(TDA) then
94+
write(*,*) 'Tamm-Dancoff approximation activated!'
95+
write(*,*)
96+
end if
97+
98+
! Initialization
99+
100+
dRPA = .false.
101+
EcRPA(:) = 0d0
102+
lambda = 1d0
103+
104+
! Memory allocation
105+
if(nC/=0) then
106+
print *, "No frozen core with CVS available !"
107+
stop
108+
end if
109+
nSCVS = (nV - nCVS)*(nO - nFC)
110+
allocate(Om(nSCVS),XpY(nSCVS,nSCVS),XmY(nSCVS,nSCVS),Aph(nSCVS,nSCVS),Bph(nSCVS,nSCVS),virtuals(nV))
111+
call non_occupied(nO,nBas,occupations,virtuals)
112+
113+
! Singlet manifold
114+
115+
if(singlet) then
116+
117+
ispin = 1
118+
119+
call CVS_phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,eHF,ERI,Aph)
120+
if(.not.TDA) call CVS_phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,ERI,Bph)
121+
122+
call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
123+
call print_excitation_energies('phRPAx@RHF','singlet',nSCVS,Om)
124+
call CVS_phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
125+
126+
end if
127+
128+
! Triplet manifold
129+
130+
if(triplet) then
131+
132+
ispin = 2
133+
134+
call CVS_phRLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,eHF,ERI,Aph)
135+
if(.not.TDA) call CVS_phRLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,lambda,ERI,Bph)
136+
137+
call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
138+
call print_excitation_energies('phRPAx@RHF','triplet',nSCVS,Om)
139+
call CVS_phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
140+
141+
end if
142+
143+
if(exchange_kernel) then
144+
145+
EcRPA(1) = 0.5d0*EcRPA(1)
146+
EcRPA(2) = 1.5d0*EcRPA(2)
147+
148+
end if
149+
150+
write(*,*)
151+
write(*,*)'-------------------------------------------------------------------------------'
152+
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPA@RHF correlation energy (singlet) = ',EcRPA(1),' au'
153+
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPA@RHF correlation energy (triplet) = ',EcRPA(2),' au'
154+
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPA@RHF correlation energy = ',sum(EcRPA),' au'
155+
write(*,'(2X,A50,F20.10,A3)') 'Tr@phRPA@RHF total energy = ',ENuc + ERHF + sum(EcRPA),' au'
156+
write(*,*)'-------------------------------------------------------------------------------'
157+
write(*,*)
158+
159+
deallocate(Om,XpY,XmY,Aph,Bph)
160+
161+
! Compute the correlation energy via the adiabatic connection
162+
163+
if(doACFDT) then
164+
print *,"No ADC implemented with CVS."
165+
end if
166+
167+
if(dotest) then
168+
169+
call dump_test_value('R','phRPA correlation energy',sum(EcRPA))
170+
171+
end if
172+
173+
end subroutine

‎src/RPA/CVS_phURPA.f90‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -137,7 +137,7 @@ subroutine CVS_phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fli
137137
if(.not.TDA) call CVS_phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nCVS,nFC,occupations_fc, virtuals,&
138138
lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
139139

140-
call phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
140+
call CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
141141
call print_excitation_energies('phRPA@UHF','spin-conserved',nSt,Om)
142142
call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, &
143143
dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
@@ -165,7 +165,7 @@ subroutine CVS_phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fli
165165
if(.not.TDA) call CVS_phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nCVS,nFC,occupations_fc,virtuals,&
166166
lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
167167

168-
call phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
168+
call CVS_phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
169169
call print_excitation_energies('phRPA@UHF','spin-flip',nSt,Om)
170170
call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
171171

‎src/RPA/CVS_phURPAx.f90‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -139,10 +139,10 @@ subroutine CVS_phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fl
139139
if(.not.TDA) call CVS_phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nCVS,nFC,occupations_fc, virtuals,&
140140
lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
141141

142-
call phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
143-
call print_excitation_energies('phRPA@UHF','spin-conserved',nSt,Om)
142+
call CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
143+
call print_excitation_energies('phRPAx@UHF','spin-conserved',nSt,Om)
144144
call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, &
145-
dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
145+
dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
146146

147147
deallocate(Aph,Bph,Om,XpY,XmY)
148148

@@ -167,8 +167,8 @@ subroutine CVS_phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fl
167167
if(.not.TDA) call CVS_phULR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nCVS,nFC,occupations_fc,virtuals,&
168168
lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Bph)
169169

170-
call phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
171-
call print_excitation_energies('phRPA@UHF','spin-flip',nSt,Om)
170+
call CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
171+
call print_excitation_energies('phRPAx@UHF','spin-flip',nSt,Om)
172172
call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
173173

174174
deallocate(Aph,Bph,Om,XpY,XmY)
@@ -196,6 +196,7 @@ subroutine CVS_phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fl
196196
write(*,*)'-------------------------------------------------------------------------------'
197197
write(*,*)
198198

199+
199200
! Compute the correlation energy via the adiabatic connection
200201

201202
if(doACFDT) then

‎src/RPA/RRPA.f90‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,11 @@ subroutine RRPA(use_gpu,dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exch
7777
if(dophRPAx) then
7878

7979
call wall_time(start_RPA)
80-
call phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
80+
if(CVS) then
81+
call CVS_phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,nCVS,FC,occupations,ENuc,ERHF,ERI,dipole_int,eHF)
82+
else
83+
call phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
84+
end if
8185
call wall_time(end_RPA)
8286

8387
t_RPA = end_RPA - start_RPA

0 commit comments

Comments
 (0)