From 0e9edb30d28120c28c8462538f8d850bb6da229a Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 28 Sep 2020 21:25:25 +0200 Subject: [PATCH] oscillator strength and dipole integrals --- examples/basis.Ne.cc-pvdz | 43 +++++----- examples/molecule.H2 | 4 +- input/basis | 9 ++ input/methods | 6 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/options | 4 +- src/QuAcK/BSE2.f90 | 30 +++---- src/QuAcK/BSE2_dynamic_perturbation.f90 | 5 +- .../BSE2_dynamic_perturbation_iterative.f90 | 6 +- src/QuAcK/Bethe_Salpeter.f90 | 15 ++-- .../Bethe_Salpeter_dynamic_perturbation.f90 | 6 +- ...alpeter_dynamic_perturbation_iterative.f90 | 6 +- src/QuAcK/CIS.f90 | 3 +- src/QuAcK/G0F2.f90 | 5 +- src/QuAcK/G0T0.f90 | 21 ++--- src/QuAcK/G0W0.f90 | 5 +- src/QuAcK/QuAcK.f90 | 85 ++++++++++++------- src/QuAcK/RPAx.f90 | 8 +- src/QuAcK/UCIS.f90 | 3 +- src/QuAcK/UG0W0.f90 | 3 +- src/QuAcK/URPAx.f90 | 7 +- src/QuAcK/UdRPA.f90 | 3 +- src/QuAcK/dRPA.f90 | 7 +- src/QuAcK/evGF2.f90 | 5 +- src/QuAcK/evGT.f90 | 19 ++--- src/QuAcK/evGW.f90 | 5 +- src/QuAcK/evUGW.f90 | 3 +- src/QuAcK/print_transition_vectors.f90 | 72 +++++++++++----- src/QuAcK/qsGW.f90 | 6 +- src/QuAcK/spatial_to_spin_MO_energy.f90 | 2 +- ...unrestricted_spatial_to_spin_MO_energy.f90 | 30 +++++++ src/utils/read_integrals.f90 | 4 + 33 files changed, 283 insertions(+), 157 deletions(-) create mode 100644 src/QuAcK/unrestricted_spatial_to_spin_MO_energy.f90 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 779d849..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.399 + H 0. 0. -0.7 + H 0. 0. 0.7 diff --git a/input/basis b/input/basis index 79a747a..fb05e68 100644 --- a/input/basis +++ b/input/basis @@ -7,3 +7,12 @@ 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 58c5d33..0d456f5 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF MOM - F T F -# MP2* MP3 MP2-F12 T F F +# MP2* MP3 MP2-F12 + F F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index fd4bfbe..7225285 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 1 0 0 0 + 2 1 1 0 0 # Znuc x y z - H 0. 0. 0. + H 0. 0. -0.7 + H 0. 0. 0.7 diff --git a/input/molecule.xyz b/input/molecule.xyz index 3dca9a4..ac9420c 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,4 @@ - 1 + 2 - H 0.0000000000 0.0000000000 0.0000000000 + H 0.0000000000 0.0000000000 -0.3704240743 + H 0.0000000000 0.0000000000 0.3704240743 diff --git a/input/options b/input/options index 0d3bf88..1d00d01 100644 --- a/input/options +++ b/input/options @@ -1,4 +1,4 @@ -# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type +# HF: maxSCF thresh DIIS n_diis guess_type ortho_type 64 0.00001 T 5 1 1 # MP: @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - F F T T + F T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/BSE2.f90 b/src/QuAcK/BSE2.f90 index 64a9e8c..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,14 +46,14 @@ 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(:,:,:,:), & + 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)) @@ -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 @@ -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 d0763ee..8ca6757 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,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 @@ -20,6 +20,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, double precision,intent(in) :: eW(nBas) double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -70,6 +71,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, 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 @@ -81,11 +84,11 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, & + 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,rho_RPA, & + 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 @@ -107,6 +110,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, 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 @@ -118,11 +123,11 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA, & + 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,rho_RPA, & + 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 diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index ba0d2dd..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,rho_RPA,OmBSE,XpY,XmY) +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,6 +18,7 @@ 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) @@ -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(*,*) diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation_iterative.f90 index 72b5e88..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,rho_RPA,OmBSE,XpY,XmY) +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,6 +18,7 @@ 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) @@ -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(*,*) 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 a01c21b..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,triplet,linearize,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,PHF,cHF,eHF,eGW) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,dipole_int,PHF,cHF,eHF,eGW) ! Perform G0W0 calculation @@ -34,6 +34,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -182,7 +183,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 2efdaa2..9c5c095 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -52,9 +52,17 @@ 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 :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) + integer :: ixyz + integer :: ispin integer :: bra1,bra2 integer :: ket1,ket2 double precision,allocatable :: ERI_MO_aaaa(:,:,:,:) @@ -225,7 +233,8 @@ program QuAcK ! Memory allocation for one- and two-electron integrals allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), & - S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas)) + S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), & + dipole_int(nBas,nBas,ncart,nspin),ERI_AO(nBas,nBas,nBas,nBas)) ! Read integrals @@ -319,7 +328,6 @@ program QuAcK write(*,*) 'AO to MO transformation... Please be patient' write(*,*) - if(doSph) then allocate(ERI_MO(nBas,nBas,nBas,nBas)) @@ -331,6 +339,15 @@ program QuAcK if(unrestricted) then + ! Read and transform dipole-related integrals + + call read_dipole_integrals(nBas,dipole_int) + do ixyz=1,ncart + do ispin=1,nspin + call AOtoMO_transform(nBas,cHF(:,:,ispin),dipole_int(:,:,ixyz,ispin)) + end do + end do + ! Memory allocation allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas)) @@ -379,7 +396,15 @@ program QuAcK ! Memory allocation allocate(ERI_MO(nBas,nBas,nBas,nBas)) - + + ! Read and transform dipole-related integrals + + ispin = 1 + call read_dipole_integrals(nBas,dipole_int) + do ixyz=1,ncart + call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz,ispin)) + end do + ! 4-index transform bra1 = 1 @@ -588,11 +613,12 @@ program QuAcK call cpu_time(start_CIS) 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,eHF) + 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,eHF) + call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF) end if call cpu_time(end_CIS) @@ -626,7 +652,7 @@ program QuAcK if(doCISD) then call cpu_time(start_CISD) - call CISD(singlet,triplet,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 @@ -645,11 +671,11 @@ program QuAcK 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,eHF) + 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,eHF) + 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) @@ -670,11 +696,11 @@ program QuAcK 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,eHF) + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF) else - call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + 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) @@ -692,8 +718,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_ppRPA) - call ppRPA(singlet,triplet, & - 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 @@ -727,7 +752,7 @@ program QuAcK call cpu_time(start_GF2) call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & - eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + 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 @@ -744,8 +769,8 @@ program QuAcK call cpu_time(start_GF2) call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, & - singlet,triplet,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 @@ -799,11 +824,11 @@ program QuAcK 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,PHF,cHF,eHF,eG0W0) + ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,PHF,cHF,eHF,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0) + linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0) end if @@ -826,13 +851,13 @@ program QuAcK call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & - ERHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,PHF,cHF,eHF,eG0W0) + ERHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,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,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0) end if call cpu_time(end_evGW) @@ -851,7 +876,7 @@ program QuAcK call cpu_time(start_qsGW) call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) + 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 @@ -869,9 +894,9 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & + 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,eHF,eG0T0) + 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 @@ -887,9 +912,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,triplet,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 @@ -898,10 +923,6 @@ program QuAcK end if -!------------------------------------------------------------------------ -! Perform evGT calculatiom -!------------------------------------------------------------------------ - !------------------------------------------------------------------------ ! Information for Monte Carlo calculations !------------------------------------------------------------------------ @@ -1006,7 +1027,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) + 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 @@ -1020,7 +1041,7 @@ program QuAcK call cpu_time(start_G0T0) 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_ERF_MO,eHF,eG0T0) + 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 37f3f25..3e523b4 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -1,4 +1,5 @@ -subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +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) @@ -24,6 +25,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR 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 @@ -62,7 +64,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR 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@HF ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -75,7 +77,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR 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@HF ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif diff --git a/src/QuAcK/UCIS.f90 b/src/QuAcK/UCIS.f90 index 9db0626..e75fc5a 100644 --- a/src/QuAcK/UCIS.f90 +++ b/src/QuAcK/UCIS.f90 @@ -1,4 +1,4 @@ -subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,eHF) +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` @@ -20,6 +20,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E 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 diff --git a/src/QuAcK/UG0W0.f90 b/src/QuAcK/UG0W0.f90 index 5e1b6fc..dd5fa5b 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, & spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS, & - ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,PHF,cHF,eHF,eGW) + ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,PHF,cHF,eHF,eGW) ! Perform unrestricted G0W0 calculation @@ -41,6 +41,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev 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 diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 index ff13fad..9a3c510 100644 --- a/src/QuAcK/URPAx.f90 +++ b/src/QuAcK/URPAx.f90 @@ -1,5 +1,5 @@ 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,e) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,e) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -28,6 +28,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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 @@ -77,7 +78,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) +! call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) deallocate(Omega_sc,XpY_sc,XmY_sc) @@ -100,7 +101,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) +! 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) diff --git a/src/QuAcK/UdRPA.f90 b/src/QuAcK/UdRPA.f90 index 0169441..5d1e7f6 100644 --- a/src/QuAcK/UdRPA.f90 +++ b/src/QuAcK/UdRPA.f90 @@ -1,5 +1,5 @@ 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,e) + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,e) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -28,6 +28,7 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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 diff --git a/src/QuAcK/dRPA.f90 b/src/QuAcK/dRPA.f90 index 3985af2..fa45771 100644 --- a/src/QuAcK/dRPA.f90 +++ b/src/QuAcK/dRPA.f90 @@ -1,4 +1,4 @@ -subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +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 @@ -24,6 +24,7 @@ subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR 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 @@ -62,7 +63,7 @@ subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR 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@HF ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -75,7 +76,7 @@ subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR 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@HF ',ispin,nS,Omega(:,ispin)) - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif 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 d7d816e..26061f8 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,6 +1,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI, & - PHF,cHF,eHF,eG0W0) + dipole_int,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -37,6 +37,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE double precision,intent(in) :: eG0W0(nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables @@ -235,7 +236,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,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 diff --git a/src/QuAcK/evUGW.f90 b/src/QuAcK/evUGW.f90 index 1ba36ac..24f8e0e 100644 --- a/src/QuAcK/evUGW.f90 +++ b/src/QuAcK/evUGW.f90 @@ -1,6 +1,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & - ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,PHF,cHF,eHF,eG0W0) + ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -46,6 +46,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE 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 diff --git a/src/QuAcK/print_transition_vectors.f90 b/src/QuAcK/print_transition_vectors.f90 index ab3e093..6986eec 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,12 +7,14 @@ 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) @@ -20,31 +22,66 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) ! Local variables 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(:,:) + + write(*,*) '------------------------' + write(*,*) ' Dipole moments (X Y Z) ' + write(*,*) '------------------------' + call matout(nS,ncart,f) + write(*,*) + + do ia=1,nS + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nS,1,os) + write(*,*) + + 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 +91,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 +102,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/qsGW.f90 b/src/QuAcK/qsGW.f90 index b1fbfe2..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, & 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,PHF,cHF,eHF) + Hc,ERI_AO_basis,ERI_MO_basis,dipole_int,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -41,6 +41,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE 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 @@ -267,7 +268,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,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 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_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