diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 9d5e8a67..d8b0a2c3 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -628,7 +628,7 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz integer :: k_pairs, kl integer :: iter2 - double precision, allocatable :: W(:,:), U(:,:), S(:,:) + double precision, allocatable :: W(:,:), U(:,:), S(:,:), overlap(:,:) double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:) double precision, allocatable :: c(:), s_(:,:), s_tmp(:,:) double precision :: diag_h_mat_elem @@ -688,16 +688,17 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz itermax = min(davidson_sze_max, sze/N_st_diag) allocate( & - W(sze_8,N_st_diag*itermax), & - U(sze_8,N_st_diag*itermax), & - S(sze_8,N_st_diag*itermax), & - h(N_st_diag*itermax,N_st_diag*itermax), & - y(N_st_diag*itermax,N_st_diag*itermax), & - s_(N_st_diag*itermax,N_st_diag*itermax), & - s_tmp(N_st_diag*itermax,N_st_diag*itermax), & + W(sze_8,N_st_diag*itermax), & + U(sze_8,N_st_diag*itermax), & + S(sze_8,N_st_diag*itermax), & + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + s_(N_st_diag*itermax,N_st_diag*itermax), & + s_tmp(N_st_diag*itermax,N_st_diag*itermax), & residual_norm(N_st_diag), & - c(N_st_diag*itermax), & - s2(N_st_diag*itermax), & + c(N_st_diag*itermax), & + s2(N_st_diag*itermax), & + overlap(N_st_diag*itermax,N_st_diag*itermax), & lambda(N_st_diag*itermax)) h = 0.d0 @@ -795,26 +796,39 @@ subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sz s2(k) = s_(k,k) + S_z2_Sz enddo - if (s2_eig) then - logical :: state_ok(N_st_diag*davidson_sze_max) - do k=1,shift2 - state_ok(k) = (dabs(s2(k)-expected_s2) < 0.6d0) + ! Compute overlap with U_in + ! ------------------------- + + integer :: coord(2), order(N_st) + overlap = -1.d0 + do k=1,N_st + do i=1,shift2 + overlap(i,k) = dabs(y(i,k)) enddo - do k=1,shift2 - if (.not. state_ok(k)) then - do l=k+1,shift2 - if (state_ok(l)) then - call dswap(shift2, y(1,k), 1, y(1,l), 1) - call dswap(1, s2(k), 1, s2(l), 1) - call dswap(1, lambda(k), 1, lambda(l), 1) - state_ok(k) = .True. - state_ok(l) = .False. - exit - endif - enddo - endif - enddo - endif + enddo + do k=1,N_st + coord = maxloc(overlap) + order( coord(2) ) = coord(1) + overlap(coord(1),coord(2)) = -1.d0 + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + overlap(k,2) = s2(k) + enddo + do k=1,N_st + l = order(k) + if (k /= l) then + lambda(k) = overlap(l,1) + s2(k) = overlap(l,2) + endif + enddo ! Express eigenvectors of h in the determinant basis diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 281b6760..0735980d 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -149,7 +149,7 @@ END_PROVIDER allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)), & eigenvalues(size(CI_electronic_energy_dressed,1))) - do mrcc_state=N_states,1,-1 + do mrcc_state=1,N_states do j=1,min(N_states,N_det) do i=1,N_det eigenvectors(i,j) = psi_coef(i,j) @@ -161,10 +161,12 @@ END_PROVIDER output_determinants,mrcc_state) CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) - enddo - do mrcc_state=N_states+1,N_states_diag - CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) - CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) + if (mrcc_state == 1) then + do mrcc_state=N_states+1,N_states_diag + CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) + CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) + enddo + endif enddo call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& N_states_diag,size(CI_eigenvectors_dressed,1)) @@ -685,7 +687,6 @@ END_PROVIDER if(is_active_exc(pp)) cycle lref = 0 AtB(pp) = 0.d0 - X(pp) = 0.d0 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 @@ -695,7 +696,6 @@ END_PROVIDER 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) = X(pp) + psi_ref_coef(II,s)*psi_ref_coef(II,s) AtB(pp) += psi_non_ref_coef(ind, s) * psi_ref_coef(II, s) * phase lref(II) = ind if(phase < 0.d0) lref(II) = -ind