From e319cb7f1d4f7761fa223988971983a362018b96 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Wed, 6 Jul 2016 11:28:52 +0200 Subject: [PATCH] ref-dependent amplitudes --- plugins/MRCC_Utils/mrcc_utils.irp.f | 59 +++++++++++++++++++++++------ plugins/mrcepa0/dressing.irp.f | 10 +++-- 2 files changed, 54 insertions(+), 15 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index bf5b7694..513e7d09 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -687,18 +687,19 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] integer, external :: searchDet, unsortedSearchDet integer(bit_kind) :: myDet(N_int, 2), myMask(N_int, 2) integer :: N, INFO, AtA_size, r1, r2 - double precision , allocatable:: AtB(:), AtA_val(:), A_val(:,:), x(:), x_new(:), A_val_mwen(:) + double precision , allocatable:: B(:), AtB(:), AtA_val(:), A_val(:,:), x(:), x_new(:), A_val_mwen(:) double precision :: t, norm, cx integer, allocatable :: A_ind(:,:), lref(:), AtA_ind(:), A_ind_mwen(:), col_shortcut(:), N_col(:) - + if(n_states /= 1) stop "n_states /= 1" + nex = hh_shortcut(hh_shortcut(0)+1)-1 print *, "TI", nex, N_det_non_ref allocate(A_ind(N_det_ref+1, nex), A_val(N_det_ref+1, nex)) allocate(AtA_ind(N_det_ref * nex), AtA_val(N_det_ref * nex)) !!!!! MAY BE TOO SMALL !!!!!!!! allocate(x(nex), AtB(nex)) allocate(A_val_mwen(nex), A_ind_mwen(nex)) - allocate(N_col(nex), col_shortcut(nex)) + allocate(N_col(nex), col_shortcut(nex), B(N_det_non_ref)) A_val = 0d0 A_ind = 0 !$OMP PARALLEL DO schedule(static,10) default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int, A_val, A_ind) & @@ -746,7 +747,7 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind) do at_row = 1, nex wk = 0 - if(mod(at_row, 1000) == 0) print *, "AtA", at_row, "/", nex + if(mod(at_row, 10000) == 0) print *, "AtA", at_row, "/", nex do i=1,N_det_ref if(A_ind(i, at_row) == 0) exit AtB(at_row) = AtB(at_row) + psi_non_ref_coef(A_ind(i, at_row), 1) * A_val(i, at_row) @@ -803,12 +804,12 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] enddo do k=0,100000 - ! df $ fg OMP PARALLEL DO default(shared) + !$OMP PARALLEL DO default(shared) do i=1,nex x_new(i) = AtB(i) enddo - ! sdf $OMP PARALLEL DO default(shared) private(cx, i) + !$OMP PARALLEL DO default(shared) private(cx, i) do a_col = 1, nex cx = 0d0 do i=col_shortcut(a_col), col_shortcut(a_col) + N_col(a_col) - 1 @@ -816,13 +817,13 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] end do x_new(a_col) += cx end do - ! sdf $OMP END PARALLEL DO + !$OMP END PARALLEL DO double precision :: norm_cas norm_cas = 0d0 do i = 1, N_det_ref norm_cas += psi_ref_coef(i,1)**2 end do - + norm = 0d0 t = 0d0 @@ -830,15 +831,37 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] t = t + X_new(j) * X_new(j) end do - x_new = x_new / sqrt(norm_cas + t) - + !t = (1d0 - norm_cas) / t + !x_new = x_new * sqrt(t) + !!!!!!!!!!!!!! + !B = 0d0 + !do i=1, nex + ! do j=1, N_det_ref + ! if(A_ind(j, i) == 0) exit + ! B(A_ind(j, i)) += A_val(j, i) * x(i) + ! end do + !end do + !t = 0d0 + !do i=1, size(B) + ! t += B(i)**2 + !end do + !print *, "NORMT", sqrt(t + norm_cas) + !x_new = x_new / sqrt(t + norm_cas) +!!!!!!!!!! + + t = (1d0 / norm_cas - 1d0) / t + x_new = x_new * sqrt(t) + do j=1, size(X) norm += (X_new(j) - X(j))**2 x(j) = x_new(j) end do - !print *, "NORM X_new", t - if(mod(k, 1000) == 0) print *, "residu ", k, norm + + if(mod(k, 50) == 0) then + print *, "residu ", k, norm, "norm t", sqrt(t) + end if + if(norm < 1d-16) exit end do print *, "CONVERGENCE : ", norm @@ -898,6 +921,18 @@ BEGIN_PROVIDER [ double precision, dIj, (hh_shortcut(hh_shortcut(0)+1)-1) ] END_PROVIDER +double precision function get_dij_index(II, i, Nint) + integer, intent(in) :: II, i, Nint + double precision, external :: get_dij + + if(dabs(psi_ref_coef(II, 1)) > 1d-1) then + get_dij_index = psi_non_ref_coef(i, 1) / psi_ref_coef(II, 1) + else + get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint) + end if +end function + + double precision function get_dij(det1, det2, Nint) use bitmasks implicit none diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 74a572f5..e4b63208 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -86,7 +86,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) integer :: mobiles(2), smallerlist logical, external :: detEq, is_generable - double precision, external :: get_dij + double precision, external :: get_dij, get_dij_index leng = max(N_det_generators, N_det_non_ref) allocate(miniList(Nint, 2, leng), tq(Nint,2,n_selected), idx_minilist(leng), hij_cache(N_det_non_ref)) @@ -225,7 +225,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe do i_state=1,N_states - dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + dIK(i_state) = get_dij_index(i_I, idx_alpha(k_sd), Nint) + !dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) + !dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state) enddo @@ -264,7 +266,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe hIl = hij_mrcc(idx_alpha(l_sd),i_I) ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) do i_state=1,N_states - dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + dka(i_state) = get_dij_index(i_I, idx_alpha(l_sd), N_int) * phase * phase2 + !dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 + !dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2 enddo endif