diff --git a/plugins/MRCC_Utils/amplitudes.irp.f b/plugins/MRCC_Utils/amplitudes.irp.f index 77360d93..095eebbe 100644 --- a/plugins/MRCC_Utils/amplitudes.irp.f +++ b/plugins/MRCC_Utils/amplitudes.irp.f @@ -135,58 +135,9 @@ END_PROVIDER 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) ] &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 @@ -195,14 +146,13 @@ 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(:), ts(:), A_val_mwen(:,:), As2_val_mwen(:,:) + double precision, allocatable :: t(:), A_val_mwen(:,:), As2_val_mwen(:,:) integer, allocatable :: A_ind_mwen(:) double precision :: sij - PROVIDE psi_non_ref mrcc_S2A_val + PROVIDE psi_non_ref 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 @@ -210,11 +160,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, ts, i, r1, r2, wk, A_ind_mwen, A_val_mwen,& + !$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen,& !$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 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) + allocate(A_val_mwen(N_states,hh_nex), As2_val_mwen(N_states,hh_nex), A_ind_mwen(hh_nex), t(N_states) ) !$OMP DO schedule(dynamic, 100) do at_roww = 1, n_exc_active ! hh_nex @@ -241,36 +191,6 @@ END_PROVIDER 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 - end if - end do - - if(a_col == at_row) then - do s=1,N_states - t(s) = t(s) + 1.d0 - enddo - end if - if(sum(abs(t)+abs(ts)) /= 0.d0) then - wk += 1 - 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 if(wk /= 0) then @@ -285,7 +205,6 @@ 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 @@ -293,7 +212,7 @@ END_PROVIDER end if end do !$OMP END DO NOWAIT - deallocate (A_ind_mwen, A_val_mwen, As2_val_mwen, t, ts) + deallocate (A_ind_mwen, A_val_mwen, As2_val_mwen, t) !$OMP END PARALLEL print *, "ATA SIZE", ata_size diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 970802de..e81fce53 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -707,23 +707,11 @@ END_PROVIDER x_new = x - 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 - s2_local = s2_init(s) !$OMP PARALLEL default(shared) private(cx, dx, i, j, a_col, a_coll) !$OMP DO @@ -732,31 +720,18 @@ END_PROVIDER enddo !$OMP END DO - !$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) 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 + (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 + x_new(a_col) = AtB(a_col) + cx * factor end do !$OMP END DO !$OMP END PARALLEL - s2(s) = s2_local - res = 0.d0 @@ -776,12 +751,11 @@ END_PROVIDER resold = res if(mod(k, 100) == 0) then - print *, "res ", k, res, s2(s) + print *, "res ", k, res end if if(res < 1d-9) exit end do - s2(s) = s2_local norm = 0.d0 do i=1,N_det_non_ref @@ -923,7 +897,6 @@ END_PROVIDER ! Avoid numerical instabilities f = min(f,2.d0) f = max(f,-2.d0) - f = 1.d0 norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s) rho_mrcc(i,s) = f @@ -944,7 +917,6 @@ 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) @@ -953,7 +925,6 @@ 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