From d3f741bacaa27dc7ed203f06c454753f32093884 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 2 Feb 2023 14:52:21 +0100 Subject: [PATCH] Changes Manus --- bin/qp_reset | 1 + external/qp2-dependencies | 2 +- src/ao_many_one_e_ints/list_grid.irp.f | 8 +- src/cipsi_tc_bi_ortho/selection.irp.f | 155 +++++++++---------------- src/csf/sigma_vector.irp.f | 2 +- src/fci_tc_bi/copy_wf.irp.f | 2 +- src/tc_bi_ortho/tc_h_eigvectors.irp.f | 10 +- src/tc_bi_ortho/tc_natorb.irp.f | 10 +- src/tc_keywords/EZFIO.cfg | 6 - 9 files changed, 71 insertions(+), 125 deletions(-) diff --git a/bin/qp_reset b/bin/qp_reset index 74dd1f78..d94ab24c 100755 --- a/bin/qp_reset +++ b/bin/qp_reset @@ -105,6 +105,7 @@ if [[ $mos -eq 1 ]] ; then echo "Warning: You will need to re-define the MO classes" fi rm --recursive --force -- ${ezfio}/mo_basis + rm --recursive --force -- ${ezfio}/bi_ortho_mos rm --recursive --force -- ${ezfio}/work/mo_ints_* fi diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 242151e0..9e5b27ce 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 242151e03d1d6bf042387226431d82d35845686a +Subproject commit 9e5b27ce5a174901765cec9db9e7b2aa6170a5de diff --git a/src/ao_many_one_e_ints/list_grid.irp.f b/src/ao_many_one_e_ints/list_grid.irp.f index ccdc33ad..d5d88007 100644 --- a/src/ao_many_one_e_ints/list_grid.irp.f +++ b/src/ao_many_one_e_ints/list_grid.irp.f @@ -8,8 +8,8 @@ thr = 1.d-11 print*,' expo_good_j_mu_1gauss = ',expo_good_j_mu_1gauss !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, r, overlap, thr,fact_gauss, alpha, center,dist,sigma,center_ij) & - !$OMP SHARED (n_points_final_grid, ao_num, ao_overlap_abs_grid,n_pts_grid_ao_prod,expo_good_j_mu_1gauss,& + !$OMP PRIVATE (ipoint, i, j, r, overlap, fact_gauss, alpha, center,dist,sigma,center_ij) & + !$OMP SHARED (n_points_final_grid, ao_num, thr, ao_overlap_abs_grid,n_pts_grid_ao_prod,expo_good_j_mu_1gauss,& !$OMP final_grid_points,ao_prod_center,ao_prod_sigma,ao_nucl) !$OMP DO do i = 1, ao_num @@ -55,5 +55,5 @@ do i = 1, ao_num list(i) = maxval(n_pts_grid_ao_prod(:,i)) enddo - max_n_pts_grid_ao_prod = maxval(list) -END_PROVIDER + max_n_pts_grid_ao_prod = maxval(list) +END_PROVIDER diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 945226de..b293946a 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -681,7 +681,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! endif !endif - ! MANU: ERREUR dans les calculs puisque < I | H | J > = 0 + ! MANU: ERREUR dans les calculs puisque < I | H | J > = 0 ! n'implique pas < I | H_TC | J > = 0 ?? !val = maxval(abs(mat(1:N_states, p1, p2))) !if( val == 0d0) cycle @@ -764,118 +764,69 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! = \sum_i c_i ! ------------------------------------------- - ! Non hermitian + ! Non hermitian ! c_alpha = /delta_E(alpha) ! e_alpha = c_alpha * ! and - ! and transpose + ! and transpose ! ------------------------------------------- - istate = 1 ! call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii) - double precision :: hmono, htwoe, hthree + double precision :: hmono, htwoe, hthree call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii) - delta_E = E0(istate) - Hii + E_shift - - double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp, error - if(debug_tc_pt2 == 1)then !! Using the old version - psi_h_alpha = 0.d0 - alpha_h_psi = 0.d0 - do iii = 1, N_det - call htilde_mu_mat_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha) - call htilde_mu_mat_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i) - psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function - alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function - enddo - else if(debug_tc_pt2 == 2)then !! debugging the new version - psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version - alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version - psi_h_alpha = 0.d0 - alpha_h_psi = 0.d0 - do iii = 1, N_det ! old version - call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) - call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) - psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function - alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function -! psi_h_alpha += i_h_alpha * 1.d0 ! left function -! alpha_h_psi += alpha_h_i * 1.d0 ! right function - enddo - if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then - error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi) - if(error.gt.1.d-2)then - print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E - endif -! if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-08 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-08)then -! call debug_det(det,N_int) -! print*,'psi_h_alpha,alpha_h_psi' -! print*,psi_h_alpha,alpha_h_psi -! print*,psi_h_alpha_tmp,alpha_h_psi_tmp -! print*,dabs(psi_h_alpha - psi_h_alpha_tmp),dabs(alpha_h_psi - alpha_h_psi_tmp) -! alpha_h_psi = 0.d0 -! psi_h_alpha = 0.d0 -! do iii = 1, N_det -! -! call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int) -! call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) -! call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) -! alpha_h_psi += alpha_h_i -! psi_h_alpha += i_h_alpha -! if(dabs(i_h_alpha).gt.1.d-10.or.dabs(alpha_h_i).gt.1.d-10)then -! call debug_det(psi_det(1,1,iii),N_int) -! print*,iii,degree,i_h_alpha,alpha_h_i -! print*,psi_h_alpha,alpha_h_psi -! print*,leigvec_tc_bi_orth(iii,1),reigvec_tc_bi_orth(iii,1) -! endif -! enddo -! stop -! endif - endif - else - psi_h_alpha = mat_l(istate, p1, p2) - alpha_h_psi = mat_r(istate, p1, p2) - endif - - !if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then - ! print*, 'E0,Hii,E_shift' - ! print*, E0(istate), Hii, E_shift - ! print*, psi_h_alpha, alpha_h_psi, delta_E - ! print*, psi_h_alpha * alpha_h_psi / delta_E - ! !if(Hii .lt. E0(istate)) then - ! ! call debug_det(det, N_int) - ! ! print*, ' |E0| < |Hii| !!!' - ! ! print*, ' E0 = ', E0(istate) - ! ! print*, ' Hii = ', Hii - ! !endif - !endif - - coef(istate) = alpha_h_psi / delta_E - e_pert(istate) = coef(istate) * psi_h_alpha - if(selection_tc == 1 )then - if(e_pert(istate).lt.0.d0)then - e_pert(istate) = 0.d0 - endif - else if(selection_tc == -1)then - if(e_pert(istate).gt.0.d0)then - e_pert(istate) = 0.d0 - endif - endif - - - !if(e_pert(istate) .gt. 1.d-15) then - ! print*, 'E0,Hii,E_shift' - ! print*, E0(istate), Hii, E_shift - ! print*, psi_h_alpha, alpha_h_psi, delta_E - ! print*, psi_h_alpha*alpha_h_psi/delta_E - !endif - -! elseif(cipsi_tc == "h_tc_2x2") then + do istate = 1,N_states + delta_E = E0(istate) - Hii + E_shift + double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp, error + if(debug_tc_pt2 == 1)then !! Using the old version + psi_h_alpha = 0.d0 + alpha_h_psi = 0.d0 + do iii = 1, N_det + call htilde_mu_mat_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha) + call htilde_mu_mat_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i) + psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function + alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function + enddo + else if(debug_tc_pt2 == 2)then !! debugging the new version + psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version + alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version + psi_h_alpha = 0.d0 + alpha_h_psi = 0.d0 + do iii = 1, N_det ! old version + call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) + call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) + psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function + alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function + enddo + if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then + error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi) + if(error.gt.1.d-2)then + print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E + endif + endif + else + psi_h_alpha = mat_l(istate, p1, p2) + alpha_h_psi = mat_r(istate, p1, p2) + endif + coef(istate) = alpha_h_psi / delta_E + e_pert(istate) = coef(istate) * psi_h_alpha +! if(selection_tc == 1 )then +! if(e_pert(istate).lt.0.d0)then +! e_pert(istate) = 0.d0 +! endif +! else if(selection_tc == -1)then +! if(e_pert(istate).gt.0.d0)then +! e_pert(istate) = 0.d0 +! endif +! endif + enddo do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1 do istate = 1, N_states - alpha_h_psi = mat(istate, p1, p2) + alpha_h_psi = mat_r(istate, p1, p2) + psi_h_alpha = mat_l(istate, p1, p2) pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate) pt2_data % variance(istate) = pt2_data % variance(istate) + dabs(e_pert(istate)) @@ -885,10 +836,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d case(5) ! Variance selection if (h0_type == 'CFG') then - w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) & + w = min(w, - psi_h_alpha * alpha_h_psi * s_weight(istate,istate)) & / c0_weight(istate) else - w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) + w = min(w, - psi_h_alpha * alpha_h_psi * s_weight(istate,istate)) endif case(6) if (h0_type == 'CFG') then diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index c6644c4c..82b720a4 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -236,7 +236,7 @@ end subroutine get_phase_qp_to_cfg ! initialization psi_coef_config = 0.d0 DetToCSFTransformationMatrix(0,:,:) = 1.d0 - do i = 2-iand(MS,1_8), NSOMOMax,2 + do i = 2-int(iand(MS,1_8),4), NSOMOMax,2 Isomo = IBSET(0_8, i) - 1_8 ! rows = Ncsfs ! cols = Ndets diff --git a/src/fci_tc_bi/copy_wf.irp.f b/src/fci_tc_bi/copy_wf.irp.f index cdb1ead8..c46f20d2 100644 --- a/src/fci_tc_bi/copy_wf.irp.f +++ b/src/fci_tc_bi/copy_wf.irp.f @@ -100,7 +100,7 @@ subroutine copy_H_apply_buffer_to_wf_tc logical :: found_duplicates call remove_duplicates_in_psi_det_tc(found_duplicates) - call bi_normalize(psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,size(psi_l_coef_bi_ortho,1),N_det,N_states) + call bi_normalize(psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,N_det,size(psi_l_coef_bi_ortho,1),N_states) SOFT_TOUCH N_det psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho end diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f index 651088be..1ab6be7a 100644 --- a/src/tc_bi_ortho/tc_h_eigvectors.irp.f +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -100,33 +100,33 @@ print*,'Computing the left-eigenvector ' vec_tmp = 0.d0 do istate = 1, N_states - vec_tmp(:,istate) = psi_l_coef_bi_ortho(:,istate) + vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate) enddo do istate = N_states+1, n_states_diag vec_tmp(istate,istate) = 1.d0 enddo call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav) do istate = 1, N_states - leigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate) + leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo print*,'Computing the right-eigenvector ' !!!! Preparing the right-eigenvector vec_tmp = 0.d0 do istate = 1, N_states - vec_tmp(:,istate) = psi_r_coef_bi_ortho(:,istate) + vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate) enddo do istate = N_states+1, n_states_diag vec_tmp(istate,istate) = 1.d0 enddo call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav) do istate = 1, N_states - reigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate) + reigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate) enddo deallocate(H_jj) endif - call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,N_det,N_det,N_states) + call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,size(reigvec_tc_bi_orth,1),N_det,N_states) print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) norm_ground_left_right_bi_orth = 0.d0 do j = 1, N_det diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index b2ccb091..33410570 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -14,7 +14,7 @@ END_DOC implicit none - integer :: i, j, k, n_real + integer :: i, j, k double precision :: thr_d, thr_nd, thr_deg, accu double precision :: accu_d, accu_nd double precision, allocatable :: dm_tmp(:,:), fock_diag(:) @@ -36,17 +36,17 @@ , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval) ! call non_hrmt_bieig( mo_num, dm_tmp& ! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& -! , n_real, natorb_tc_eigval ) +! , mo_num, natorb_tc_eigval ) accu = 0.d0 - do i = 1, n_real + do i = 1, mo_num print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i) accu += -natorb_tc_eigval(i) enddo print *, ' accu = ', accu dm_tmp = 0.d0 - do i = 1, n_real + do i = 1, mo_num accu = 0.d0 do k = 1, mo_num accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i) @@ -54,7 +54,7 @@ accu = 1.d0/dsqrt(dabs(accu)) natorb_tc_reigvec_mo(:,i) *= accu natorb_tc_leigvec_mo(:,i) *= accu - do j = 1, n_real + do j = 1, mo_num do k = 1, mo_num dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j) enddo diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 24404b13..dee13b25 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -94,12 +94,6 @@ doc: Maximum number of SCF iterations interface: ezfio,provider,ocaml default: 100 -[selection_tc] -type: integer -doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative -interface: ezfio,provider,ocaml -default: 0 - [j1b_pen] type: double precision doc: exponents of the 1-body Jastrow