From b8bf488a9a1e0e11fffe26673f23e438c8307cf2 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 22 Sep 2020 23:08:47 +0200 Subject: [PATCH] spin conserved and spin flip --- input/options | 4 +- src/QuAcK/AOtoMO_integral_transform.f90 | 16 ++-- src/QuAcK/QuAcK.f90 | 94 +++++++++++++++-------- src/QuAcK/UG0W0.f90 | 53 ++++++------- src/QuAcK/qsGW.f90 | 2 +- src/QuAcK/read_options.f90 | 26 ++++--- src/QuAcK/unrestricted_Bethe_Salpeter.f90 | 60 ++++++++------- 7 files changed, 143 insertions(+), 112 deletions(-) diff --git a/input/options b/input/options index a5abf91..7be6358 100644 --- a/input/options +++ b/input/options @@ -4,8 +4,8 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 -# spin: singlet triplet TDA - T T F +# spin: singlet triplet spin_conserved spinf_flip TDA + T T T F F # 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 diff --git a/src/QuAcK/AOtoMO_integral_transform.f90 b/src/QuAcK/AOtoMO_integral_transform.f90 index 23e4198..dabeee8 100644 --- a/src/QuAcK/AOtoMO_integral_transform.f90 +++ b/src/QuAcK/AOtoMO_integral_transform.f90 @@ -1,15 +1,15 @@ -subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis) +subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI_MO_basis) ! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm -! bra and ket are the spin of (bra|ket) +! bra and ket are the spin of (bra1 bra2|ket1 ket2) implicit none include 'parameters.h' ! Input variables - integer,intent(in) :: bra - integer,intent(in) :: ket + integer,intent(in) :: bra1,bra2 + integer,intent(in) :: ket1,ket2 integer,intent(in) :: nBas double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin) @@ -35,7 +35,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis) do la=1,nBas do nu=1,nBas do mu=1,nBas - scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket) + scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket2) enddo enddo enddo @@ -49,7 +49,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis) do nu=1,nBas do i=1,nBas do mu=1,nBas - ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra)*scr(mu,nu,la,l) + ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra1)*scr(mu,nu,la,l) enddo enddo enddo @@ -63,7 +63,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis) do la=1,nBas do nu=1,nBas do i=1,nBas - scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra) + scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra2) enddo enddo enddo @@ -77,7 +77,7 @@ subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis) do j=1,nBas do i=1,nBas do nu=1,nBas - ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket)*scr(i,nu,k,l) + ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l) enddo ! print*,i,k,j,l,ERI_MO_basis(i,j,k,l) enddo diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 0bb56a2..6f62482 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -55,11 +55,12 @@ program QuAcK double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:) double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) - integer :: bra - integer :: ket + integer :: bra1,bra2 + integer :: ket1,ket2 double precision,allocatable :: ERI_MO_aaaa(:,:,:,:) double precision,allocatable :: ERI_MO_aabb(:,:,:,:) double precision,allocatable :: ERI_MO_bbbb(:,:,:,:) + double precision,allocatable :: ERI_MO_abab(:,:,:,:) double precision,allocatable :: ERI_ERF_AO(:,:,:,:) double precision,allocatable :: ERI_ERF_MO(:,:,:,:) double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:) @@ -101,8 +102,10 @@ program QuAcK double precision :: thresh_CC logical :: DIIS_CC - logical :: singlet_manifold - logical :: triplet_manifold + logical :: singlet + logical :: triplet + logical :: spin_conserved + logical :: spin_flip logical :: TDA integer :: maxSCF_GF,n_diis_GF,renormGF @@ -156,7 +159,7 @@ program QuAcK call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - singlet_manifold,triplet_manifold,TDA, & + singlet,triplet,spin_conserved,spin_flip,TDA, & maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & COHSEX,SOSEX,TDA_W,G0W,GW0, & @@ -334,21 +337,42 @@ program QuAcK ! 4-index transform for (aa|aa) block - bra = 1 - ket = 1 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aaaa) + bra1 = 1 + bra2 = 1 + ket1 = 1 + ket2 = 1 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aaaa) ! 4-index transform for (aa|bb) block - bra = 1 - ket = 2 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aabb) + bra1 = 1 + bra2 = 1 + ket1 = 2 + ket2 = 2 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_aabb) ! 4-index transform for (bb|bb) block - bra = 2 - ket = 2 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bbbb) + bra1 = 2 + bra2 = 2 + ket1 = 2 + ket2 = 2 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_bbbb) + + if(spin_flip) then + + allocate(ERI_MO_abab(nBas,nBas,nBas,nBas)) + + ! 4-index transform for (ab|ab) block + + bra1 = 1 + bra2 = 2 + ket1 = 1 + ket2 = 2 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO_abab) + + end if + else @@ -358,9 +382,11 @@ program QuAcK ! 4-index transform - bra = 1 - ket = 1 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO) + bra1 = 1 + bra2 = 1 + ket1 = 1 + ket2 = 1 + call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO) end if @@ -560,7 +586,7 @@ program QuAcK if(doCIS) then call cpu_time(start_CIS) - call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF) + call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF) call cpu_time(end_CIS) t_CIS = end_CIS - start_CIS @@ -576,7 +602,7 @@ program QuAcK if(doCID) then call cpu_time(start_CID) -! call CID(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF) +! call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF) call cpu_time(end_CID) t_CID = end_CID - start_CID @@ -592,7 +618,7 @@ program QuAcK if(doCISD) then call cpu_time(start_CISD) - call CISD(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ERI_MO,eHF) + call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,eHF) call cpu_time(end_CISD) t_CISD = end_CISD - start_CISD @@ -608,7 +634,7 @@ program QuAcK if(doRPA) then call cpu_time(start_RPA) - call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, & + call RPA(doACFDT,exchange_kernel,singlet,triplet,0d0, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPA) @@ -625,7 +651,7 @@ program QuAcK if(doRPAx) then call cpu_time(start_RPAx) - call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, & + call RPAx(doACFDT,exchange_kernel,singlet,triplet,0d0, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPAx) @@ -642,7 +668,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_ppRPA) - call ppRPA(singlet_manifold,triplet_manifold, & + call ppRPA(singlet,triplet, & nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_ppRPA) @@ -659,7 +685,7 @@ program QuAcK ! if(doADC) then ! call cpu_time(start_ADC) -! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & +! call ADC(singlet,triplet,maxSCF_GF,thresh_GF,n_diis_GF, & ! nBas,nC,nO,nV,nR,eHF,ERI_MO) ! call cpu_time(end_ADC) @@ -676,7 +702,7 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGF, & + call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) @@ -694,7 +720,7 @@ program QuAcK call cpu_time(start_GF2) call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, & - singlet_manifold,triplet_manifold,linGF, & + singlet,triplet,linGF, & eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) @@ -748,12 +774,12 @@ program QuAcK if(unrestricted) then call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & - singlet_manifold,triplet_manifold,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, & + 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,PHF,cHF,eHF,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & + dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0) end if @@ -774,7 +800,7 @@ program QuAcK call cpu_time(start_evGW) 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_manifold,triplet_manifold,eta_GW, & + BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) @@ -792,7 +818,7 @@ program QuAcK call cpu_time(start_qsGW) call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & - BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & + BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) call cpu_time(end_qsGW) @@ -812,7 +838,7 @@ program QuAcK call cpu_time(start_G0T0) call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & + dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_G0T0) @@ -830,7 +856,7 @@ program QuAcK call cpu_time(start_evGT) call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & + BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_evGT) @@ -947,7 +973,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & + dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -961,7 +987,7 @@ program QuAcK call cpu_time(start_G0T0) call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + singlet,triplet,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) call cpu_time(end_G0T0) diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index 0e0b38e..07de370 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -45,21 +45,13 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev logical :: print_W = .true. integer :: ispin - integer :: iblock - integer :: bra - integer :: ket - integer :: nSa - integer :: nSb - integer :: nSt - double precision :: EcRPA - double precision :: EcBSE + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) - double precision,allocatable :: Omega(:) - double precision,allocatable :: XpY(:,:) - double precision,allocatable :: XmY(:,:) - double precision,allocatable :: rho(:,:,:,:) + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Omega_sc(:),XpY_sc(:,:),XmY_sc(:,:),rho_sc(:,:,:,:) double precision,allocatable :: eGWlin(:,:) @@ -103,12 +95,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Memory allocation - nSa = nS(1) - nSb = nS(2) - nSt = nSa + nSb + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb - allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega(nSt),XpY(nSt,nSt),XmY(nSt,nSt), & - rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin)) + allocate(SigC(nBas,nspin),Z(nBas,nspin),Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc), & + rho_sc(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin)) !-------------------! ! Compute screening ! @@ -118,28 +110,28 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ispin = 1 - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, & - eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY) + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) - if(print_W) call print_excitation('RPA@UHF',5,nSt,Omega) + if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc) !----------------------! ! Excitation densities ! !----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_sc,rho_sc) !---------------------! ! Compute self-energy ! !---------------------! - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,SigC) + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,SigC) !--------------------------------! ! Compute renormalization factor ! !--------------------------------! - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,Z) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,eHF,Omega_sc,rho_sc,Z) !-----------------------------------! ! Solve the quasi-particle equation ! @@ -171,23 +163,26 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Compute the RPA correlation energy - call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY) + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(ispin) + write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(ispin) write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Free memory + + deallocate(Omega_sc,XpY_sc,XmY_sc,rho_sc) + ! Perform BSE calculation if(BSE) then call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb, & - eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcRPA,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 8191c4b..dad60a5 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -149,7 +149,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! AO to MO transformation of two-electron integrals - call AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis) + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO_basis,ERI_MO_basis) ! Compute linear response diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index e5e8e14..c9e72f3 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,6 +1,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - singlet_manifold,triplet_manifold,TDA, & + singlet,triplet,spin_conserved,spin_flip,TDA, & maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & COHSEX,SOSEX,TDA_W,G0W,GW0, & @@ -26,8 +26,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: DIIS_CC integer,intent(out) :: n_diis_CC - logical,intent(out) :: singlet_manifold - logical,intent(out) :: triplet_manifold + logical,intent(out) :: singlet + logical,intent(out) :: triplet + logical,intent(out) :: spin_conserved + logical,intent(out) :: spin_flip logical,intent(out) :: TDA integer,intent(out) :: maxSCF_GF @@ -113,16 +115,20 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t ! Read excited state options - singlet_manifold = .false. - triplet_manifold = .false. - TDA = .false. + singlet = .false. + triplet = .false. + spin_conserved = .false. + spin_flip = .false. + TDA = .false. read(1,*) - read(1,*) answer1,answer2,answer3 + read(1,*) answer1,answer2,answer3,answer4,answer5 - if(answer1 == 'T') singlet_manifold = .true. - if(answer2 == 'T') triplet_manifold = .true. - if(answer3 == 'T') TDA = .true. + if(answer1 == 'T') singlet = .true. + if(answer2 == 'T') triplet = .true. + if(answer3 == 'T') spin_conserved = .true. + if(answer4 == 'T') spin_flip = .true. + if(answer5 == 'T') TDA = .true. ! Read Green function options diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 index a75c9a7..ef85215 100644 --- a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 +++ b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 @@ -1,6 +1,5 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & - nBas,nC,nO,nV,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb, & - eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) + nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,eW,eGW,EcRPA,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -23,36 +22,32 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, integer,intent(in) :: nO(nspin) integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) - integer,intent(in) :: nSa - integer,intent(in) :: nSb - integer,intent(in) :: nSt + integer,intent(in) :: nS(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) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) - double precision :: OmRPA(nSt) - double precision :: XpY_RPA(nSt,nSt) - double precision :: XmY_RPA(nSt,nSt) - double precision :: rho_RPA(nBas,nBas,nSt,nspin) ! Local variables integer :: ispin integer :: isp_W - double precision,allocatable :: OmBSE(:) - double precision,allocatable :: XpY_BSE(:,:) - double precision,allocatable :: XmY_BSE(:,:) + integer :: nS_aa,nS_bb,nS_sc + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: OmRPA_sc(:) + double precision,allocatable :: XpY_RPA_sc(:,:) + double precision,allocatable :: XmY_RPA_sc(:,:) + double precision,allocatable :: rho_RPA_sc(:,:,:,:) + double precision,allocatable :: OmBSE_sc(:) + double precision,allocatable :: XpY_BSE_sc(:,:) + double precision,allocatable :: XmY_BSE_sc(:,:) ! Output variables - double precision,intent(out) :: EcRPA - double precision,intent(out) :: EcBSE - -! Memory allocation - - allocate(OmBSE(nSt),XpY_BSE(nSt,nSt),XmY_BSE(nSt,nSt)) + double precision,intent(out) :: EcRPA(nspin) + double precision,intent(out) :: EcBSE(nspin) !----------------------------! ! Spin-conserved excitations ! @@ -62,23 +57,32 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, ispin = 1 isp_W = 1 - EcBSE = 0d0 + EcBSE(ispin) = 0d0 - ! Compute spin-conserved RPA screening + ! Memory allocation - call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, & - eW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(OmRPA_sc(nS_sc),XpY_RPA_sc(nS_sc,nS_sc),XmY_RPA_sc(nS_sc,nS_sc),rho_RPA_sc(nBas,nBas,nS_sc,nspin)) + allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc)) - call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + ! Compute spin-conserved RPA screening - ! Compute BSE excitation energies + call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & + eW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA_sc,EcRPA(ispin),OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc) - OmBSE(:) = OmRPA(:) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA_sc,rho_RPA_sc) - call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, & - eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA,EcBSE,OmBSE,XpY_BSE,XmY_BSE) + ! Compute spin-conserved BSE excitation energies - call print_excitation('BSE@UG0W0',5,nSt,OmBSE) + OmBSE_sc(:) = OmRPA_sc(:) + + call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA_sc,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc) + + call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc) !------------------------------------------------- ! Compute the dynamical screening at the BSE level