From 3433e3f99e1c572c2de1ade193f87b848440cdd7 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 6 Oct 2020 14:24:54 +0200 Subject: [PATCH] sf-BSE --- input/methods | 4 +- input/options | 4 +- src/CI/UCIS.f90 | 14 +-- src/LR/linear_response.f90 | 92 ++++++++++--------- src/LR/print_excitation.f90 | 2 +- .../print_unrestricted_transition_vectors.f90 | 2 +- src/LR/unrestricted_S2_expval.f90 | 4 +- src/LR/unrestricted_linear_response.f90 | 70 +++++++------- src/MBPT/G0W0.f90 | 7 +- src/MBPT/UG0W0.f90 | 14 +-- src/MBPT/evGW.f90 | 7 +- src/MBPT/evUGW.f90 | 7 +- src/MBPT/print_qsGW.f90 | 4 +- src/MBPT/qsGW.f90 | 2 +- src/MBPT/unrestricted_Bethe_Salpeter.f90 | 16 ++-- ...ed_Bethe_Salpeter_dynamic_perturbation.f90 | 4 +- src/QuAcK/QuAcK.f90 | 20 ++-- src/QuAcK/read_options.f90 | 3 +- 18 files changed, 137 insertions(+), 139 deletions(-) diff --git a/input/methods b/input/methods index f7bba22..d16be27 100644 --- a/input/methods +++ b/input/methods @@ -7,13 +7,13 @@ # drCCD rCCD lCCD pCCD F F F F # CIS* CIS(D) CID CISD - T F F F + F F F F # RPA* RPAx* ppRPA F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 5459e80..e075a66 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + T T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T F T F + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index e9596e5..19bace6 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -38,9 +38,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E integer :: nS_ab,nS_ba,nS_sf double precision,allocatable :: A_sf(:,:) - double precision,allocatable :: Z_sf(:,:) double precision,allocatable :: Omega_sf(:) - integer ,allocatable :: order(:) ! Hello world @@ -80,9 +78,10 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E endif call diagonalize_matrix(nS_sc,A_sc,Omega_sc) + A_sc(:,:) = transpose(A_sc) call print_excitation('UCIS ',5,nS_sc,Omega_sc) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & - cHF,S,Omega_sc,transpose(A_sc),transpose(A_sc)) + cHF,S,Omega_sc,A_sc,A_sc) if(dump_trans) then print*,'Spin-conserved CIS transition vectors' @@ -108,7 +107,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) nS_sf = nS_ab + nS_ba - allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf),Z_sf(nS_sf,nS_sf)) + allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf)) call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf) @@ -120,13 +119,10 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E endif call diagonalize_matrix(nS_sf,A_sf,Omega_sf) -! allocate(order(nS_sf)) -! call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf) -! call quick_sort(Omega_sf,order(:),nS_sf) - + A_sf(:,:) = transpose(A_sf) call print_excitation('UCIS ',6,nS_sf,Omega_sf) call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & - cHF,S,Omega_sf,transpose(A_sf),transpose(A_sf)) + cHF,S,Omega_sf,A_sf,A_sf) if(dump_trans) then print*,'Spin-flip CIS transition vectors' diff --git a/src/LR/linear_response.f90 b/src/LR/linear_response.f90 index 427b4b8..fa1ae3f 100644 --- a/src/LR/linear_response.f90 +++ b/src/LR/linear_response.f90 @@ -19,7 +19,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Local variables - integer :: ia double precision :: trace_matrix double precision,allocatable :: A(:,:) double precision,allocatable :: B(:,:) @@ -48,55 +47,62 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Tamm-Dancoff approximation - B = 0d0 - if(.not. TDA) then + if(TDA) then + + B(:,:) = 0d0 + XpY(:,:) = A(:,:) + XmY(:,:) = 0d0 + call diagonalize_matrix(nS,XpY,Omega) + XpY(:,:) = transpose(XpY(:,:)) + + else call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B) + ! Build A + B and A - B matrices + + ApB = A + B + AmB = A - B + + ! Diagonalize linear response matrix + + call diagonalize_matrix(nS,AmB,Omega) + + if(minval(Omega) < 0d0) & + call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') + +! do ia=1,nS +! if(Omega(ia) < 0d0) Omega(ia) = 0d0 +! end do + + call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq) + call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv) + + Z = matmul(AmBSq,matmul(ApB,AmBSq)) + + call diagonalize_matrix(nS,Z,Omega) + + if(minval(Omega) < 0d0) & + call print_warning('You may have instabilities in linear response: negative excitations!!') + + ! do ia=1,nS + ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 + ! end do + + Omega = sqrt(Omega) + + XpY = matmul(transpose(Z),AmBSq) + call DA(nS,1d0/sqrt(Omega),XpY) + + XmY = matmul(transpose(Z),AmBIv) + call DA(nS,1d0*sqrt(Omega),XmY) + end if -! Build A + B and A - B matrices + ! Compute the RPA correlation energy - ApB = A + B - AmB = A - B - -! Diagonalize linear response matrix - - call diagonalize_matrix(nS,AmB,Omega) - - if(minval(Omega) < 0d0) & - call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - -! do ia=1,nS -! if(Omega(ia) < 0d0) Omega(ia) = 0d0 -! end do - - call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq) - call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv) - - Z = matmul(AmBSq,matmul(ApB,AmBSq)) - - call diagonalize_matrix(nS,Z,Omega) - - if(minval(Omega) < 0d0) & - call print_warning('You may have instabilities in linear response: negative excitations!!') - -! do ia=1,nS -! if(Omega(ia) < 0d0) Omega(ia) = 0d0 -! end do - - Omega = sqrt(Omega) - - XpY = matmul(transpose(Z),AmBSq) - call DA(nS,1d0/sqrt(Omega),XpY) - - XmY = matmul(transpose(Z),AmBIv) - call DA(nS,1d0*sqrt(Omega),XmY) - -! Compute the RPA correlation energy - - EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) end subroutine linear_response diff --git a/src/LR/print_excitation.f90 b/src/LR/print_excitation.f90 index 26ee587..1720300 100644 --- a/src/LR/print_excitation.f90 +++ b/src/LR/print_excitation.f90 @@ -14,7 +14,7 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Local variables character*14 :: spin_manifold - integer,parameter :: maxS = 10 + integer,parameter :: maxS = 20 integer :: ia if(ispin == 1) spin_manifold = 'singlet' diff --git a/src/LR/print_unrestricted_transition_vectors.f90 b/src/LR/print_unrestricted_transition_vectors.f90 index 515248c..7816192 100644 --- a/src/LR/print_unrestricted_transition_vectors.f90 +++ b/src/LR/print_unrestricted_transition_vectors.f90 @@ -29,7 +29,7 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n ! Local variables integer :: ia,jb,j,b - integer :: maxS = 10 + integer :: maxS = 20 double precision,parameter :: thres_vec = 0.1d0 double precision,allocatable :: X(:) double precision,allocatable :: Y(:) diff --git a/src/LR/unrestricted_S2_expval.f90 b/src/LR/unrestricted_S2_expval.f90 index 6dc7018..c62bf2a 100644 --- a/src/LR/unrestricted_S2_expval.f90 +++ b/src/LR/unrestricted_S2_expval.f90 @@ -62,7 +62,7 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S !-------------------------! S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) - S2_gs = S2_exact + nO(2) - sum(OO(:,:)**2) + S2_gs = S2_exact + dble(nO(2)) - sum(OO(:,:)**2) !------------------------------------------! ! for spin-conserved-excited states ! @@ -152,7 +152,7 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S Yat(:,:) = transpose(Ya(:,:)) Ybt(:,:) = transpose(Yb(:,:)) - S2(m) = S2_gs & + S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0 & + trace_matrix(nV(1),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) & - trace_matrix(nO(2),matmul(Xb,matmul(VO,matmul(VOt,Xbt)))) & diff --git a/src/LR/unrestricted_linear_response.f90 b/src/LR/unrestricted_linear_response.f90 index 32b2d6f..2571f6a 100644 --- a/src/LR/unrestricted_linear_response.f90 +++ b/src/LR/unrestricted_linear_response.f90 @@ -33,7 +33,6 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Local variables - integer :: ia double precision,external :: trace_matrix double precision,allocatable :: A(:,:) double precision,allocatable :: B(:,:) @@ -65,8 +64,15 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, ! Tamm-Dancoff approximation - B = 0d0 - if(.not. TDA) then + if(TDA) then + + B(:,:) = 0d0 + XpY(:,:) = A(:,:) + XmY(:,:) = 0d0 + call diagonalize_matrix(nSt,XpY,Omega) + XpY(:,:) = transpose(XpY(:,:)) + + else call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & ERI_aaaa,ERI_aabb,ERI_bbbb,B) @@ -75,45 +81,45 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,B) - end if + ! Build A + B and A - B matrices -! Build A + B and A - B matrices + ApB = A + B + AmB = A - B - ApB = A + B - AmB = A - B + ! Diagonalize linear response matrix -! Diagonalize linear response matrix + call diagonalize_matrix(nSt,AmB,Omega) - call diagonalize_matrix(nSt,AmB,Omega) + if(minval(Omega) < 0d0) & + call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - if(minval(Omega) < 0d0) & - call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') + ! do ia=1,nSt + ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 + ! end do -! do ia=1,nSt -! if(Omega(ia) < 0d0) Omega(ia) = 0d0 -! end do - - call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq) - call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv) - - Z = matmul(AmBSq,matmul(ApB,AmBSq)) - - call diagonalize_matrix(nSt,Z,Omega) - - if(minval(Omega) < 0d0) & - call print_warning('You may have instabilities in linear response: negative excitations!!') + call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq) + call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv) -! do ia=1,nSt -! if(Omega(ia) < 0d0) Omega(ia) = 0d0 -! end do + Z = matmul(AmBSq,matmul(ApB,AmBSq)) + + call diagonalize_matrix(nSt,Z,Omega) - Omega = sqrt(Omega) + if(minval(Omega) < 0d0) & + call print_warning('You may have instabilities in linear response: negative excitations!!') + + ! do ia=1,nSt + ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 + ! end do - XpY = matmul(transpose(Z),AmBSq) - call DA(nSt,1d0/sqrt(Omega),XpY) + Omega = sqrt(Omega) + + XpY = matmul(transpose(Z),AmBSq) + call DA(nSt,1d0/sqrt(Omega),XpY) + + XmY = matmul(transpose(Z),AmBIv) + call DA(nSt,1d0*sqrt(Omega),XmY) - XmY = matmul(transpose(Z),AmBIv) - call DA(nSt,1d0*sqrt(Omega),XmY) + end if ! Compute the RPA correlation energy diff --git a/src/MBPT/G0W0.f90 b/src/MBPT/G0W0.f90 index 058a0c8..a69b76b 100644 --- a/src/MBPT/G0W0.f90 +++ b/src/MBPT/G0W0.f90 @@ -1,6 +1,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,dipole_int,PHF,cHF,eHF,eGW) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eGW) ! Perform G0W0 calculation @@ -29,12 +29,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nBas) ! Local variables diff --git a/src/MBPT/UG0W0.f90 b/src/MBPT/UG0W0.f90 index 843bf25..ad33756 100644 --- a/src/MBPT/UG0W0.f90 +++ b/src/MBPT/UG0W0.f90 @@ -1,5 +1,5 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb, & + linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW) ! Perform unrestricted G0W0 calculation @@ -33,15 +33,15 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev integer,intent(in) :: nS(nspin) double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF - double precision,intent(in) :: eHF(nBas,nspin) - double precision,intent(in) :: cHF(nBas,nBas,nspin) - double precision,intent(in) :: PHF(nBas,nBas,nspin) - double precision,intent(in) :: Hc(nBas,nBas,nspin) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: cHF(nBas,nBas,nspin) + double precision,intent(in) :: PHF(nBas,nBas,nspin) ! Local variables @@ -180,8 +180,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & - ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,eGW,EcBSE) + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/MBPT/evGW.f90 b/src/MBPT/evGW.f90 index 26061f8..d471d70 100644 --- a/src/MBPT/evGW.f90 +++ b/src/MBPT/evGW.f90 @@ -1,6 +1,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI, & - dipole_int,PHF,cHF,eHF,eG0W0) + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI, & + dipole_int,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -32,10 +32,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eG0W0(nBas) - double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) diff --git a/src/MBPT/evUGW.f90 b/src/MBPT/evUGW.f90 index 3323394..f322982 100644 --- a/src/MBPT/evUGW.f90 +++ b/src/MBPT/evUGW.f90 @@ -1,6 +1,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & - ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0) + ERHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -39,9 +39,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin) - double precision,intent(in) :: PHF(nBas,nBas,nspin) double precision,intent(in) :: eG0W0(nBas,nspin) - double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) @@ -255,7 +254,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE if(BSE) then call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & - ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eGW,eGW,EcBSE) + S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE) ! if(exchange_kernel) then diff --git a/src/MBPT/print_qsGW.f90 b/src/MBPT/print_qsGW.f90 index 9f3cf62..d57fc00 100644 --- a/src/MBPT/print_qsGW.f90 +++ b/src/MBPT/print_qsGW.f90 @@ -1,4 +1,4 @@ -subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW) +subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW) ! Print one-electron energies and other stuff for qsGW @@ -11,7 +11,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig integer,intent(in) :: nBas,nO,nSCF double precision,intent(in) :: ENuc,EcRPA,EcGM,Conv,thresh double precision,intent(in) :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas) double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) double precision,intent(in) :: Z(nBas),SigC(nBas,nBas) diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90 index 2bdfaf3..bb594ef 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/MBPT/qsGW.f90 @@ -223,7 +223,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Print results ! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW) ! Increment diff --git a/src/MBPT/unrestricted_Bethe_Salpeter.f90 b/src/MBPT/unrestricted_Bethe_Salpeter.f90 index ac8bf1c..22aeb6a 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter.f90 @@ -1,6 +1,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int_aa,dipole_int_bb,eW,eGW,EcBSE) + nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -24,6 +24,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: cW(nBas,nBas,nspin) double precision,intent(in) :: eW(nBas,nspin) double precision,intent(in) :: eGW(nBas,nspin) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) @@ -96,9 +98,9 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin), & OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) - call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) -! call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & -! OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) + call print_excitation('BSE@UGW ',5,nS_sc,OmBSE_sc) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + cW,S,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -136,7 +138,9 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin), & OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) - call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf) + call print_excitation('BSE@UGW ',6,nS_sf,OmBSE_sf) + call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & + cW,S,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) !------------------------------------------------- ! Compute the dynamical screening at the BSE level diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 index 03660ee..4bb3cd7 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 @@ -1,6 +1,6 @@ subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nS_sc,eGW, & - ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & - OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) ! Compute dynamical effects via perturbation theory for BSE diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 58de284..f2504bf 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -56,7 +56,6 @@ program QuAcK double precision,allocatable :: T(:,:) double precision,allocatable :: V(:,:) double precision,allocatable :: Hc(:,:) - double precision,allocatable :: H(:,:) double precision,allocatable :: X(:,:) double precision,allocatable :: dipole_int(:,:,:) double precision,allocatable :: dipole_int_aa(:,:,:) @@ -64,7 +63,6 @@ program QuAcK double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz - integer :: ispin integer :: bra1,bra2 integer :: ket1,ket2 double precision,allocatable :: ERI_MO_aaaa(:,:,:,:) @@ -234,7 +232,7 @@ program QuAcK ! Memory allocation for one- and two-electron integrals allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), & - S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), & + S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), & dipole_int(nBas,nBas,ncart)) ! Read integrals @@ -322,10 +320,6 @@ program QuAcK ! AO to MO integral transform for post-HF methods !------------------------------------------------------------------------ -! Compute Hartree Hamiltonian in the MO basis - - call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO,H) - call cpu_time(start_AOtoMO) write(*,*) @@ -811,12 +805,12 @@ program QuAcK if(unrestricted) then call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0) + linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0W0) end if @@ -839,14 +833,14 @@ program QuAcK call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & - EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & - PHF,cHF,eHF,eG0W0) + EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & + cHF,eHF,eG0W0) else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0W0) end if call cpu_time(end_evGW) @@ -1016,7 +1010,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,dipole_int,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int,eHF,eG0W0) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0 diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index a667a3c..1abed8a 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -71,8 +71,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4,answer5, & - answer6,answer7,answer8,answer9,answer10 + character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7 ! Open file with method specification