From 2f1c7c5ce90f1b107992b5b7c0e461cf8a3d1a63 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 18 Oct 2016 23:07:03 +0200 Subject: [PATCH] Small changes in MRCC --- plugins/MRCC_Utils/mrcc_dress.irp.f | 2 +- plugins/MRCC_Utils/mrcc_utils.irp.f | 323 +++++++++++++++++++--------- plugins/mrcepa0/EZFIO.cfg | 2 +- plugins/mrcepa0/dressing.irp.f | 37 ++-- 4 files changed, 241 insertions(+), 123 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index e6d0fb81..5c2f5efc 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -271,7 +271,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge !delta_ii_(i_state,i_I) = 0.d0 do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0 * dIa_hla(i_state,k_sd) enddo endif enddo diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 14885153..f1c4b4a3 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -685,7 +685,7 @@ END_PROVIDER - do s = 1, N_states + do s=1, N_states A_val = 0d0 A_ind = 0 @@ -698,61 +698,61 @@ END_PROVIDER !$OMP PARALLEL default(none) shared(psi_non_ref, hh_exists, pp_exists, N_int, A_val, A_ind)& !$OMP shared(s, hh_shortcut, psi_ref_coef, N_det_non_ref, psi_non_ref_sorted, psi_non_ref_sorted_idx, psi_ref, N_det_ref)& - !$OMP shared(active, active_hh_idx, active_pp_idx, nactive)& + !$OMP shared(active, active_hh_idx, active_pp_idx, nactive) & !$OMP private(lref, pp, II, ok, myMask, myDet, ind, phase, wk, ppp, hh) - allocate(lref(N_det_non_ref)) - !$OMP DO schedule(static,10) - do ppp=1,nactive - pp = active_pp_idx(ppp) - hh = active_hh_idx(ppp) - lref = 0 - 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) then - call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int) - if (phase > 0.d0) then - lref(psi_non_ref_sorted_idx(ind)) = II - else - lref(psi_non_ref_sorted_idx(ind)) = -II - endif - end if - end do - wk = 0 - do i=1, N_det_non_ref - if(lref(i) > 0) then - wk += 1 - A_val(wk, ppp) = psi_ref_coef(lref(i), s) - A_ind(wk, ppp) = i - else if(lref(i) < 0) then - wk += 1 - A_val(wk, ppp) = -psi_ref_coef(-lref(i), s) - A_ind(wk, ppp) = i - end if - end do - A_ind(0,ppp) = wk + allocate(lref(N_det_non_ref)) + !$OMP DO schedule(static,10) + do ppp=1,nactive + pp = active_pp_idx(ppp) + hh = active_hh_idx(ppp) + lref = 0 + 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) then + call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int) + if (phase > 0.d0) then + lref(psi_non_ref_sorted_idx(ind)) = II + else + lref(psi_non_ref_sorted_idx(ind)) = -II + endif + end if end do + wk = 0 + do i=1, N_det_non_ref + if(lref(i) > 0) then + wk += 1 + A_val(wk, ppp) = psi_ref_coef(lref(i), s) + A_ind(wk, ppp) = i + else if(lref(i) < 0) then + wk += 1 + A_val(wk, ppp) = -psi_ref_coef(-lref(i), s) + A_ind(wk, ppp) = i + end if + end do + A_ind(0,ppp) = wk + end do !$OMP END DO deallocate(lref) - !$OMP END PARALLEL - - + !$OMP END PARALLEL + + print *, 'Done building A_val, A_ind' AtA_size = 0 col_shortcut = 0 N_col = 0 - integer :: a_coll, at_roww + integer :: a_coll, at_roww !$OMP PARALLEL default(none) shared(k, psi_non_ref_coef, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref)& !$OMP private(at_row, a_col, t, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)& !$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind, s, nactive, active_pp_idx) allocate(A_val_mwen(nex), A_ind_mwen(nex)) - + !$OMP DO schedule(dynamic, 100) do at_roww = 1, nactive ! nex at_row = active_pp_idx(at_roww) @@ -762,8 +762,8 @@ END_PROVIDER j = active_pp_idx(i) AtB(at_row) = AtB(at_row) + psi_non_ref_coef(A_ind(i, at_roww), s) * A_val(i, at_roww) end do - - do a_coll = 1, nactive + + do a_coll = 1, nactive a_col = active_pp_idx(a_coll) t = 0d0 r1 = 1 @@ -795,12 +795,12 @@ END_PROVIDER col_shortcut(at_roww) = AtA_size+1 N_col(at_roww) = wk if (AtA_size+wk > size(AtA_ind,1)) then - print *, AtA_size+wk , size(AtA_ind,1) - stop 'too small' + print *, AtA_size+wk , size(AtA_ind,1) + stop 'too small' endif do i=1,wk - AtA_ind(AtA_size+i) = A_ind_mwen(i) - AtA_val(AtA_size+i) = A_val_mwen(i) + AtA_ind(AtA_size+i) = A_ind_mwen(i) + AtA_val(AtA_size+i) = A_val_mwen(i) enddo AtA_size += wk !$OMP END CRITICAL @@ -822,41 +822,41 @@ END_PROVIDER rho_mrcc_init = 0d0 allocate(lref(N_det_ref)) - !$OMP PARALLEL DO default(shared) schedule(static, 1) & + !$OMP PARALLEL DO default(shared) schedule(static, 1) & !$OMP private(lref, hh, pp, II, myMask, myDet, ok, ind, phase) do hh = 1, hh_shortcut(0) - do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 - if(active(pp)) cycle - lref = 0 - 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 - ind = psi_non_ref_sorted_idx(ind) - call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int) - X(pp) += psi_ref_coef(II,s)**2 - AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase - lref(II) = ind - if(phase < 0d0) lref(II) = -ind + do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1 + if(active(pp)) cycle + lref = 0 + 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 + ind = psi_non_ref_sorted_idx(ind) + call get_phase(myDet(1,1), psi_ref(1,1,II), phase, N_int) + X(pp) += psi_ref_coef(II,s)**2 + AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase + lref(II) = ind + if(phase < 0d0) lref(II) = -ind + end do + X(pp) = AtB(pp) / X(pp) + do II=1,N_det_ref + if(lref(II) > 0) then + rho_mrcc_init(lref(II),s) = psi_ref_coef(II,s) * X(pp) + else if(lref(II) < 0) then + rho_mrcc_init(-lref(II),s) = -psi_ref_coef(II,s) * X(pp) + end if + end do end do - X(pp) = AtB(pp) / X(pp) - do II=1,N_det_ref - if(lref(II) > 0) then - rho_mrcc_init(lref(II),s) = psi_ref_coef(II,s) * X(pp) - else if(lref(II) < 0) then - rho_mrcc_init(-lref(II),s) = -psi_ref_coef(II,s) * X(pp) - end if - end do - end do end do !$OMP END PARALLEL DO x_new = x - - double precision :: factor, resold + + double precision :: factor, resold factor = 1.d0 resold = huge(1.d0) do k=0,100000 @@ -882,10 +882,10 @@ END_PROVIDER !$OMP END PARALLEL res = 0.d0 - + if (res < resold) then - do a_coll=1,nactive ! nex + do a_coll=1,nactive ! nex a_col = active_pp_idx(a_coll) do j=1,N_det_non_ref i = A_ind(j,a_coll) @@ -894,60 +894,172 @@ END_PROVIDER enddo res = res + (X_new(a_col) - X(a_col))*(X_new(a_col) - X(a_col)) X(a_col) = X_new(a_col) - end do - factor = 1.d0 + end do + factor = 1.d0 else factor = -factor * 0.5d0 endif resold = res - - if(mod(k, 5) == 0) then + + if(mod(k, 100) == 0) then print *, "res ", k, res end if - if(res < 1d-12) exit + 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) - enddo - ! Norm now contains the norm of A.X - - do i=1,N_det_ref - norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) - enddo - ! Norm now contains the norm of Psi + A.X - - print *, k, "res : ", res, "norm : ", sqrt(norm) - - !dIj_unique(:size(X), s) = X(:) + 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 + + do i=1,N_det_ref + norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) + enddo + ! Norm now contains the norm of Psi + A.X + + print *, k, "res : ", res, "norm : ", sqrt(norm) + +!--------------- +! double precision :: e_0, overlap +! double precision, allocatable :: u_0(:) +! integer(bit_kind), allocatable :: keys_tmp(:,:,:) +! allocate (u_0(N_det), keys_tmp(N_int,2,N_det) ) +! k=0 +! overlap = 0.d0 +! do i=1,N_det_ref +! k = k+1 +! u_0(k) = psi_ref_coef(i,1) +! keys_tmp(:,:,k) = psi_ref(:,:,i) +! overlap += u_0(k)*psi_ref_coef(i,1) +! enddo +! norm = 0.d0 +! do i=1,N_det_non_ref +! k = k+1 +! u_0(k) = psi_non_ref_coef(i,1) +! keys_tmp(:,:,k) = psi_non_ref(:,:,i) +! overlap += u_0(k)*psi_non_ref_coef(i,1) +! enddo +! +! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det) +! print *, 'Energy of |Psi_CASSD> : ', e_0 + nuclear_repulsion, overlap +! +! k=0 +! overlap = 0.d0 +! do i=1,N_det_ref +! k = k+1 +! u_0(k) = psi_ref_coef(i,1) +! keys_tmp(:,:,k) = psi_ref(:,:,i) +! overlap += u_0(k)*psi_ref_coef(i,1) +! enddo +! norm = 0.d0 +! do i=1,N_det_non_ref +! k = k+1 +! ! f is such that f.\tilde{c_i} = c_i +! f = psi_non_ref_coef(i,1) / rho_mrcc(i,1) +! +! ! Avoid numerical instabilities +! f = min(f,2.d0) +! f = max(f,-2.d0) +! +! f = 1.d0 +! +! u_0(k) = rho_mrcc(i,1)*f +! keys_tmp(:,:,k) = psi_non_ref(:,:,i) +! norm += u_0(k)**2 +! overlap += u_0(k)*psi_non_ref_coef(i,1) +! enddo +! +! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det) +! print *, 'Energy of |(1+T)Psi_0> : ', e_0 + nuclear_repulsion, overlap +! +! f = 1.d0/norm +! norm = 1.d0 +! do i=1,N_det_ref +! norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s) +! enddo +! f = dsqrt(f*norm) +! overlap = norm +! do i=1,N_det_non_ref +! u_0(k) = rho_mrcc(i,1)*f +! overlap += u_0(k)*psi_non_ref_coef(i,1) +! enddo +! +! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det) +! print *, 'Energy of |(1+T)Psi_0> (normalized) : ', e_0 + nuclear_repulsion, overlap +! +! k=0 +! overlap = 0.d0 +! do i=1,N_det_ref +! k = k+1 +! u_0(k) = psi_ref_coef(i,1) +! keys_tmp(:,:,k) = psi_ref(:,:,i) +! overlap += u_0(k)*psi_ref_coef(i,1) +! enddo +! norm = 0.d0 +! do i=1,N_det_non_ref +! k = k+1 +! ! f is such that f.\tilde{c_i} = c_i +! f = psi_non_ref_coef(i,1) / rho_mrcc(i,1) +! +! ! Avoid numerical instabilities +! f = min(f,2.d0) +! f = max(f,-2.d0) +! +! u_0(k) = rho_mrcc(i,1)*f +! keys_tmp(:,:,k) = psi_non_ref(:,:,i) +! norm += u_0(k)**2 +! overlap += u_0(k)*psi_non_ref_coef(i,1) +! enddo +! +! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det) +! print *, 'Energy of |(1+T)Psi_0> (mu_i): ', e_0 + nuclear_repulsion, overlap +! +! f = 1.d0/norm +! norm = 1.d0 +! do i=1,N_det_ref +! norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s) +! enddo +! overlap = norm +! f = dsqrt(f*norm) +! do i=1,N_det_non_ref +! u_0(k) = rho_mrcc(i,1)*f +! overlap += u_0(k)*psi_non_ref_coef(i,1) +! enddo +! +! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det) +! print *, 'Energy of |(1+T)Psi_0> (normalized mu_i) : ', e_0 + nuclear_repulsion, overlap +! +! deallocate(u_0, keys_tmp) +! +!--------------- norm = 0.d0 - double precision :: f + double precision :: f do i=1,N_det_non_ref if (rho_mrcc(i,s) == 0.d0) then rho_mrcc(i,s) = 1.d-32 endif - + ! f is such that f.\tilde{c_i} = c_i f = psi_non_ref_coef(i,s) / rho_mrcc(i,s) - + ! Avoid numerical instabilities f = min(f,2.d0) f = max(f,-2.d0) - + norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s) 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 = 1.d0 do i=1,N_det_ref norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s) @@ -955,22 +1067,23 @@ END_PROVIDER ! norm now contains f = dsqrt(f*norm) ! f normalises T.Psi_0 such that (1+T)|Psi> is normalized - + norm = norm*f print *, 'norm of |T Psi_0> = ', dsqrt(norm) - + do i=1,N_det_ref norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s) enddo - + 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 - dIj_unique(:size(X), s) = X(:) + dIj_unique(:size(X), s) = X(:) end do + END_PROVIDER diff --git a/plugins/mrcepa0/EZFIO.cfg b/plugins/mrcepa0/EZFIO.cfg index d792390d..61f3392f 100644 --- a/plugins/mrcepa0/EZFIO.cfg +++ b/plugins/mrcepa0/EZFIO.cfg @@ -23,7 +23,7 @@ interface: ezfio type: Threshold doc: Threshold on the convergence of the dressed CI energy interface: ezfio,provider,ocaml -default: 1.e-4 +default: 5.e-5 [n_it_max_dressed_ci] type: Strictly_positive_int diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 8df7e91a..4f355f2b 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -299,7 +299,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe delta_ii_(i_state,i_I) = 0.d0 do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0*dIa_hla(i_state,k_sd) enddo endif enddo @@ -554,7 +554,7 @@ END_PROVIDER do k=1,N_det_non_ref call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) - call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) +! call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) !print *, Hjk * get_dij(psi_ref(1,1,i), psi_non_ref(1,1,k), N_int), Hki * get_dij(psi_ref(1,1,j), psi_non_ref(1,1,k), N_int) @@ -647,7 +647,7 @@ end function integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref) logical :: ok double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) - double precision :: contrib, HIIi, HJk, wall + double precision :: contrib, contrib2, HIIi, HJk, wall integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ integer(bit_kind) :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2) integer(bit_kind),allocatable :: sortRef(:,:,:) @@ -677,7 +677,7 @@ end function delta_mrcepa0_ij(:,:,:) = 0d0 !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii) & - !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib) & + !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2) & !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) & !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) @@ -720,16 +720,18 @@ end function !$OMP ATOMIC notf = notf+1 - call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) - !contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k)) +! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) - !$OMP ATOMIC - delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then + contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) !$OMP ATOMIC - delta_mrcepa0_ii(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) + delta_mrcepa0_ii(J,i_state) -= contrib2 + else + contrib = contrib * 0.5d0 end if + !$OMP ATOMIC + delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib end do kloop end do @@ -753,7 +755,7 @@ END_PROVIDER integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_ logical :: ok double precision :: phase_Ji, phase_Ik, phase_Ii - double precision :: contrib, delta_IJk, HJk, HIk, HIl + double precision :: contrib, contrib2, delta_IJk, HJk, HIk, HIl integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2) integer, allocatable :: idx_sorted_bit(:) @@ -778,7 +780,7 @@ END_PROVIDER !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii) & !$OMP private(i, J, k, degree, degree2, l, deg, ni) & !$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_) & - !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) & + !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib2, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) & !$OMP private(det_tmp, det_tmp2, II, blok) & !$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) & !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb) @@ -827,13 +829,16 @@ END_PROVIDER delta_IJk = HJk * HIk * lambda_mrcc(i_state, k) call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) if(ok) cycle - contrib = delta_IJk * HIl * lambda_mrcc(i_state,l) + contrib = delta_IJk * HIl * lambda_mrcc(i_state,l) + if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then + contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) + !$OMP ATOMIC + delta_sub_ii(II,i_state) -= contrib2 + else + contrib = contrib * 0.5d0 + endif !$OMP ATOMIC delta_sub_ij(II, i, i_state) += contrib - if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then - !$OMP ATOMIC - delta_sub_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) - endif end do end do end do