LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity

Changes Manus

This commit is contained in:
Anthony Scemama 2023-02-02 14:52:21 +01:00
parent f4de97a71d
commit d3f741baca
9 changed files with 71 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
bin/qp_reset
View File

@ -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

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit 242151e03d1d6bf042387226431d82d35845686a
Subproject commit 9e5b27ce5a174901765cec9db9e7b2aa6170a5de

8
src/ao_many_one_e_ints/list_grid.irp.f
View File

@ -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

155
src/cipsi_tc_bi_ortho/selection.irp.f
View File

@ -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
! <alpha|H|psi_0> = \sum_i c_i <alpha|H|i>
! -------------------------------------------
! Non hermitian
! Non hermitian
! c_alpha = <alpha|H(j)|psi_0>/delta_E(alpha)
! e_alpha = c_alpha * <psi_0|H(j)|alpha>
! <alpha|H|psi_0> and <psi_0|H|alpha>
! <det|H(j)|psi_0> and transpose
! <det|H(j)|psi_0> 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

2
src/csf/sigma_vector.irp.f
View File

@ -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

2
src/fci_tc_bi/copy_wf.irp.f
View File

@ -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

10
src/tc_bi_ortho/tc_h_eigvectors.irp.f
View File

@ -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

10
src/tc_bi_ortho/tc_natorb.irp.f
View File

@ -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

6
src/tc_keywords/EZFIO.cfg
View File

@ -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