4
1
mirror of https://github.com/pfloos/quack synced 2024-12-23 12:56:38 +01:00

looking for a bug

This commit is contained in:
Pierre-Francois Loos 2020-01-14 22:56:20 +01:00
parent c2ac210e4c
commit 43a4ed34c2
10 changed files with 133 additions and 47 deletions

View File

@ -1,9 +1,26 @@
1 3 1 3
S 3 1.00 S 3 1.00
38.3600000 0.0238090 0.9910616896D+02 0.1543289673D+00
5.7700000 0.1548910 0.1805231239D+02 0.5353281423D+00
1.2400000 0.4699870 0.4885660238D+01 0.4446345422D+00
S 1 1.00 S 3 1.00
0.2976000 1.0000000 0.3780455879D+01 -0.9996722919D-01
P 1 1.00 0.8784966449D+00 0.3995128261D+00
1.2750000 1.0000000 0.2857143744D+00 0.7001154689D+00
P 3 1.00
0.3780455879D+01 0.1559162750D+00
0.8784966449D+00 0.6076837186D+00
0.2857143744D+00 0.3919573931D+00
2 3
S 3 1.00
0.9910616896D+02 0.1543289673D+00
0.1805231239D+02 0.5353281423D+00
0.4885660238D+01 0.4446345422D+00
S 3 1.00
0.3780455879D+01 -0.9996722919D-01
0.8784966449D+00 0.3995128261D+00
0.2857143744D+00 0.7001154689D+00
P 3 1.00
0.3780455879D+01 0.1559162750D+00
0.8784966449D+00 0.6076837186D+00
0.2857143744D+00 0.3919573931D+00

View File

@ -5,11 +5,11 @@
# CCD CCSD CCSD(T) ringCCD ladderCCD # CCD CCSD CCSD(T) ringCCD ladderCCD
F F F F F F F F F F
# CIS RPA RPAx ppRPA ADC # CIS RPA RPAx ppRPA ADC
F T T F F F T F F F
# GF2 GF3 # GF2 GF3
F F F F
# G0W0 evGW qsGW # G0W0 evGW qsGW
T T T F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F F F F
# MCMP2 # MCMP2

View File

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

View File

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

View File

@ -1,9 +1,26 @@
1 3 1 3
S 3 1.00 S 3 1.00
38.3600000 0.0238090 0.9910616896D+02 0.1543289673D+00
5.7700000 0.1548910 0.1805231239D+02 0.5353281423D+00
1.2400000 0.4699870 0.4885660238D+01 0.4446345422D+00
S 1 1.00 S 3 1.00
0.2976000 1.0000000 0.3780455879D+01 -0.9996722919D-01
P 1 1.00 0.8784966449D+00 0.3995128261D+00
1.2750000 1.0000000 0.2857143744D+00 0.7001154689D+00
P 3 1.00
0.3780455879D+01 0.1559162750D+00
0.8784966449D+00 0.6076837186D+00
0.2857143744D+00 0.3919573931D+00
2 3
S 3 1.00
0.9910616896D+02 0.1543289673D+00
0.1805231239D+02 0.5353281423D+00
0.4885660238D+01 0.4446345422D+00
S 3 1.00
0.3780455879D+01 -0.9996722919D-01
0.8784966449D+00 0.3995128261D+00
0.2857143744D+00 0.7001154689D+00
P 3 1.00
0.3780455879D+01 0.1559162750D+00
0.8784966449D+00 0.6076837186D+00
0.2857143744D+00 0.3919573931D+00

View File

@ -34,6 +34,7 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
integer :: ispin integer :: ispin
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM double precision :: EcGM
double precision,allocatable :: SigC(:) double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
@ -113,6 +114,15 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
@ -139,7 +149,16 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani
end if end if
call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho) nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if end if

View File

@ -31,6 +31,7 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
double precision :: rho double precision :: rho
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision :: EcAC(nspin)
! Hello world ! Hello world
@ -43,6 +44,7 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
! Initialization ! Initialization
EcRPA(:) = 0d0 EcRPA(:) = 0d0
EcAC(:) = 0d0
! Memory allocation ! Memory allocation
@ -91,7 +93,18 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu
write(*,*) write(*,*)
call ACFDT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, & call ACFDT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho) nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@RPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@RPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@RPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if end if

View File

@ -45,6 +45,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
double precision :: Conv double precision :: Conv
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM double precision :: EcGM
double precision :: alpha double precision :: alpha
double precision,allocatable :: error_diis(:,:) double precision,allocatable :: error_diis(:,:)
@ -218,6 +219,15 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
@ -244,7 +254,16 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
end if end if
call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho) nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if end if

View File

@ -51,19 +51,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
ApB = A + B ApB = A + B
AmB = A - B AmB = A - B
! print*,'A' ! Diagonalize linear response matrix
! call matout(nS,nS,A)
! print*,'B'
! call matout(nS,nS,B)
! print*,'A+B'
! call matout(nS,nS,ApB)
! print*,'A-B'
! call matout(nS,nS,AmB)
! Diagonalize TD-HF matrix
call diagonalize_matrix(nS,AmB,Omega) call diagonalize_matrix(nS,AmB,Omega)
@ -71,10 +59,8 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
call print_warning('You may have instabilities in linear response!!') call print_warning('You may have instabilities in linear response!!')
call ADAt(nS,AmB,sqrt(Omega),AmBSq) call ADAt(nS,AmB,sqrt(Omega),AmBSq)
Z = matmul(AmBSq,matmul(ApB,AmBSq))
! print*,'Z' Z = matmul(AmBSq,matmul(ApB,AmBSq))
! call matout(nS,nS,Z)
call diagonalize_matrix(nS,Z,Omega) call diagonalize_matrix(nS,Z,Omega)
@ -85,15 +71,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
XpY = matmul(transpose(Z),AmBSq) XpY = matmul(transpose(Z),AmBSq)
call DA(nS,1d0/sqrt(abs(Omega)),XpY) call DA(nS,1d0/sqrt(abs(Omega)),XpY)
call ADAt(nS,AmB,1d0/sqrt(Omega),AmBSq)
XmY = matmul(transpose(Z),AmBSq) XmY = matmul(transpose(Z),AmBSq)
call DA(nS,sqrt(abs(Omega)),XmY) call DA(nS,sqrt(abs(Omega)),XmY)
! print*,'X+Y'
! call matout(nS,nS,XpY)
! print*,'RPA excitations'
! call matout(nS,1,Omega)
! Compute the RPA correlation energy ! Compute the RPA correlation energy
EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A))

View File

@ -45,6 +45,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
double precision :: EqsGW double precision :: EqsGW
double precision :: EcRPA(nspin) double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin) double precision :: EcBSE(nspin)
double precision :: EcAC(nspin)
double precision :: EcGM double precision :: EcGM
double precision :: Conv double precision :: Conv
double precision :: rcond double precision :: rcond
@ -248,6 +249,15 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A40,F15.6)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
write(*,*) write(*,*)
write(*,*)'-------------------------------------------------------------------------------' write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1) write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
@ -274,7 +284,16 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW
end if end if
call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho) nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A40,F15.6)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
end if end if