From 1acf76b15650c26f3737155fef99a43b702e9b93 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 30 Oct 2019 10:31:05 +0100 Subject: [PATCH] remove ERI copy in qsGW --- input/basis | 19 +++++-------- input/molecule | 5 ++-- input/options | 2 +- input/weight | 19 +++++-------- src/IntPak/Makefile | 4 +-- src/QuAcK/Makefile | 4 +-- src/QuAcK/QuAcK.f90 | 2 +- src/QuAcK/linear_response_pp.f90 | 33 ++++++++++++++++++---- src/QuAcK/qsGW.f90 | 14 ++++++---- src/QuAcK/soG0T0.f90 | 26 +++++++++--------- src/QuAcK/sort_ppRPA.f90 | 12 ++++---- src/utils/Makefile | 2 +- src/utils/wrap_lapack.f90 | 47 ++++++++++++++------------------ 13 files changed, 99 insertions(+), 90 deletions(-) diff --git a/input/basis b/input/basis index a18b478..b9ca7b5 100644 --- a/input/basis +++ b/input/basis @@ -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 diff --git a/input/molecule b/input/molecule index 81c624a..c78e87e 100644 --- a/input/molecule +++ b/input/molecule @@ -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 diff --git a/input/options b/input/options index d8f69d6..6786f5f 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/input/weight b/input/weight index a18b478..b9ca7b5 100644 --- a/input/weight +++ b/input/weight @@ -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 diff --git a/src/IntPak/Makefile b/src/IntPak/Makefile index 6b0ec1d..ee03122 100644 --- a/src/IntPak/Makefile +++ b/src/IntPak/Makefile @@ -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) diff --git a/src/QuAcK/Makefile b/src/QuAcK/Makefile index 5183c73..8ba0980 100644 --- a/src/QuAcK/Makefile +++ b/src/QuAcK/Makefile @@ -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) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 116bb21..c4b6976 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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 diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index a5533f8..ca8a956 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -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(:,:)) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 86e88a4..b12c070 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -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 diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index 9997967..01a1602 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -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 diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index 4308944..c26ea66 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/src/QuAcK/sort_ppRPA.f90 @@ -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 diff --git a/src/utils/Makefile b/src/utils/Makefile index d30e9a9..55227cc 100644 --- a/src/utils/Makefile +++ b/src/utils/Makefile @@ -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) diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 1f11f3f..724e248 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -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