Add transistion vectors in unrestricted MOM rpas

lburth · lburth · commit 1380207673c1 · 2026-02-27T17:52:28.000+01:00
diff --git a/src/LR/CVS_phULR_oscillator_strength.f90 b/src/LR/CVS_phULR_oscillator_strength.f90
@@ -0,0 +1,88 @@
+subroutine CVS_phULR_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nCVS,nFC,maxS,occupations,virtuals,dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os)
+
+! Compute linear response
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nS(nspin)
+  integer,intent(in)            :: nSa
+  integer,intent(in)            :: nSb
+  integer,intent(in)            :: nSt
+  integer,intent(in)            :: maxS
+  integer,intent(in)            :: nCVS(nspin),nFC(nspin)
+  integer,intent(in)            :: occupations(maxval(nO),nspin)
+  integer,intent(in)            :: virtuals(nBas - minval(nO),nspin)
+  double precision              :: dipole_int_aa(nBas,nBas,ncart)
+  double precision              :: dipole_int_bb(nBas,nBas,ncart)
+  double precision,intent(in)   :: Omega(nSt)
+  double precision,intent(in)   :: XpY(nSt,nSt)
+  double precision,intent(in)   :: XmY(nSt,nSt)
+
+! Local variables
+
+  integer                       :: ia,jb,i,j,a,b
+  integer                       :: ixyz
+
+  double precision,allocatable  :: f(:,:)
+
+! Output variables
+
+  double precision              :: os(maxS)
+
+! Memory allocation
+
+  allocate(f(maxS,ncart))
+
+! Initialization
+   
+  f(:,:) = 0d0
+
+! Compute dipole moments and oscillator strengths
+
+  do ia=1,maxS
+    do ixyz=1,ncart
+
+      jb = 0
+      do j=1,nO(1) - nFC(1)
+        do b=nCVS(1) +1,nBas-nO(1)
+          jb = jb + 1
+          f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(occupations(j,1),virtuals(b,1),ixyz)*XpY(ia,jb)
+        end do
+      end do
+
+      jb = 0
+      do j=1,nO(2) - nFC(2)
+        do b=nCVS(2)+1,nBas-nO(2)
+          jb = jb + 1
+          f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(occupations(j,2),virtuals(b,2),ixyz)*XpY(ia,nSa+jb)
+        end do
+      end do
+
+    end do
+  end do
+
+  do ia=1,maxS
+    os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
+  end do
+    
+  write(*,*) '---------------------------------------------------------------'
+  write(*,*) '                Transition dipole moment (au)                  '
+  write(*,*) '---------------------------------------------------------------'
+  write(*,'(A3,5A12)') '#','X','Y','Z','dip. str.','osc. str.'
+  write(*,*) '---------------------------------------------------------------'
+  do ia=1,maxS
+    write(*,'(I3,5F12.6)') ia,(f(ia,ixyz),ixyz=1,ncart),sum(f(ia,:)**2),os(ia)
+  end do
+  write(*,*) '---------------------------------------------------------------'
+  write(*,*)
+
+end subroutine 
diff --git a/src/LR/CVS_phULR_transition_vectors.f90 b/src/LR/CVS_phULR_transition_vectors.f90
@@ -1,5 +1,5 @@
 subroutine CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,&
-                nCVS,occupations,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
+                nCVS,nFC,occupations,virtuals,dipole_int_aa,dipole_int_bb,c,S,Om,XpY,XmY)
 
 ! Print transition vectors for linear response calculation
 
@@ -19,6 +19,7 @@ subroutine CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,&
   integer,intent(in)            :: nSb
   integer,intent(in)            :: nSt
   integer,intent(in)            :: nCVS(nspin)
+  integer,intent(in)            :: nFC(nspin)
   integer,intent(in)            :: occupations(maxval(nO),nspin)
   integer,intent(in)            :: virtuals(nBas - minval(nO),nspin)
   double precision              :: dipole_int_aa(nBas,nBas,ncart)
@@ -32,142 +33,141 @@ subroutine CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,&
 ! Local variables
 
   integer                       :: ia,jb,j,b
-  integer                       :: maxS = 30
+  integer                       :: maxS = 10
   double precision,parameter    :: thres_vec = 0.1d0
   double precision,allocatable  :: X(:)
   double precision,allocatable  :: Y(:)
   double precision,allocatable  :: os(:)
   double precision,allocatable  :: S2(:)
 
 ! Memory allocation
-print *, "Transistion vectors not yet implemented for CVS !"
-!  maxS = min(nSt,maxS)
-!  allocate(X(nSt),Y(nSt),os(maxS),S2(maxS))
-!
-!! Compute oscillator strengths
-!
-!  os(:) = 0d0
-!  if(ispin == 1) call phULR_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS, & 
-!                                                dipole_int_aa,dipole_int_bb,Om,XpY,XmY,os)
-!
+  maxS = min(nSt,maxS)
+  allocate(X(nSt),Y(nSt),os(maxS),S2(maxS))
+
+! Compute oscillator strengths
+
+  os(:) = 0d0
+  if(ispin == 1) call CVS_phULR_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nCVS,nFC,maxS, & 
+                                                occupations,virtuals,dipole_int_aa,dipole_int_bb,Om,XpY,XmY,os)
+
 !! Compute <S**2>
 !
 !  call S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Om,XpY,XmY,S2)
-!
-!! Print details about spin-conserved excitations
-!
-!  if(ispin == 1) then
-!
-!    do ia=1,maxS
-!
-!      X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
-!      Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
-!
-!      print*,'-------------------------------------------------------------'
-!      write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & 
-!              ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV','  f = ',os(ia),'  <S**2> = ',S2(ia)
-!      print*,'-------------------------------------------------------------'
-!
-!      ! Spin-up transitions
-!
-!      jb = 0
-!      do j=nC(1)+1,nO(1)
-!        do b=nO(1)+1,nBas-nR(1)
-!          jb = jb + 1
-!          if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)
-!        end do
-!      end do
-!   
-!      jb = 0
-!      do j=nC(1)+1,nO(1)
-!        do b=nO(1)+1,nBas-nR(1)
-!          jb = jb + 1
-!          if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)
-!        end do
-!      end do
-!
-!      ! Spin-down transitions
-!
-!      jb = 0
-!      do j=nC(2)+1,nO(2)
-!        do b=nO(2)+1,nBas-nR(2)
-!          jb = jb + 1
-!          if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(nSa+jb)
-!        end do
-!      end do
-!   
-!      jb = 0
-!      do j=nC(2)+1,nO(2)
-!        do b=nO(2)+1,nBas-nR(2)
-!          jb = jb + 1
-!          if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(nSa+jb)
-!        end do
-!      end do
-!     write(*,*)
-!
-!    end do
-!
-!  end if
-!
-!! Print details about spin-flip excitations
-!
-!  if(ispin == 2) then
-!
-!    do ia=1,maxS
-!
-!      X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
-!      Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
-!
-!
-!      print*,'-------------------------------------------------------------'
-!      write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & 
-!              ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV','  f = ',os(ia),'  <S**2> = ',S2(ia)
-!      print*,'-------------------------------------------------------------'
-!
-!      ! Spin-up transitions
-!
-!      jb = 0
-!      do j=nC(1)+1,nO(1)
-!        do b=nO(2)+1,nBas-nR(2)
-!          jb = jb + 1
-!          if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'B = ',X(jb)
-!        end do
-!      end do
-!   
-!      jb = 0
-!      do j=nC(1)+1,nO(1)
-!        do b=nO(2)+1,nBas-nR(2)
-!          jb = jb + 1
-!          if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'B = ',Y(jb)
-!        end do
-!      end do
-!
-!      ! Spin-down transitions
-!
-!      jb = 0
-!      do j=nC(2)+1,nO(2)
-!        do b=nO(1)+1,nBas-nR(1)
-!          jb = jb + 1
-!          if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'B = ',X(nSa+jb)
-!        end do
-!      end do
-!   
-!      jb = 0
-!      do j=nC(2)+1,nO(2)
-!        do b=nO(1)+1,nBas-nR(1)
-!          jb = jb + 1
-!          if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'B = ',Y(nSa+jb)
-!        end do
-!      end do
-!     write(*,*)
-!
-!    end do
-!
-!  end if
-!
-!! Thomas-Reiche-Kuhn sum rule
-!
-!  write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
-!  write(*,*)
-!
+
+! Print details about spin-conserved excitations
+
+  if(ispin == 1) then
+
+    do ia=1,maxS
+
+      X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
+      Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
+
+      print*,'-------------------------------------------------------------'
+      write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & 
+              ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV','  f = ',os(ia)!,'  <S**2> = ',S2(ia)
+      print*,'-------------------------------------------------------------'
+
+      ! Spin-up transitions
+
+      jb = 0
+      do j=1,nO(1)-nFC(1)
+        do b=nCVS(1)+1,nBas-nO(1)
+          jb = jb + 1
+          if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,1),'A -> ',virtuals(b,1),'A = ',X(jb)
+        end do
+      end do
+   
+      jb = 0
+      do j=1,nO(1)-nFC(1)
+        do b=nCVS(1)+1,nBas-nO(1)
+          jb = jb + 1
+          if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,1),'A <- ',virtuals(b,1),'A = ',Y(jb)
+        end do
+      end do
+
+      ! Spin-down transitions
+
+      jb = 0
+      do j=1,nO(2) - nFC(2)
+        do b=nCVS(2)+1,nBas-nO(2)
+          jb = jb + 1
+          if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,2),'B -> ',virtuals(b,2),'B = ',X(nSa+jb)
+        end do
+      end do
+   
+      jb = 0
+      do j=1,nO(2)-nFC(2)
+        do b=nCVS(2)+1,nBas-nO(2)
+          jb = jb + 1
+          if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,2),'B <- ',virtuals(b,2),'B = ',Y(nSa+jb)
+        end do
+      end do
+     write(*,*)
+
+    end do
+
+  end if
+
+! Print details about spin-flip excitations
+
+  if(ispin == 2) then
+
+    do ia=1,maxS
+
+      X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
+      Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
+
+
+      print*,'-------------------------------------------------------------'
+      write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & 
+              ' Excitation n. ',ia,': ',Om(ia)*HaToeV,' eV','  f = ',os(ia)!,'  <S**2> = ',S2(ia)
+      print*,'-------------------------------------------------------------'
+
+      ! Spin-up transitions
+
+      jb = 0
+      do j=1,nO(1)-nFC(1)
+        do b=nCVS(2)+1,nBas-nO(2)
+          jb = jb + 1
+          if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,1),'A -> ',virtuals(b,2),'B = ',X(jb)
+        end do
+      end do
+   
+      jb = 0
+      do j=1,nO(1)-nFC(1)
+        do b=nCVS(2)+1,nBas-nO(2)
+          jb = jb + 1
+          if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,1),'A <- ',virtuals(b,2),'B = ',Y(jb)
+        end do
+      end do
+
+      ! Spin-down transitions
+
+      jb = 0
+      do j=1,nO(2)-nFC(2)
+        do b=nCVS(1)+1,nBas-nO(1)
+          jb = jb + 1
+          if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,2),'A -> ',virtuals(b,1),'B = ',X(nSa+jb)
+        end do
+      end do
+   
+      jb = 0
+      do j=1,nO(2)-nFC(2)
+        do b=nCVS(1)+1,nBas-nO(1)
+          jb = jb + 1
+          if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') occupations(j,2),'A <- ',virtuals(b,1),'B = ',Y(nSa+jb)
+        end do
+      end do
+     write(*,*)
+
+    end do
+
+  end if
+
+! Thomas-Reiche-Kuhn sum rule
+
+  write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:))
+  write(*,*)
+
 end subroutine 
diff --git a/src/LR/complex_phLR_transition_vectors.f90 b/src/LR/complex_phLR_transition_vectors.f90
@@ -26,7 +26,7 @@ subroutine complex_phLR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,&
 ! Local variables
 
   integer                       :: ia,jb,j,b
-  integer                       :: maxS = 30
+  integer                       :: maxS = 10
   double precision,parameter    :: thres_vec = 0.1d0
   complex*16,allocatable        :: X(:)
   complex*16,allocatable        :: Y(:)
diff --git a/src/LR/complex_phULR_transition_vectors.f90 b/src/LR/complex_phULR_transition_vectors.f90
@@ -33,7 +33,7 @@ subroutine complex_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nS
 ! Local variables
 
   integer                       :: ia,jb,j,b
-  integer                       :: maxS = 30
+  integer                       :: maxS = 10
   double precision,parameter    :: thres_vec = 0.1d0
   complex*16,allocatable        :: X(:)
   complex*16,allocatable        :: Y(:)
diff --git a/src/RPA/CVS_phURPA.f90 b/src/RPA/CVS_phURPA.f90
diff --git a/src/RPA/CVS_phURPAx.f90 b/src/RPA/CVS_phURPAx.f90
