MRCC eigenfunction of S2

plugins/MRCC_Utils/amplitudes.irp.f

@@ -137,6 +137,7 @@ END_PROVIDER
 
  BEGIN_PROVIDER [ integer, mrcc_AtA_ind, (N_det_ref * n_exc_active) ]
 &BEGIN_PROVIDER [ double precision, mrcc_AtA_val, (N_states, N_det_ref * n_exc_active) ]
+&BEGIN_PROVIDER [ double precision, mrcc_AtS2A_val, (N_states, N_det_ref * n_exc_active) ]
 &BEGIN_PROVIDER [ integer, mrcc_col_shortcut, (n_exc_active) ]
 &BEGIN_PROVIDER [ integer, mrcc_N_col, (n_exc_active) ]
  implicit none
@@ -145,11 +146,15 @@ END_PROVIDER
  END_DOC
  integer                        :: AtA_size, i,k
  integer                        :: at_roww, at_row, wk, a_coll, a_col, r1, r2, s
- double precision, allocatable  :: t(:), A_val_mwen(:,:)
+ double precision, allocatable  :: t(:), ts(:), A_val_mwen(:,:), As2_val_mwen(:,:)
  integer, allocatable           :: A_ind_mwen(:)
+ integer(bit_kind)              :: det1(N_int,2), det2(N_int,2)
+ double precision               :: sij
+ PROVIDE psi_non_ref S_z2_Sz S_z
 
  mrcc_AtA_ind(:) = 0
  mrcc_AtA_val(:,:) = 0.d0
+ mrcc_AtS2A_val(:,:) = 0.d0
  mrcc_col_shortcut(:) = 0
  mrcc_N_col(:) = 0
  AtA_size = 0
@@ -157,9 +162,11 @@ END_PROVIDER
 
   !$OMP PARALLEL default(none) shared(k, active_excitation_to_determinants_idx,&
       !$OMP   active_excitation_to_determinants_val, hh_nex)         &
-      !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
+      !$OMP private(at_row, a_col, t, ts, i, r1, r2, wk, A_ind_mwen, A_val_mwen,&
-      !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtA_size, mrcc_AtA_val, mrcc_AtA_ind, n_exc_active, active_pp_idx)
+      !$OMP   det1,det2,As2_val_mwen, a_coll, at_roww,sij) &
-  allocate(A_val_mwen(N_states,hh_nex), A_ind_mwen(hh_nex), t(N_states))
+      !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtA_size, mrcc_AtA_val, mrcc_AtA_ind, &
+      !$OMP n_exc_active, active_pp_idx,psi_non_ref,N_int,S_z2_Sz, mrcc_AtS2A_val)
+  allocate(A_val_mwen(N_states,hh_nex), As2_val_mwen(N_states,hh_nex), A_ind_mwen(hh_nex), t(N_states), ts(N_states))
 
   !$OMP DO schedule(dynamic, 100)
   do at_roww = 1, n_exc_active ! hh_nex
@@ -170,6 +177,7 @@ END_PROVIDER
     do a_coll = 1, n_exc_active
       a_col = active_pp_idx(a_coll)
       t(:) = 0d0
+      ts(:) = 0d0
       r1 = 1
       r2 = 1
       do while ((active_excitation_to_determinants_idx(r1, at_roww) /= 0).and.(active_excitation_to_determinants_idx(r2, a_coll) /= 0))
@@ -178,8 +186,12 @@ END_PROVIDER
         else if(active_excitation_to_determinants_idx(r1, at_roww) < active_excitation_to_determinants_idx(r2, a_coll)) then
           r1 = r1+1
         else
+          det1(:,:) = psi_non_ref(:,:, active_excitation_to_determinants_idx(r1,at_roww))
+          det2(:,:) = psi_non_ref(:,:, active_excitation_to_determinants_idx(r2,a_coll))
+          call get_s2(det1, det2,N_int,sij)
           do s=1,N_states
             t(s) = t(s) - active_excitation_to_determinants_val(s,r1, at_roww) * active_excitation_to_determinants_val(s,r2, a_coll)
+            ts(s) = ts(s) - active_excitation_to_determinants_val(s,r1, at_roww) * active_excitation_to_determinants_val(s,r2, a_coll) * sij
           enddo
           r1 = r1+1
           r2 = r2+1
@@ -189,6 +201,7 @@ END_PROVIDER
       if(a_col == at_row) then
         do s=1,N_states
           t(s) = t(s) + 1.d0
+          ts(s) = ts(s) + S_z2_Sz
         enddo
       end if
       if(sum(abs(t)) /= 0.d0) then
@@ -196,6 +209,7 @@ END_PROVIDER
         A_ind_mwen(wk) = a_col
         do s=1,N_states
           A_val_mwen(s,wk) = t(s)
+          As2_val_mwen(s,wk) = ts(s)
         enddo
       end if
     end do
@@ -212,6 +226,7 @@ END_PROVIDER
         mrcc_AtA_ind(AtA_size+i) = A_ind_mwen(i)
         do s=1,N_states
           mrcc_AtA_val(s,AtA_size+i) = A_val_mwen(s,i)
+          mrcc_AtS2A_val(s,AtA_size+i) = As2_val_mwen(s,i)
         enddo
       enddo
       AtA_size += wk
@@ -219,7 +234,7 @@ END_PROVIDER
     end if
   end do
   !$OMP END DO NOWAIT
-  deallocate (A_ind_mwen, A_val_mwen, t)
+  deallocate (A_ind_mwen, A_val_mwen, As2_val_mwen, t, ts)
   !$OMP END PARALLEL
 
   print *, "ATA SIZE", ata_size

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -626,8 +626,6 @@ END_PROVIDER
   double precision               :: norm, cx, res
   integer, allocatable           :: lref(:), A_ind_mwen(:)
   double precision               :: phase
-!  double precision , allocatable :: mrcc_AtA_val(:,:)
-!  integer, allocatable           :: mrcc_AtA_ind(:), col_shortcut(:), , mrcc_N_col(:)
   
   
   double precision, allocatable :: rho_mrcc_init(:,:)
@@ -635,13 +633,10 @@ END_PROVIDER
   
   print *, "TI", hh_nex, N_det_non_ref
 
-  
-  
   allocate(rho_mrcc_init(N_det_non_ref, N_states))
-  
+  allocate(x_new(hh_nex))
   allocate(x(hh_nex), AtB(hh_nex))
   x = 0d0
-  allocate(x_new(hh_nex))
         
 
   do s=1,N_states
@@ -712,29 +707,38 @@ END_PROVIDER
     
     x_new = x
     
+    double precision               :: s2(N_states), s2_local, dx
     double precision               :: factor, resold
     factor = 1.d0
     resold = huge(1.d0)
     do k=0,100000
-      !$OMP PARALLEL default(shared) private(cx, i, j, a_col, a_coll)
+      !$OMP PARALLEL default(shared) private(cx, dx, i, j, a_col, a_coll, s2_local)
       
       !$OMP DO
       do i=1,N_det_non_ref
-        rho_mrcc(i,s) = rho_mrcc_init(i,s) ! 0d0
+        rho_mrcc(i,s) = rho_mrcc_init(i,s) 
       enddo
       !$OMP END DO
       
+      s2(s) = 0.d0
       !$OMP DO
-      do a_coll = 1, n_exc_active !: hh_nex
+      do a_coll = 1, n_exc_active
         a_col = active_pp_idx(a_coll)
-        cx = 0d0
+        cx = 0.d0
+        dx = 0.d0
         do i=mrcc_col_shortcut(a_coll), mrcc_col_shortcut(a_coll) + mrcc_N_col(a_coll) - 1
           cx = cx + x(mrcc_AtA_ind(i)) * mrcc_AtA_val(s,i)
+          dx = dx + x(mrcc_AtA_ind(i)) * mrcc_AtS2A_val(s,i)
+          s2_local = s2_local + X(a_col)*X(mrcc_AtA_ind(i))*mrcc_AtS2A_val(s,i)
         end do
-        x_new(a_col) = AtB(a_col) + cx * factor
+        x_new(a_col) = AtB(a_col) + (cx+dx) * factor
       end do
       !$OMP END DO
 
+      !$OMP CRITICAL
+      s2(s) = s2(s) + s2_local
+      !$OMP END CRITICAL
+      
       !$OMP END PARALLEL
       
       
@@ -756,14 +760,13 @@ END_PROVIDER
       resold = res
       
       if(mod(k, 100) == 0) then
-        print *, "res ", k, res, factor
+        print *, "res ", k, res, s2(s)
       end if
       
       if(res < 1d-9) exit
     end do
     
     
-    
     norm = 0.d0
     do i=1,N_det_non_ref
       norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)