Added Phi_S and Natural transition MOs with thibaudetienne

plugins/analyze_wf/attachment.irp.f

@@ -0,0 +1,64 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_detachment, (mo_tot_num,mo_tot_num,2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_attachment, (mo_tot_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Detachment and attachment density matrices in MO basis
+ END_DOC
+ integer :: i,j, k, istate
+ double precision :: km(mo_tot_num), kp(mo_tot_num)
+
+ one_body_dm_mo_detachment = 0.d0
+ one_body_dm_mo_attachment = 0.d0
+
+ do istate=2,N_states
+
+   km(:) = 0.d0
+   kp(:) = 0.d0
+   do i=1,mo_tot_num
+      if (one_body_dm_mo_diff_eigvalues(i,istate) < 0) then
+        km(i) = -one_body_dm_mo_diff_eigvalues(i,istate)
+      else
+        kp(i) =  one_body_dm_mo_diff_eigvalues(i,istate)
+      endif
+   enddo
+
+   ! Attachment
+  do k=1,mo_tot_num
+     do j=1,mo_tot_num
+       do i=1,mo_tot_num
+         one_body_dm_mo_detachment(i,j,istate) = one_body_dm_mo_detachment(i,j,istate) + &
+           one_body_dm_mo_diff_eigvectors(i,k,istate)*km(k)* &
+           one_body_dm_mo_diff_eigvectors(j,k,istate)
+         one_body_dm_mo_attachment(i,j,istate) = one_body_dm_mo_attachment(i,j,istate) + &
+           one_body_dm_mo_diff_eigvectors(i,k,istate)*kp(k)* &
+           one_body_dm_mo_diff_eigvectors(j,k,istate)
+       enddo
+    enddo
+  enddo
+
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, one_body_dm_ao_detachment, (ao_num,ao_num,2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_ao_attachment, (ao_num,ao_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Detachment and attachment density matrices in AO basis
+ END_DOC
+ integer :: istate
+ do istate=2,N_states
+   call mo_to_ao_no_overlap(                                         &
+       one_body_dm_mo_attachment(1,1,istate),                        &
+       size(one_body_dm_mo_attachment,1),                            &
+       one_body_dm_ao_attachment(1,1,istate),                        &
+       size(one_body_dm_ao_attachment,1) )
+   call mo_to_ao_no_overlap(                                         &
+       one_body_dm_mo_detachment(1,1,istate),                        &
+       size(one_body_dm_mo_detachment,1),                            &
+       one_body_dm_ao_detachment(1,1,istate),                        &
+       size(one_body_dm_ao_detachment,1) )
+ enddo
+
+END_PROVIDER
+

plugins/analyze_wf/dump_nto.irp.f

@@ -0,0 +1,19 @@
+program dump_nto
+  implicit none
+  integer :: i,j, istate
+
+  print *,  'Phi_S'
+  do i=2,N_states
+    print *,  i, Phi_S(i)
+  enddo
+
+  do istate=2,N_states
+    do j=1,mo_tot_num
+      print *,  'MO: ', j, 'State:', istate, 'Eig:', one_body_dm_mo_diff_eigvalues(j,istate)
+      do i=1,ao_num
+        print *,  i, transition_natorb(i,j,istate)
+      enddo
+    enddo
+  enddo
+
+end

plugins/analyze_wf/dump_one_body_mos.irp.f

@@ -0,0 +1,38 @@
+program dump_one_body_mos
+  implicit none
+  BEGIN_DOC
+! Output density matrices of all the states
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run()
+end
+
+subroutine run
+  implicit none
+  integer                        :: istate
+  integer, parameter             :: iunit = 66
+  character*(64)                 :: filename, fmt
+  integer                        :: i,j,k
+
+  write(fmt,'(''('',I4.4,''(X,E20.14))'')') mo_tot_num
+  do istate=1,N_states
+    write(filename,'(''state.'',I2.2)') istate
+    open(unit=iunit, form='formatted', file=filename)
+    write(iunit,*) mo_tot_num
+    do j=1,mo_tot_num
+      write(iunit,fmt) one_body_dm_mo_alpha(1:mo_tot_num,j,istate) + one_body_dm_mo_beta(1:mo_tot_num,j,istate)
+    enddo
+  enddo
+
+
+  call run2()
+end
+subroutine run2
+  integer :: i,j, istate
+  print *,  'Phi_S'
+  do i=2,N_states
+    print *,  i, Phi_S(i)
+  enddo
+end
+

plugins/analyze_wf/phi_s.irp.f

@@ -0,0 +1,128 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvalues,  (mo_tot_num, 2:N_states) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvectors, (mo_tot_num, mo_tot_num, 2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Eigenvalues and eigenvectors of one_body_dm_mo_diff
+ END_DOC
+ integer                        :: i,j,istate
+ integer                        :: liwork, lwork, n, info
+ integer, allocatable           :: iwork(:)
+ double precision, allocatable  :: work(:)
+ 
+ 
+ one_body_dm_mo_diff_eigvectors(1:mo_tot_num, 1:mo_tot_num, 2:N_states) =&
+     one_body_dm_mo_diff(1:mo_tot_num, 1:mo_tot_num, 2:N_states)
+ 
+ n = mo_tot_num
+ lwork = 1+6*n + 2*n*n
+ liwork = 3 + 5*n
+ 
+ allocate(work(lwork), iwork(liwork))
+ 
+ lwork=-1
+ liwork=-1
+ istate=2
+ 
+ call dsyevd( 'V', 'U', mo_tot_num,                                  &
+     one_body_dm_mo_diff_eigvectors(1,1,istate),                     &
+     size(one_body_dm_mo_diff_eigvectors,1),                         &
+     one_body_dm_mo_diff_eigvalues(1,istate),                        &
+     work, lwork, iwork, liwork, info)
+ 
+ 
+ if (info /= 0) then
+   print *,  irp_here//' DSYEVD failed : ', info
+   stop 1
+ endif
+ lwork = int(work(1))
+ liwork = iwork(1)
+ deallocate(iwork,work)
+ 
+ allocate(work(lwork), iwork(liwork))
+ 
+ do istate=2,N_states
+   call dsyevd( 'V', 'U', mo_tot_num,                                &
+       one_body_dm_mo_diff_eigvectors(1,1,istate),                   &
+       size(one_body_dm_mo_diff_eigvectors,1),                       &
+       one_body_dm_mo_diff_eigvalues(1,istate),                      &
+       work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+       print *,  irp_here//' DSYEVD failed : ', info
+       stop 1
+   endif
+
+ enddo
+
+ deallocate(iwork,work)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, transition_natorb, (ao_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Natural transition molecular orbitals
+ END_DOC
+
+ integer :: istate
+
+ do istate=2,N_states
+   call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
+       mo_coef, size(mo_coef,1),                                     &
+       one_body_dm_mo_diff_eigvectors(1,1,istate),                   &
+       size(one_body_dm_mo_diff_eigvectors,1), 0.d0,                 &
+       transition_natorb(1,1,istate), size(transition_natorb,1))
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, phi_s, (2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! 
+ END_DOC
+
+ integer :: i,istate
+ double precision, allocatable :: T(:,:), A(:,:), D(:,:)
+ double precision :: trace, norm
+ allocate(T(ao_num,ao_num), A(ao_num,ao_num), D(ao_num,ao_num))
+
+  do istate=2,N_states
+
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        S_half, size(S_half,1),                                      &
+        one_body_dm_ao_attachment(1,1,istate), size(one_body_dm_ao_attachment,1), 0.d0,&
+        T, size(T,1))
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        T, size(T,1),                                                &
+        S_half, size(S_half,1), 0.d0,                                &
+        A, size(A,1))
+!
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        S_half, size(S_half,1),                                      &
+        one_body_dm_ao_detachment(1,1,istate), size(one_body_dm_ao_detachment,1), 0.d0,&
+        T, size(T,1))
+    call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
+        T, size(T,1),                                                &
+        S_half, size(S_half,1), 0.d0,                                &
+        D, size(D,1))
+
+    trace = 0.d0
+    do i=1,ao_num
+      trace = trace + A(i,i)
+    enddo
+    norm = 0.d0
+    do i=1,ao_num
+      norm = norm + D(i,i)
+    enddo
+    norm = 0.5d0*(norm + trace)
+
+    trace = 0.d0
+    do i=1,mo_tot_num
+      trace = trace + dsqrt(A(i,i)*D(i,i))
+    enddo
+    phi_s(istate) = trace/norm 
+  enddo
+
+END_PROVIDER
+

src/AO_Basis/ao_overlap.irp.f

@@ -194,3 +194,39 @@ BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
   
 
 END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num)  ]
+ implicit none
+ BEGIN_DOC
+ ! S^{1/2}
+ END_DOC
+
+  integer :: i,j,k
+  double precision, allocatable  :: U(:,:)
+  double precision, allocatable  :: Vt(:,:)
+  double precision, allocatable  :: D(:)
+  
+  allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num))
+
+  call svd(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num)
+
+  do i=1,ao_num
+    D(i) = dsqrt(D(i))
+    do j=1,ao_num
+      S_half(j,i) = 0.d0
+    enddo
+  enddo
+
+  do k=1,ao_num
+      do j=1,ao_num
+        do i=1,ao_num
+          S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)
+        enddo
+      enddo
+  enddo
+  
+  deallocate(U,Vt,D)
+
+END_PROVIDER
+

src/Determinants/density_matrix.irp.f

@@ -15,6 +15,32 @@
    enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff, (mo_tot_num,mo_tot_num,2:N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Difference of the one-body density matrix with respect to the ground state
+ END_DOC
+ integer :: i,j, istate
+
+ do istate=2,N_states
+  do j=1,mo_tot_num
+    do i=1,mo_tot_num
+      one_body_dm_mo_diff(i,j,istate) = &
+        one_body_dm_mo_alpha(i,j,istate) - one_body_dm_mo_alpha(i,j,1) + &
+        one_body_dm_mo_beta (i,j,istate) - one_body_dm_mo_beta (i,j,1) 
+    enddo
+  enddo
+  double precision :: trace
+  trace = 0.d0
+  do i=1,mo_tot_num
+    trace += one_body_dm_mo_diff(i,i,istate)
+  enddo
+  print *, irp_here, trace
+ enddo
+ 
+END_PROVIDER
+
+
  BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ]
  implicit none 
  integer :: i,j,ispin,istate

src/Utils/LinearAlgebra.irp.f

@@ -69,7 +69,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   double precision, allocatable  :: U(:,:)
   double precision, allocatable  :: Vt(:,:)
   double precision, allocatable  :: D(:)
-  double precision, allocatable  :: S_half(:,:)
+  double precision, allocatable  :: S(:,:)
   !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
   integer                        :: info, i, j
   
@@ -77,7 +77,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
     return
   endif
 
-  allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n))
+  allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n))
 
   call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
 
@@ -103,13 +103,13 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
 
 
   !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
+  !$OMP SHARED(S,U,D,Vt,n,C,m) &
   !$OMP PRIVATE(i,j)
 
   !$OMP DO
   do j=1,n
     do i=1,n
-      S_half(i,j) = U(i,j)*D(j)
+      S(i,j) = U(i,j)*D(j)
     enddo
     do i=1,n
       U(i,j) = C(i,j)
@@ -119,8 +119,8 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   
   !$OMP END PARALLEL
 
-  call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1))
+  call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1))
-  deallocate (U, Vt, D, S_half)
+  deallocate (U, Vt, D, S)
   
 end
 
@@ -210,19 +210,19 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   double precision, allocatable  :: U(:,:)
   double precision, allocatable  :: Vt(:,:)
   double precision, allocatable  :: D(:)
-  double precision, allocatable  :: S_half(:,:)
+  double precision, allocatable  :: S(:,:)
   integer                        :: info, i, j, k
   
   if (n < 2) then
     return
   endif
 
-  allocate(U(ldc,n),Vt(lda,n),S_half(lda,n),D(n))
+  allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n))
 
   call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
 
   !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
+  !$OMP SHARED(S,U,D,Vt,n,C,m) &
   !$OMP PRIVATE(i,j,k)
 
   !$OMP DO
@@ -233,7 +233,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
       D(i) = 1.d0/dsqrt(D(i))
     endif
     do j=1,n
-      S_half(j,i) = 0.d0
+      S(j,i) = 0.d0
     enddo
   enddo
   !$OMP END DO
@@ -243,7 +243,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
       !$OMP DO
       do j=1,n
         do i=1,n
-          S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)
+          S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j)
         enddo
       enddo
       !$OMP END DO NOWAIT
@@ -261,9 +261,9 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   
   !$OMP END PARALLEL
 
-  call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1))
+  call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1))
   
-  deallocate(U,Vt,S_half,D)
+  deallocate(U,Vt,S,D)
 end