diff --git a/plugins/MRCC_Utils/amplitudes.irp.f b/plugins/MRCC_Utils/amplitudes.irp.f index 0e6a4cf4..1dcf2a2b 100644 --- a/plugins/MRCC_Utils/amplitudes.irp.f +++ b/plugins/MRCC_Utils/amplitudes.irp.f @@ -23,33 +23,39 @@ allocate(pathTo(N_det_non_ref)) pathTo(:) = 0 - is_active_exc(:) = .false. + is_active_exc(:) = .True. n_exc_active = 0 - do hh = 1, hh_shortcut(0) - do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 - do II = 1, N_det_ref +! do hh = 1, hh_shortcut(0) +! do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 +! do II = 1, N_det_ref +! +! call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int) +! if(.not. ok) cycle +! +! call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int) +! if(.not. ok) cycle +! +! ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int) +! if(ind == -1) cycle +! +! logical, external :: is_a_two_holes_two_particles +! if (is_a_two_holes_two_particles(myDet)) then +! is_active_exc(pp) = .False. +! endif - call apply_hole_local(psi_ref(1,1,II), hh_exists(1, hh), myMask, ok, N_int) - if(.not. ok) cycle +! ind = psi_non_ref_sorted_idx(ind) +! if(pathTo(ind) == 0) then +! pathTo(ind) = pp +! else +! is_active_exc(pp) = .true. +! is_active_exc(pathTo(ind)) = .true. +! end if - call apply_particle_local(myMask, pp_exists(1, pp), myDet, ok, N_int) - if(.not. ok) cycle +! end do +! end do +! end do - ind = searchDet(psi_non_ref_sorted(1,1,1), myDet(1,1), N_det_non_ref, N_int) - if(ind == -1) cycle - - ind = psi_non_ref_sorted_idx(ind) - if(pathTo(ind) == 0) then - pathTo(ind) = pp - else - is_active_exc(pp) = .true. - is_active_exc(pathTo(ind)) = .true. - end if - 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 @@ -66,6 +72,32 @@ END_PROVIDER +BEGIN_PROVIDER [ logical, has_a_unique_parent, (N_det_non_ref) ] + implicit none + BEGIN_DOC + ! True if the determinant in the non-reference has a unique parent + END_DOC + integer :: i,j,n + integer :: degree + do j=1,N_det_non_ref + has_a_unique_parent(j) = .True. + n=0 + do i=1,N_det_ref + call get_excitation_degree(psi_ref(1,1,i), psi_non_ref(1,1,j), degree, N_int) + if (degree < 2) then + n = n+1 + if (n > 1) then + has_a_unique_parent(j) = .False. + exit + endif + endif + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ integer, n_exc_active_sze ] implicit none BEGIN_DOC @@ -95,7 +127,7 @@ END_PROVIDER !$OMP active_excitation_to_determinants_val, active_excitation_to_determinants_idx)& !$OMP shared(hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, & !$OMP psi_non_ref_sorted_idx, psi_ref, N_det_ref, N_states)& - !$OMP shared(is_active_exc, active_hh_idx, active_pp_idx, n_exc_active)& + !$OMP shared(active_hh_idx, active_pp_idx, n_exc_active)& !$OMP private(lref, pp, II, ok, myMask, myDet, ind, phase, wk, ppp, hh, s) allocate(lref(N_det_non_ref)) !$OMP DO schedule(dynamic) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 0958b8d1..41435688 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -693,12 +693,12 @@ END_PROVIDER double precision :: phase - double precision, allocatable :: rho_mrcc_init(:) + double precision, allocatable :: rho_mrcc_inact(:) integer :: a_coll, at_roww print *, "TI", hh_nex, N_det_non_ref - allocate(rho_mrcc_init(N_det_non_ref)) + allocate(rho_mrcc_inact(N_det_non_ref)) allocate(x_new(hh_nex)) allocate(x(hh_nex), AtB(hh_nex)) @@ -710,7 +710,7 @@ END_PROVIDER !$OMP private(at_row, a_col, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)& !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtB, mrcc_AtA_val, mrcc_AtA_ind, s, n_exc_active, active_pp_idx) - !$OMP DO schedule(dynamic, 100) + !$OMP DO schedule(static, 100) do at_roww = 1, n_exc_active ! hh_nex at_row = active_pp_idx(at_roww) do i=1,active_excitation_to_determinants_idx(0,at_roww) @@ -729,7 +729,7 @@ END_PROVIDER X(a_col) = AtB(a_col) end do - rho_mrcc_init = 0d0 + rho_mrcc_inact(:) = 0d0 allocate(lref(N_det_ref)) do hh = 1, hh_shortcut(0) @@ -753,19 +753,15 @@ END_PROVIDER X(pp) = AtB(pp) do II=1,N_det_ref if(lref(II) > 0) then - rho_mrcc_init(lref(II)) = psi_ref_coef(II,s) * X(pp) + rho_mrcc_inact(lref(II)) = psi_ref_coef(II,s) * X(pp) else if(lref(II) < 0) then - rho_mrcc_init(-lref(II)) = -psi_ref_coef(II,s) * X(pp) + rho_mrcc_inact(-lref(II)) = -psi_ref_coef(II,s) * X(pp) end if end do end do end do deallocate(lref) - do i=1,N_det_non_ref - rho_mrcc(i,s) = rho_mrcc_init(i) - enddo - x_new = x double precision :: factor, resold @@ -793,7 +789,10 @@ END_PROVIDER print *, k, res, 1.d0 - res/resold endif - if ( (res < 1d-10).or.(res/resold > 0.99d0) ) then + if ( res < 1d-10 ) then + exit + endif + if ( (res/resold > 0.99d0) ) then exit endif resold = res @@ -802,38 +801,60 @@ END_PROVIDER dIj_unique(1:size(X), s) = X(1:size(X)) print *, k, res, 1.d0 - res/resold - enddo - do s=1,N_states + do i=1,N_det_non_ref + rho_mrcc(i,s) = 0.d0 + enddo 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) if (i==0) exit + if (rho_mrcc_inact(i) /= 0.d0) then + call debug_det(psi_non_ref(1,1,i),N_int) + stop + endif rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * dIj_unique(a_col,s) enddo end do - norm = 0.d0 - do i=1,N_det_non_ref - norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s) - enddo - ! Norm now contains the norm of A.X - + double precision :: norm2_ref, norm2_inact, a, b, c, Delta + ! Psi = Psi_ref + Psi_inactive + f*Psi_active + ! Find f to normalize Psi + + norm2_ref = 0.d0 do i=1,N_det_ref - norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) + norm2_ref = norm2_ref + psi_ref_coef(i,s)*psi_ref_coef(i,s) enddo - ! Norm now contains the norm of Psi + A.X - + + a = 0.d0 + do i=1,N_det_non_ref + a = a + rho_mrcc(i,s)*rho_mrcc(i,s) + enddo + + norm = a + norm2_ref print *, "norm : ", sqrt(norm) - enddo + + norm = sqrt((1.d0-norm2_ref)/a) + + ! Renormalize Psi+A.X + do i=1,N_det_non_ref + rho_mrcc(i,s) = rho_mrcc(i,s) * norm + enddo + +!norm = norm2_ref +!do i=1,N_det_non_ref +! norm = norm + rho_mrcc(i,s)**2 +!enddo +!print *, 'check', norm +!stop + - do s=1,N_states norm = 0.d0 double precision :: f, g, gmax - gmax = 1.d0*maxval(dabs(psi_non_ref_coef(:,s))) + gmax = maxval(dabs(psi_non_ref_coef(:,s))) do i=1,N_det_non_ref if (lambda_type == 2) then f = 1.d0 @@ -845,41 +866,22 @@ END_PROVIDER f = psi_non_ref_coef(i,s) / rho_mrcc(i,s) ! Avoid numerical instabilities -! g = 1.d0+dabs(gmax / psi_non_ref_coef(i,s) ) g = 2.d0+100.d0*exp(-20.d0*dabs(psi_non_ref_coef(i,s)/gmax)) f = min(f, g) f = max(f,-g) + endif - norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s) + norm = norm + (rho_mrcc(i,s)*f)**2 rho_mrcc(i,s) = f enddo - ! norm now contains the norm of |T.Psi_0> - ! rho_mrcc now contains the f factors - - f = 1.d0/norm - ! f now contains 1/ - - norm = 0.d0 - do i=1,N_det_non_ref - norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s) - enddo - ! norm now contains - f = dsqrt(f*norm) - ! f normalises T.Psi_0 such that (1+T)|Psi> is normalized + ! rho_mrcc now contains the mu_i factors print *, 'norm of |T Psi_0> = ', dsqrt(norm) - norm = norm*f - if (dsqrt(norm) > 1.d0) then + if (norm > 1.d0) then stop 'Error : Norm of the SD larger than the norm of the reference.' endif - do i=1,N_det_non_ref - rho_mrcc(i,s) = rho_mrcc(i,s) * f - enddo - ! rho_mrcc now contains the product of the scaling factors and the - ! normalization constant - end do END_PROVIDER diff --git a/plugins/analyze_wf/analyze_wf.irp.f b/plugins/analyze_wf/analyze_wf.irp.f index 061ff940..7d005a05 100644 --- a/plugins/analyze_wf/analyze_wf.irp.f +++ b/plugins/analyze_wf/analyze_wf.irp.f @@ -30,13 +30,13 @@ subroutine run if (occupation(i) > 1.999d0) then class(0,1) += 1 class( class(0,1), 1) = i - else if (occupation(i) > 1.95d0) then + else if (occupation(i) > 1.97d0) then class(0,2) += 1 class( class(0,2), 2) = i else if (occupation(i) < 0.001d0) then class(0,5) += 1 class( class(0,5), 5) = i - else if (occupation(i) < 0.01d0) then + else if (occupation(i) < 0.03d0) then class(0,4) += 1 class( class(0,4), 4) = i else diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 2fd40838..2820750f 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -38,7 +38,7 @@ use bitmasks do p=hh_shortcut(h), hh_shortcut(h+1)-1 call apply_particle_local(mask, pp_exists(1, p), buf(1,1,n), ok, N_int) if(ok) n = n + 1 - if(n > N_det_non_ref) stop "MRCC..." + if(n > N_det_non_ref) stop "Buffer too small in MRCC..." end do n = n - 1 diff --git a/tests/bats/mrcepa0.bats b/tests/bats/mrcepa0.bats index 2420955c..9a62885e 100644 --- a/tests/bats/mrcepa0.bats +++ b/tests/bats/mrcepa0.bats @@ -32,7 +32,7 @@ source $QP_ROOT/tests/bats/common.bats.sh ezfio set mrcepa0 n_it_max_dressed_ci 3 qp_run $EXE $INPUT energy="$(ezfio get mrcepa0 energy_pt2)" - eq $energy -76.2381673136696 2.e-4 + eq $energy -76.2381754078899 1.e-4 } @test "MRSC2 H2O cc-pVDZ" {