transition vectors for restricted rpa on mom

lburth · lburth · commit 0c012dcc57fd · 2026-02-27T18:53:00.000+01:00
diff --git a/src/HF/MOM_RHF.f90 b/src/HF/MOM_RHF.f90
@@ -243,6 +243,7 @@ subroutine MOM_RHF(dotest,doaordm,maxSCF,thresh,max_diis,guess_type,level_shift,
 
   end if
 
+  call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
   call print_MOM_RHF(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,ERHF,dipole,occupations)
 
 ! Print the 1-RDM and 2-RDM in AO basis
diff --git a/src/LR/CVS_phLR_oscillator_strength.f90 b/src/LR/CVS_phLR_oscillator_strength.f90
@@ -0,0 +1,73 @@
+subroutine CVS_phLR_oscillator_strength(nOrb,nC,nO,nV,nR,nS,nCVS,nFC,maxS,occupations,virtuals,dipole_int,Om,XpY,XmY,os)
+
+! Compute oscillator strength from a ph linear response calculation
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nOrb
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  integer,intent(in)            :: nCVS,nFC
+  integer,intent(in)            :: occupations(nO-nFC),virtuals(nOrb-nO)
+  integer,intent(in)            :: maxS
+  double precision              :: dipole_int(nOrb,nOrb,ncart)
+  double precision,intent(in)   :: Om(nS)
+  double precision,intent(in)   :: XpY(nS,nS)
+  double precision,intent(in)   :: XmY(nS,nS)
+  
+! Local variables
+
+  integer                       :: m,jb,i,j,a,b
+  integer                       :: ixyz
+
+  double precision,allocatable  :: f(:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: os(nS)
+
+! Memory allocation
+
+  allocate(f(maxS,ncart))
+
+! Initialization
+   
+  f(:,:) = 0d0
+
+! Compute dipole moments and oscillator strengths
+
+  do m=1,maxS
+    do ixyz=1,ncart
+      jb = 0
+      do j=1,nO-nFC
+        do b=nCVS+1,nOrb-nO
+          jb = jb + 1
+          f(m,ixyz) = f(m,ixyz) + dipole_int(occupations(j),virtuals(b),ixyz)*XpY(m,jb)
+        end do
+      end do
+    end do
+  end do
+  f(:,:) = sqrt(2d0)*f(:,:)
+
+  do m=1,maxS
+    os(m) = 2d0/3d0*Om(m)*sum(f(m,:)**2)
+  end do
+ 
+  write(*,*) '---------------------------------------------------------------'
+  write(*,*) '                Transition dipole moment (au)                  '
+  write(*,*) '---------------------------------------------------------------'
+  write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.'
+  write(*,*) '---------------------------------------------------------------'
+  do m=1,maxS
+    write(*,'(I3,5F12.6)') m,(f(m,ixyz),ixyz=1,ncart),sum(f(m,:)**2),os(m)
+  end do
+  write(*,*) '---------------------------------------------------------------'
+  write(*,*)
+
+end subroutine 
diff --git a/src/LR/CVS_phLR_transition_vectors.f90 b/src/LR/CVS_phLR_transition_vectors.f90
@@ -1,4 +1,4 @@
-subroutine CVS_phLR_transition_vectors(spin_allowed,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY)
+subroutine CVS_phLR_transition_vectors(spin_allowed,nOrb,nC,nO,nV,nR,nS,nCVS,nFC,occupations,virtuals,dipole_int,Om,XpY,XmY)
 
 ! Print transition vectors for linear response calculation
 
@@ -14,6 +14,8 @@ subroutine CVS_phLR_transition_vectors(spin_allowed,nOrb,nC,nO,nV,nR,nS,dipole_i
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
+  integer,intent(in)            :: nCVS,nFC
+  integer,intent(in)            :: occupations(nO-nFC),virtuals(nOrb-nO)
   double precision              :: dipole_int(nOrb,nOrb,ncart)
   double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: XpY(nS,nS)
@@ -22,66 +24,65 @@ subroutine CVS_phLR_transition_vectors(spin_allowed,nOrb,nC,nO,nV,nR,nS,dipole_i
 ! Local variables
 
   integer                       :: ia,jb,j,b
-  integer                       :: maxS = 30
+  integer                       :: maxS = 10
   double precision              :: S2
   double precision,parameter    :: thres_vec = 0.1d0
   double precision,allocatable  :: X(:)
   double precision,allocatable  :: Y(:)
   double precision,allocatable  :: os(:)
 
-  print *, "Transistion vectors for this method not implemented yet, sry..."
-!! Memory allocation
-!
-!  maxS = min(nS,maxS)
-!  allocate(X(nS),Y(nS),os(maxS))
-!
-!! Compute oscillator strengths
-!
-!  os(:) = 0d0
-!  if(spin_allowed) call phLR_oscillator_strength(nOrb,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os)
-!
-!! Print details about excitations
-!
-!  do ia=1,maxS
-!
-!    X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
-!    Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
-!
-!    ! <S**2> values
-!
-!    if(spin_allowed) then 
-!      S2 = 0d0
-!    else
-!      S2 = 2d0
-!    end if
-!
-!    print*,'-------------------------------------------------------------'
-!    write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
-!            ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV','  f = ',os(ia),'  <S**2> = ',S2
-!    print*,'-------------------------------------------------------------'
-!
-!    jb = 0
-!    do j=nC+1,nO
-!      do b=nO+1,nOrb-nR
-!        jb = jb + 1
-!        if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0)
-!      end do
-!    end do
-! 
-!    jb = 0
-!    do j=nC+1,nO
-!      do b=nO+1,nOrb-nR
-!        jb = jb + 1
-!        if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0)
-!      end do
-!    end do
-!   write(*,*)
-!
-!  end do
-!
-!! Thomas-Reiche-Kuhn sum rule
-!
-!  write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
-!  write(*,*)
+! Memory allocation
+
+  maxS = min(nS,maxS)
+  allocate(X(nS),Y(nS),os(maxS))
+
+! Compute oscillator strengths
+
+  os(:) = 0d0
+  if(spin_allowed) call CVS_phLR_oscillator_strength(nOrb,nC,nO,nV,nR,nS,nCVS,nFC,maxS,occupations,virtuals,dipole_int,Om,XpY,XmY,os)
+
+! Print details about excitations
+
+  do ia=1,maxS
+
+    X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
+    Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
+
+    ! <S**2> values
+
+    if(spin_allowed) then 
+      S2 = 0d0
+    else
+      S2 = 2d0
+    end if
+
+    print*,'-------------------------------------------------------------'
+    write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
+            ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV','  f = ',os(ia),'  <S**2> = ',S2
+    print*,'-------------------------------------------------------------'
+
+    jb = 0
+    do j=1,nO-nFC
+      do b=nCVS+1,nOrb-nO
+        jb = jb + 1
+        if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') occupations(j),' -> ',virtuals(b),' = ',X(jb)/sqrt(2d0)
+      end do
+    end do
+ 
+    jb = 0
+    do j=1,nO-nFC
+      do b=nCVS+1,nOrb-nO
+        jb = jb + 1
+        if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') occupations(j),' <- ',virtuals(b),' = ',Y(jb)/sqrt(2d0)
+      end do
+    end do
+   write(*,*)
+
+  end do
+
+! Thomas-Reiche-Kuhn sum rule
+
+  write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
+  write(*,*)
 
 end subroutine 
diff --git a/src/RPA/CVS_phRRPA.f90 b/src/RPA/CVS_phRRPA.f90
@@ -121,7 +121,7 @@ subroutine CVS_phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC
 
     call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPA@RHF','singlet',nSCVS,Om)
-    call CVS_phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
+    call CVS_phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,dipole_int,Om,XpY,XmY)
 
   end if
 
@@ -136,7 +136,7 @@ subroutine CVS_phRRPA(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC
 
     call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPA@RHF','triplet',nSCVS,Om)
-    call CVS_phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
+    call CVS_phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,dipole_int,Om,XpY,XmY)
 
   end if
 
diff --git a/src/RPA/CVS_phRRPAx.f90 b/src/RPA/CVS_phRRPAx.f90
@@ -121,7 +121,7 @@ subroutine CVS_phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,n
 
     call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPAx@RHF','singlet',nSCVS,Om)
-    call CVS_phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
+    call CVS_phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,dipole_int,Om,XpY,XmY)
 
   end if
 
@@ -136,7 +136,7 @@ subroutine CVS_phRRPAx(dotest,TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,n
 
     call CVS_phRLR(TDA,nSCVS,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPAx@RHF','triplet',nSCVS,Om)
-    call CVS_phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nSCVS,dipole_int,Om,XpY,XmY)
+    call CVS_phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nSCVS,nCVS,nFC,occupations_fc,virtuals,dipole_int,Om,XpY,XmY)
 
   end if
 
diff --git a/src/RPA/CVS_phURPA.f90 b/src/RPA/CVS_phURPA.f90
@@ -140,7 +140,7 @@ subroutine CVS_phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fli
     call CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPA@UHF','spin-conserved',nSt,Om)
     call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nCVS,nFC, &
-                               occupations,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
+                               occupations_fc,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
 
@@ -168,7 +168,7 @@ subroutine CVS_phURPA(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fli
     call CVS_phULR(TDA,nSa,nSa,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPA@UHF','spin-flip',nSt,Om)
     call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nCVS,nFC, &
-                               occupations,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
+                               occupations_fc,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
 
diff --git a/src/RPA/CVS_phURPAx.f90 b/src/RPA/CVS_phURPAx.f90
@@ -140,7 +140,7 @@ subroutine CVS_phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fl
     call CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPAx@UHF','spin-conserved',nSt,Om)
     call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nCVS,nFC, &
-                                      occupations,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
+                                      occupations_fc,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
 
@@ -168,7 +168,7 @@ subroutine CVS_phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fl
     call CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA(ispin),Om,XpY,XmY)
     call print_excitation_energies('phRPAx@UHF','spin-flip',nSt,Om)
     call CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nCVS,nFC, &
-                                      occupations,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
+                                      occupations_fc,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
     deallocate(Aph,Bph,Om,XpY,XmY)
 
