diff --git a/examples/basis.Ne.cc-pvdz b/examples/basis.Ne.cc-pvdz index f19a2d0..ae7d362 100644 --- a/examples/basis.Ne.cc-pvdz +++ b/examples/basis.Ne.cc-pvdz @@ -1,30 +1,29 @@ 1 6 S 8 - 1 17880.0000000 0.0007380 - 2 2683.0000000 0.0056770 - 3 611.5000000 0.0288830 - 4 173.5000000 0.1085400 - 5 56.6400000 0.2909070 - 6 20.4200000 0.4483240 - 7 7.8100000 0.2580260 - 8 1.6530000 0.0150630 + 1 17880.0000000 0.0007380 + 2 2683.0000000 0.0056770 + 3 611.5000000 0.0288830 + 4 173.5000000 0.1085400 + 5 56.6400000 0.2909070 + 6 20.4200000 0.4483240 + 7 7.8100000 0.2580260 + 8 1.6530000 0.0150630 S 8 - 1 17880.0000000 -0.0001720 - 2 2683.0000000 -0.0013570 - 3 611.5000000 -0.0067370 - 4 173.5000000 -0.0276630 - 5 56.6400000 -0.0762080 - 6 20.4200000 -0.1752270 - 7 7.8100000 -0.1070380 - 8 1.6530000 0.5670500 + 1 17880.0000000 -0.0001720 + 2 2683.0000000 -0.0013570 + 3 611.5000000 -0.0067370 + 4 173.5000000 -0.0276630 + 5 56.6400000 -0.0762080 + 6 20.4200000 -0.1752270 + 7 7.8100000 -0.1070380 + 8 1.6530000 0.5670500 S 1 - 1 0.4869000 1.0000000 + 1 0.4869000 1.0000000 P 3 - 1 28.3900000 0.0460870 - 2 6.2700000 0.2401810 - 3 1.6950000 0.5087440 + 1 28.3900000 0.0460870 + 2 6.2700000 0.2401810 + 3 1.6950000 0.5087440 P 1 - 1 0.4317000 1.0000000 + 1 0.4317000 1.0000000 D 1 1 2.2020000 1.0000000 - diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 81c624a..7225285 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd 2 1 1 0 0 # Znuc x y z - H 0. 0. 0. - H 0. 0. 1.4 + H 0. 0. -0.7 + H 0. 0. 0.7 diff --git a/input/basis b/input/basis index fb05e68..79a747a 100644 --- a/input/basis +++ b/input/basis @@ -7,12 +7,3 @@ S 1 1 0.1220000 1.0000000 P 1 1 0.7270000 1.0000000 -2 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -P 1 - 1 0.7270000 1.0000000 diff --git a/input/methods b/input/methods index ba17592..b1f2358 100644 --- a/input/methods +++ b/input/methods @@ -1,20 +1,21 @@ # RHF UHF MOM F T F -# MP2 MP3 MP2-F12 +# MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD F F F F -# CIS CID CISD +# CIS* CID CISD F F F -# RPA RPAx ppRPA +# RPA* RPAx* ppRPA F F F # G0F2 evGF2 G0F3 evGF3 F F F F -# G0W0 evGW qsGW +# G0W0* evGW* qsGW T F F # G0T0 evGT qsGT F F F # MCMP2 F +# * unrestricted version available diff --git a/input/molecule b/input/molecule index 81c624a..e2a4fd3 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 2 4 3 0 0 # Znuc x y z - H 0. 0. 0. - H 0. 0. 1.4 + C 0. 0. -0.16245872 + H 0. 0. 1.93436816 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..7a4f218 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + C 0.0000000000 0.0000000000 -0.0859694585 + H 0.0000000000 0.0000000000 1.0236236215 diff --git a/input/options b/input/options index daaa7f4..d9b23c7 100644 --- a/input/options +++ b/input/options @@ -1,11 +1,11 @@ -# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type - 64 0.0000001 T 5 1 1 +# HF: maxSCF thresh DIIS n_diis guess_type ortho_type + 64 0.00001 T 5 1 1 # MP: # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 -# spin: singlet triplet TDA - T T F +# spin: singlet triplet spin_conserved spin_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 @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T F F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/ACFDT.f90 b/src/QuAcK/ACFDT.f90 index 3e5b6a5..cd8843b 100644 --- a/src/QuAcK/ACFDT.f90 +++ b/src/QuAcK/ACFDT.f90 @@ -1,5 +1,4 @@ -subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eW,e,Omega,XpY,XmY,rho,EcAC) +subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,e,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem @@ -15,8 +14,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: BSE - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS @@ -24,11 +23,6 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision :: Omega(nS,nspin) - double precision :: XpY(nS,nS,nspin) - double precision :: XmY(nS,nS,nspin) - double precision :: rho(nBas,nBas,nS,nspin) - ! Local variables integer :: ispin @@ -37,6 +31,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl double precision :: lambda double precision,allocatable :: Ec(:,:) + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) + ! Output variables double precision,intent(out) :: EcAC(nspin) @@ -44,6 +48,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl ! Memory allocation allocate(Ec(nAC,nspin)) + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Antisymmetrized kernel version @@ -58,12 +64,20 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl EcAC(:) = 0d0 Ec(:,:) = 0d0 +! Compute (singlet) RPA screening + + isp_W = 1 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + ! Singlet manifold - if(singlet_manifold) then + if(singlet) then ispin = 1 - isp_W = 1 write(*,*) '--------------' write(*,*) 'Singlet states' @@ -80,17 +94,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(doXBS) then - call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & + rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -98,6 +112,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) write(*,*) '-----------------------------------------------------------------------------------' @@ -107,7 +123,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl ! Triplet manifold - if(triplet_manifold) then + if(triplet) then ispin = 2 isp_W = 1 @@ -127,17 +143,16 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(doXBS) then - call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & + rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) + call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) @@ -145,6 +160,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) + if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) + write(*,*) '-----------------------------------------------------------------------------------' write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) write(*,*) '-----------------------------------------------------------------------------------' 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/BSE2.f90 b/src/QuAcK/BSE2.f90 index cff9aef..123cfea 100644 --- a/src/QuAcK/BSE2.f90 +++ b/src/QuAcK/BSE2.f90 @@ -1,5 +1,4 @@ -subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & - eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,EcBSE) +subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -12,8 +11,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -25,6 +24,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -46,15 +46,15 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & ! Singlet manifold !------------------- - if(singlet_manifold) then + if(singlet) then ispin = 1 EcBSE(ispin) = 0d0 ! Compute BSE2 excitation energies - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & - rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI, & + OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -63,12 +63,12 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & if(evDyn) then - call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & + OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) else - call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & + OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if @@ -80,7 +80,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & ! Triplet manifold !------------------- - if(triplet_manifold) then + if(triplet) then ispin = 2 EcBSE(ispin) = 0d0 @@ -88,7 +88,7 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & ! Compute BSE2 excitation energies call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & - rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -97,12 +97,12 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & if(evDyn) then - call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & + OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) else - call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & + OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if diff --git a/src/QuAcK/BSE2_dynamic_perturbation.f90 b/src/QuAcK/BSE2_dynamic_perturbation.f90 index 30883d8..94c7c59 100644 --- a/src/QuAcK/BSE2_dynamic_perturbation.f90 +++ b/src/QuAcK/BSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,OmBSE,XpY,XmY) +subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -18,6 +18,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF, integer,intent(in) :: nS double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) double precision,intent(in) :: OmBSE(nS) @@ -56,7 +57,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF, ! Print main components of transition vectors - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) gapGF = eGF(nO+1) - eGF(nO) diff --git a/src/QuAcK/BSE2_dynamic_perturbation_iterative.f90 b/src/QuAcK/BSE2_dynamic_perturbation_iterative.f90 index 26d5905..eedfad1 100644 --- a/src/QuAcK/BSE2_dynamic_perturbation_iterative.f90 +++ b/src/QuAcK/BSE2_dynamic_perturbation_iterative.f90 @@ -1,4 +1,5 @@ -subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,OmBSE,XpY,XmY) +subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & + eHF,eGF,OmBSE,XpY,XmY) ! Compute self-consistently the dynamical effects via perturbation theory for BSE2 @@ -18,6 +19,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n integer,intent(in) :: nS double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) double precision,intent(in) :: OmBSE(nS) @@ -58,7 +60,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n ! Print main components of transition vectors - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) if(dTDA) then write(*,*) diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index d01f978..8ca6757 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) +subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -13,60 +12,67 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eW(nBas) double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - - double precision :: OmRPA(nS,nspin) - double precision :: XpY_RPA(nS,nS,nspin) - double precision :: XmY_RPA(nS,nS,nspin) - double precision :: rho_RPA(nBas,nBas,nS,nspin) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables integer :: ispin integer :: isp_W + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) ! Output variables - double precision,intent(out) :: EcRPA(nspin) double precision,intent(out) :: EcBSE(nspin) ! Memory allocation + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) +!--------------------------------- +! Compute (singlet) RPA screening +!--------------------------------- + + isp_W = 1 + EcRPA = 0d0 + + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + !------------------- ! Singlet manifold !------------------- - if(singlet_manifold) then + if(singlet) then ispin = 1 - isp_W = 1 EcBSE(ispin) = 0d0 - ! Compute (singlet) RPA screening - - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) - ! Compute BSE excitation energies - OmBSE(:,ispin) = OmRPA(:,ispin) - - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -78,12 +84,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) end if end if @@ -94,25 +100,18 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man ! Triplet manifold !------------------- - if(triplet_manifold) then + if(triplet) then ispin = 2 - isp_W = 1 EcBSE(ispin) = 0d0 - ! Compute (singlet) RPA screening - - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) - ! Compute BSE excitation energies - OmBSE(:,ispin) = OmRPA(:,ispin) - - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) !------------------------------------------------- ! Compute the dynamical screening at the BSE level @@ -124,12 +123,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) end if end if diff --git a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 index 25e7dea..2d82750 100644 --- a/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & - Ap,Am,Bp,Bm) +subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -13,8 +12,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -60,8 +59,8 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS - chi_A = chi_A + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - chi_B = chi_B + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) + chi_A = chi_A + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) + chi_B = chi_B + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) enddo @@ -80,28 +79,28 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) enddo diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 index 61dfa01..634323f 100644 --- a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn) +subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -48,7 +48,7 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om chi = 0d0 do kc=1,maxS - chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) enddo @@ -58,10 +58,10 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om do kc=1,maxS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) enddo diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 index c323753..7e000ee 100644 --- a/src/QuAcK/Bethe_Salpeter_B_matrix.f90 +++ b/src/QuAcK/Bethe_Salpeter_B_matrix.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr) +subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,B_lr) ! Compute the extra term for Bethe-Salpeter equation for linear response @@ -11,8 +11,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) ! Local variables @@ -35,8 +35,8 @@ subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho, chi = 0d0 do kc=1,nS - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps + eps = OmRPA(kc)**2 + eta**2 + chi = chi + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps enddo B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 index 14cd120..c7550da 100644 --- a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho, & +subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE, & ZAp,ZAm,ZBp,ZBm) ! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices @@ -13,8 +13,8 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -63,28 +63,28 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, do kc=1,maxS eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 + chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 + chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 + chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 + chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 enddo diff --git a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 index eb71cb4..e15677d 100644 --- a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,ZA_dyn) +subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -12,8 +12,8 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O double precision,intent(in) :: lambda double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -49,10 +49,10 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O do kc=1,maxS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 enddo diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index d01e611..197f2b2 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -1,4 +1,5 @@ -subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, & + OmRPA,rho_RPA,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -17,11 +18,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, integer,intent(in) :: nS double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -55,7 +57,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, ! Print main components of transition vectors - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) if(dTDA) then write(*,*) @@ -84,11 +86,11 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, ! Resonant part of the BSE correction for dynamical TDA - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,Ap_dyn) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn) ! Renormalization factor of the resonant parts for dynamical TDA - call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,ZAp_dyn) + call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) @@ -97,12 +99,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, ! Resonant and anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), & Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) ! Renormalization factor of the resonant and anti-resonant parts - call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & + call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), & ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index 9435c09..c25da2a 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 @@ -1,4 +1,5 @@ -subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, & + OmRPA,rho_RPA,OmBSE,XpY,XmY) ! Compute self-consistently the dynamical effects via perturbation theory for BSE @@ -17,11 +18,12 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, integer,intent(in) :: nS double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -54,7 +56,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Print main components of transition vectors - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) if(dTDA) then write(*,*) @@ -97,8 +99,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Resonant part of the BSE correction - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & - Ap_dyn(:,:)) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) @@ -106,8 +107,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmOld(ia),rho(:,:,:), & - Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia), & + Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & diff --git a/src/QuAcK/CIS.f90 b/src/QuAcK/CIS.f90 index a18660f..ac4923c 100644 --- a/src/QuAcK/CIS.f90 +++ b/src/QuAcK/CIS.f90 @@ -1,5 +1,5 @@ subroutine CIS(singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ERI,eHF) + nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) ! Perform configuration interaction single calculation` @@ -13,6 +13,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, & integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90 index 88a41ce..8204ad1 100644 --- a/src/QuAcK/G0F2.f90 +++ b/src/QuAcK/G0F2.f90 @@ -1,5 +1,5 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) + linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a one-shot second-order Green function calculation @@ -27,6 +27,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -115,7 +116,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & if(BSE) then - call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) + call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) end if diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 231c842..c7ed983 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -1,5 +1,5 @@ -subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) +subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eG0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -17,8 +17,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta @@ -32,6 +32,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -212,10 +213,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing if(BSE) then - allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) - - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eG0T0,EcBSE) if(exchange_kernel) then @@ -249,8 +248,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,EcAC) if(exchange_kernel) then @@ -270,6 +269,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,sing end if + + allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) end if end subroutine G0T0 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index bd76850..058a0c8 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,6 +1,6 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW) +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) ! Perform G0W0 calculation @@ -21,8 +21,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta @@ -33,24 +33,23 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, 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) :: H(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables logical :: print_W = .true. integer :: ispin - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) double precision,allocatable :: eGWlin(:) @@ -68,27 +67,35 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, ! Initialization - EcRPA(:) = 0d0 + EcRPA = 0d0 ! SOSEX correction - if(SOSEX) write(*,*) 'SOSEX correction activated!' - write(*,*) + if(SOSEX) then + write(*,*) 'SOSEX correction activated but BUG!' + stop + end if ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if ! Spin manifold @@ -96,34 +103,43 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, ! Memory allocation - allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & - rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin),eGWlin(nBas)) + allocate(SigC(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGWlin(nBas)) -! Compute linear response +!-------------------! +! Compute screening ! +!-------------------! - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & + eHF,ERI,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - if(print_W) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) -! Compute correlation part of the self-energy +!--------------------------! +! Compute spectral weights ! +!--------------------------! - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) +!------------------------! +! Compute GW self-energy ! +!------------------------! - call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) -! Compute renormalization factor +!--------------------------------! +! Compute renormalization factor ! +!--------------------------------! - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) -! Solve the quasi-particle equation +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! eGWlin(:) = eHF(:) + Z(:)*SigC(:) + ! Linearized or graphical solution? + if(linearize) then write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' @@ -131,51 +147,48 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, eGW(:) = eGWlin(:) -! Find all the roots of the QP equation if necessary - -! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) - else - - ! Find graphical solution of the QP equation - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin,eGW) + + ! Find all the roots of the QP equation if necessary + + ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) end if -! Dump results - - call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) - ! Compute the RPA correlation energy - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) +!--------------! +! Dump results ! +!--------------! + + call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + +! Deallocate memory + + deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) ! Plot stuff -! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) +! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) ! Perform BSE calculation if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) if(exchange_kernel) then - EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(1) + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) end if @@ -204,15 +217,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcAC) - - if(exchange_kernel) then - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) - - end if + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d7f4fbc..4a906c9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -4,7 +4,7 @@ program QuAcK include 'parameters.h' logical :: doSph - logical :: unrestricted + logical :: unrestricted = .false. logical :: doRHF,doUHF,doMOM logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT @@ -52,14 +52,25 @@ program QuAcK integer :: TrialType double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:) - double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:) + double precision,allocatable :: S(:,:) + 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(:,:,:) + double precision,allocatable :: dipole_int_bb(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) - integer :: bra - integer :: ket - double precision,allocatable :: ERI_MO_aa(:,:,:,:) - double precision,allocatable :: ERI_MO_ab(:,:,:,:) - double precision,allocatable :: ERI_MO_bb(:,:,:,:) + integer :: ixyz + integer :: ispin + 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 +112,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 +169,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, & @@ -316,7 +329,6 @@ program QuAcK write(*,*) 'AO to MO transformation... Please be patient' write(*,*) - if(doSph) then allocate(ERI_MO(nBas,nBas,nBas,nBas)) @@ -328,39 +340,81 @@ program QuAcK if(unrestricted) then + ! Read and transform dipole-related integrals + + allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart)) + + call read_dipole_integrals(nBas,dipole_int_aa) + call read_dipole_integrals(nBas,dipole_int_bb) + do ixyz=1,ncart + call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz)) + call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz)) + end do + ! Memory allocation - allocate(ERI_MO_aa(nBas,nBas,nBas,nBas),ERI_MO_ab(nBas,nBas,nBas,nBas),ERI_MO_bb(nBas,nBas,nBas,nBas)) + allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas)) ! 4-index transform for (aa|aa) block - bra = 1 - ket = 1 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aa) + 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 (bb|bb) block - - bra = 1 - ket = 2 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_ab) - ! 4-index transform for (aa|bb) block - bra = 2 - ket = 2 - call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bb) + 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 + + 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 ! Memory allocation allocate(ERI_MO(nBas,nBas,nBas,nBas)) - + + ! Read and transform dipole-related integrals + + allocate(dipole_int(nBas,nBas,ncart)) + call read_dipole_integrals(nBas,dipole_int) + do ixyz=1,ncart + call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz)) + end do + ! 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 @@ -382,7 +436,7 @@ program QuAcK if(unrestricted) then - call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,ENuc,EUHF,eHF,EcMP2) + call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ENuc,EUHF,eHF,EcMP2) else @@ -560,7 +614,16 @@ 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) + if(unrestricted) then + + call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_MO_aaaa,ERI_MO_aabb, & + ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF) + + else + + call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF) + + end if call cpu_time(end_CIS) t_CIS = end_CIS - start_CIS @@ -576,7 +639,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 +655,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,8 +671,16 @@ program QuAcK if(doRPA) then call cpu_time(start_RPA) - call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + if(unrestricted) then + + call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF) + + else + + call dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + + end if call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -625,8 +696,16 @@ program QuAcK if(doRPAx) then call cpu_time(start_RPAx) - call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + if(unrestricted) then + + call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) + + else + + call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) + + end if call cpu_time(end_RPAx) t_RPAx = end_RPAx - start_RPAx @@ -642,8 +721,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_ppRPA) - call ppRPA(singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + call ppRPA(singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_ppRPA) t_ppRPA = end_ppRPA - start_ppRPA @@ -659,7 +737,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,8 +754,8 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGF, & - eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & + eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -694,8 +772,8 @@ 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, & - eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & + ERI_MO,dipole_int,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -747,14 +825,13 @@ program QuAcK call cpu_time(start_G0W0) 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, & - ENuc,EUHF,Hc,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,PHF,cHF,eHF,eG0W0) + 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,ERI_MO_abab, & + 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_manifold,triplet_manifold,linGW,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0) + 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) end if @@ -773,9 +850,19 @@ program QuAcK if(doevGW) then 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, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) + if(unrestricted) then + + 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,ERI_MO_abab,dipole_int_aa,dipole_int_bb, & + PHF,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) + end if call cpu_time(end_evGW) t_evGW = end_evGW - start_evGW @@ -792,8 +879,8 @@ 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, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) + 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,dipole_int,PHF,cHF,eHF) call cpu_time(end_qsGW) t_qsGW = end_qsGW - start_qsGW @@ -811,9 +898,9 @@ program QuAcK if(doG0T0) then 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, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) + call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & + dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0T0) call cpu_time(end_G0T0) t_G0T0 = end_G0T0 - start_G0T0 @@ -829,9 +916,9 @@ program QuAcK if(doevGT) then 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, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) + call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & + BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0T0) call cpu_time(end_evGT) t_evGT = end_evGT - start_evGT @@ -840,10 +927,6 @@ program QuAcK end if -!------------------------------------------------------------------------ -! Perform evGT calculatiom -!------------------------------------------------------------------------ - !------------------------------------------------------------------------ ! Information for Monte Carlo calculations !------------------------------------------------------------------------ @@ -947,8 +1030,8 @@ 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, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) + 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) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0 @@ -961,8 +1044,8 @@ 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, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) + singlet,triplet,linGW,eta_GW, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int,eHF,eG0T0) call cpu_time(end_G0T0) t_G0T0 = end_G0T0 - start_G0T0 diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 index 8afd56d..3e523b4 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -1,5 +1,5 @@ -subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & + ERI,dipole_int,eHF) ! Perform random phase approximation calculation with exchange (aka TDHF) @@ -9,11 +9,12 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Input variables + logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel - logical,intent(in) :: singlet_manifold + logical,intent(in) :: singlet double precision,intent(in) :: eta - logical,intent(in) :: triplet_manifold + logical,intent(in) :: triplet integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -22,8 +23,9 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: e(nBas) + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -55,28 +57,27 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Singlet manifold - if(singlet_manifold) then + if(singlet) then ispin = 1 - call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - + call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif ! Triplet manifold - if(triplet_manifold) then + if(triplet) then ispin = 2 - call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -105,15 +106,8 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & write(*,*) '-------------------------------------------------------' write(*,*) - call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) - - if(exchange_kernel) then - - EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(2) - - end if + call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/SigmaC.f90 b/src/QuAcK/SigmaC.f90 index 32a07ca..10eddc2 100644 --- a/src/QuAcK/SigmaC.f90 +++ b/src/QuAcK/SigmaC.f90 @@ -1,4 +1,4 @@ -double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) +double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Compute diagonal of the correlation part of the self-energy @@ -7,7 +7,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Input variables - integer,intent(in) :: x + integer,intent(in) :: p double precision,intent(in) :: w double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -22,7 +22,7 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Local variables - integer :: i,j,a,b,p,jb + integer :: i,a,jb double precision :: eps ! Initialize @@ -32,26 +32,18 @@ double precision function SigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) ! Occupied part of the correlation self-energy do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = w - e(i) + Omega(jb) - SigmaC = SigmaC + 2d0*rho(x,i,jb)**2*eps/(eps**2 + eta**2) - enddo + do jb=1,nS + eps = w - e(i) + Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2) enddo enddo ! Virtual part of the correlation self-energy do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = w - e(a) - Omega(jb) - SigmaC = SigmaC + 2d0*rho(x,a,jb)**2*eps/(eps**2 + eta**2) - enddo + do jb=1,nS + eps = w - e(a) - Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2) enddo enddo diff --git a/src/QuAcK/UCIS.f90 b/src/QuAcK/UCIS.f90 new file mode 100644 index 0000000..e75fc5a --- /dev/null +++ b/src/QuAcK/UCIS.f90 @@ -0,0 +1,128 @@ +subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,eHF) + +! Perform configuration interaction single calculation` + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + double precision,intent(in) :: eHF(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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart,nspin) + +! Local variables + + logical :: dump_matrix = .false. + logical :: dump_trans = .false. + integer :: ispin + double precision :: lambda + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: A_sc(:,:) + double precision,allocatable :: Omega_sc(:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: A_sf(:,:) + double precision,allocatable :: Omega_sf(:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Configuration Interaction Singles |' + write(*,*)'************************************************' + write(*,*) + +! Adiabatic connection scaling + + lambda = 1d0 + +!----------------------------! +! Spin-conserved transitions ! +!----------------------------! + + if(spin_conserved) then + + ispin = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(A_sc(nS_sc,nS_sc),Omega_sc(nS_sc)) + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sc) + + if(dump_matrix) then + print*,'CIS matrix (spin-conserved transitions)' + call matout(nS_sc,nS_sc,A_sc) + write(*,*) + endif + + call diagonalize_matrix(nS_sc,A_sc,Omega_sc) + call print_excitation('UCIS ',5,nS_sc,Omega_sc) + + if(dump_trans) then + print*,'Spin-conserved CIS transition vectors' + call matout(nS_sc,nS_sc,A_sc) + write(*,*) + endif + + deallocate(A_sc,Omega_sc) + + endif + +!-----------------------! +! Spin-flip transitions ! +!-----------------------! + + if(spin_flip) then + + ispin = 2 + + ! Memory allocation + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + 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)) + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf) + + if(dump_matrix) then + print*,'CIS matrix (spin-conserved transitions)' + call matout(nS_sf,nS_sf,A_sf) + write(*,*) + endif + + call diagonalize_matrix(nS_sf,A_sf,Omega_sf) + call print_excitation('UCIS ',6,nS_sf,Omega_sf) + + if(dump_trans) then + print*,'Spin-flip CIS transition vectors' + call matout(nS_sf,nS_sf,A_sf) + write(*,*) + endif + + deallocate(A_sf,Omega_sf) + + endif + +end subroutine UCIS diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index b71122b..53e3e41 100644 --- a/src/QuAcK/UG0W0.f90 +++ b/src/QuAcK/UG0W0.f90 @@ -1,6 +1,6 @@ -subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, & - singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS, & - ENuc,EUHF,Hc,ERI_aa,ERI_ab,ERI_bb,PHF,cHF,eHF,eGW) +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,ERI_abab, & + dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW) ! Perform unrestricted G0W0 calculation @@ -20,8 +20,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip logical,intent(in) :: linearize double precision,intent(in) :: eta @@ -37,31 +37,28 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev 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) :: ERI_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,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) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) ! Local variables logical :: print_W = .true. + integer :: is integer :: ispin - integer :: iblock - integer :: bra - integer :: ket - integer :: nSa - integer :: nSb - integer :: nSt - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision,allocatable :: SigC(:,:) double precision,allocatable :: Z(:,:) - double precision,allocatable :: Omega(:) - double precision,allocatable :: XpY_a(:,:) - double precision,allocatable :: XpY_b(:,:) - double precision,allocatable :: XmY_a(:,:) - double precision,allocatable :: XmY_b(:,:) - double precision,allocatable :: rho(:,:,:,:) + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) double precision,allocatable :: eGWlin(:,:) @@ -80,77 +77,73 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Initialization - EcRPA(:) = 0d0 + EcRPA = 0d0 ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if ! 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_a(nSa,nSa),XpY_b(nSb,nSb),XmY_a(nSa,nSa),XmY_b(nSb,nSb), & - rho(nBas,nBas,nSt,nspin),eGWlin(nBas,nspin)) + allocate(SigC(nBas,nspin),Z(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc), & + rho_RPA(nBas,nBas,nS_sc,nspin),eGWlin(nBas,nspin)) -! Compute linear response +!-------------------! +! Compute screening ! +!-------------------! -!---------------------------------------------- -! alpha-alpha block -!---------------------------------------------- +! Spin-conserving transition - ispin = 1 - iblock = 3 + ispin = 1 - call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSa,1d0, & - eHF(:,ispin),ERI_aa,rho(:,:,1:nSa,ispin),EcRPA(ispin),Omega(1:nSa),XpY_a,XmY_a) + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - if(print_W) call print_excitation('RPA@HF alpha',iblock,nSa,Omega(1:nSa)) + if(print_W) call print_excitation('RPA@UHF ',5,nS_sc,OmRPA) -!---------------------------------------------- -! alpha-beta block -!---------------------------------------------- - - ispin = 2 - iblock = 3 - - call linear_response(iblock,.true.,TDA_W,.false.,eta,nBas,nC(ispin),nO(ispin),nV(ispin),nR(ispin),nSb,1d0, & - eHF(:,ispin),ERI_bb,rho(:,:,nSa+1:nSt,ispin),EcRPA(ispin),Omega(nSa+1:nSt),XpY_b,XmY_b) - - if(print_W) call print_excitation('RPA@HF beta ',iblock,nSb,Omega(nSa+1:nSt)) - -!---------------------------------------------- -! Excitation densities for alpha self-energy -!---------------------------------------------- +!----------------------! +! Excitation densities ! +!----------------------! - call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho) + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) -!---------------------- -! Compute self-energy -!---------------------- +!---------------------! +! 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_sc,eHF,OmRPA,rho_RPA,SigC) -! Compute renormalization factor +!--------------------------------! +! Compute renormalization factor ! +!--------------------------------! -! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & -! Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! -! Solve the quasi-particle equation - Z(:,:) = 1d0 eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:) if(linearize) then @@ -164,39 +157,32 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! Find graphical solution of the QP equation -! do is=1,nspin -! call QP_graph(nBas,nC(:,is),nO(:,is),nV(:,is),nR(:,is),nS(:,is),eta,eHF(:,is),Omega(:,is), & -! rho(:,:,:,ispin),eGWlin(:,is),eGW(:,is)) -! end do + do is=1,nspin + call unrestricted_QP_graph(nBas,nC(is),nO(is),nV(is),nR(is),nS_sc,eta,eHF(:,is),OmRPA, & + rho_RPA(:,:,:,is),eGWlin(:,is),eGW(:,is)) + end do end if +! Compute RPA correlation energy + + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + ! Dump results - do ispin=1,nspin - call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA(ispin)) - end do + call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA) -! Compute the RPA correlation energy +! Free memory -! call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & -! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) + deallocate(OmRPA,XpY_RPA,XmY_RPA,rho_RPA) ! Perform BSE calculation -! if(BSE) then + if(BSE) then -! call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & -! nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + 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,ERI_abab,dipole_int_aa,dipole_int_bb,eHF,eGW,EcBSE) ! if(exchange_kernel) then ! @@ -251,6 +237,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev ! end if -! end if + end if end subroutine UG0W0 diff --git a/src/QuAcK/UHF.f90 b/src/QuAcK/UHF.f90 index 28fc3c8..34eb9ef 100644 --- a/src/QuAcK/UHF.f90 +++ b/src/QuAcK/UHF.f90 @@ -166,7 +166,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH n_diis = min(n_diis+1,max_diis) do ispin=1,nspin - call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) + if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin), & + err(:,:,ispin),F(:,:,ispin)) end do ! Reset DIIS if required @@ -232,6 +233,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH ! Compute final UHF energy - call print_UHF(nBas,nO(:),e(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),EUHF) + call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) end subroutine UHF diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 new file mode 100644 index 0000000..d7cf08c --- /dev/null +++ b/src/QuAcK/URPAx.f90 @@ -0,0 +1,158 @@ +subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e) + +! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + double precision,intent(in) :: eta + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: e(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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Omega_sc(:) + double precision,allocatable :: XpY_sc(:,:) + double precision,allocatable :: XmY_sc(:,:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: Omega_sf(:) + double precision,allocatable :: XpY_sf(:,:) + double precision,allocatable :: XmY_sf(:,:) + + double precision :: rho_sc,rho_sf + double precision :: EcRPAx(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'*********************************************************************' + write(*,*)'| Unrestricted random phase approximation calculation with exchange |' + write(*,*)'*********************************************************************' + write(*,*) + +! Initialization + + EcRPAx(:) = 0d0 + EcAC(:) = 0d0 + +! Spin-conserved transitions + + if(spin_conserved) then + + ispin = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + + call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) + call print_excitation('URPAx ',5,nS_sc,Omega_sc) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + Omega_sc,XpY_sc,XmY_sc) + + deallocate(Omega_sc,XpY_sc,XmY_sc) + + endif + +! Spin-flip transitions + + if(spin_flip) then + + ispin = 2 + + ! Memory allocation + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + + call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) + call print_excitation('URPAx ',6,nS_sf,Omega_sf) +! call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + deallocate(Omega_sf,XpY_sf,XmY_sf) + + endif + +! if(exchange_kernel) then + +! EcRPAx(1) = 0.5d0*EcRPAx(1) +! EcRPAx(2) = 1.5d0*EcRPAx(2) + +! end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-conserved) =',EcRPAx(1) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy (spin-flip) =',EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx correlation energy =',EcRPAx(1) + EcRPAx(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPAx total energy =',ENuc + EUHF + EcRPAx(1) + EcRPAx(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '-------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of RPAx correlation energy' +! write(*,*) '-------------------------------------------------------' +! write(*,*) + +! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & +! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + +! if(exchange_kernel) then + +! EcAC(1) = 0.5d0*EcAC(1) +! EcAC(2) = 1.5d0*EcAC(2) + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@RPAx correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@RPAx total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine URPAx diff --git a/src/QuAcK/USigmaC.f90 b/src/QuAcK/USigmaC.f90 new file mode 100644 index 0000000..c5e0817 --- /dev/null +++ b/src/QuAcK/USigmaC.f90 @@ -0,0 +1,48 @@ +double precision function USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,jb + double precision :: eps + +! Initialize + + USigmaC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + USigmaC = uSigmaC + rho(p,i,jb)**2*eps/(eps**2 + eta**2) + end do + end do + + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + USigmaC = USigmaC + rho(p,a,jb)**2*eps/(eps**2 + eta**2) + end do + end do + +end function USigmaC diff --git a/src/QuAcK/UdRPA.f90 b/src/QuAcK/UdRPA.f90 new file mode 100644 index 0000000..fd1aa0e --- /dev/null +++ b/src/QuAcK/UdRPA.f90 @@ -0,0 +1,156 @@ +subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e) + +! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + double precision,intent(in) :: eta + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + double precision,intent(in) :: e(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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Omega_sc(:) + double precision,allocatable :: XpY_sc(:,:) + double precision,allocatable :: XmY_sc(:,:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: Omega_sf(:) + double precision,allocatable :: XpY_sf(:,:) + double precision,allocatable :: XmY_sf(:,:) + + double precision :: rho_sc,rho_sf + double precision :: EcRPA(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'**************************************************************' + write(*,*)'| Unrestricted direct random phase approximation calculation |' + write(*,*)'**************************************************************' + write(*,*) + +! Initialization + + EcRPA(:) = 0d0 + EcAC(:) = 0d0 + +! Spin-conserved transitions + + if(spin_conserved) then + + ispin = 1 + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) + + call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) + call print_excitation('URPA ',5,nS_sc,Omega_sc) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + Omega_sc,XpY_sc,XmY_sc) + + + endif + +! Spin-flip transitions + + if(spin_flip) then + + ispin = 2 + + ! Memory allocation + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) + + call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) + call print_excitation('URPA ',6,nS_sf,Omega_sf) +! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + + + endif + +! if(exchange_kernel) then + +! EcRPA(1) = 0.5d0*EcRPA(1) +! EcRPA(2) = 1.5d0*EcRPA(2) + +! end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-conserved) =',EcRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy (spin-flip) =',EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@URPA total energy =',ENuc + EUHF + EcRPA(1) + EcRPA(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '-------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of RPA correlation energy' +! write(*,*) '-------------------------------------------------------' +! write(*,*) + +! call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & +! nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + +! if(exchange_kernel) then + +! EcAC(1) = 0.5d0*EcAC(1) +! EcAC(2) = 1.5d0*EcAC(2) + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + EUHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +end subroutine UdRPA diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/dRPA.f90 similarity index 72% rename from src/QuAcK/RPA.f90 rename to src/QuAcK/dRPA.f90 index 32d9587..fa45771 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/dRPA.f90 @@ -1,5 +1,4 @@ -subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a direct random phase approximation calculation @@ -9,10 +8,11 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Input variables + logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -22,8 +22,9 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: e(nBas) + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -40,7 +41,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & write(*,*) write(*,*)'***********************************************' - write(*,*)'| random-phase approximation calculation |' + write(*,*)'| Random-phase approximation calculation |' write(*,*)'***********************************************' write(*,*) @@ -55,32 +56,34 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & ! Singlet manifold - if(singlet_manifold) then + if(singlet) then ispin = 1 - call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif ! Triplet manifold - if(triplet_manifold) then + if(triplet) then ispin = 2 - call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif if(exchange_kernel) then EcRPA(1) = 0.5d0*EcRPA(1) - EcRPA(2) = 1.5d0*EcRPA(1) + EcRPA(2) = 1.5d0*EcRPA(2) end if @@ -103,13 +106,13 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & write(*,*) '------------------------------------------------------' write(*,*) - call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,e,e,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) if(exchange_kernel) then EcAC(1) = 0.5d0*EcAC(1) - EcAC(2) = 1.5d0*EcAC(1) + EcAC(2) = 1.5d0*EcAC(2) end if @@ -125,4 +128,4 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & end if -end subroutine RPA +end subroutine dRPA diff --git a/src/QuAcK/dUSigmaC.f90 b/src/QuAcK/dUSigmaC.f90 new file mode 100644 index 0000000..418e131 --- /dev/null +++ b/src/QuAcK/dUSigmaC.f90 @@ -0,0 +1,50 @@ +double precision function dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) + +! Compute the derivative of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: p + double precision,intent(in) :: w + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,jb + double precision :: eps + +! Initialize + + dUSigmaC = 0d0 + +! Occupied part of the correlation self-energy + + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + dUSigmaC = dUSigmaC + rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + +! Virtual part of the correlation self-energy + + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + dUSigmaC = dUSigmaC + rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + +end function dUSigmaC diff --git a/src/QuAcK/evGF2.f90 b/src/QuAcK/evGF2.f90 index d8b6411..168a762 100644 --- a/src/QuAcK/evGF2.f90 +++ b/src/QuAcK/evGF2.f90 @@ -1,5 +1,5 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) + linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -30,6 +30,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -168,7 +169,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold if(BSE) then - call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) + call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) end if diff --git a/src/QuAcK/evGT.f90 b/src/QuAcK/evGT.f90 index 1e8c030..c923a20 100644 --- a/src/QuAcK/evGT.f90 +++ b/src/QuAcK/evGT.f90 @@ -1,6 +1,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & - eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) + BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eG0T0) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -21,8 +21,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -35,6 +35,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eG0T0(nBas) @@ -260,10 +261,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & if(BSE) then - allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) - - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGT,eGT,EcRPA,EcBSE) if(exchange_kernel) then @@ -297,8 +296,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,EcAC) if(exchange_kernel) then diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 625c1fa..26061f8 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,5 +1,6 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0, & - singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) +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) ! Perform self-consistent eigenvalue-only GW calculation @@ -26,8 +27,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: evDyn logical,intent(in) :: G0W logical,intent(in) :: GW0 - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eHF(nBas) @@ -35,8 +36,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eG0W0(nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: H(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -46,7 +47,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE integer :: n_diis double precision :: rcond double precision :: Conv - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM @@ -57,11 +58,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) double precision,allocatable :: SigC(:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) ! Hello world @@ -73,23 +73,45 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! SOSEX correction - if(SOSEX) write(*,*) 'SOSEX correction activated!' - write(*,*) + if(SOSEX) then + write(*,*) 'SOSEX correction activated but BUG!' + stop + end if ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! GW0 + + if(GW0) then + write(*,*) 'GW0 scheme activated!' + write(*,*) + end if + +! G0W + + if(G0W) then + write(*,*) 'G0W scheme activated!' + write(*,*) + end if ! Linear mixing @@ -98,10 +120,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Memory allocation - allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), & - XpY(nS,nS,nspin),XmY(nS,nS,nspin), & - rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & - error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Initialization @@ -121,36 +141,30 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE do while(Conv > thresh .and. nSCF <= maxSCF) - ! Compute linear response + ! Compute screening if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & + rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) endif +! Compute spectral weights + + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + ! Compute correlation part of the self-energy - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) - - ! Correlation self-energy - if(G0W) then - call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else - call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) endif @@ -164,8 +178,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE ! Print results -! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW) + call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) ! Linear mixing or DIIS extrapolation @@ -215,23 +228,22 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE endif -! Dump the RPA correlation energy +! Deallocate memory - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) + deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis) ! Perform BSE calculation if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -258,8 +270,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/evUGW.f90 b/src/QuAcK/evUGW.f90 new file mode 100644 index 0000000..322e13e --- /dev/null +++ b/src/QuAcK/evUGW.f90 @@ -0,0 +1,308 @@ +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,ERI_abab,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0) + +! Perform self-consistent eigenvalue-only GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: COHSEX + logical,intent(in) :: BSE + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: G0W + logical,intent(in) :: GW0 + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + double precision,intent(in) :: eta + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + + 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) :: 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) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + logical :: linear_mixing + integer :: is + integer :: ispin + integer :: nSCF + integer :: n_diis + double precision :: rcond(nspin) + double precision :: Conv + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcGM + double precision :: alpha + double precision,allocatable :: error_diis(:,:,:) + double precision,allocatable :: e_diis(:,:,:) + double precision,allocatable :: eGW(:,:) + double precision,allocatable :: eOld(:,:) + double precision,allocatable :: Z(:,:) + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: SigC(:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent evGW calculation |' + write(*,*)'************************************************' + write(*,*) + +! COHSEX approximation + + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! GW0 + + if(GW0) then + write(*,*) 'GW0 scheme activated!' + write(*,*) + end if + +! G0W + + if(G0W) then + write(*,*) 'G0W scheme activated!' + write(*,*) + end if + +! Linear mixing + + linear_mixing = .false. + alpha = 0.2d0 + +! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(eGW(nBas,nspin),eOld(nBas,nspin),Z(nBas,nspin),SigC(nBas,nspin),OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc), & + XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin),error_diis(nBas,max_diis,nspin),e_diis(nBas,max_diis,nspin)) + +! Initialization + + nSCF = 0 + ispin = 1 + n_diis = 0 + Conv = 1d0 + e_diis(:,:,:) = 0d0 + error_diis(:,:,:) = 0d0 + eGW(:,:) = eG0W0(:,:) + eOld(:,:) = eGW(:,:) + Z(:,:) = 1d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Compute screening + + if(.not. GW0 .or. nSCF == 0) then + + call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + + + endif + + !----------------------! + ! Excitation densities ! + !----------------------! + + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + !------------------------------------------------! + ! Compute self-energy and renormalization factor ! + !------------------------------------------------! + + if(G0W) then + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + else + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + + endif + + !-----------------------------------! + ! Solve the quasi-particle equation ! + !-----------------------------------! + + eGW(:,:) = eHF(:,:) + SigC(:,:) + + ! Convergence criteria + + Conv = maxval(abs(eGW(:,:) - eOld(:,:))) + + ! Print results + + call print_evUGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA) + + ! Linear mixing or DIIS extrapolation + + if(linear_mixing) then + + eGW(:,:) = alpha*eGW(:,:) + (1d0 - alpha)*eOld(:,:) + + else + + n_diis = min(n_diis+1,max_diis) + do is=1,nspin + call DIIS_extrapolation(rcond(ispin),nBas,nBas,n_diis,error_diis(:,:,is), & + e_diis(:,:,is),eGW(:,is)-eOld(:,is),eGW(:,is)) + end do + +! Reset DIIS if required + + if(minval(rcond(:)) < 1d-15) n_diis = 0 + + endif + + ! Save quasiparticles energy for next cycle + + eOld(:,:) = eGW(:,:) + + ! Increment + + nSCF = nSCF + 1 + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Plot stuff + +! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Deallocate memory + + deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis) + +! 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,nS, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,eGW,eGW,EcBSE) + +! if(exchange_kernel) then + +! EcBSE(1) = 0.5d0*EcBSE(1) +! EcBSE(2) = 1.5d0*EcBSE(2) + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of BSE correlation energy' +! write(*,*) '------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + + endif + +end subroutine evUGW diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90 index c3ea157..e443d52 100644 --- a/src/QuAcK/linear_response.f90 +++ b/src/QuAcK/linear_response.f90 @@ -1,4 +1,4 @@ -subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY) +subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY) ! Compute linear response @@ -12,8 +12,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas) - double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: Omega_RPA(nS) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS) ! Local variables @@ -42,7 +44,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) - if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A) + if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,A) ! Tamm-Dancoff approximation @@ -51,7 +53,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E 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,rho,B) + if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B) end if diff --git a/src/QuAcK/plot_GW.f90 b/src/QuAcK/plot_GW.f90 index 11c7f0c..fedc40c 100644 --- a/src/QuAcK/plot_GW.f90 +++ b/src/QuAcK/plot_GW.f90 @@ -1,4 +1,4 @@ -subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox) +subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho) ! Dump several GW quantities for external plotting @@ -8,7 +8,7 @@ subroutine plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega,rho,rhox) ! Input variables integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS) + double precision,intent(in) :: eHF(nBas),eGW(nBas),Omega(nS),rho(nBas,nBas,nS) ! Local variables diff --git a/src/QuAcK/print_G0W0.f90 b/src/QuAcK/print_G0W0.f90 index 95d1051..782643b 100644 --- a/src/QuAcK/print_G0W0.f90 +++ b/src/QuAcK/print_G0W0.f90 @@ -40,11 +40,11 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA -! write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA -! write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM -! write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM -! write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM + write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM + write(*,*)'-------------------------------------------------------------------------------' write(*,*) end subroutine print_G0W0 diff --git a/src/QuAcK/print_UG0W0.f90 b/src/QuAcK/print_UG0W0.f90 index b73e886..0e8092e 100644 --- a/src/QuAcK/print_UG0W0.f90 +++ b/src/QuAcK/print_UG0W0.f90 @@ -1,44 +1,71 @@ -subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA) +subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA) ! Print one-electron energies and other stuff for G0W0 implicit none include 'parameters.h' - integer,intent(in) :: nBas,nO + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) double precision,intent(in) :: ENuc double precision,intent(in) :: EHF double precision,intent(in) :: EcRPA - double precision,intent(in) :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: SigC(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGW(nBas,nspin) - integer :: x,HOMO,LUMO - double precision :: Gap + integer :: p + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) ! HOMO and LUMO - HOMO = nO - LUMO = HOMO + 1 - Gap = eGW(LUMO)-eGW(HOMO) + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGW(nO(ispin),ispin) + LUMO(ispin) = eGW(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1,ispin) + Gap(ispin) = 0d0 + end if + end do ! Dump results - write(*,*)'-------------------------------------------------------------------------------' - write(*,*)' One-shot G0W0 calculation' - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & - '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|' - write(*,*)'-------------------------------------------------------------------------------' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*)' Unrestricted one-shot G0W0 calculation (eV)' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' + write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & + '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' - do x=1,nBas - write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' enddo - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',eGW(HOMO)*HaToeV - write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV - write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV - write(*,*)'-------------------------------------------------------------------------------' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',maxval(HOMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',minval(LUMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' write(*,*) end subroutine print_UG0W0 diff --git a/src/QuAcK/print_UHF.f90 b/src/QuAcK/print_UHF.f90 index 6c52e4d..640d023 100644 --- a/src/QuAcK/print_UHF.f90 +++ b/src/QuAcK/print_UHF.f90 @@ -16,18 +16,24 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: EUHF - integer :: HOMO(nspin) - integer :: LUMO(nspin) + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) double precision :: Gap(nspin) ! HOMO and LUMO - HOMO(:) = nO(:) - - LUMO(:) = HOMO(:) + 1 - - Gap(1) = e(LUMO(1),1) - e(HOMO(1),1) - Gap(2) = e(LUMO(2),2) - e(HOMO(2),2) + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = e(nO(ispin),ispin) + LUMO(ispin) = e(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1,ispin) + Gap(ispin) = 0d0 + end if + end do ! Dump results @@ -62,12 +68,12 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' write(*,'(A40,1X,F16.10,A3)') ' UHF energy: ',EUHF + ENuc,' au' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',e(HOMO(1),1)*HatoeV,' eV' - write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',e(LUMO(1),1)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',HOMO(1)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',LUMO(1)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',e(HOMO(2),2)*HatoeV,' eV' - write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',e(LUMO(2),2)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',HOMO(2)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' write(*,*) @@ -78,6 +84,7 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) write(*,'(A50)') 'UHF spin-up orbital coefficients ' write(*,'(A50)') '-----------------------------------------' call matout(nBas,nBas,c(:,:,1)) + write(*,*) write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') 'UHF spin-down orbital coefficients ' write(*,'(A50)') '-----------------------------------------' diff --git a/src/QuAcK/print_evGW.f90 b/src/QuAcK/print_evGW.f90 index 0bf95b5..8db2956 100644 --- a/src/QuAcK/print_evGW.f90 +++ b/src/QuAcK/print_evGW.f90 @@ -1,4 +1,4 @@ -subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) +subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA,EcGM) ! Print one-electron energies and other stuff for evGW @@ -8,9 +8,15 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) integer,intent(in) :: nBas,nO,nSCF double precision,intent(in) :: ENuc double precision,intent(in) :: EHF - double precision,intent(in) :: Conv,e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas) + double precision,intent(in) :: Conv + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: SigC(nBas) + double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: EcRPA + double precision,intent(in) :: EcGM - integer :: x,HOMO,LUMO + integer :: p,HOMO,LUMO double precision :: Gap ! HOMO and LUMO @@ -32,9 +38,9 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do x=1,nBas + do p=1,nBas write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & - '|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|' + '|',p,'|',e(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' enddo write(*,*)'-------------------------------------------------------------------------------' @@ -45,11 +51,11 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW) write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA -! write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA -! write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM -! write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM -! write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM + write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM + write(*,*)'-------------------------------------------------------------------------------' write(*,*) end subroutine print_evGW diff --git a/src/QuAcK/print_evUGW.f90 b/src/QuAcK/print_evUGW.f90 new file mode 100644 index 0000000..984d88c --- /dev/null +++ b/src/QuAcK/print_evUGW.f90 @@ -0,0 +1,81 @@ +subroutine print_evUGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigC,Z,eGW,EcRPA) + +! Print one-electron energies and other stuff for evGW + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nSCF + double precision,intent(in) :: ENuc + double precision,intent(in) :: EHF + double precision,intent(in) :: EcRPA + double precision,intent(in) :: Conv + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: SigC(nBas,nspin) + double precision,intent(in) :: Z(nBas,nspin) + double precision,intent(in) :: eGW(nBas,nspin) + + integer :: p + integer :: ispin + double precision :: HOMO(nspin) + double precision :: LUMO(nspin) + double precision :: Gap(nspin) + +! HOMO and LUMO + + do ispin=1,nspin + if(nO(ispin) > 0) then + HOMO(ispin) = eGW(nO(ispin),ispin) + LUMO(ispin) = eGW(nO(ispin)+1,ispin) + Gap(ispin) = LUMO(ispin) - HOMO(ispin) + else + HOMO(ispin) = 0d0 + LUMO(ispin) = e(1,ispin) + Gap(ispin) = 0d0 + end if + end do + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent evG',nSCF,'W',nSCF,' calculation' + endif + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') & + '|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|' + write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') & + '|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|' + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + + do p=1,nBas + write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') & + '|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & + Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'evGW HOMO energy (eV):',maxval(HOMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',minval(LUMO(:))*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA + write(*,*)'-------------------------------------------------------------------------------& + -------------------------------------------------' + write(*,*) + +end subroutine print_evUGW diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index ea8fa8b..c2bbdb4 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -7,24 +7,26 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Input variables - character*12,intent(in) :: method + character*12,intent(in) :: method integer,intent(in) :: ispin,nS double precision,intent(in) :: Omega(nS) ! Local variables - character*7 :: spin_manifold + character*14 :: spin_manifold integer,parameter :: maxS = 32 integer :: ia if(ispin == 1) spin_manifold = 'singlet' if(ispin == 2) spin_manifold = 'triplet' - if(ispin == 3) spin_manifold = 'alp-bet' - if(ispin == 4) spin_manifold = 'alp-alp' + if(ispin == 3) spin_manifold = 'alpha-beta' + if(ispin == 4) spin_manifold = 'alpha-alpha' + if(ispin == 5) spin_manifold = 'spin-conserved' + if(ispin == 6) spin_manifold = 'spin-flip' write(*,*) write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A1,1X,A14,A14,A7,A9,A15)')'|',method,' calculation: ',spin_manifold,' manifold',' |' + write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold' write(*,*)'-------------------------------------------------------------' write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' diff --git a/src/QuAcK/print_qsGW.f90 b/src/QuAcK/print_qsGW.f90 index b61bdd3..9f3cf62 100644 --- a/src/QuAcK/print_qsGW.f90 +++ b/src/QuAcK/print_qsGW.f90 @@ -70,8 +70,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig write(*,'(2X,A30,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM write(*,'(2X,A30,F15.6)') 'qsGW exchange energy =',Ex write(*,'(2X,A30,F15.6)') 'qsGW correlation energy =',Ec -! write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA -! write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM + write(*,'(2X,A30,F15.6)') 'RPA@qsGW correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@qsGW correlation energy =',EcGM write(*,*)'-------------------------------------------' write(*,*) @@ -95,7 +95,6 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig write(*,'(A32,1X,F16.10)') ' Electronic energy ',EqsGW write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc write(*,'(A32,1X,F16.10)') ' qsGW energy ',ENuc + EqsGW -! write(*,'(A32,1X,F16.10)') ' RPA corr. energy ',EcRPA write(*,'(A50)') '---------------------------------------' write(*,*) diff --git a/src/QuAcK/print_transition_vectors.f90 b/src/QuAcK/print_transition_vectors.f90 index ab3e093..3c30621 100644 --- a/src/QuAcK/print_transition_vectors.f90 +++ b/src/QuAcK/print_transition_vectors.f90 @@ -1,4 +1,4 @@ -subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) +subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY) ! Print transition vectors for linear response calculation @@ -7,44 +7,86 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) ! 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) :: nS + double precision :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) ! Local variables + logical :: debug = .false. integer :: ia,jb,i,j,a,b + integer :: ixyz integer,parameter :: maxS = 10 - double precision,parameter :: thres_vec = 0.1d0 double precision :: norm + double precision,parameter :: thres_vec = 0.1d0 double precision,allocatable :: X(:) double precision,allocatable :: Y(:) + double precision,allocatable :: f(:,:) + double precision,allocatable :: os(:) ! Memory allocation - allocate(X(nS),Y(nS)) + allocate(X(nS),Y(nS),f(nS,ncart),os(nS)) + +! Compute dipole moments and oscillator strengths + + + f(:,:) = 0d0 + if(spin_allowed) then + + do ia=1,nS + do ixyz=1,ncart + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz)*XpY(ia,jb) + end do + end do + end do + end do + f(:,:) = sqrt(2d0)*f(:,:) + + do ia=1,nS + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + if(debug) then + + write(*,*) '------------------------' + write(*,*) ' Dipole moments (X Y Z) ' + write(*,*) '------------------------' + call matout(nS,ncart,f) + write(*,*) + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nS,1,os) + write(*,*) + + end if + + end if + +! Print details about excitations - write(*,*) do ia=1,min(nS,maxS) X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - - norm = 0d0 - do jb=1,nS - norm = norm + X(jb)*X(jb) - end do - norm = sqrt(norm) - print*,'--------------------------------' - write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV' - print*,'--------------------------------' + print*,'---------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')' + print*,'---------------------------------------------' jb = 0 do j=nC+1,nO @@ -54,12 +96,6 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) end do end do - norm = 0d0 - do jb=1,nS - norm = norm + Y(jb)*Y(jb) - end do - norm = sqrt(norm) - jb = 0 do j=nC+1,nO do b=nO+1,nBas-nR @@ -71,4 +107,7 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) end do + write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) + write(*,*) + end subroutine print_transition_vectors diff --git a/src/QuAcK/print_unrestricted_transition_vectors.f90 b/src/QuAcK/print_unrestricted_transition_vectors.f90 new file mode 100644 index 0000000..256f83c --- /dev/null +++ b/src/QuAcK/print_unrestricted_transition_vectors.f90 @@ -0,0 +1,148 @@ +subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, & + Omega,XpY,XmY) + +! Print transition vectors for linear response calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: spin_allowed + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + double precision :: dipole_int_aa(nBas,nBas,ncart) + double precision :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: XpY(nSt,nSt) + double precision,intent(in) :: XmY(nSt,nSt) + +! Local variables + + logical :: debug = .false. + integer :: ia,jb,i,j,a,b + integer :: ixyz + integer :: ispin + integer,parameter :: maxS = 10 + double precision :: norm + double precision,parameter :: thres_vec = 0.1d0 + double precision,allocatable :: X(:) + double precision,allocatable :: Y(:) + double precision,allocatable :: f(:,:) + double precision,allocatable :: os(:) + +! Memory allocation + + allocate(X(nSt),Y(nSt),f(nSt,ncart),os(nSt)) + +! Compute dipole moments and oscillator strengths + + + f(:,:) = 0d0 + if(spin_allowed) then + + do ia=1,nSt + do ixyz=1,ncart + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(j,b,ixyz)*XpY(ia,jb) + end do + end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(j,b,ixyz)*XpY(ia,nSa+jb) + end do + end do + + end do + end do + + do ia=1,nSt + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + if(debug) then + + write(*,*) '----------------' + write(*,*) ' Dipole moments ' + write(*,*) '----------------' + call matout(nSt,ncart,f(:,:)) + write(*,*) + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nSt,1,os(:)) + write(*,*) + + end if + + end if + +! Print details about excitations + + do ia=1,min(nSt,maxS) + + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + print*,'---------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')' + print*,'---------------------------------------------' + + ! Spin-up transitions + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)/sqrt(2d0) + end do + end do + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)/sqrt(2d0) + end do + end do + + ! Spin-down transitions + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(jb)/sqrt(2d0) + end do + end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(jb)/sqrt(2d0) + end do + end do + write(*,*) + + end do + + write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) + write(*,*) + +end subroutine print_unrestricted_transition_vectors diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 8191c4b..2bdfaf3 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,6 +1,6 @@ -subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) +subroutine qsGW(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,S,X,T,V, & + Hc,ERI_AO_basis,ERI_MO_basis,dipole_int,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -25,8 +25,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, logical,intent(in) :: evDyn logical,intent(in) :: G0W logical,intent(in) :: GW0 - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc @@ -41,6 +41,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_MO_basis(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -49,7 +50,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, integer :: ispin integer :: n_diis double precision :: EqsGW - double precision :: EcRPA(nspin) + double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) double precision :: EcGM @@ -58,11 +59,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, double precision,external :: trace_matrix double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGW(:) @@ -96,29 +96,37 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! SOSEX correction - if(SOSEX) write(*,*) 'SOSEX correction activated!' - write(*,*) + if(SOSEX) then + write(*,*) 'SOSEX correction activated but BUG!' + stop + end if ! COHSEX approximation - if(COHSEX) write(*,*) 'COHSEX approximation activated!' - write(*,*) + if(COHSEX) then + write(*,*) 'COHSEX approximation activated!' + write(*,*) + end if ! TDA for W - if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if ! TDA - if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if ! Memory allocation allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & - Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & + OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -149,36 +157,30 @@ 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 if(.not. GW0 .or. nSCF == 0) then call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) endif ! Compute correlation part of the self-energy - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) - - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY_RPA,rho_RPA) if(G0W) then - call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else - call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & - Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) endif @@ -221,7 +223,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! 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(ispin),EcGM,EqsGW) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW) ! Increment @@ -258,23 +260,23 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, endif -! Dump RPA correlation energy +! Deallocate memory - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis) ! Perform BSE calculation if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,dipole_int, & + eGW,eGW,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -301,8 +303,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, end if - call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcAC) + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' 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/renormalization_factor.f90 b/src/QuAcK/renormalization_factor.f90 index 070c1e5..10f2849 100644 --- a/src/QuAcK/renormalization_factor.f90 +++ b/src/QuAcK/renormalization_factor.f90 @@ -1,4 +1,4 @@ -subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) +subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) ! Compute renormalization factor for GW @@ -8,13 +8,11 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r ! Input variables logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables @@ -35,42 +33,8 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r Z(:) = 1d0 return - - end if - -! Occupied part of the correlation self-energy - - do x=nC+1,nBas-nR - do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) - Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - end do - -! Virtual part of the correlation self-energy - - do x=nC+1,nBas-nR - do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) - Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - end do - - ! SOSEX correction - - if(SOSEX) then + + else ! Occupied part of the correlation self-energy @@ -81,7 +45,7 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r do b=nO+1,nBas-nR jb = jb + 1 eps = e(x) - e(i) + Omega(jb) - Z(x) = Z(x) - (rho(x,i,jb)/eps)*(rhox(x,i,jb)/eps) + Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2 end do end do end do @@ -96,13 +60,13 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r do b=nO+1,nBas-nR jb = jb + 1 eps = e(x) - e(a) - Omega(jb) - Z(x) = Z(x) - (rho(x,a,jb)/eps)*(rhox(x,a,jb)/eps) + Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2 end do end do end do end do - endif + end if ! Compute renormalization factor from derivative of SigC diff --git a/src/QuAcK/self_energy_correlation.f90 b/src/QuAcK/self_energy_correlation.f90 index 5baccea..e47a9b4 100644 --- a/src/QuAcK/self_energy_correlation.f90 +++ b/src/QuAcK/self_energy_correlation.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) +subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Compute correlation part of the self-energy @@ -8,13 +8,11 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, ! Input variables logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables @@ -30,7 +28,9 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, SigC = 0d0 -! COHSEX static approximation +!-----------------------------! +! COHSEX static approximation ! +!-----------------------------! if(COHSEX) then @@ -65,6 +65,10 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, else +!----------------! +! GW self-energy ! +!----------------! + ! Occupied part of the correlation self-energy do x=nC+1,nBas-nR @@ -91,36 +95,6 @@ subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega, enddo enddo - if(SOSEX) then - - ! SOSEX: occupied part of the correlation self-energy - - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(x) - e(i) + Omega(jb) - SigC(x,y) = SigC(x,y) - rho(x,i,jb)*rhox(y,i,jb)*eps/(eps**2 + eta**2) - enddo - enddo - enddo - enddo - - ! SOSEX: virtual part of the correlation self-energy - - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(x) - e(a) - Omega(jb) - SigC(x,y) = SigC(x,y) - rho(x,a,jb)*rhox(y,a,jb)*eps/(eps**2 + eta**2) - enddo - enddo - enddo - enddo - - endif - endif end subroutine self_energy_correlation diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90 index 03afb7a..59aab3d 100644 --- a/src/QuAcK/self_energy_correlation_diag.f90 +++ b/src/QuAcK/self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) +subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Compute diagonal of the correlation part of the self-energy @@ -8,7 +8,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O ! Input variables logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -19,7 +18,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) - double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables @@ -66,46 +64,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O EcGM = 0d0 do i=nC+1,nO - EcGM = EcGM + 0.5d0*SigC(i) - end do - -!----------------------------- -! SOSEX self-energy *BUG* -!----------------------------- - - elseif(SOSEX) then - - ! SOSEX: occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - SigC(p) = SigC(p) - rho(p,i,jb)*rhox(p,i,jb)*eps/(eps**2 + eta**2) - end do - end do - end do - - ! SOSEX: virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - SigC(p) = SigC(p) - rho(p,a,jb)*rhox(p,a,jb)*eps/(eps**2 + eta**2) - end do - end do - end do - - ! GM correlation energy - - do i=nC+1,nO - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2) - end do - end do + EcGM = EcGM - SigC(i) end do !----------------------------- @@ -143,7 +102,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,O do a=nO+1,nBas-nR do jb=1,nS eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) + EcGM = EcGM - 4d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) end do end do end do diff --git a/src/QuAcK/spatial_to_spin_MO_energy.f90 b/src/QuAcK/spatial_to_spin_MO_energy.f90 index 688dc1d..0a35a8c 100644 --- a/src/QuAcK/spatial_to_spin_MO_energy.f90 +++ b/src/QuAcK/spatial_to_spin_MO_energy.f90 @@ -1,6 +1,6 @@ subroutine spatial_to_spin_MO_energy(nBas,e,nBas2,se) -! Convert ERIs from spatial to spin orbitals +! Convert MO energies from spatial to spin orbitals implicit none diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 new file mode 100644 index 0000000..1b7566a --- /dev/null +++ b/src/QuAcK/unrestricted_Bethe_Salpeter.f90 @@ -0,0 +1,171 @@ +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,ERI_abab, & + dipole_int_aa,dipole_int_bb,eW,eGW,EcBSE) + +! Compute the Bethe-Salpeter excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + 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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + integer :: isp_W + + integer :: nS_aa,nS_bb,nS_sc + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) + + double precision,allocatable :: OmBSE_sc(:) + double precision,allocatable :: XpY_BSE_sc(:,:) + double precision,allocatable :: XmY_BSE_sc(:,:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: OmBSE_sf(:) + double precision,allocatable :: XpY_BSE_sf(:,:) + double precision,allocatable :: XmY_BSE_sf(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE(nspin) + + ! Memory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(OmRPA(nS_sc),XpY_RPA(nS_sc,nS_sc),XmY_RPA(nS_sc,nS_sc),rho_RPA(nBas,nBas,nS_sc,nspin)) + +!--------------------------! +! Spin-conserved screening ! +!--------------------------! + + isp_W = 1 + EcRPA = 0d0 + + ! Compute spin-conserved RPA screening + + call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & + eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + + call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) + + +!----------------------------! +! Spin-conserved excitations ! +!----------------------------! + + if(spin_conserved) then + + ispin = 1 + EcBSE(ispin) = 0d0 + + allocate(OmBSE_sc(nS_sc),XpY_BSE_sc(nS_sc,nS_sc),XmY_BSE_sc(nS_sc,nS_sc)) + + ! Compute spin-conserved BSE excitation energies + + 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,ERI_abab,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(.true.,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) + + !------------------------------------------------- + ! Compute the dynamical screening at the BSE level + !------------------------------------------------- + +! if(dBSE) then + +! ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) +! +! if(evDyn) then +! +! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & +! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) +! else +! +! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & +! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) +! end if + +! end if + + end if + +!-----------------------! +! Spin-flip excitations ! +!-----------------------! + + if(spin_flip) then + + ispin = 2 + EcBSE(ispin) = 0d0 + + ! Memory allocation + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + allocate(OmBSE_sf(nS_sf),XpY_BSE_sf(nS_sf,nS_sf),XmY_BSE_sf(nS_sf,nS_sf)) + + call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sc,1d0, & + eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcBSE(ispin), & + OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) + + call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf) + + !------------------------------------------------- + ! Compute the dynamical screening at the BSE level + !------------------------------------------------- + +! if(dBSE) then + +! ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) + +! if(evDyn) then +! +! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & +! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) +! else +! +! call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & +! XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) +! end if + +! end if + + end if + +end subroutine unrestricted_Bethe_Salpeter diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 new file mode 100644 index 0000000..1b1e415 --- /dev/null +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90 @@ -0,0 +1,152 @@ +subroutine unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,A_lr) + +! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + 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_sc + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + 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) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS_sc) + double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: A_lr(nSt,nSt) + +!--------------------------------------------------! +! Build BSE matrix for spin-conserving transitions ! +!--------------------------------------------------! + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,1)*rho(a,b,kc,1)*Omega(kc)/eps + enddo + + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_aaaa(i,b,j,a) + 2d0*lambda*chi + + enddo + enddo + enddo + enddo + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,2)*Omega(kc)/eps + enddo + + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,b,j,a) + 2d0*lambda*chi + + enddo + enddo + enddo + enddo + + end if + +!--------------------------------------------! +! Build BSE matrix for spin-flip transitions ! +!--------------------------------------------! + + if(ispin == 2) then + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,1)*rho(a,b,kc,2)*Omega(kc)/eps + enddo + + A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI_abab(i,b,j,a) + 2d0*lambda*chi + + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,j,kc,2)*rho(a,b,kc,1)*Omega(kc)/eps + enddo + + A_lr(nSa+ia,nSa+jb) = A_lr(nSa+ia,nSa+jb) - lambda*ERI_abab(b,i,a,j) + 2d0*lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine unrestricted_Bethe_Salpeter_A_matrix diff --git a/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 new file mode 100644 index 0000000..2f559ee --- /dev/null +++ b/src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90 @@ -0,0 +1,154 @@ +subroutine unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega,rho,B_lr) + +! Compute the extra term for Bethe-Salpeter equation for linear response + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + 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_sc + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + 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) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Omega(nS_sc) + double precision,intent(in) :: rho(nBas,nBas,nS_sc,nspin) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: B_lr(nSt,nSt) + +!--------------------------------------------------! +! Build BSE matrix for spin-conserving transitions ! +!--------------------------------------------------! + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,1)*rho(a,j,kc,1)*Omega(kc)/eps + enddo + + B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI_aaaa(i,j,b,a) + 2d0*lambda*chi + + enddo + enddo + enddo + enddo + + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,2)*Omega(kc)/eps + enddo + + B_lr(nSa+ia,nSa+jb) = B_lr(nSa+ia,nSa+jb) - lambda*ERI_bbbb(i,j,b,a) + 2d0*lambda*chi + + enddo + enddo + enddo + enddo + + end if + + +!--------------------------------------------! +! Build BSE matrix for spin-flip transitions ! +!--------------------------------------------! + + if(ispin == 2) then + + ! abba block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,1)*rho(a,j,kc,2)*Omega(kc)/eps + enddo + + B_lr(ia,nSa+jb) = B_lr(ia,nSa+jb) - lambda*ERI_abab(i,a,b,j) + 2d0*lambda*chi + + end do + end do + end do + end do + + ! baab block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + eps = Omega(kc)**2 + eta**2 + chi = chi + rho(i,b,kc,2)*rho(a,j,kc,1)*Omega(kc)/eps + enddo + + B_lr(nSa+ia,jb) = B_lr(nSa+ia,jb) - lambda*ERI_abab(b,j,i,a) + 2d0*lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine unrestricted_Bethe_Salpeter_B_matrix diff --git a/src/QuAcK/unrestricted_QP_graph.f90 b/src/QuAcK/unrestricted_QP_graph.f90 new file mode 100644 index 0000000..96585a9 --- /dev/null +++ b/src/QuAcK/unrestricted_QP_graph.f90 @@ -0,0 +1,83 @@ +subroutine unrestricted_QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) + +! Compute the graphical solution of the QP equation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + + double precision,intent(in) :: eGWlin(nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 10 + double precision,parameter :: thresh = 1d-6 + double precision,external :: USigmaC,dUSigmaC + double precision :: sig,dsig + double precision :: f,df + double precision :: w + +! Output variables + + double precision,intent(out) :: eGW(nBas) + +! Run Newton's algorithm to find the root + + do p=nC+1,nBas-nR + + write(*,*) '-----------------' + write(*,'(A10,I3)') 'Orbital ',p + write(*,*) '-----------------' + + w = eGWlin(p) + nIt = 0 + f = 1d0 + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + sig = USigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + dsig = dUSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + f = w - eHF(p) - sig + df = 1d0 - dsig + + w = w - f/df + + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig + + + end do + + if(nIt == maxIt) then + + write(*,*) 'Newton root search has not converged!' + eGW(p) = eGWlin(p) + + else + + eGW(p) = w + + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV + write(*,*) + + end if + + end do + +end subroutine unrestricted_QP_graph diff --git a/src/QuAcK/unrestricted_excitation_density.f90 b/src/QuAcK/unrestricted_excitation_density.f90 index a705278..5822d65 100644 --- a/src/QuAcK/unrestricted_excitation_density.f90 +++ b/src/QuAcK/unrestricted_excitation_density.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY_a,XpY_b,rho) +subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho) ! Compute excitation densities for unrestricted reference @@ -14,11 +14,10 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ integer,intent(in) :: nSa integer,intent(in) :: nSb integer,intent(in) :: nSt - double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas) - double precision,intent(in) :: XpY_a(nSa,nSa) - double precision,intent(in) :: XpY_b(nSb,nSb) + 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) :: XpY(nSt,nSt) ! Local variables @@ -32,34 +31,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ rho(:,:,:,:) = 0d0 -!------------- -! alpha block -!------------- +!-------------! +! alpha block ! +!-------------! do p=nC(1)+1,nBas-nR(1) do q=nC(1)+1,nBas-nR(1) ! Same-spin contribution - do ia=1,nSa + do ia=1,nSt jb = 0 do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY_a(ia,jb) + rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aaaa(p,j,q,b)*XpY(ia,jb) enddo enddo enddo ! Opposite-spin contribution - do ia=1,nSb - jb = 0 + do ia=1,nSt + jb = nSa do j=nC(2)+1,nO(2) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - rho(p,q,nSa+ia,1) = rho(p,q,nSa+ia,1) + ERI_ab(p,j,q,b)*XpY_b(ia,jb) + rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aabb(p,j,q,b)*XpY(ia,jb) enddo enddo @@ -68,34 +67,34 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ enddo enddo -!------------ -! Beta block -!------------ +!------------! +! Beta block ! +!------------! do p=nC(2)+1,nBas-nR(2) do q=nC(2)+1,nBas-nR(2) - ! Same-spin contribution - do ia=1,nSb - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - - rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bb(p,j,q,b)*XpY_b(ia,jb) - - enddo - enddo - enddo - ! Opposite-spin contribution - do ia=1,nSa + do ia=1,nSt jb = 0 do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - rho(p,q,nSb+ia,2) = rho(p,q,nSb+ia,2) + ERI_ab(j,p,b,q)*XpY_a(ia,jb) + rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(j,p,b,q)*XpY(ia,jb) + + enddo + enddo + enddo + + ! Same-spin contribution + do ia=1,nSt + jb = nSa + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bbbb(p,j,q,b)*XpY(ia,jb) enddo enddo diff --git a/src/QuAcK/unrestricted_linear_response.f90 b/src/QuAcK/unrestricted_linear_response.f90 new file mode 100644 index 0000000..e816162 --- /dev/null +++ b/src/QuAcK/unrestricted_linear_response.f90 @@ -0,0 +1,123 @@ +subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,EcRPA,Omega,XpY,XmY) + +! Compute linear response for unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dRPA + logical,intent(in) :: TDA + logical,intent(in) :: BSE + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + 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_sc + double precision,intent(in) :: lambda + double precision,intent(in) :: e(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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: OmRPA(nS_sc) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) + +! Local variables + + integer :: ia + double precision :: trace_matrix + double precision,allocatable :: A(:,:) + double precision,allocatable :: B(:,:) + double precision,allocatable :: ApB(:,:) + double precision,allocatable :: AmB(:,:) + double precision,allocatable :: AmBSq(:,:) + double precision,allocatable :: AmBIv(:,:) + double precision,allocatable :: Z(:,:) + +! Output variables + + double precision,intent(out) :: EcRPA + double precision,intent(out) :: Omega(nSt) + double precision,intent(out) :: XpY(nSt,nSt) + double precision,intent(out) :: XmY(nSt,nSt) + +! Memory allocation + + allocate(A(nSt,nSt),B(nSt,nSt),ApB(nSt,nSt),AmB(nSt,nSt),AmBSq(nSt,nSt),AmBIv(nSt,nSt),Z(nSt,nSt)) + +! Build A and B matrices + + call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A) + + if(BSE) & + call unrestricted_Bethe_Salpeter_A_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,A) + +! Tamm-Dancoff approximation + + B = 0d0 + if(.not. TDA) then + + call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B) + + if(BSE) & + 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,ERI_abab,OmRPA,rho_RPA,B) + + end if + +! Build A + B and A - B matrices + + ApB = A + B + AmB = A - B + +! Diagonalize linear response matrix + + 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!!') + + 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!!') + + do ia=1,nSt + if(Omega(ia) < 0d0) Omega(ia) = 0d0 + end do + + 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) + +! Compute the RPA correlation energy + + EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nSt,A)) + +end subroutine unrestricted_linear_response diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 new file mode 100644 index 0000000..2f2d1b8 --- /dev/null +++ b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 @@ -0,0 +1,175 @@ +subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_lr) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + 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 + double precision,intent(in) :: lambda + double precision,intent(in) :: e(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,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,ia,jb + +! Output variables + + double precision,intent(out) :: A_lr(nSt,nSt) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +!----------------------------------------------- +! Build A matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + A_lr(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + lambda*ERI_aaaa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,b,j,a) + + end do + end do + end do + end do + + ! aabb block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + A_lr(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j) + + end do + end do + end do + end do + + ! bbaa block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + A_lr(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a) + + end do + end do + end do + end do + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + A_lr(nSa+ia,nSa+jb) = (e(a,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + lambda*ERI_bbbb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,b,j,a) + + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build A matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + A_lr(:,:) = 0d0 + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) + + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a) + + end do + end do + end do + end do + + end if + + +end subroutine unrestricted_linear_response_A_matrix diff --git a/src/QuAcK/unrestricted_linear_response_B_matrix.f90 b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 new file mode 100644 index 0000000..5cf001c --- /dev/null +++ b/src/QuAcK/unrestricted_linear_response_B_matrix.f90 @@ -0,0 +1,170 @@ +subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B_lr) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + 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 + double precision,intent(in) :: lambda + 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) :: ERI_abab(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,ia,jb + +! Output variables + + double precision,intent(out) :: B_lr(nSt,nSt) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +!----------------------------------------------- +! Build B matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + B_lr(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a) + + end do + end do + end do + end do + + ! aabb block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + B_lr(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b) + + end do + end do + end do + end do + + ! bbaa block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + B_lr(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a) + + end do + end do + end do + end do + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + B_lr(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a) + + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build B matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + B_lr(:,:) = 0d0 + + ! abba block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + B_lr(ia,nSa+jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(i,a,b,j) + + end do + end do + end do + end do + + ! baab block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + B_lr(nSa+ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a) + + end do + end do + end do + end do + + end if + + +end subroutine unrestricted_linear_response_B_matrix diff --git a/src/QuAcK/unrestricted_renormalization_factor.f90 b/src/QuAcK/unrestricted_renormalization_factor.f90 new file mode 100644 index 0000000..f175d93 --- /dev/null +++ b/src/QuAcK/unrestricted_renormalization_factor.f90 @@ -0,0 +1,90 @@ +subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z) + +! Compute the renormalization factor in the unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,j,a,b,p,q,jb + double precision :: eps + +! Output variables + + double precision,intent(out) :: Z(nBas,nspin) + +! Initialize + + Z(:,:) = 0d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the renormalization factor + + do p=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do jb=1,nSt + eps = e(p,1) - e(i,1) + Omega(jb) + Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(p,1) - e(a,1) - Omega(jb) + Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do jb=1,nSt + eps = e(p,2) - e(i,2) + Omega(jb) + Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(p,2) - e(a,2) - Omega(jb) + Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + +! Final rescaling + + Z(:,:) = 1d0/(1d0 + Z(:,:)) + +end subroutine unrestricted_renormalization_factor diff --git a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 b/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 index 3c91767..fa4084c 100644 --- a/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 +++ b/src/QuAcK/unrestricted_self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,SigC) +subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC) ! Compute diagonal of the correlation part of the self-energy @@ -13,8 +13,6 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS 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 double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: Omega(nSt) @@ -33,9 +31,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS SigC(:,:) = 0d0 -!-------------- -! Spin-up part -!-------------- +!--------------! +! Spin-up part ! +!--------------! ! Occupied part of the correlation self-energy @@ -59,9 +57,9 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSa,nS end do end do -!---------------- -! Spin-down part -!---------------- +!----------------! +! Spin-down part ! +!----------------! ! Occupied part of the correlation self-energy diff --git a/src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 b/src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 new file mode 100644 index 0000000..b48fbe4 --- /dev/null +++ b/src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 @@ -0,0 +1,30 @@ +subroutine unrestricted_spatial_to_spin_MO_energy(nBas,e,nBas2,se) + +! Convert MO energies from unrestricted spatial to spin orbitals + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nBas2 + double precision,intent(in) :: e(nBas,nspin) + +! Local variables + + integer :: p + +! Output variables + + double precision,intent(out) :: se(nBas2) + + do p=1,nBas2,2 + se(p) = e(p,1) + enddo + + do p=2,nBas2,2 + se(p) = e(p,2) + enddo + +end subroutine unrestricted_spatial_to_spin_MO_energy diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index 719df20..02c91a7 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -37,6 +37,10 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) open(unit=10,file='int/Nuc.dat') open(unit=11,file='int/ERI.dat') + open(unit=21,file='int/x.dat') + open(unit=22,file='int/y.dat') + open(unit=23,file='int/z.dat') + ! Read overlap integrals S(:,:) = 0d0