remove ERI copy in qsGW

This commit is contained in:
Pierre-Francois Loos 2019-10-30 10:31:05 +01:00
parent 72abe8a4f6
commit 1acf76b156
13 changed files with 99 additions and 90 deletions

View File

@ -1,14 +1,9 @@
1 2
1 3
S 3 1.00
18.7311370 0.03349460
2.8253937 0.23472695
0.6401217 0.81375733
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
S 1 1.00
0.1612778 1.0000000
2 2
S 3 1.00
18.7311370 0.03349460
2.8253937 0.23472695
0.6401217 0.81375733
S 1 1.00
0.1612778 1.0000000
0.2976000 1.0000000
P 1 1.00
1.2750000 1.0000000

View File

@ -1,5 +1,4 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
1 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.4
He 0.0 0.0 0.0

View File

@ -1,5 +1,5 @@
# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 5 2 1
64 0.000000001 T 5 1 1
# MP:
# CC: maxSCF thresh DIIS n_diis

View File

@ -1,14 +1,9 @@
1 2
1 3
S 3 1.00
18.7311370 0.03349460
2.8253937 0.23472695
0.6401217 0.81375733
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
S 1 1.00
0.1612778 1.0000000
2 2
S 3 1.00
18.7311370 0.03349460
2.8253937 0.23472695
0.6401217 0.81375733
S 1 1.00
0.1612778 1.0000000
0.2976000 1.0000000
P 1 1.00
1.2750000 1.0000000

View File

@ -9,8 +9,8 @@ else
FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O3
endif
LIBS = ~/Dropbox/quack/lib/*.a
#LIBS = -lblas -llapack
LIBS = ../../lib/*.a
LIBS += -lblas -llapack
SRCF90 = $(wildcard *.f90)

View File

@ -14,8 +14,8 @@ ifeq ($(PROFIL),1)
FC += -p -fno-inline
endif
LIBS = ~/Dropbox/quack/lib/*.a
#LIBS = -lblas -llapack
LIBS = ../../lib/*.a
LIBS += -lblas -llapack
SRCF90 = $(wildcard *.f90)

View File

@ -516,7 +516,7 @@ program QuAcK
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
call cpu_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW

View File

@ -16,7 +16,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
! Local variables
integer :: p,q
integer :: ab,cd,ij,kl
double precision :: trace_matrix
double precision,allocatable :: B(:,:)
double precision,allocatable :: C(:,:)
@ -64,11 +64,32 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO)
M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO))
! do p=1,nOO+nVV
! do q=1,nOO+nVV
! write(42,*) p,q,M(p,q)
! end do
! end do
open(unit=42,file='B.dat')
open(unit=43,file='C.dat')
open(unit=44,file='D.dat')
do ab=1,nVV
do ij=1,nOO
write(42,*) ab,ij,B(ab,ij)
end do
end do
do ab=1,nVV
do cd=1,nVV
write(43,*) ab,cd,C(ab,cd)
end do
end do
do ij=1,nOO
do kl=1,nOO
write(44,*) ij,kl,D(ij,kl)
end do
end do
close(42)
close(43)
close(44)
! print*, 'pp-RPA matrix'
! call matout(nOO+nVV,nOO+nVV,M(:,:))

View File

@ -1,5 +1,5 @@
subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
! Compute linear response
@ -32,6 +32,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO_basis(nBas,nBas,nBas,nBas)
! Local variables
@ -65,7 +66,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
double precision,allocatable :: SigCp(:,:)
double precision,allocatable :: SigCm(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: ERI_MO_basis(:,:,:,:)
double precision,allocatable :: error(:,:)
! Hello world
@ -76,6 +76,11 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
write(*,*)'************************************************'
write(*,*)
! Warning
write(*,*) '!! ERIs in MO basis will be overwritten in qsGW !!'
write(*,*)
! Stuff
nBasSq = nBas*nBas
@ -90,7 +95,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(COHSEX) write(*,*) 'COHSEX approximation activated!'
write(*,*)
! Switch off exchange for G0W0
! Switch off exchange for qsGW
dRPA = .true.
@ -98,9 +103,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
allocate(e(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), &
ERI_MO_basis(nBas,nBas,nBas,nBas),error(nBas,nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Initialization

View File

@ -83,39 +83,39 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2, &
EcRPA)
! call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
! call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:))
! Compute excitation densities for the T-matrix
! call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), &
! X1(:,:),Y1(:,:),rho1(:,:,:), &
! X2(:,:),Y2(:,:),rho2(:,:,:))
call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), &
X1(:,:),Y1(:,:),rho1(:,:,:), &
X2(:,:),Y2(:,:),rho2(:,:,:))
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
! call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
! Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
! SigT(:))
call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
SigT(:))
! Compute renormalization factor for T-matrix self-energy
! call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
! Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
! Z(:))
call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
Z(:))
!----------------------------------------------
! Solve the quasi-particle equation
!----------------------------------------------
! eG0T0(:) = seHF(:) + Z(:)*SigT(:)
eG0T0(:) = seHF(:) + Z(:)*SigT(:)
!----------------------------------------------
! Dump results
!----------------------------------------------
! call print_G0T0(nBas2,nO2,seHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA)
call print_G0T0(nBas2,nO2,seHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA)
end subroutine soG0T0

View File

@ -69,12 +69,12 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
call matout(nOO,1,Omega2(:))
write(*,*)
! write(*,*) 'X1.X1 - Y1.Y1'
! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
! write(*,*) 'X2.X2 - Y2.Y2'
! call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2))
! write(*,*) 'X1.X2 - Y1.Y2'
! call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2))
write(*,*) 'X1t.X1 - Y1t.Y1'
call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
write(*,*) 'X2t.X2 - Y2t.Y2'
call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2))
write(*,*) 'X1t.X2 - Y1t.Y2'
call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2))
end subroutine sort_ppRPA

View File

@ -13,7 +13,7 @@ FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O3
endif
LIBS = ../../lib/*.a
#LIBS = -lblas -llapack
LIBS += -lblas -llapack
SRCF90 = $(wildcard *.f90)

View File

@ -47,45 +47,40 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
! Local variables
integer :: lwork,info
double precision,allocatable :: work(:),WI(:),VL(:,:),tmp1(:,:),tmp2(:,:)
double precision,allocatable :: work(:),WI(:),VL(:,:)
! Memory allocation
lwork = 4*N
allocate(WI(N),VL(N,N),work(lwork),tmp1(N,N),tmp2(N,N))
tmp1 = A
allocate(WI(N),VL(N,N),work(lwork))
! tmp1 = A
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
if(info /= 0) then
print*,'Problem in diagonalize_general_matrix (dgeev)!!'
endif
call matout(N,1,WI)
! tmp2 = 0d0
! do i=1,N
! do j=1,N
! do k=1,N
! tmp2(i,j) = tmp2(i,j) + tmp1(i,k)*vr(k,j)
! end do
! end do
! end do
tmp2 = 0d0
do i=1,N
do j=1,N
do k=1,N
tmp2(i,j) = tmp2(i,j) + vl(k,i)*tmp1(k,j)
end do
end do
end do
! print*,'tmp2'
! call matout(N,N,tmp2)
print*,'tmp2'
call matout(N,N,tmp2)
! tmp1 = 0d0
! do i=1,N
! do j=1,N
! tmp1(i,j) = wr(j)*vr(i,j)
! end do
! end do
tmp1 = 0d0
do i=1,N
do j=1,N
tmp1(i,j) = wr(i)*vl(i,j)
end do
end do
print*,'tmp1'
call matout(N,N,tmp1)
print*,'coucou'
print*,maxval(tmp1-tmp2),minval(tmp1-tmp2)
! print*,'tmp1'
! call matout(N,N,tmp1)
end subroutine diagonalize_general_matrix