diff --git a/plugins/MRCC_Utils/amplitudes.irp.f b/plugins/MRCC_Utils/amplitudes.irp.f index 2694aa75..77360d93 100644 --- a/plugins/MRCC_Utils/amplitudes.irp.f +++ b/plugins/MRCC_Utils/amplitudes.irp.f @@ -49,7 +49,7 @@ end do end do end do - +!is_active_exc=.true. do hh = 1, hh_shortcut(0) do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 if(is_active_exc(pp)) then @@ -133,6 +133,55 @@ END_PROVIDER deallocate(lref) !$OMP END PARALLEL +END_PROVIDER + + BEGIN_PROVIDER [ integer, mrcc_S2A_ind, (0:N_det_ref*mo_tot_num, n_exc_active) ] +&BEGIN_PROVIDER [ double precision, mrcc_S2A_val, (N_states, N_det_ref*mo_tot_num, n_exc_active) ] + implicit none + BEGIN_DOC + ! A is active_excitation_to_determinants in S^2. + END_DOC + integer :: a_coll, a_col + integer :: i,j,idx,s + double precision :: sij + double precision, allocatable :: tmp(:,:) + logical, allocatable :: ok(:) + mrcc_S2A_val = 0.d0 + print *, 'Computing S2A' + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_coll,idx,i,sij,s,tmp,ok,j) + allocate(tmp(N_states,N_det_non_ref), ok(N_det_non_ref)) + !$OMP DO + do a_coll=1, n_exc_active + tmp = 0.d0 + ok = .False. + do idx=1, active_excitation_to_determinants_idx(0,a_coll) + i = active_excitation_to_determinants_idx(idx,a_coll) + do j=1,N_det_non_ref + call get_s2(psi_non_ref(1,1,i), psi_non_ref(1,1,j), N_int, sij) + if (sij /= 0.d0) then + do s=1,N_states + tmp(s,j) = tmp(s,j) + sij*active_excitation_to_determinants_val(s,idx,a_coll) + enddo + ok(j) = .True. + endif + enddo + enddo + idx = 0 + do j=1,N_det_non_ref + if (ok(j)) then + idx = idx+1 + mrcc_S2A_ind(idx,a_coll) = j + do s=1,N_states + mrcc_S2A_val(s,idx,a_coll) = tmp(s,j) + enddo + endif + enddo + mrcc_S2A_ind(0,a_coll) = idx + enddo + !$OMP END DO + deallocate(tmp,ok) + !$OMP END PARALLEL + END_PROVIDER BEGIN_PROVIDER [ integer, mrcc_AtA_ind, (N_det_ref * n_exc_active) ] @@ -148,9 +197,8 @@ END_PROVIDER integer :: at_roww, at_row, wk, a_coll, a_col, r1, r2, s 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 + PROVIDE psi_non_ref mrcc_S2A_val mrcc_AtA_ind(:) = 0 mrcc_AtA_val(:,:) = 0.d0 @@ -163,9 +211,9 @@ 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, ts, i, r1, r2, wk, A_ind_mwen, A_val_mwen,& - !$OMP det1,det2,As2_val_mwen, a_coll, at_roww,sij) & - !$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) + !$OMP As2_val_mwen, a_coll, at_roww,sij) & + !$OMP shared(N_states,mrcc_col_shortcut, mrcc_S2A_val, mrcc_S2A_ind, mrcc_N_col, AtA_size, mrcc_AtA_val, mrcc_AtA_ind, & + !$OMP n_exc_active, active_pp_idx,psi_non_ref,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) @@ -177,7 +225,6 @@ 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)) @@ -186,12 +233,25 @@ 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 + end if + end do + + ts(:) = 0d0 + r1 = 1 + r2 = 1 + do while ((active_excitation_to_determinants_idx(r1, at_roww) /= 0).and.(mrcc_S2A_ind(r2, a_coll) /= 0)) + if(active_excitation_to_determinants_idx(r1, at_roww) > mrcc_S2A_ind(r2, a_coll)) then + r2 = r2+1 + else if(active_excitation_to_determinants_idx(r1, at_roww) < mrcc_S2A_ind(r2, a_coll)) then + r1 = r1+1 + else + do s=1,N_states + ts(s) = ts(s) + active_excitation_to_determinants_val(s,r1, at_roww) * mrcc_S2A_val(s,r2, a_coll) enddo r1 = r1+1 r2 = r2+1 @@ -201,10 +261,9 @@ 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 + if(sum(abs(t)+abs(ts)) /= 0.d0) then wk += 1 A_ind_mwen(wk) = a_col do s=1,N_states diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index f864dd08..970802de 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -707,12 +707,24 @@ END_PROVIDER x_new = x - double precision :: s2(N_states), s2_local, dx + double precision :: s2(N_states), s2_local, dx, s2_init(N_states) double precision :: factor, resold factor = 1.d0 resold = huge(1.d0) + + s2_init(s) = S_z2_Sz + do hh = 1, hh_shortcut(0) + do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 + if(is_active_exc(pp)) cycle + do i=mrcc_col_shortcut(a_coll), mrcc_col_shortcut(a_coll) + mrcc_N_col(a_coll) - 1 + s2_init(s) = s2_init(s) + X(pp)*X(mrcc_AtA_ind(i))*mrcc_AtS2A_val(s,i) + end do + enddo + end do + do k=0,100000 - !$OMP PARALLEL default(shared) private(cx, dx, i, j, a_col, a_coll, s2_local) + s2_local = s2_init(s) + !$OMP PARALLEL default(shared) private(cx, dx, i, j, a_col, a_coll) !$OMP DO do i=1,N_det_non_ref @@ -720,7 +732,15 @@ END_PROVIDER enddo !$OMP END DO - s2(s) = 0.d0 + !$OMP DO REDUCTION(+:s2_local) + do a_coll = 1, n_exc_active + a_col = active_pp_idx(a_coll) + do i=mrcc_col_shortcut(a_coll), mrcc_col_shortcut(a_coll) + mrcc_N_col(a_coll) - 1 + s2_local = s2_local + X(a_col)*X(mrcc_AtA_ind(i))*mrcc_AtS2A_val(s,i) + end do + end do + !$OMP END DO + !$OMP DO do a_coll = 1, n_exc_active a_col = active_pp_idx(a_coll) @@ -728,23 +748,19 @@ END_PROVIDER 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) + dx = dx + (s2_local-expected_s2)*x(mrcc_AtA_ind(i)) * mrcc_AtS2A_val(s,i) end do 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 + s2(s) = s2_local res = 0.d0 - do a_coll=1,n_exc_active ! hh_nex + do a_coll=1,n_exc_active a_col = active_pp_idx(a_coll) do j=1,N_det_non_ref i = active_excitation_to_determinants_idx(j,a_coll) @@ -765,7 +781,7 @@ END_PROVIDER if(res < 1d-9) exit end do - + s2(s) = s2_local norm = 0.d0 do i=1,N_det_non_ref @@ -928,6 +944,7 @@ END_PROVIDER norm = norm*f print *, 'norm of |T Psi_0> = ', dsqrt(norm) + print *, 'S^2 |T Psi_0> = ', s2(s) do i=1,N_det_ref norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) @@ -936,6 +953,7 @@ END_PROVIDER do i=1,N_det_non_ref rho_mrcc(i,s) = rho_mrcc(i,s) * f enddo +rho_mrcc = 1.d0 ! rho_mrcc now contains the product of the scaling factors and the ! normalization constant