mirror of
https://github.com/pfloos/quack
synced 2024-06-18 11:15:30 +02:00
ppRPA done
This commit is contained in:
parent
f8a9321cf0
commit
0c6f1ee7d5
43
input/basis
43
input/basis
|
|
@ -1,16 +1,37 @@
|
|||
1 6
|
||||
S 4 1.00
|
||||
234.0000000 0.0025870
|
||||
35.1600000 0.0195330
|
||||
7.9890000 0.0909980
|
||||
2.2120000 0.2720500
|
||||
1 10
|
||||
S 8 1.00
|
||||
24350.0000000 0.0005020
|
||||
3650.0000000 0.0038810
|
||||
829.6000000 0.0199970
|
||||
234.0000000 0.0784180
|
||||
75.6100000 0.2296760
|
||||
26.7300000 0.4327220
|
||||
9.9270000 0.3506420
|
||||
1.1020000 -0.0076450
|
||||
S 8 1.00
|
||||
24350.0000000 -0.0001180
|
||||
3650.0000000 -0.0009150
|
||||
829.6000000 -0.0047370
|
||||
234.0000000 -0.0192330
|
||||
75.6100000 -0.0603690
|
||||
26.7300000 -0.1425080
|
||||
9.9270000 -0.1777100
|
||||
1.1020000 0.6058360
|
||||
S 1 1.00
|
||||
0.6669000 1.0000000
|
||||
2.8360000 1.0000000
|
||||
S 1 1.00
|
||||
0.2089000 1.0000000
|
||||
0.3782000 1.0000000
|
||||
P 3 1.00
|
||||
54.7000000 0.0171510
|
||||
12.4300000 0.1076560
|
||||
3.6790000 0.3216810
|
||||
P 1 1.00
|
||||
3.0440000 1.0000000
|
||||
1.1430000 1.0000000
|
||||
P 1 1.00
|
||||
0.7580000 1.0000000
|
||||
0.3300000 1.0000000
|
||||
D 1 1.00
|
||||
1.9650000 1.0000000
|
||||
4.0140000 1.0000000
|
||||
D 1 1.00
|
||||
1.0960000 1.0000000
|
||||
F 1 1.00
|
||||
2.5440000 1.0000000
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
# nAt nEla nElb nCore nRyd
|
||||
1 1 1 0 0
|
||||
1 5 5 0 0
|
||||
# Znuc x y z
|
||||
He 0.0 0.0 0.0
|
||||
Ne 0.0 0.0 0.0
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
# MP:
|
||||
|
||||
# CC: maxSCF thresh DIIS n_diis
|
||||
64 0.00001 F 1
|
||||
64 0.0000001 F 1
|
||||
# CIS/TDHF/BSE: singlet triplet
|
||||
T T
|
||||
# GF: maxSCF thresh DIIS n_diis renormalization
|
||||
|
|
|
|||
43
input/weight
43
input/weight
|
|
@ -1,16 +1,37 @@
|
|||
1 6
|
||||
S 4 1.00
|
||||
234.0000000 0.0025870
|
||||
35.1600000 0.0195330
|
||||
7.9890000 0.0909980
|
||||
2.2120000 0.2720500
|
||||
1 10
|
||||
S 8 1.00
|
||||
24350.0000000 0.0005020
|
||||
3650.0000000 0.0038810
|
||||
829.6000000 0.0199970
|
||||
234.0000000 0.0784180
|
||||
75.6100000 0.2296760
|
||||
26.7300000 0.4327220
|
||||
9.9270000 0.3506420
|
||||
1.1020000 -0.0076450
|
||||
S 8 1.00
|
||||
24350.0000000 -0.0001180
|
||||
3650.0000000 -0.0009150
|
||||
829.6000000 -0.0047370
|
||||
234.0000000 -0.0192330
|
||||
75.6100000 -0.0603690
|
||||
26.7300000 -0.1425080
|
||||
9.9270000 -0.1777100
|
||||
1.1020000 0.6058360
|
||||
S 1 1.00
|
||||
0.6669000 1.0000000
|
||||
2.8360000 1.0000000
|
||||
S 1 1.00
|
||||
0.2089000 1.0000000
|
||||
0.3782000 1.0000000
|
||||
P 3 1.00
|
||||
54.7000000 0.0171510
|
||||
12.4300000 0.1076560
|
||||
3.6790000 0.3216810
|
||||
P 1 1.00
|
||||
3.0440000 1.0000000
|
||||
1.1430000 1.0000000
|
||||
P 1 1.00
|
||||
0.7580000 1.0000000
|
||||
0.3300000 1.0000000
|
||||
D 1 1.00
|
||||
1.9650000 1.0000000
|
||||
4.0140000 1.0000000
|
||||
D 1 1.00
|
||||
1.0960000 1.0000000
|
||||
F 1 1.00
|
||||
2.5440000 1.0000000
|
||||
|
|
|
|||
|
|
@ -345,8 +345,8 @@ program QuAcK
|
|||
|
||||
call cpu_time(start_CCD)
|
||||
! call ring_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
! call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
! call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
|
||||
call cpu_time(end_CCD)
|
||||
|
||||
t_CCD = end_CCD - start_CCD
|
||||
|
|
|
|||
59
src/QuAcK/form_ladder_r.f90
Normal file
59
src/QuAcK/form_ladder_r.f90
Normal file
|
|
@ -0,0 +1,59 @@
|
|||
subroutine form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
! Form residuals for ladder CCD
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: OOOO(nO,nO,nO,nO)
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: VVVV(nV,nV,nV,nV)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,k,l
|
||||
integer :: a,b,c,d
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: r2(nO,nO,nV,nV)
|
||||
|
||||
r2(:,:,:,:) = 0d0
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(k,l,a,b)*OOOO(i,j,k,l)
|
||||
end do
|
||||
end do
|
||||
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*VVVV(c,d,a,b)*t2(i,j,c,d)
|
||||
end do
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + 0.25d0*t2(i,j,c,d)*OOVV(k,l,c,d)*t2(k,l,a,b)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine form_ladder_r
|
||||
55
src/QuAcK/form_ring_r.f90
Normal file
55
src/QuAcK/form_ring_r.f90
Normal file
|
|
@ -0,0 +1,55 @@
|
|||
subroutine form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
|
||||
! Form residuals for ring CCD
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nO,nV
|
||||
double precision,intent(in) :: t2(nO,nO,nV,nV)
|
||||
double precision,intent(in) :: OVVO(nO,nV,nV,nO)
|
||||
double precision,intent(in) :: OOVV(nO,nO,nV,nV)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: i,j,k,l
|
||||
integer :: a,b,c,d
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: r2(nO,nO,nV,nV)
|
||||
|
||||
r2(:,:,:,:) = 0d0
|
||||
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
|
||||
do k=1,nO
|
||||
do c=1,nV
|
||||
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,a,k)*t2(k,j,c,b) + OVVO(j,c,b,k)*t2(i,k,a,c)
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
do c=1,nV
|
||||
do d=1,nV
|
||||
|
||||
r2(i,j,a,b) = r2(i,j,a,b) + t2(i,k,a,c)*OOVV(k,l,c,d)*t2(l,j,d,b)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine form_ring_r
|
||||
193
src/QuAcK/ladder_CCD.f90
Normal file
193
src/QuAcK/ladder_CCD.f90
Normal file
|
|
@ -0,0 +1,193 @@
|
|||
subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Ladder CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nEl
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nBas2
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
double precision :: ECCD,EcCCD
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVOV(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: X1(:,:,:,:)
|
||||
double precision,allocatable :: X2(:,:)
|
||||
double precision,allocatable :: X3(:,:)
|
||||
double precision,allocatable :: X4(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: u(:,:,:,:)
|
||||
double precision,allocatable :: v(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'**************************************'
|
||||
write(*,*)'| ladder-CCD calculation |'
|
||||
write(*,*)'**************************************'
|
||||
write(*,*)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*nEl
|
||||
nV = nBas2 - nO
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = seHF(1:nO)
|
||||
eV(:) = seHF(nO+1:nBas2)
|
||||
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2)
|
||||
OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO ,nO+1:nBas2)
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t2(nO,nO,nV,nV))
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
EcMP4 = 0d0
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV))
|
||||
allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| ladder-CCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:,:,:)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
|
||||
if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.))
|
||||
|
||||
! Dump results
|
||||
|
||||
ECCD = ERHF + EcCCD
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||
'|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
! Moller-Plesset energies
|
||||
|
||||
write(*,*)
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
|
||||
write(*,*)
|
||||
|
||||
end subroutine ladder_CCD
|
||||
|
|
@ -70,23 +70,25 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
|
|||
|
||||
Z(:,:) = M(:,:)
|
||||
call diagonalize_general_matrix(nOO+nVV,M(:,:),Omega(:),Z(:,:))
|
||||
! call diagonalize_matrix(nOO+nVV,Z(:,:),Omega(:))
|
||||
|
||||
write(*,*) 'pp-RPA excitation energies'
|
||||
call matout(nOO+nVV,1,Omega(:))
|
||||
write(*,*)
|
||||
! write(*,*) 'pp-RPA excitation energies'
|
||||
! call matout(nOO+nVV,1,Omega(:))
|
||||
! write(*,*)
|
||||
|
||||
! call matout(nOO+nVV,nOO+nVV,matmul(transpose(Z(:,:)),matmul(W(:,:),Z(:,:))))
|
||||
! write(*,*)
|
||||
|
||||
! Split the various quantities in p-p and h-h parts
|
||||
|
||||
Omega1(:) = Omega(nOO+1:nOO+nVV)
|
||||
Omega2(:) = Omega(1:nOO)
|
||||
call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
|
||||
! Omega1(:) = Omega(nOO+1:nOO+nVV)
|
||||
! Omega2(:) = Omega(1:nOO)
|
||||
|
||||
X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV)
|
||||
Y1(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV)
|
||||
X2(:,:) = Z(nOO+1:nOO+nVV, 1:nOO )
|
||||
Y2(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV)
|
||||
|
||||
if(minval(Omega1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!')
|
||||
if(maxval(Omega2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!')
|
||||
! X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV)
|
||||
! Y1(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV)
|
||||
! X2(:,:) = Z(nOO+1:nOO+nVV, 1:nOO )
|
||||
! Y2(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV)
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
|
||||
|
|
|
|||
|
|
@ -110,9 +110,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
|
|||
write(*,*)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (singlet) =',Ec_ppRPA(1)
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (triplet) =',Ec_ppRPA(2)
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy =',Ec_ppRPA(1) + Ec_ppRPA(2)
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + Ec_ppRPA(2)
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (triplet) =',3d0*Ec_ppRPA(2)
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
|
||||
write(*,'(2X,A40,F15.6)') 'pp-RPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
|
|
|
|||
193
src/QuAcK/ring_CCD.f90
Normal file
193
src/QuAcK/ring_CCD.f90
Normal file
|
|
@ -0,0 +1,193 @@
|
|||
subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
|
||||
|
||||
! Ring CCD module
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
double precision,intent(in) :: thresh
|
||||
|
||||
integer,intent(in) :: nBas,nEl
|
||||
double precision,intent(in) :: ENuc,ERHF
|
||||
double precision,intent(in) :: eHF(nBas)
|
||||
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nBas2
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
double precision :: EcMP2,EcMP3,EcMP4
|
||||
double precision :: ECCD,EcCCD
|
||||
double precision,allocatable :: seHF(:)
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
double precision,allocatable :: delta_OOVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: X1(:,:,:,:)
|
||||
double precision,allocatable :: X2(:,:)
|
||||
double precision,allocatable :: X3(:,:)
|
||||
double precision,allocatable :: X4(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: u(:,:,:,:)
|
||||
double precision,allocatable :: v(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: r2(:,:,:,:)
|
||||
double precision,allocatable :: t2(:,:,:,:)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'**************************************'
|
||||
write(*,*)'| ring-CCD calculation |'
|
||||
write(*,*)'**************************************'
|
||||
write(*,*)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas2 = 2*nBas
|
||||
|
||||
allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF)
|
||||
call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas2,nBas2,nBas2,nBas2))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas2,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Define occupied and virtual spaces
|
||||
|
||||
nO = 2*nEl
|
||||
nV = nBas2 - nO
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
allocate(delta_OOVV(nO,nO,nV,nV))
|
||||
|
||||
eO(:) = seHF(1:nO)
|
||||
eV(:) = seHF(nO+1:nBas2)
|
||||
|
||||
call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV)
|
||||
|
||||
deallocate(seHF)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVVV(nV,nV,nV,nV))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2, 1:nO )
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! MP2 guess amplitudes
|
||||
|
||||
allocate(t2(nO,nO,nV,nV))
|
||||
|
||||
t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
EcMP4 = 0d0
|
||||
|
||||
! Initialization
|
||||
|
||||
allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV))
|
||||
allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV))
|
||||
|
||||
Conv = 1d0
|
||||
nSCF = 0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,*)'| ring-CCD calculation |'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Compute residual
|
||||
|
||||
call form_ring_r(nO,nV,OVVO,OOVV,t2,r2)
|
||||
|
||||
r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:)
|
||||
|
||||
! Check convergence
|
||||
|
||||
Conv = maxval(abs(r2(:,:,:,:)))
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:)
|
||||
|
||||
! Compute correlation energy
|
||||
|
||||
EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.))
|
||||
|
||||
if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.))
|
||||
|
||||
! Dump results
|
||||
|
||||
ECCD = ERHF + EcCCD
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||
'|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|'
|
||||
|
||||
enddo
|
||||
write(*,*)'----------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
endif
|
||||
|
||||
! Moller-Plesset energies
|
||||
|
||||
write(*,*)
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3
|
||||
write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4
|
||||
write(*,*)
|
||||
|
||||
end subroutine ring_CCD
|
||||
88
src/QuAcK/sort_ppRPA.f90
Normal file
88
src/QuAcK/sort_ppRPA.f90
Normal file
|
|
@ -0,0 +1,88 @@
|
|||
subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
|
||||
|
||||
! Compute the metric matrix for pp-RPA
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nOO
|
||||
integer,intent(in) :: nVV
|
||||
double precision,intent(in) :: Omega(nOO+nVV)
|
||||
double precision,intent(in) :: Z(nOO+nVV,nOO+nVV)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: pq,ab,ij
|
||||
double precision,allocatable :: s(:,:)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: Omega1(nVV)
|
||||
double precision,intent(out) :: X1(nVV,nVV)
|
||||
double precision,intent(out) :: Y1(nOO,nVV)
|
||||
double precision,intent(out) :: Omega2(nOO)
|
||||
double precision,intent(out) :: X2(nVV,nOO)
|
||||
double precision,intent(out) :: Y2(nOO,nOO)
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(s(nOO+nVV,nOO+nVV))
|
||||
|
||||
! Initializatiom
|
||||
|
||||
Omega1(:) = 0d0
|
||||
Omega2(:) = 0d0
|
||||
|
||||
ab = 0
|
||||
ij = 0
|
||||
do pq=1,nOO+nVV
|
||||
|
||||
if(Omega(pq) > 0d0) then
|
||||
|
||||
ab = ab + 1
|
||||
Omega1(ab) = Omega(pq)
|
||||
X1(1:nVV,ab) = Z( 1: nVV,pq)
|
||||
Y1(1:nOO,ab) = Z(nVV+1:nOO+nVV,pq)
|
||||
|
||||
else
|
||||
|
||||
ij = ij + 1
|
||||
Omega2(ij) = Omega(pq)
|
||||
X2(1:nVV,ij) = Z( 1: nVV,pq)
|
||||
Y2(1:nOO,ij) = Z(nVV+1:nOO+nVV,pq)
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
|
||||
if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
|
||||
|
||||
write(*,*) 'pp-RPA positive excitation energies'
|
||||
call matout(nVV,1,Omega1(:))
|
||||
write(*,*)
|
||||
|
||||
write(*,*) 'pp-RPA negative excitation energies'
|
||||
call matout(nOO,1,Omega2(:))
|
||||
write(*,*)
|
||||
|
||||
! Check eigenvector signatures
|
||||
|
||||
! s( 1: nVV, 1: nVV) = matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1)
|
||||
! s(nVV+1:nOO+nVV,nVV+1:nOO+nVV) = matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2)
|
||||
|
||||
! write(*,*) 'Signatures pp'
|
||||
! do ab=1,nVV
|
||||
! write(*,'(I6,F10.6)') ab,s(ab,ab)
|
||||
! end do
|
||||
! write(*,*)
|
||||
|
||||
! write(*,*) 'Signatures hh'
|
||||
! do ij=1,nOO
|
||||
! write(*,'(I6,F10.6)') ij,s(ij,ij)
|
||||
! end do
|
||||
! write(*,*)
|
||||
|
||||
end subroutine sort_ppRPA
|
||||
Loading…
Reference in New Issue
Block a user