diff --git a/input/methods b/input/methods index 3a1c5a4..61f2d19 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - F T F + T F F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) @@ -9,7 +9,7 @@ # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA - F F F + F T F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW diff --git a/src/LR/oscillator_strength.f90 b/src/LR/oscillator_strength.f90 new file mode 100644 index 0000000..5524169 --- /dev/null +++ b/src/LR/oscillator_strength.f90 @@ -0,0 +1,76 @@ +subroutine oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: XpY(nS,nS) + double precision,intent(in) :: XmY(nS,nS) + +! Local variables + + logical :: debug = .false. + integer :: ia,jb,i,j,a,b + integer :: ixyz + + double precision,allocatable :: f(:,:) + +! Output variables + + double precision :: os(nS) + +! Memory allocation + + allocate(f(nS,ncart)) + +! Initialization + + f(:,:) = 0d0 + +! Compute dipole moments and oscillator strengths + + do ia=1,nS + do ixyz=1,ncart + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz)*XpY(ia,jb) + end do + end do + end do + end do + f(:,:) = sqrt(2d0)*f(:,:) + + do ia=1,nS + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + if(debug) then + + write(*,*) '------------------------' + write(*,*) ' Dipole moments (X Y Z) ' + write(*,*) '------------------------' + call matout(nS,ncart,f) + write(*,*) + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nS,1,os) + write(*,*) + + end if + +end subroutine oscillator_strength diff --git a/src/LR/print_transition_vectors.f90 b/src/LR/print_transition_vectors.f90 index ef29043..c07e4a0 100644 --- a/src/LR/print_transition_vectors.f90 +++ b/src/LR/print_transition_vectors.f90 @@ -25,61 +25,21 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, integer :: ia,jb,i,j,a,b integer :: ixyz integer,parameter :: maxS = 10 - double precision :: norm + double precision :: S2 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),f(nS,ncart),os(nS)) + allocate(X(nS),Y(nS),os(nS)) -! Initialization - - f(:,:) = 0d0 - os(:) = 0d0 +! Compute oscillator strengths -! Compute dipole moments and oscillator strengths + os(:) = 0d0 + if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,dipole_int,Omega,XpY,XmY,os) - if(spin_allowed) then - - do ia=1,nS - do ixyz=1,ncart - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz)*XpY(ia,jb) - end do - end do - end do - end do - f(:,:) = sqrt(2d0)*f(:,:) - - do ia=1,nS - os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) - end do - - if(debug) then - - write(*,*) '------------------------' - write(*,*) ' Dipole moments (X Y Z) ' - write(*,*) '------------------------' - call matout(nS,ncart,f) - write(*,*) - - write(*,*) '----------------------' - write(*,*) ' Oscillator strengths ' - write(*,*) '----------------------' - call matout(nS,1,os) - write(*,*) - - end if - - end if - ! Print details about excitations do ia=1,min(nS,maxS) @@ -87,9 +47,10 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - print*,'---------------------------------------------' - write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')' - print*,'---------------------------------------------' + print*,'-------------------------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & + ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2 + print*,'-------------------------------------------------------------' jb = 0 do j=nC+1,nO @@ -108,8 +69,6 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, end do write(*,*) - print*,' = ',2d0*sum(X(:)**2 + Y(:)**2) - end do ! Thomas-Reiche-Kuhn sum rule diff --git a/src/LR/print_unrestricted_transition_vectors.f90 b/src/LR/print_unrestricted_transition_vectors.f90 index ab9edfc..11b5fb5 100644 --- a/src/LR/print_unrestricted_transition_vectors.f90 +++ b/src/LR/print_unrestricted_transition_vectors.f90 @@ -40,62 +40,14 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n ! Memory allocation - allocate(X(nSt),Y(nSt),f(nSt,ncart),os(nSt)) + allocate(X(nSt),Y(nSt),os(nSt)) -! Initialization - - f(:,:) = 0d0 - os(:) = 0d0 +! Compute oscillator strengths -! Compute dipole moments and oscillator strengths + os(:) = 0d0 + if(spin_allowed) call unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt, & + dipole_int_aa,dipole_int_bb,Omega,XpY,XmY,os) - - if(spin_allowed) then - - do ia=1,nSt - do ixyz=1,ncart - - jb = 0 - do j=nC(1)+1,nO(1) - do b=nO(1)+1,nBas-nR(1) - jb = jb + 1 - f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(j,b,ixyz)*XpY(ia,jb) - end do - end do - - jb = 0 - do j=nC(2)+1,nO(2) - do b=nO(2)+1,nBas-nR(2) - jb = jb + 1 - f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(j,b,ixyz)*XpY(ia,nSa+jb) - end do - end do - - end do - end do - - do ia=1,nSt - os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) - end do - - if(debug) then - - write(*,*) '----------------' - write(*,*) ' Dipole moments ' - write(*,*) '----------------' - call matout(nSt,ncart,f(:,:)) - write(*,*) - - write(*,*) '----------------------' - write(*,*) ' Oscillator strengths ' - write(*,*) '----------------------' - call matout(nSt,1,os(:)) - write(*,*) - - end if - - end if - ! Print details about excitations do ia=1,min(nSt,maxS) diff --git a/src/LR/unrestricted_oscillator_strength.f90 b/src/LR/unrestricted_oscillator_strength.f90 new file mode 100644 index 0000000..4a4f8b4 --- /dev/null +++ b/src/LR/unrestricted_oscillator_strength.f90 @@ -0,0 +1,90 @@ +subroutine unrestricted_oscillator_strength(nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,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 + double precision :: dipole_int_aa(nBas,nBas,ncart) + double precision :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: XpY(nSt,nSt) + double precision,intent(in) :: XmY(nSt,nSt) + +! Local variables + + logical :: debug = .false. + integer :: ia,jb,i,j,a,b + integer :: ixyz + + double precision,allocatable :: f(:,:) + +! Output variables + + double precision :: os(nSt) + +! Memory allocation + + allocate(f(nSt,ncart)) + +! Initialization + + f(:,:) = 0d0 + +! Compute dipole moments and oscillator strengths + + do ia=1,nSt + do ixyz=1,ncart + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(j,b,ixyz)*XpY(ia,jb) + end do + end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(j,b,ixyz)*XpY(ia,nSa+jb) + end do + end do + + end do + end do + + do ia=1,nSt + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + if(debug) then + + write(*,*) '------------------------' + write(*,*) ' Dipole moments (X Y Z) ' + write(*,*) '------------------------' + call matout(nS,ncart,f) + write(*,*) + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nS,1,os) + write(*,*) + + end if + +end subroutine unrestricted_oscillator_strength