diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 374b003d..043bc5ca 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -5,10 +5,10 @@ use bitmasks END_PROVIDER + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ] &BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] &BEGIN_PROVIDER [ integer, lambda_mrcc_kept, (0:psi_det_size) ] -&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states, N_det_non_ref) ] implicit none BEGIN_DOC ! cm/ or perturbative 1/Delta_E(m) @@ -16,7 +16,7 @@ END_PROVIDER integer :: i,k double precision :: ihpsi_current(N_states) integer :: i_pert_count - double precision :: hii, E2(N_states), E2var(N_states) + double precision :: hii, lambda_pert integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3 i_pert_count = 0 @@ -26,8 +26,6 @@ END_PROVIDER lambda_mrcc_pt2(0) = 0 lambda_mrcc_kept(0) = 0 - E2 = 0.d0 - E2var = 0.d0 do i=1,N_det_non_ref call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& size(psi_ref_coef,1), N_states,ihpsi_current) @@ -36,51 +34,24 @@ END_PROVIDER if (ihpsi_current(k) == 0.d0) then ihpsi_current(k) = 1.d-32 endif - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - E2(k) += ihpsi_current(k)*ihpsi_current(k) / (psi_ref_energy_diagonalized(k)-hii) - E2var(k) += ihpsi_current(k) * psi_non_ref_coef(i,k) - enddo - enddo - - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& - size(psi_ref_coef,1), N_states,ihpsi_current) - call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 - endif - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) * E2var(k)/E2(k) -! lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) ) -! if ( lambda_pert / lambda_mrcc(k,i) < 0.5d0) then -! print *, ' xxx ', dabs(lambda_pert - lambda_mrcc(k,i)) * ihpsi_current(k) -! if ( dabs( (lambda_pert - lambda_mrcc(k,i)) * ihpsi_current(k) ) > 1.d-3) then -! print *, 'lambda mrcc', lambda_mrcc(k,i) -! print *, 'lambda pert', lambda_pert -! print *, 'coef: ', psi_non_ref_coef(i,k) -! call debug_det(psi_non_ref(1,1,i), N_int) -! call i_H_j(psi_ref(1,1,1),psi_non_ref(1,1,i),N_int,hii) -! print *, hii -! call i_H_j(psi_ref(1,1,2),psi_non_ref(1,1,i),N_int,hii) -! print *, hii -! print *, '---' +! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) + lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) ) + lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) + if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then ! Ignore lamdba -! i_pert_count += 1 -! lambda_mrcc(k,i) = 0.d0 -! lambda_mrcc(k,i) = lambda_pert * E2var(k)/E2(k) -! if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then -! N_lambda_mrcc_pt2 += 1 -! lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i -! endif -! else -! ! Keep lamdba -! if (lambda_mrcc_kept(N_lambda_mrcc_pt3) /= i) then -! N_lambda_mrcc_pt3 += 1 -! lambda_mrcc_kept(N_lambda_mrcc_pt3) = i -! endif -! endif + i_pert_count += 1 + lambda_mrcc(k,i) = 0.d0 + if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then + N_lambda_mrcc_pt2 += 1 + lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i + endif + else + ! Keep lamdba + if (lambda_mrcc_kept(N_lambda_mrcc_pt3) /= i) then + N_lambda_mrcc_pt3 += 1 + lambda_mrcc_kept(N_lambda_mrcc_pt3) = i + endif + endif enddo enddo lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 @@ -92,6 +63,65 @@ END_PROVIDER END_PROVIDER +! BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states, N_det_non_ref) ] +!&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] +!&BEGIN_PROVIDER [ integer, lambda_mrcc_kept, (0:psi_det_size) ] +!&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states, N_det_non_ref) ] +! implicit none +! BEGIN_DOC +! ! cm/ or perturbative 1/Delta_E(m) +! END_DOC +! integer :: i,k +! double precision :: ihpsi_current(N_states) +! integer :: i_pert_count +! double precision :: hii, E2(N_states), E2var(N_states) +! integer :: N_lambda_mrcc_pt2, N_lambda_mrcc_pt3 +! +! i_pert_count = 0 +! lambda_mrcc = 0.d0 +! N_lambda_mrcc_pt2 = 0 +! N_lambda_mrcc_pt3 = 0 +! lambda_mrcc_pt2(0) = 0 +! lambda_mrcc_kept(0) = 0 +! +! E2 = 0.d0 +! E2var = 0.d0 +! do i=1,N_det_non_ref +! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& +! size(psi_ref_coef,1), N_states,ihpsi_current) +! call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) +! do k=1,N_states +! if (ihpsi_current(k) == 0.d0) then +! ihpsi_current(k) = 1.d-32 +! endif +! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) +! lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) +! E2(k) += ihpsi_current(k)*ihpsi_current(k) / (psi_ref_energy_diagonalized(k)-hii) +! E2var(k) += ihpsi_current(k) * psi_non_ref_coef(i,k) +! enddo +! enddo +! +! do i=1,N_det_non_ref +! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& +! size(psi_ref_coef,1), N_states,ihpsi_current) +! call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) +! do k=1,N_states +! if (ihpsi_current(k) == 0.d0) then +! ihpsi_current(k) = 1.d-32 +! endif +! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) +! lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) * E2var(k)/E2(k) +! enddo +! enddo +! lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 +! lambda_mrcc_kept(0) = N_lambda_mrcc_pt3 +! print*,'N_det_non_ref = ',N_det_non_ref +! print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) +! print*,'lambda max = ',maxval(dabs(lambda_mrcc)) + print*,'Number of ignored determinants = ',i_pert_count + +END_PROVIDER + BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] @@ -814,8 +844,8 @@ END_PROVIDER f = psi_non_ref_coef(i,s) / rho_mrcc(i,s) ! Avoid numerical instabilities -! f = min(f,2.d0) -! f = max(f,-2.d0) + f = min(f,2.d0) + f = max(f,-2.d0) endif norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s) @@ -869,45 +899,64 @@ BEGIN_PROVIDER [ double precision, dij, (N_det_ref, N_det_non_ref, N_states) ] dij(j, i, s) = get_dij_index(j, i, s, N_int) end do end do - end do - print *, "done computing amplitudes" -END_PROVIDER +! end do +! print *, "done computing amplitudes" +!END_PROVIDER -double precision function f_fit(x) - implicit none - double precision :: x - f_fit = 0.d0 - return - if (x < 0.d0) then - f_fit = 0.d0 - else if (x < 1.d0) then - f_fit = 1.d0/0.367879441171442 * ( x**2 * exp(-x**2)) - else - f_fit = 1.d0 - endif -end +!double precision function f_fit(x) +! implicit none +! double precision :: x +! f_fit = 0.d0 +! return +! if (x < 0.d0) then +! f_fit = 0.d0 +! else if (x < 1.d0) then +! f_fit = 1.d0/0.367879441171442 * ( x**2 * exp(-x**2)) +! else +! f_fit = 1.d0 +! endif +!end +! +!double precision function get_dij_index(II, i, s, Nint) +! integer, intent(in) :: II, i, s, Nint +! double precision, external :: get_dij +! double precision :: HIi, phase, c, a, b, d +! +! call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) +! call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) +! +! a = lambda_pert(s,i) +! b = lambda_mrcc(s,i) +! c = f_fit(a/b) +! +! d = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase* rho_mrcc(i,s) +! +! c = f_fit(a*HIi/d) +! +! get_dij_index = HIi * a * c + (1.d0 - c) * d +! get_dij_index = d +! return +! +! if(lambda_type == 0) then +! call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) +! get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase +! get_dij_index = get_dij_index * rho_mrcc(i,s) +! else if(lambda_type == 1) then +! call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) +! get_dij_index = HIi * lambda_mrcc(s, i) +! else if(lambda_type == 2) then +! call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) +! get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase +! get_dij_index = get_dij_index * rho_mrcc(i,s) +! end if +!end function double precision function get_dij_index(II, i, s, Nint) integer, intent(in) :: II, i, s, Nint double precision, external :: get_dij - double precision :: HIi, phase, c, a, b, d - - call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) - call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) - - a = lambda_pert(s,i) -! b = lambda_mrcc(s,i) -! c = f_fit(a/b) - - d = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase* rho_mrcc(i,s) - - c = f_fit(a*HIi/d) - -! get_dij_index = HIi * a * c + (1.d0 - c) * d - get_dij_index = d - return + double precision :: HIi, phase if(lambda_type == 0) then call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)