oscillator strength and dipole integrals

examples/basis.Ne.cc-pvdz

@@ -1,30 +1,29 @@
 1   6
 S   8
-  1  17880.0000000              0.0007380        
+  1  17880.0000000              0.0007380
-  2   2683.0000000              0.0056770        
+  2   2683.0000000              0.0056770
-  3    611.5000000              0.0288830        
+  3    611.5000000              0.0288830
-  4    173.5000000              0.1085400        
+  4    173.5000000              0.1085400
-  5     56.6400000              0.2909070        
+  5     56.6400000              0.2909070
-  6     20.4200000              0.4483240        
+  6     20.4200000              0.4483240
-  7      7.8100000              0.2580260        
+  7      7.8100000              0.2580260
-  8      1.6530000              0.0150630        
+  8      1.6530000              0.0150630
 S   8
-  1  17880.0000000             -0.0001720        
+  1  17880.0000000             -0.0001720
-  2   2683.0000000             -0.0013570        
+  2   2683.0000000             -0.0013570
-  3    611.5000000             -0.0067370        
+  3    611.5000000             -0.0067370
-  4    173.5000000             -0.0276630        
+  4    173.5000000             -0.0276630
-  5     56.6400000             -0.0762080        
+  5     56.6400000             -0.0762080
-  6     20.4200000             -0.1752270        
+  6     20.4200000             -0.1752270
-  7      7.8100000             -0.1070380        
+  7      7.8100000             -0.1070380
-  8      1.6530000              0.5670500        
+  8      1.6530000              0.5670500
 S   1
-  1      0.4869000              1.0000000        
+  1      0.4869000              1.0000000
 P   3
-  1     28.3900000              0.0460870        
+  1     28.3900000              0.0460870
-  2      6.2700000              0.2401810        
+  2      6.2700000              0.2401810
-  3      1.6950000              0.5087440        
+  3      1.6950000              0.5087440
 P   1
-  1      0.4317000              1.0000000        
+  1      0.4317000              1.0000000
 D   1
   1      2.2020000              1.0000000
-

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.        -0.7
-  H     0.           0.         1.399
+  H     0.           0.         0.7

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

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

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

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

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

src/QuAcK/BSE2.f90

@@ -1,5 +1,4 @@
-subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & 
+subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE)
-                eta,nBas,nC,nO,nV,nR,nS,ERI,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)            :: singlet
-  logical,intent(in)            :: triplet_manifold
+  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, & 
+      call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                               ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+                                               OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
      else
  
-      call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & 
+      call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                     ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+                                     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, & 
+      call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & 
-                                               ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+                                               OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
      else
  
-      call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & 
+      call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, &
-                                     ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+                                     OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
  
       end if
 

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) 
 

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(*,*)

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
 

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(*,*)

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(*,*)

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
 

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
 

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, & 
+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,eHF,eG0T0)
+                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)            :: singlet
-  logical,intent(in)            :: triplet_manifold
+  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,triplet,eta, &
-
+                        nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eG0T0,EcBSE)
-    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)
 
     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, &
+      call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, &
-                 nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcAC)
+                 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

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
  

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, & 
+    call ppRPA(singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF)
-               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,                     & 
+               singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & 
-               eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
+               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,            &
+    call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, &
-              BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & 
+              BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,        & 
-              nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0)
+              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

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
 

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
 

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
 

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)
 

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
 

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
 

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
 

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,      & 
+                BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,      & 
-                eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0)
+                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)            :: singlet
-  logical,intent(in)            :: triplet_manifold
+  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,triplet,eta, &
-
+                        nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGT,eGT,EcRPA,EcBSE)
-    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)
 
     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, &
+      call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta, &
-                 nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcAC)
+                 nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,EcAC)
 
       if(exchange_kernel) then
 

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
 

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
 

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*,'--------------------------------'
+    print*,'---------------------------------------------'
-    write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV'
+    write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')'
-    print*,'--------------------------------'
+    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

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
 

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
 

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

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