S2 OK in TC

src/tc_bi_ortho/dav_h_tc_s2.irp.f

@@ -1,7 +1,7 @@
 
 ! ---
 
-subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, energies, sze, N_st, N_st_diag_in, converged, hcalc)
+subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, s2_out,energies, sze, N_st, N_st_diag_in, converged, hcalc)
 
   use mmap_module
 
@@ -30,6 +30,7 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, energies, sze, N_st, N_st_dia
   logical,          intent(inout) :: converged
   double precision, intent(inout) :: u_in(sze,N_st_diag_in)
   double precision, intent(out)   :: energies(N_st)
+  double precision,  intent(inout) :: s2_out(N_st)
   external                           hcalc
 
   character*(16384)               :: write_buffer
@@ -528,6 +529,7 @@ subroutine davidson_hs2_nonsym_b1space(u_in, H_jj, energies, sze, N_st, N_st_dia
 
   do k = 1, N_st
     energies(k) = lambda(k)
+    s2_out(k) = s2(k)
   enddo
   write_buffer = '====='
   do i = 1, N_st

src/tc_bi_ortho/h_tc_s2_u0.irp.f

@@ -1,3 +1,33 @@
+
+subroutine get_H_tc_s2_l0_r0(l_0,r_0,N_st,sze,energies, s2)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Computes $e_0 = \langle l_0 | H | r_0\rangle$.
+  !
+  ! Computes $s_0 = \langle l_0 | S^2 | r_0\rangle$.
+  !
+  ! Assumes that the determinants are in psi_det
+  !
+  ! istart, iend, ishift, istep are used in ZMQ parallelization.
+  END_DOC
+  integer, intent(in)              :: N_st,sze
+  double precision, intent(in)     :: l_0(sze,N_st), r_0(sze,N_st)
+  double precision, intent(out)    :: energies(N_st), s2(N_st)
+  logical           :: do_right 
+  integer :: istate
+  double precision, allocatable :: s_0(:,:), v_0(:,:)
+  double precision :: u_dot_v, norm
+  allocate(s_0(sze,N_st), v_0(sze,N_st))
+  do_right = .True.
+  call H_tc_s2_u_0_opt(v_0,s_0,r_0,N_st,sze)
+  do istate = 1, N_st
+   norm = u_dot_v(l_0(1,istate),r_0(1,istate),sze)
+   energies(istate) = u_dot_v(l_0(1,istate),v_0(1,istate),sze)/norm
+   s2(istate) = u_dot_v(l_0(1,istate),s_0(1,istate),sze)/norm
+  enddo
+end
+
 subroutine H_tc_s2_u_0_opt(v_0,s_0,u_0,N_st,sze)
   use bitmasks
   implicit none

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -35,6 +35,7 @@ end
 &BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)]
 &BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)]
 &BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)]
+&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth, (N_states)]
 &BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ]
 
   BEGIN_DOC
@@ -46,64 +47,153 @@ end
   logical                       :: converged, dagger
   integer                       :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
   double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
+  double precision, allocatable :: s2_values_tmp(:), H_prime(:,:), expect_e(:)
+  double precision, parameter :: alpha = 0.1d0
+  integer                        :: i_good_state,i_other_state, i_state
+  integer, allocatable           :: index_good_state_array(:)
+  logical, allocatable           :: good_state_array(:)
+  double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
+  integer, allocatable :: iorder(:)
 
   PROVIDE N_det N_int
 
    if(n_det.le.N_det_max_full)then
-    allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
+    allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det),expect_e(N_det))
-    call non_hrmt_real_diag(N_det,htilde_matrix_elmt_bi_ortho,& 
+    allocate (H_prime(N_det,N_det),s2_values_tmp(N_det))
+    H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det)
+    if(s2_eig)then
+       H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det)
+       do j=1,N_det
+         H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
+       enddo
+    endif
+    call non_hrmt_real_diag(N_det,H_prime,& 
          leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& 
          n_real_tc_bi_orth_eigval_right,eigval_right_tmp)
-    double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
+!    do i = 1, N_det
-    integer, allocatable :: iorder(:)
+!     call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp(1,i),reigvec_tc_bi_orth_tmp(1,i),1,N_det,expect_e(i), s2_values_tmp(i))
-    allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
+!    enddo
-    do i = 1,N_det
+    call get_H_tc_s2_l0_r0(leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,N_det,N_det,expect_e, s2_values_tmp)
-     iorder(i) = i
+    allocate(index_good_state_array(N_det),good_state_array(N_det))
-     coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
+    i_state = 0
-    enddo
+    good_state_array = .False.
-    call dsort(coef_hf_r,iorder,N_det)
+    if(s2_eig)then
-    igood_r = iorder(1)
+     if (only_expected_s2) then
-    print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1)
+       do j=1,N_det
-    do i = 1,N_det
+         ! Select at least n_states states with S^2 values closed to "expected_s2"
-     iorder(i) = i
+!         print*,'s2_values_tmp(j) = ',s2_values_tmp(j),eigval_right_tmp(j),expect_e(j)
-     coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
+         if(dabs(s2_values_tmp(j)-expected_s2).le.0.5d0)then
-    enddo
+           i_state +=1
-    call dsort(coef_hf_l,iorder,N_det)
+           index_good_state_array(i_state) = j
-    igood_l = iorder(1)
+           good_state_array(j) = .True.
-    print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1)
+         endif
-
+         if(i_state.eq.N_states) then
-    if(igood_r.ne.igood_l.and.igood_r.ne.1)then
+           exit
-     print *,''
+         endif
-     print *,'Warning, the left and right eigenvectors are "not the same" '
-     print *,'Warning, the ground state is not dominated by HF...'
-     print *,'State with largest RIGHT coefficient of HF ',igood_r
-     print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r)
-     print *,'State with largest LEFT  coefficient of HF ',igood_l
-     print *,'coef of HF in LEFT  eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
-    endif
-    if(state_following_tc)then
-     print *,'Following the states with the largest coef on HF'
-     print *,'igood_r,igood_l',igood_r,igood_l
-     i= igood_r
-     eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
-     do j = 1, N_det
-       reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
-!       print*,reigvec_tc_bi_orth(j,1)
-     enddo
-     i= igood_l
-     eigval_left_tc_bi_orth(1)  = eigval_right_tmp(i)
-     do j = 1, N_det
-       leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
-     enddo
-    else 
-     do i = 1, N_states
-       eigval_right_tc_bi_orth(i) = eigval_right_tmp(i)
-       eigval_left_tc_bi_orth(i)  = eigval_right_tmp(i)
-       do j = 1, N_det
-         reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i)
-         leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i)
        enddo
-     enddo
+     else
+       do j=1,N_det
+         index_good_state_array(j) = j
+         good_state_array(j) = .True.
+       enddo
+     endif
+     if(i_state .ne.0)then
+       ! Fill the first "i_state" states that have a correct S^2 value
+       do j = 1, i_state
+         do i=1,N_det
+           reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,index_good_state_array(j))
+           leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,index_good_state_array(j))
+         enddo
+         eigval_right_tc_bi_orth(j) = expect_e(index_good_state_array(j))
+         eigval_left_tc_bi_orth(j)  = expect_e(index_good_state_array(j))
+         s2_eigvec_tc_bi_orth(j)    = s2_values_tmp(index_good_state_array(j))
+       enddo
+       i_other_state = 0
+       do j = 1, N_det
+         if(good_state_array(j))cycle
+         i_other_state +=1
+         if(i_state+i_other_state.gt.n_states)then
+           exit
+         endif
+         do i=1,N_det
+           reigvec_tc_bi_orth(i,i_state+i_other_state) = reigvec_tc_bi_orth_tmp(i,j)
+           leigvec_tc_bi_orth(i,i_state+i_other_state) = leigvec_tc_bi_orth_tmp(i,j)
+         enddo
+         eigval_right_tc_bi_orth(i_state+i_other_state) = eigval_right_tmp(j)
+         eigval_left_tc_bi_orth (i_state+i_other_state) = eigval_right_tmp(j)
+         s2_eigvec_tc_bi_orth(i_state+i_other_state)    = s2_values_tmp(i_state+i_other_state)
+       enddo
+     else ! istate == 0
+         print*,''
+         print*,'!!!!!!!!   WARNING  !!!!!!!!!'
+         print*,'  Within the ',N_det,'determinants selected'
+         print*,'  and the ',N_states_diag,'states requested'
+         print*,'  We did not find only states with S^2 values close to ',expected_s2
+         print*,'  We will then set the first N_states eigenvectors of the H matrix'
+         print*,'  as the CI_eigenvectors'
+         print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
+         print*,''
+         do j=1,min(N_states_diag,N_det)
+           do i=1,N_det
+             leigvec_tc_bi_orth(i,j) = leigvec_tc_bi_orth_tmp(i,j)
+             reigvec_tc_bi_orth(i,j) = reigvec_tc_bi_orth_tmp(i,j)
+           enddo
+           eigval_right_tc_bi_orth(j) = eigval_right_tmp(j)
+           eigval_left_tc_bi_orth (j) = eigval_right_tmp(j)
+           s2_eigvec_tc_bi_orth(j)    = s2_values_tmp(j)
+         enddo
+     endif ! istate .ne. 0
+
+    else ! s2_eig
+       allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
+       do i = 1,N_det
+        iorder(i) = i
+        coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
+       enddo
+       call dsort(coef_hf_r,iorder,N_det)
+       igood_r = iorder(1)
+       print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1)
+       do i = 1,N_det
+        iorder(i) = i
+        coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
+       enddo
+       call dsort(coef_hf_l,iorder,N_det)
+       igood_l = iorder(1)
+       print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1)
+       
+       if(igood_r.ne.igood_l.and.igood_r.ne.1)then
+        print *,''
+        print *,'Warning, the left and right eigenvectors are "not the same" '
+        print *,'Warning, the ground state is not dominated by HF...'
+        print *,'State with largest RIGHT coefficient of HF ',igood_r
+        print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r)
+        print *,'State with largest LEFT  coefficient of HF ',igood_l
+        print *,'coef of HF in LEFT  eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
+       endif
+       if(state_following_tc)then
+        print *,'Following the states with the largest coef on HF'
+        print *,'igood_r,igood_l',igood_r,igood_l
+        i= igood_r
+        eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
+        do j = 1, N_det
+          reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
+!          print*,reigvec_tc_bi_orth(j,1)
+        enddo
+        i= igood_l
+        eigval_left_tc_bi_orth(1)  = eigval_right_tmp(i)
+        do j = 1, N_det
+          leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
+        enddo
+       else 
+        do i = 1, N_states
+          eigval_right_tc_bi_orth(i) = eigval_right_tmp(i)
+          eigval_left_tc_bi_orth(i)  = eigval_right_tmp(i)
+          do j = 1, N_det
+            reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i)
+            leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i)
+          enddo
+        enddo
+       endif
     endif
    else 
     double precision, allocatable :: H_jj(:),vec_tmp(:,:)
@@ -111,6 +201,8 @@ end
     external                         htcdag_bi_ortho_calc_tdav
     external                         H_tc_u_0_opt
     external                         H_tc_dagger_u_0_opt
+    external                         H_tc_s2_dagger_u_0_opt
+    external                         H_tc_s2_u_0_opt
     allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
     do i = 1, N_det
       call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
@@ -125,7 +217,8 @@ end
      vec_tmp(istate,istate) = 1.d0
     enddo
 !    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav)
-    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt)
+!    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_dagger_u_0_opt)
+    call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
     do istate = 1, N_states
      leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
     enddo
@@ -140,7 +233,8 @@ end
      vec_tmp(istate,istate) = 1.d0
     enddo
 !    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
-    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt)
+!    call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_u_0_opt)
+    call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2_eigvec_tc_bi_orth, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
     do istate = 1, N_states
      reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
     enddo
@@ -154,6 +248,7 @@ end
     norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)
    enddo
    print*,'norm l/r = ',norm_ground_left_right_bi_orth
+   print*,'<S2>     = ',s2_eigvec_tc_bi_orth(1)
 
 END_PROVIDER 
 

src/tc_bi_ortho/test_s2_tc.irp.f

@@ -84,12 +84,12 @@ end
 
 subroutine routine_test_s2_davidson
  implicit none
- double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) 
+ double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) , s2(:)
  integer :: i,istate
  logical :: converged 
  external H_tc_s2_dagger_u_0_opt
  external H_tc_s2_u_0_opt
- allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag))
+ allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag), s2(n_states_diag))
  do i = 1, N_det
    call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
  enddo
@@ -105,8 +105,7 @@ subroutine routine_test_s2_davidson
  do istate = 1, N_states
   leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
  enddo
- call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
+ call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_dagger_u_0_opt)
- print*,'energies = ',energies
  double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
  integer :: sze,N_st
  logical           :: do_right 
@@ -122,6 +121,8 @@ subroutine routine_test_s2_davidson
   accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
   accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
  enddo
+ print*,'energies = ',energies
+ print*,'s2       = ',s2
  print*,'accu_e_0',accu_e_0
  print*,'accu_s_0',accu_s_0
 
@@ -137,8 +138,7 @@ subroutine routine_test_s2_davidson
  do istate = 1, N_states
   leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
  enddo
- call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_u_0_opt)
+ call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, converged, H_tc_s2_u_0_opt)
- print*,'energies = ',energies
  sze = N_det
  N_st = 1
  do_right = .True.
@@ -151,6 +151,8 @@ subroutine routine_test_s2_davidson
   accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
   accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
  enddo
+ print*,'energies = ',energies
+ print*,'s2       = ',s2
  print*,'accu_e_0',accu_e_0
  print*,'accu_s_0',accu_s_0