From 43a4ed34c26035b6495d8c0dde1c0a94e9ee95d6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 14 Jan 2020 22:56:20 +0100 Subject: [PATCH] looking for a bug --- input/basis | 33 +++++++++++++++++++++++++-------- input/methods | 4 ++-- input/molecule | 5 +++-- input/options | 2 +- input/weight | 33 +++++++++++++++++++++++++-------- src/QuAcK/G0W0.f90 | 21 ++++++++++++++++++++- src/QuAcK/RPA.f90 | 15 ++++++++++++++- src/QuAcK/evGW.f90 | 21 ++++++++++++++++++++- src/QuAcK/linear_response.f90 | 25 +++---------------------- src/QuAcK/qsGW.f90 | 21 ++++++++++++++++++++- 10 files changed, 133 insertions(+), 47 deletions(-) diff --git a/input/basis b/input/basis index b9ca7b5..b246175 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,26 @@ -1 3 +1 3 S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 -S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 + 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 +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 diff --git a/input/methods b/input/methods index 4d309d6..3c25e69 100644 --- a/input/methods +++ b/input/methods @@ -5,11 +5,11 @@ # CCD CCSD CCSD(T) ringCCD ladderCCD F F F F F # CIS RPA RPAx ppRPA ADC - F T T F F + F T F F F # GF2 GF3 F F # G0W0 evGW qsGW - T T T + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index c78e87e..6905a9c 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 2 7 7 0 0 # Znuc x y z - He 0.0 0.0 0.0 + N 0. 0. 0. + N 0. 0. 2. diff --git a/input/options b/input/options index 07e4f0b..ea27add 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 1 1 + 64 0.0000001 T 5 2 1 # MP: # CC: maxSCF thresh DIIS n_diis diff --git a/input/weight b/input/weight index b9ca7b5..b246175 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,26 @@ -1 3 +1 3 S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 -S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 + 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 +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 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index b9d4a4e..cf9a4e4 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -34,6 +34,7 @@ subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_mani integer :: ispin double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) double precision :: EcGM double precision,allocatable :: SigC(:) 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, & 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(*,'(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 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 diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90 index bf9bac9..c7f256d 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/RPA.f90 @@ -31,6 +31,7 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu double precision :: rho double precision :: EcRPA(nspin) + double precision :: EcAC(nspin) ! Hello world @@ -43,6 +44,7 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu ! Initialization EcRPA(:) = 0d0 + EcAC(:) = 0d0 ! Memory allocation @@ -91,7 +93,18 @@ subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENu write(*,*) 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 diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 91c39fa..0648467 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -45,6 +45,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW double precision :: Conv double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) double precision :: EcGM double precision :: alpha 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, & 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(*,'(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 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 diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90 index 4a9e66c..a3bb787 100644 --- a/src/QuAcK/linear_response.f90 +++ b/src/QuAcK/linear_response.f90 @@ -51,19 +51,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r ApB = A + B AmB = A - B -! print*,'A' -! 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 +! Diagonalize linear response matrix 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 ADAt(nS,AmB,sqrt(Omega),AmBSq) - Z = matmul(AmBSq,matmul(ApB,AmBSq)) -! print*,'Z' -! call matout(nS,nS,Z) + Z = matmul(AmBSq,matmul(ApB,AmBSq)) 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) call DA(nS,1d0/sqrt(abs(Omega)),XpY) + call ADAt(nS,AmB,1d0/sqrt(Omega),AmBSq) XmY = matmul(transpose(Z),AmBSq) 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 EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 375039d..74e5e93 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -45,6 +45,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW double precision :: EqsGW double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) double precision :: EcGM double precision :: Conv 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, & 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(*,'(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 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