From 3acc640eee6bd78311cad35f6fd4389c8721a13a Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 17 Aug 2022 15:52:13 +0200 Subject: [PATCH] excitation vectors for ppRPA and ppBSE --- input/options | 2 +- src/GW/Bethe_Salpeter_pp.f90 | 4 + src/GW/static_screening_WC_pp.f90 | 8 +- src/GW/static_screening_WD_pp.f90 | 8 +- src/LR/print_transition_vectors_pp.f90 | 183 +++++++++++++++++++++++++ src/QuAcK/QuAcK.f90 | 2 +- src/RPA/ppRPA.f90 | 67 ++++----- 7 files changed, 233 insertions(+), 41 deletions(-) create mode 100644 src/LR/print_transition_vectors_pp.f90 diff --git a/input/options b/input/options index 5a29752..62bd5b2 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + T T F T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 3 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg diff --git a/src/GW/Bethe_Salpeter_pp.f90 b/src/GW/Bethe_Salpeter_pp.f90 index 32eac69..5e46305 100644 --- a/src/GW/Bethe_Salpeter_pp.f90 +++ b/src/GW/Bethe_Salpeter_pp.f90 @@ -97,6 +97,8 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) + call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + deallocate(Omega1,X1,Y1,Omega2,X2,Y2,WB,WC,WD) end if @@ -129,6 +131,8 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) + call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + deallocate(Omega1,X1,Y1,Omega2,X2,Y2,WB,WC,WD) end if diff --git a/src/GW/static_screening_WC_pp.f90 b/src/GW/static_screening_WC_pp.f90 index 2c41605..7422f76 100644 --- a/src/GW/static_screening_WC_pp.f90 +++ b/src/GW/static_screening_WC_pp.f90 @@ -44,11 +44,11 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ab = 0 do a=nO+1,nBas-nR - do b=a,nBas-nR + do b=a,nBas-nR ab = ab + 1 cd = 0 do c=nO+1,nBas-nR - do d=c,nBas-nR + do d=c,nBas-nR cd = cd + 1 chi = 0d0 @@ -75,11 +75,11 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ab = 0 do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do b=a+1,nBas-nR ab = ab + 1 cd = 0 do c=nO+1,nBas-nR - do d=c+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 chi = 0d0 diff --git a/src/GW/static_screening_WD_pp.f90 b/src/GW/static_screening_WD_pp.f90 index 4dfdf9e..ae4ffef 100644 --- a/src/GW/static_screening_WD_pp.f90 +++ b/src/GW/static_screening_WD_pp.f90 @@ -44,11 +44,11 @@ subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ij = 0 do i=nC+1,nO - do j=i,nO + do j=i,nO ij = ij + 1 kl = 0 do k=nC+1,nO - do l=k,nO + do l=k,nO kl = kl + 1 chi = 0d0 @@ -74,11 +74,11 @@ subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ij = 0 do i=nC+1,nO - do j=i+1,nO + do j=i+1,nO ij = ij + 1 kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=k+1,nO kl = kl + 1 chi = 0d0 diff --git a/src/LR/print_transition_vectors_pp.f90 b/src/LR/print_transition_vectors_pp.f90 new file mode 100644 index 0000000..620ed28 --- /dev/null +++ b/src/LR/print_transition_vectors_pp.f90 @@ -0,0 +1,183 @@ +subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + +! Print transition vectors for p-p linear response calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: spin_allowed + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision :: dipole_int(nBas,nBas,ncart) + double precision,intent(out) :: Omega1(nVV) + double precision,intent(out) :: X1(nVV,nVV) + double precision,intent(out) :: Y1(nOO,nVV) + double precision,intent(out) :: Omega2(nOO) + double precision,intent(out) :: X2(nVV,nOO) + double precision,intent(out) :: Y2(nOO,nOO) + +! Local variables + + integer :: a,b,c,d,ab,cd + integer :: i,j,k,l,ij,kl + integer :: maxOO + integer :: maxVV = 10 + double precision :: S2 + double precision,parameter :: thres_vec = 0.1d0 + double precision,allocatable :: osOO(:) + double precision,allocatable :: osVV(:) + +! Memory allocation + + maxOO = min(nOO,maxOO) + maxVV = min(nVV,maxVV) + + allocate(osOO(maxOO),osVV(maxVV)) + +! Compute oscillator strengths + + osOO(:) = 0d0 + osVV(:) = 0d0 + +! if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os) + +!-----------------------------------------------! +! Print details about excitations for pp sector ! +!-----------------------------------------------! + + do ab=1,maxVV + + ! values + + if(spin_allowed) then + S2 = 0d0 + else + S2 = 2d0 + end if + + print*,'-------------------------------------------------------------' + write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & + ' p-p excitation n. ',ab,': ',Omega1(ab)*HaToeV,' eV',' f = ',osVV(ab),' = ',S2 + print*,'-------------------------------------------------------------' + + if(spin_allowed) then + + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + if(abs(X1(cd,ab)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') c,' -- ',d,' = ',X1(cd,ab)/sqrt(2d0) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + if(abs(Y1(kl,ab)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') k,' -- ',l,' = ',Y1(kl,ab)/sqrt(2d0) + end do + end do + + else + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + if(abs(X1(cd,ab)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') c,' -- ',d,' = ',X1(cd,ab)/sqrt(2d0) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + if(abs(Y1(kl,ab)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') k,' -- ',l,' = ',Y1(kl,ab)/sqrt(2d0) + end do + end do + + end if + + write(*,*) + + end do + +! Thomas-Reiche-Kuhn sum rule + + write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(osVV(:)) + write(*,*) + +!-----------------------------------------------! +! Print details about excitations for hh sector ! +!-----------------------------------------------! + + do ij=1,maxOO + + ! values + + if(spin_allowed) then + S2 = 0d0 + else + S2 = 2d0 + end if + + print*,'-------------------------------------------------------------' + write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & + ' h-h excitation n. ',ij,': ',Omega2(ij)*HaToeV,' eV',' f = ',osOO(ij),' = ',S2 + print*,'-------------------------------------------------------------' + + if(spin_allowed) then + + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + if(abs(X2(cd,ij)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') c,' -- ',d,' = ',X2(cd,ij)/sqrt(2d0) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + if(abs(Y2(kl,ij)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') k,' -- ',l,' = ',Y2(kl,ij)/sqrt(2d0) + end do + end do + + else + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + if(abs(X2(cd,ij)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') c,' -- ',d,' = ',X2(cd,ij)/sqrt(2d0) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + if(abs(Y2(kl,ij)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') k,' -- ',l,' = ',Y2(kl,ij)/sqrt(2d0) + end do + end do + + end if + + write(*,*) + + end do + +! Thomas-Reiche-Kuhn sum rule + + write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(osOO(:)) + write(*,*) + +end subroutine print_transition_vectors_pp diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index b83fabb..f46fc9d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -845,7 +845,7 @@ program QuAcK else - call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) end if diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index 3e0b18a..2e6894f 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) +subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipole_int,e) ! Perform pp-RPA calculation @@ -20,19 +20,20 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) double precision,intent(in) :: ERHF double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables integer :: ispin integer :: nS - integer :: nOOs,nOOt - integer :: nVVs,nVVt - double precision,allocatable :: Omega1s(:),Omega1t(:) - double precision,allocatable :: X1s(:,:),X1t(:,:) - double precision,allocatable :: Y1s(:,:),Y1t(:,:) - double precision,allocatable :: Omega2s(:),Omega2t(:) - double precision,allocatable :: X2s(:,:),X2t(:,:) - double precision,allocatable :: Y2s(:,:),Y2t(:,:) + integer :: nOO + integer :: nVV + double precision,allocatable :: Omega1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: Omega2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) double precision :: Ec_ppRPA(nspin) double precision :: EcAC(nspin) @@ -54,31 +55,26 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) nS = nO*nV - nOOs = nO*(nO+1)/2 - nVVs = nV*(nV+1)/2 - - nOOt = nO*(nO-1)/2 - nVVt = nV*(nV-1)/2 - - ! Memory allocation - - allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & - Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) - - allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & - Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) - ! Singlet manifold if(singlet) then ispin = 1 - call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin)) + nOO = nO*(nO+1)/2 + nVV = nV*(nV+1)/2 - call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s) - call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s) + allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) + + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & + Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin)) + + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) + + call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + + deallocate(Omega1,X1,Y1,Omega2,X2,Y2) endif @@ -88,11 +84,20 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ispin = 2 - call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin)) + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 - call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t) - call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t) + allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) + + call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & + Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin)) + + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) + + call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) + + deallocate(Omega1,X1,Y1,Omega2,X2,Y2) endif