PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-11 22:03:55 +02:00
Code Releases Activity

Merge branch 'eginer-dev' into dev

Browse Source
This commit is contained in:
Anthony Scemama 2023-01-26 14:06:19 +01:00
commit 7a92c7aedd
19 changed files with 1248 additions and 301 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0
Subproject commit 242151e03d1d6bf042387226431d82d35845686a

23
src/ao_one_e_ints/EZFIO.cfg
View File

@ -106,26 +106,3 @@ interface: ezfio,provider,ocaml
default: 1.e-15
ezfio_name: threshold_ao
[n_pts_charge]
type: integer
doc: Number of point charges to be added to the potential
interface: ezfio
default: 0
[pts_charge_z]
type: double precision
doc: Charge associated to each point charge
interface: ezfio
size: (ao_one_e_ints.n_pts_charge)
[pts_charge_coord]
type: double precision
doc: Coordinate of each point charge.
interface: ezfio
size: (ao_one_e_ints.n_pts_charge,3)
[point_charges]
type: logical
doc: If |true|, point charges (see ao_one_e_ints/write_pt_charges.py) are added to the one-electron potential
interface: ezfio,provider,ocaml
default: False

108
src/ao_one_e_ints/pot_pt_charges.irp.f Normal file
View File

@ -0,0 +1,108 @@
BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)]
BEGIN_DOC
! Point charge-electron interaction, in the |AO| basis set.
!
! :math:`\langle \chi_i | -\sum_charge charge * \frac{1}{|r-R_charge|} | \chi_j \rangle`
!
! Notice the minus sign convention as it is supposed to be for electrons.
END_DOC
implicit none
integer :: num_A, num_B, power_A(3), power_B(3)
integer :: i, j, k, l, n_pt_in, m
double precision :: alpha, beta
double precision :: A_center(3),B_center(3),C_center(3)
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
ao_integrals_pt_chrg = 0.d0
! if (read_ao_integrals_pt_chrg) then
!
! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg)
! print *, 'AO N-e integrals read from disk'
!
! else
! if(use_cosgtos) then
! !print *, " use_cosgtos for ao_integrals_pt_chrg ?", use_cosgtos
!
! do j = 1, ao_num
! do i = 1, ao_num
! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cosgtos(i,j)
! enddo
! enddo
!
! else
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,&
!$OMP num_A,num_B,Z,c,c1,n_pt_in) &
!$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,&
!$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z)
n_pt_in = n_pt_max_integrals
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3)= ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3)= ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
double precision :: c, c1
c = 0.d0
do k = 1, n_pts_charge
double precision :: Z
Z = pts_charge_z(k)
C_center(1:3) = pts_charge_coord(k,1:3)
c1 = NAI_pol_mult( A_center, B_center, power_A, power_B &
, alpha, beta, C_center, n_pt_in )
c = c - Z * c1
enddo
ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
! endif
! IF(do_pseudo) THEN
! ao_integrals_pt_chrg += ao_pseudo_integrals
! ENDIF
! endif
! if (write_ao_integrals_pt_chrg) then
! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg)
! print *, 'AO N-e integrals written to disk'
! endif
END_PROVIDER

2
src/cas_based_on_top/two_body_dens_rout.irp.f
View File

@ -132,7 +132,7 @@ end
subroutine give_n2_cas(r1,r2,istate,n2_psi)
implicit none
BEGIN_DOC
! returns mu(r), f_psi, n2_psi for a general cas wave function
! returns n2_psi for a general cas wave function
END_DOC
integer, intent(in) :: istate
double precision, intent(in) :: r1(3),r2(3)

612
src/cipsi_tc_bi_ortho/get_d.irp.f
View File

File diff suppressed because it is too large Load Diff

2
src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f
View File

@ -130,7 +130,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted_tc
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc psi_det_sorted_tc
PROVIDE psi_det_hii selection_weight pseudo_sym
PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
PROVIDE excitation_beta_max excitation_alpha_max excitation_max

2
src/cipsi_tc_bi_ortho/run_selection_slave.irp.f
View File

@ -23,7 +23,7 @@ subroutine run_selection_slave(thread, iproc, energy)
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym
PROVIDE psi_selectors_coef_transp psi_det_sorted_tc weight_selection
PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc weight_selection
call pt2_alloc(pt2_data,N_states)

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

@ -81,7 +81,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_tc_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp_tc
PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
PROVIDE banned_excitation
@ -511,7 +511,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
integer(bit_kind) :: phasemask(N_int,2)
PROVIDE psi_selectors_coef_transp psi_det_sorted_tc
PROVIDE psi_selectors_coef_transp_tc psi_det_sorted_tc
PROVIDE psi_selectors_rcoef_bi_orth_transp psi_selectors_lcoef_bi_orth_transp
@ -564,29 +564,30 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp (1, interesting(i)) )
call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_m, mat_p, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
, psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
do j=2,N_int
perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
end do
! call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) )
! call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_m, mat_p, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
! , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
!perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
!perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
!do j=2,N_int
! perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
! perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
!end do
!call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
!call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
!call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
!if(nt == 4) then
! call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
!elseif(nt == 3) then
! call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat , mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
!else
! call get_d0 (det(1,1,i), phasemask, bannedOrb, banned, mat , mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
!endif
call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
if(nt == 4) then
call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i)))
elseif(nt == 3) then
call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i)))
else
call get_d0 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i)))
endif
elseif(nt == 4) then
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
@ -775,17 +776,63 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! call get_excitation_degree( HF_bitmask, det, degree, N_int)
! psi_h_alpha = mat_m(istate, p1, p2)
! alpha_h_psi = mat_p(istate, p1, p2)
double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp
psi_h_alpha_tmp = mat_m(istate, p1, p2)
alpha_h_psi_tmp = mat_p(istate, p1, p2)
!
psi_h_alpha = 0.d0
alpha_h_psi = 0.d0
do iii = 1, N_det
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)
psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1)
alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1)
! psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1)
! alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1)
psi_h_alpha += i_h_alpha * 1.d0
alpha_h_psi += alpha_h_i * 1.d0
enddo
! print*,'---',p1,p2
! call debug_det(det,N_int)
! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi
! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp
! if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-10 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-10)then
! if(dabs(psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d+10)then
if(dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d-10)then
! print*,'---'
! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi
! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp
call debug_det(det,N_int)
print*,dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp),psi_h_alpha *alpha_h_psi,psi_h_alpha_tmp*alpha_h_psi_tmp
print*,'-- Good '
print*, psi_h_alpha, alpha_h_psi
print*,'-- bad '
print*,psi_h_alpha_tmp,alpha_h_psi_tmp
print*,'-- details good'
double precision :: accu_1, accu_2
accu_1 = 0.d0
accu_2 = 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)
print*,iii,degree,i_h_alpha,alpha_h_i
accu_1 += i_h_alpha
accu_2 += alpha_h_i
print*,accu_1,accu_2
enddo
! if(dabs(psi_h_alpha*alpha_h_psi).gt.1.d-10)then
! print*,p1,p2
! print*,det(1,1), det(1,2)
! call debug_det(det,N_int)
! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi
! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp
! print*, dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp),&
! psi_h_alpha *alpha_h_psi,psi_h_alpha_tmp*alpha_h_psi_tmp
stop
endif
! endif
! stop
! endif
!if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then
! print*, 'E0,Hii,E_shift'

8
src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
View File

@ -44,9 +44,10 @@ subroutine run_stochastic_cipsi
pt2_data % overlap= 0.d0
pt2_data % variance = huge(1.e0)
if (s2_eig) then
call make_s2_eigenfunction
endif
!!!! WARNING !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION
! if (s2_eig) then
! call make_s2_eigenfunction
! endif
print_pt2 = .False.
call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
! call routine_save_right
@ -89,6 +90,7 @@ subroutine run_stochastic_cipsi
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
stop
N_iter += 1

34
src/determinants/fock_diag.irp.f
View File

@ -33,59 +33,59 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint)
! Occupied MOs
do ii=1,elec_alpha_num
i = occ(ii,1)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i)
E0 = E0 + mo_one_e_integrals(i,i)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_one_e(i,i)
E0 = E0 + mo_bi_ortho_tc_one_e(i,i)
do jj=1,elec_alpha_num
j = occ(jj,1)
if (i==j) cycle
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j)
E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j)
E0 = E0 + 0.5d0*mo_bi_ortho_tc_two_e_jj_anti(i,j)
enddo
do jj=1,elec_beta_num
j = occ(jj,2)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j)
E0 = E0 + mo_two_e_integrals_jj(i,j)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj(i,j)
E0 = E0 + mo_bi_ortho_tc_two_e_jj(i,j)
enddo
enddo
do ii=1,elec_beta_num
i = occ(ii,2)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i)
E0 = E0 + mo_one_e_integrals(i,i)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_one_e(i,i)
E0 = E0 + mo_bi_ortho_tc_one_e(i,i)
do jj=1,elec_beta_num
j = occ(jj,2)
if (i==j) cycle
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j)
E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j)
E0 = E0 + 0.5d0*mo_bi_ortho_tc_two_e_jj_anti(i,j)
enddo
do jj=1,elec_alpha_num
j = occ(jj,1)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj(i,j)
enddo
enddo
! Virtual MOs
do i=1,mo_num
if (fock_diag_tmp(1,i) /= 0.d0) cycle
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_one_e(i,i)
do jj=1,elec_alpha_num
j = occ(jj,1)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j)
enddo
do jj=1,elec_beta_num
j = occ(jj,2)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j)
fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj(i,j)
enddo
enddo
do i=1,mo_num
if (fock_diag_tmp(2,i) /= 0.d0) cycle
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_one_e(i,i)
do jj=1,elec_beta_num
j = occ(jj,2)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j)
enddo
do jj=1,elec_alpha_num
j = occ(jj,1)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j)
fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj(i,j)
enddo
enddo

10
src/fci_tc_bi/selectors.irp.f
View File

@ -32,6 +32,7 @@ END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef_tc, (psi_selectors_size,2,N_states) ]
implicit none
BEGIN_DOC
! Determinants on which we apply <i|H|psi> for perturbation.
@ -47,12 +48,17 @@ END_PROVIDER
do k=1,N_states
do i=1,N_det_selectors
psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
! psi_selectors_coef_tc(i,1,k) = psi_r_coef_sorted_bi_ortho(i,k)
! psi_selectors_coef_tc(i,2,k) = psi_l_coef_sorted_bi_ortho(i,k)
psi_selectors_coef_tc(i,1,k) = 1.d0
psi_selectors_coef_tc(i,2,k) = 1.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_selectors_size) ]
BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_selectors_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp_tc, (N_states,2,psi_selectors_size) ]
implicit none
BEGIN_DOC
! Transposed psi_selectors
@ -62,6 +68,8 @@ BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_sele
do i=1,N_det_selectors
do k=1,N_states
psi_selectors_coef_transp(k,i) = psi_selectors_coef(i,k)
psi_selectors_coef_transp_tc(k,1,i) = psi_selectors_coef_tc(i,1,k)
psi_selectors_coef_transp_tc(k,2,i) = psi_selectors_coef_tc(i,2,k)
enddo
enddo
END_PROVIDER

24
src/nuclei/EZFIO.cfg
View File

@ -37,3 +37,27 @@ type: logical
doc: If true, the calculation uses periodic boundary conditions
interface: ezfio, provider, ocaml
default: false
[n_pts_charge]
type: integer
doc: Number of point charges to be added to the potential
interface: ezfio
default: 0
[pts_charge_z]
type: double precision
doc: Charge associated to each point charge
interface: ezfio
size: (nuclei.n_pts_charge)
[pts_charge_coord]
type: double precision
doc: Coordinate of each point charge.
interface: ezfio
size: (nuclei.n_pts_charge,3)
[point_charges]
type: logical
doc: If |true|, point charges (see nuclei/write_pt_charges.py) are added to the one-electron potential
interface: ezfio,provider,ocaml
default: False

3
src/nuclei/nuclei.irp.f
View File

@ -205,6 +205,9 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
enddo
enddo
nuclear_repulsion *= 0.5d0
if(point_charges)then
nuclear_repulsion += pt_chrg_nuclei_repulsion + pt_chrg_repulsion
endif
end if
call write_time(6)

209
src/nuclei/point_charges.irp.f Normal file
View File

@ -0,0 +1,209 @@
! ---
BEGIN_PROVIDER [ integer, n_pts_charge ]
implicit none
BEGIN_DOC
! Number of point charges to be added to the potential
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_nuclei_n_pts_charge(has)
if (has) then
write(6,'(A)') '.. >>>>> [ IO READ: n_pts_charge ] <<<<< ..'
call ezfio_get_nuclei_n_pts_charge(n_pts_charge)
else
print *, 'nuclei/n_pts_charge not found in EZFIO file'
stop 1
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( n_pts_charge, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read n_pts_charge with MPI'
endif
IRP_ENDIF
call write_time(6)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pts_charge_z, (n_pts_charge) ]
BEGIN_DOC
! Charge associated to each point charge.
END_DOC
implicit none
logical :: exists
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_nuclei_pts_charge_z(exists)
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read pts_charge_z with MPI'
endif
IRP_ENDIF
if (exists) then
if (mpi_master) then
write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_z ] <<<<< ..'
call ezfio_get_nuclei_pts_charge_z(pts_charge_z)
IRP_IF MPI
call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read pts_charge_z with MPI'
endif
IRP_ENDIF
endif
else
integer :: i
do i = 1, n_pts_charge
pts_charge_z(i) = 0.d0
enddo
endif
print*,'Point charges '
do i = 1, n_pts_charge
print*,'i,pts_charge_z(i)',i,pts_charge_z(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, pts_charge_coord, (n_pts_charge,3) ]
BEGIN_DOC
! Coordinates of each point charge.
END_DOC
implicit none
logical :: exists
PROVIDE ezfio_filename
if (mpi_master) then
call ezfio_has_nuclei_pts_charge_coord(exists)
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read pts_charge_coord with MPI'
endif
IRP_ENDIF
if (exists) then
if (mpi_master) then
double precision, allocatable :: buffer(:,:)
allocate (buffer(n_pts_charge,3))
write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_coord ] <<<<< ..'
call ezfio_get_nuclei_pts_charge_coord(buffer)
integer :: i,j
do i=1,3
do j=1,n_pts_charge
pts_charge_coord(j,i) = buffer(j,i)
enddo
enddo
deallocate(buffer)
IRP_IF MPI
call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read pts_charge_coord with MPI'
endif
IRP_ENDIF
endif
else
do i = 1, n_pts_charge
pts_charge_coord(i,:) = 0.d0
enddo
endif
print*,'Coordinates for the point charges '
do i = 1, n_pts_charge
write(*,'(I3,X,3(F16.8,X))') i,pts_charge_coord(i,1:3)
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, pt_chrg_repulsion]
implicit none
BEGIN_DOC
! repulsion between the point charges
END_DOC
integer :: i,j
double precision :: Z_A, z_B,A_center(3), B_center(3), dist
pt_chrg_repulsion = 0.d0
do i = 1, n_pts_charge
Z_A = pts_charge_z(i)
A_center(1:3) = pts_charge_coord(i,1:3)
do j = i+1, n_pts_charge
Z_B = pts_charge_z(j)
B_center(1:3) = pts_charge_coord(j,1:3)
dist = (A_center(1)-B_center(1))**2 + (A_center(2)-B_center(2))**2 + (A_center(3)-B_center(3))**2
dist = dsqrt(dist)
pt_chrg_repulsion += Z_A*Z_B/dist
enddo
enddo
print*,'Repulsion of point charges '
print*,'pt_chrg_repulsion = ',pt_chrg_repulsion
END_PROVIDER
BEGIN_PROVIDER [ double precision, pt_chrg_nuclei_repulsion]
implicit none
BEGIN_DOC
! repulsion between the point charges and the nuclei
END_DOC
integer :: i,j
double precision :: Z_A, z_B,A_center(3), B_center(3), dist
pt_chrg_nuclei_repulsion = 0.d0
do i = 1, n_pts_charge
Z_A = pts_charge_z(i)
A_center(1:3) = pts_charge_coord(i,1:3)
do j = 1, nucl_num
Z_B = nucl_charge(j)
B_center(1:3) = nucl_coord(j,1:3)
dist = (A_center(1)-B_center(1))**2 + (A_center(2)-B_center(2))**2 + (A_center(3)-B_center(3))**2
dist = dsqrt(dist)
pt_chrg_nuclei_repulsion += Z_A*Z_B/dist
enddo
enddo
print*,'Repulsion between point charges and nuclei'
print*,'pt_chrg_nuclei_repulsion = ',pt_chrg_nuclei_repulsion
END_PROVIDER

4
src/ao_one_e_ints/write_pt_charges.py → src/nuclei/write_pt_charges.py
View File

@ -13,11 +13,11 @@ def zip_in_ezfio(ezfio,tmp):
cmdzip="gzip -c "+tmp+" > "+tmpzip
os.system(cmdzip)
os.system("rm "+tmp)
cmdmv="mv "+tmpzip+" "+EZFIO+"/ao_one_e_ints/"+tmpzip
cmdmv="mv "+tmpzip+" "+EZFIO+"/nuclei/"+tmpzip
os.system(cmdmv)
def mv_in_ezfio(ezfio,tmp):
cmdmv="mv "+tmp+" "+EZFIO+"/ao_one_e_ints/"+tmp
cmdmv="mv "+tmp+" "+EZFIO+"/nuclei/"+tmp
os.system(cmdmv)

42
src/tc_bi_ortho/slater_tc_opt.irp.f
View File

@ -42,3 +42,45 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
end
! ---
subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
BEGIN_DOC
!
! <key_j | H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the detail of the matrix element WITHOUT ANY CONTRIBUTION FROM THE THREE ELECTRON TERMS
!! WARNING !!
!
! Non hermitian !!
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: htot
integer :: degree
htot = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.gt.2) return
if(degree == 0)then
call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,htot)
else if (degree == 1)then
call single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint,key_j, key_i , htot)
else if(degree == 2)then
call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
endif
if(degree==0) then
htot += nuclear_repulsion
endif
end
! ---

194
src/tc_bi_ortho/slater_tc_opt_diag.irp.f
View File

@ -277,3 +277,197 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
end
subroutine diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, det_in,htot)
implicit none
BEGIN_DOC
! Computes $\langle i|H|i \rangle$. WITHOUT ANY CONTRIBUTIONS FROM 3E TERMS
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_in(Nint,2)
double precision, intent(out) :: htot
double precision :: hmono,htwoe
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: particle(Nint,2)
integer :: i, nexc(2), ispin
integer :: occ_particle(Nint*bit_kind_size,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer(bit_kind) :: det_tmp(Nint,2)
integer :: na, nb
ASSERT (Nint > 0)
ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num)
nexc(1) = 0
nexc(2) = 0
do i=1,Nint
hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1))
hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2))
particle(i,1) = iand(hole(i,1),det_in(i,1))
particle(i,2) = iand(hole(i,2),det_in(i,2))
hole(i,1) = iand(hole(i,1),ref_bitmask(i,1))
hole(i,2) = iand(hole(i,2),ref_bitmask(i,2))
nexc(1) = nexc(1) + popcnt(hole(i,1))
nexc(2) = nexc(2) + popcnt(hole(i,2))
enddo
if (nexc(1)+nexc(2) == 0) then
hmono = ref_tc_energy_1e
htwoe = ref_tc_energy_2e
htot = ref_tc_energy_tot
return
endif
!call debug_det(det_in,Nint)
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha
ASSERT (tmp(2) == nexc(2)) ! Number of particle beta
!DIR$ FORCEINLINE
call bitstring_to_list_ab(hole, occ_hole, tmp, Nint)
ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
ASSERT (tmp(2) == nexc(2)) ! Number of holes beta
det_tmp = ref_bitmask
hmono = ref_tc_energy_1e
htwoe = ref_tc_energy_2e
do ispin=1,2
na = elec_num_tab(ispin)
nb = elec_num_tab(iand(ispin,1)+1)
do i=1,nexc(ispin)
!DIR$ FORCEINLINE
call ac_tc_operator_no_3e( occ_particle(i,ispin), ispin, det_tmp, hmono,htwoe, Nint,na,nb)
!DIR$ FORCEINLINE
call a_tc_operator_no_3e ( occ_hole (i,ispin), ispin, det_tmp, hmono,htwoe, Nint,na,nb)
enddo
enddo
htot = hmono+htwoe
end
subroutine ac_tc_operator_no_3e(iorb,ispin,key,hmono,htwoe,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Routine that computes one- and two-body energy corresponding
!
! to the ADDITION of an electron in an orbital 'iorb' of spin 'ispin'
!
! onto a determinant 'key'.
!
! in output, the determinant key is changed by the ADDITION of that electron
!
! and the quantities hmono,htwoe are INCREMENTED
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hmono,htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i,jj,mm,j,m
double precision :: direct_int, exchange_int
if (iorb < 1) then
print *, irp_here, ': iorb < 1'
print *, iorb, mo_num
stop -1
endif
if (iorb > mo_num) then
print *, irp_here, ': iorb > mo_num'
print *, iorb, mo_num
stop -1
endif
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
ASSERT (tmp(1) == elec_alpha_num)
ASSERT (tmp(2) == elec_beta_num)
k = shiftr(iorb-1,bit_kind_shift)+1
ASSERT (k >0)
l = iorb - shiftl(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
key(k,ispin) = ibset(key(k,ispin),l)
other_spin = iand(ispin,1)+1
hmono = hmono + mo_bi_ortho_tc_one_e(iorb,iorb)
! Same spin
do i=1,na
htwoe = htwoe + mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
enddo
na = na+1
end
subroutine a_tc_operator_no_3e(iorb,ispin,key,hmono,htwoe,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Routine that computes one- and two-body energy corresponding
!
! to the REMOVAL of an electron in an orbital 'iorb' of spin 'ispin'
!
! onto a determinant 'key'.
!
! in output, the determinant key is changed by the REMOVAL of that electron
!
! and the quantities hmono,htwoe are INCREMENTED
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hmono,htwoe
double precision :: direct_int, exchange_int
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i,jj,mm,j,m
integer :: tmp(2)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
k = shiftr(iorb-1,bit_kind_shift)+1
ASSERT (k>0)
l = iorb - shiftl(k-1,bit_kind_shift)-1
key(k,ispin) = ibclr(key(k,ispin),l)
other_spin = iand(ispin,1)+1
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
na = na-1
hmono = hmono - mo_bi_ortho_tc_one_e(iorb,iorb)
! Same spin
do i=1,na
htwoe= htwoe- mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
enddo
end

55
src/tc_bi_ortho/slater_tc_opt_double.irp.f
View File

@ -419,3 +419,58 @@ subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
enddo
end
subroutine double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
BEGIN_DOC
! <key_j | H_tilde | key_i> for double excitation ONLY FOR ONE- AND TWO-BODY TERMS
!!
!! WARNING !!
!
! Non hermitian !!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: htot
double precision :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
double precision :: get_mo_two_e_integral_tc_int,phase
call get_excitation_degree(key_i, key_j, degree, Nint)
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
if(degree.ne.2)then
return
endif
integer :: degree_i,degree_j
call get_excitation_degree(ref_bitmask,key_i,degree_i,N_int)
call get_excitation_degree(ref_bitmask,key_j,degree_j,N_int)
call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(s1.ne.s2)then
! opposite spin two-body
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
else
! same spin two-body
! direct terms
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
! exchange terms
htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1)
endif
htwoe *= phase
htot = htwoe
end

112
src/tc_bi_ortho/slater_tc_opt_single.irp.f
View File

@ -458,3 +458,115 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)
END_PROVIDER
subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
BEGIN_DOC
! <key_j | H_tilde | key_i> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS
!!
!! WARNING !!
!
! Non hermitian !!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: htot
double precision :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
double precision :: get_mo_two_e_integral_tc_int, phase
double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13
integer :: other_spin(2)
integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2)
other_spin(1) = 2
other_spin(2) = 1
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.ne.1)then
return
endif
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
call get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,htot)
end
subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,hmono,htwoe,htot)
use bitmasks
implicit none
integer,intent(in) :: h,p,spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2)
double precision, intent(out) :: hmono,htwoe,htot
integer(bit_kind) :: differences(N_int,2)
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: partcl(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i
double precision :: buffer_c(mo_num),buffer_x(mo_num)
do i=1, mo_num
buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
enddo
do i = 1, N_int
differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
partcl(i,1) = iand(differences(i,1),key_i(i,1))
partcl(i,2) = iand(differences(i,2),key_i(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
hmono = mo_bi_ortho_tc_one_e(p,h)
htwoe = fock_op_2_e_tc_closed_shell(p,h)
! holes :: direct terms
do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,1)
htwoe -= buffer_c(i)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
htwoe -= buffer_c(i)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
htwoe += buffer_x(i)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
htwoe += buffer_c(i)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
htwoe += buffer_c(i)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
htwoe -= buffer_x(i)
enddo
htwoe = htwoe * phase
hmono = hmono * phase
htot = htwoe + hmono
end