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

it works with same left and right coefficients for 2 dets in WF

Browse Source
This commit is contained in:
eginer 2023-01-24 13:49:04 +01:00
parent 1e1db71df6
commit af60a02919
11 changed files with 890 additions and 271 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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)

97
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,57 @@ 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'
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
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'

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

@ -89,6 +89,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

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