PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-19 19:52:20 +02:00
Code Releases Activity

does not compile but working on it

Browse Source
This commit is contained in:
eginer 2024-05-06 18:30:05 +02:00
parent 7335e29784
commit 109a956f0d
13 changed files with 769 additions and 606 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

391
plugins/local/slater_tc/h_mat_triple.irp.f
View File

@ -1,391 +0,0 @@
subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
use bitmasks
BEGIN_DOC
! returns the array, for each spin, of holes/particles between key_i and key_j
!
! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2)
double precision, intent(out) :: phase
integer :: ispin,k,i,pos
integer(bit_kind) :: key_hole, key_particle
integer(bit_kind) :: xorvec(N_int_max,2)
holes_array = -1
particles_array = -1
degree_array = 0
do i = 1, N_int
xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
degree_array(1) += popcnt(xorvec(i,1))
degree_array(2) += popcnt(xorvec(i,2))
enddo
degree_array(1) = shiftr(degree_array(1),1)
degree_array(2) = shiftr(degree_array(2),1)
do ispin = 1, 2
k = 1
!!! GETTING THE HOLES
do i = 1, N_int
key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
do ispin = 1, 2
k = 1
!!! GETTING THE PARTICLES
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general '
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
integer :: h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree_array(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine get_holes_general(key_i, key_j,Nint, holes_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of holes between key_i and key_j
!
! with the following convention: a_{hole}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2)
integer(bit_kind) :: key_hole
integer :: ispin,k,i,pos
holes_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_particles_general(key_i, key_j,Nint,particles_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of particles between key_i and key_j
!
! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: particles_array(100,2)
integer(bit_kind) :: key_particle
integer :: ispin,k,i,pos
particles_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'Those are the two determinants'
call debug_det(key_i, N_int)
call debug_det(key_j, N_int)
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
implicit none
integer, intent(in) :: degree(2), Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
integer, intent(in) :: holes_array(100,2),particles_array(100,2)
double precision, intent(out) :: phase
integer :: i,ispin,h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine H_tc_s2_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
BEGIN_DOC
! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
use bitmasks
implicit none
integer, intent(in) :: N_st,sze
double precision, intent(in) :: u_0(sze,N_st)
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze)
integer :: i,j,degree,ist
double precision :: hmono, htwoe, hthree, htot
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree .ne. 3)cycle
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot)
do ist = 1, N_st
v_0(i,ist) += htot * u_0(j,ist)
enddo
enddo
enddo
end
subroutine H_tc_s2_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze)
BEGIN_DOC
! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
use bitmasks
implicit none
integer, intent(in) :: N_st,sze
double precision, intent(in) :: u_0(sze,N_st)
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze)
integer :: i,j,degree,ist
double precision :: hmono, htwoe, hthree, htot
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
!$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0) &
!$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot)
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree .ne. 3)cycle
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot)
do ist = 1, N_st
v_0(i,ist) += htot * u_0(j,ist)
enddo
enddo
enddo
!$OMP END PARALLEL DO
end
! ---
subroutine H_tc_s2_dagger_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
BEGIN_DOC
! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
use bitmasks
implicit none
integer, intent(in) :: N_st,sze
double precision, intent(in) :: u_0(sze,N_st)
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze)
integer :: i,j,degree,ist
double precision :: hmono, htwoe, hthree, htot
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree .ne. 3)cycle
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot)
do ist = 1, N_st
v_0(i,ist) += htot * u_0(j,ist)
enddo
enddo
enddo
end
subroutine H_tc_s2_dagger_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze)
BEGIN_DOC
! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
use bitmasks
implicit none
integer, intent(in) :: N_st,sze
double precision, intent(in) :: u_0(sze,N_st)
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze)
integer :: i,j,degree,ist
double precision :: hmono, htwoe, hthree, htot
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
!$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0) &
!$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot)
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree .ne. 3)cycle
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot)
do ist = 1, N_st
v_0(i,ist) += htot * u_0(j,ist)
enddo
enddo
enddo
!$OMP END PARALLEL DO
end
! ---
subroutine triple_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
use bitmasks
BEGIN_DOC
! <key_j | H_tilde | key_i> for triple excitation
!!
!! WARNING !!
!
! Genuine triple excitations of the same spin are not yet implemented
END_DOC
implicit none
integer(bit_kind), intent(in) :: key_j(N_int,2),key_i(N_int,2)
integer, intent(in) :: Nint
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: degree
integer :: h1, p1, h2, p2, s1, s2, h3, p3, s3
integer :: holes_array(100,2),particles_array(100,2),degree_array(2)
double precision :: phase,sym_3_e_int_from_6_idx_tensor
hmono = 0.d0
htwoe = 0.d0
hthree = 0.d0
htot = 0.d0
call get_excitation_general(key_j, key_i, Nint,degree_array,holes_array, particles_array,phase)
degree = degree_array(1) + degree_array(2)
if(degree .ne. 3)return
if(degree_array(1)==3.or.degree_array(2)==3)then
if(degree_array(1) == 3)then
h1 = holes_array(1,1)
h2 = holes_array(2,1)
h3 = holes_array(3,1)
p1 = particles_array(1,1)
p2 = particles_array(2,1)
p3 = particles_array(3,1)
else
h1 = holes_array(1,2)
h2 = holes_array(2,2)
h3 = holes_array(3,2)
p1 = particles_array(1,2)
p2 = particles_array(2,2)
p3 = particles_array(3,2)
endif
hthree = sym_3_e_int_from_6_idx_tensor(p3, p2, p1, h3, h2, h1)
else
if(degree_array(1) == 2.and.degree_array(2) == 1)then ! double alpha + single beta
h1 = holes_array(1,1)
h2 = holes_array(2,1)
h3 = holes_array(1,2)
p1 = particles_array(1,1)
p2 = particles_array(2,1)
p3 = particles_array(1,2)
else if(degree_array(2) == 2 .and. degree_array(1) == 1)then ! double beta + single alpha
h1 = holes_array(1,2)
h2 = holes_array(2,2)
h3 = holes_array(1,1)
p1 = particles_array(1,2)
p2 = particles_array(2,2)
p3 = particles_array(1,1)
else
print*,'PB !!'
stop
endif
hthree = three_body_ints_bi_ort(p3,p2,p1,h3,h2,h1) - three_body_ints_bi_ort(p3,p2,p1,h3,h1,h2)
endif
hthree *= phase
htot = hthree
end

311
plugins/local/slater_tc/slater_tc_opt_diag.irp.f
View File

@ -19,13 +19,13 @@
PROVIDE HF_bitmask
PROVIDE mo_l_coef mo_r_coef
call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot)
call diag_htc_bi_orth_2e_brute(N_int, HF_bitmask, hmono, htwoe, htot)
ref_tc_energy_1e = hmono
ref_tc_energy_2e = htwoe
if(three_body_h_tc) then
call diag_htilde_three_body_ints_bi_ort_slow(N_int, HF_bitmask, hthree)
call diag_htc_bi_orth_3e_brute(N_int, HF_bitmask, hthree)
ref_tc_energy_3e = hthree
else
ref_tc_energy_3e = 0.d0
@ -524,3 +524,310 @@ end
! ---
subroutine diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot)
BEGIN_DOC
!
! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
double precision :: get_mo_two_e_integral_tc_int
integer(bit_kind) :: key_i_core(Nint,2)
PROVIDE mo_bi_ortho_tc_two_e
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
do ispin = 1, 2
do i = 1, Ne(ispin)
ii = occ(i,ispin)
hmono += mo_bi_ortho_tc_one_e(ii,ii)
enddo
enddo
! alpha/beta two-body
ispin = 1
jspin = 2
do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
ii = occ(i,ispin)
do j = 1, Ne(jspin) ! electron 2
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
enddo
enddo
! alpha/alpha two-body
do i = 1, Ne(ispin)
ii = occ(i,ispin)
do j = i+1, Ne(ispin)
jj = occ(j,ispin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
! beta/beta two-body
do i = 1, Ne(jspin)
ii = occ(i,jspin)
do j = i+1, Ne(jspin)
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
htot = hmono + htwoe
end
! ---
subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
BEGIN_DOC
! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hthree
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
integer(bit_kind) :: key_i_core(Nint,2)
double precision :: direct_int, exchange_int, ref
double precision, external :: sym_3_e_int_from_6_idx_tensor
double precision, external :: three_e_diag_parrallel_spin
PROVIDE mo_l_coef mo_r_coef
if(core_tc_op) then
do i = 1, Nint
key_i_core(i,1) = xor(key_i(i,1), core_bitmask(i,1))
key_i_core(i,2) = xor(key_i(i,2), core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
endif
hthree = 0.d0
if((Ne(1)+Ne(2)) .ge. 3) then
! alpha/alpha/beta three-body
do i = 1, Ne(1)
ii = occ(i,1)
do j = i+1, Ne(1)
jj = occ(j,1)
do m = 1, Ne(2)
mm = occ(m,2)
!direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) !uses the 6-idx tensor
!exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) !uses the 6-idx tensor
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) !uses 3-idx tensor
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) !uses 3-idx tensor
hthree += direct_int - exchange_int
enddo
enddo
enddo
! beta/beta/alpha three-body
do i = 1, Ne(2)
ii = occ(i,2)
do j = i+1, Ne(2)
jj = occ(j,2)
do m = 1, Ne(1)
mm = occ(m,1)
!direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) !uses the 6-idx tensor
!exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) !uses the 6-idx tensor
direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii)
exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii)
hthree += direct_int - exchange_int
enddo
enddo
enddo
! alpha/alpha/alpha three-body
do i = 1, Ne(1)
ii = occ(i,1) ! 1
do j = i+1, Ne(1)
jj = occ(j,1) ! 2
do m = j+1, Ne(1)
mm = occ(m,1) ! 3
!hthree += sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) !uses the 6 idx tensor
hthree += three_e_diag_parrallel_spin(mm,jj,ii) !uses only 3-idx tensors
enddo
enddo
enddo
! beta/beta/beta three-body
do i = 1, Ne(2)
ii = occ(i,2) ! 1
do j = i+1, Ne(2)
jj = occ(j,2) ! 2
do m = j+1, Ne(2)
mm = occ(m,2) ! 3
!hthree += sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) !uses the 6 idx tensor
hthree += three_e_diag_parrallel_spin(mm,jj,ii) !uses only 3-idx tensors
enddo
enddo
enddo
endif
end
BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS
!
! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form <mji|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin
three_e_diag_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_diag_parrallel_spin_prov ...'
integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of <mjk|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0, three_e_single_parrallel_spin
three_e_single_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_single_parrallel_spin_prov ...'
integral = three_e_single_parrallel_spin(1,1,1,1)
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_double_parrallel_spin_prov(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0, three_e_double_parrallel_spin
three_e_double_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_double_parrallel_spin_prov ...'
call wall_time(wall0)
integral = three_e_double_parrallel_spin(1,1,1,1,1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0
END_PROVIDER

140
plugins/local/slater_tc/symmetrized_3_e_int_prov.irp.f
View File

@ -1,140 +0,0 @@
BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS
!
! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form <mji|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin
three_e_diag_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_diag_parrallel_spin_prov ...'
integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of <mjk|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0, three_e_single_parrallel_spin
three_e_single_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_single_parrallel_spin_prov ...'
integral = three_e_single_parrallel_spin(1,1,1,1)
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_double_parrallel_spin_prov(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0, three_e_double_parrallel_spin
three_e_double_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_double_parrallel_spin_prov ...'
call wall_time(wall0)
integral = three_e_double_parrallel_spin(1,1,1,1,1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0
END_PROVIDER

59
plugins/local/slater_tc_no_opt/.gitignore vendored Normal file
View File

@ -0,0 +1,59 @@
IRPF90_temp/
IRPF90_man/
build.ninja
irpf90.make
ezfio_interface.irp.f
irpf90_entities
tags
Makefile
ao_basis
ao_one_e_ints
ao_two_e_erf_ints
ao_two_e_ints
aux_quantities
becke_numerical_grid
bitmask
cis
cisd
cipsi
davidson
davidson_dressed
davidson_undressed
density_for_dft
determinants
dft_keywords
dft_utils_in_r
dft_utils_one_e
dft_utils_two_body
dressing
dummy
electrons
ezfio_files
fci
generators_cas
generators_full
hartree_fock
iterations
kohn_sham
kohn_sham_rs
mo_basis
mo_guess
mo_one_e_ints
mo_two_e_erf_ints
mo_two_e_ints
mpi
mrpt_utils
nuclei
perturbation
pseudo
psiref_cas
psiref_utils
scf_utils
selectors_cassd
selectors_full
selectors_utils
single_ref_method
slave
tools
utils
zmq

8
plugins/local/slater_tc_no_opt/NEED Normal file
View File

@ -0,0 +1,8 @@
determinants
normal_order_old
bi_ort_ints
bi_ortho_mos
tc_keywords
non_hermit_dav
dav_general_mat
tc_scf

4
plugins/local/slater_tc_no_opt/README.rst Normal file
View File

@ -0,0 +1,4 @@
================
slater_tc_no_opt
================

0
plugins/local/slater_tc/h_biortho.irp.f → plugins/local/slater_tc_no_opt/h_biortho.irp.f
View File

193
plugins/local/slater_tc_no_opt/h_mat_triple.irp.f Normal file
View File

@ -0,0 +1,193 @@
subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
use bitmasks
BEGIN_DOC
! returns the array, for each spin, of holes/particles between key_i and key_j
!
! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2)
double precision, intent(out) :: phase
integer :: ispin,k,i,pos
integer(bit_kind) :: key_hole, key_particle
integer(bit_kind) :: xorvec(N_int_max,2)
holes_array = -1
particles_array = -1
degree_array = 0
do i = 1, N_int
xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
degree_array(1) += popcnt(xorvec(i,1))
degree_array(2) += popcnt(xorvec(i,2))
enddo
degree_array(1) = shiftr(degree_array(1),1)
degree_array(2) = shiftr(degree_array(2),1)
do ispin = 1, 2
k = 1
!!! GETTING THE HOLES
do i = 1, N_int
key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
do ispin = 1, 2
k = 1
!!! GETTING THE PARTICLES
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general '
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
integer :: h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree_array(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine get_holes_general(key_i, key_j,Nint, holes_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of holes between key_i and key_j
!
! with the following convention: a_{hole}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2)
integer(bit_kind) :: key_hole
integer :: ispin,k,i,pos
holes_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_particles_general(key_i, key_j,Nint,particles_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of particles between key_i and key_j
!
! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: particles_array(100,2)
integer(bit_kind) :: key_particle
integer :: ispin,k,i,pos
particles_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'Those are the two determinants'
call debug_det(key_i, N_int)
call debug_det(key_j, N_int)
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
implicit none
integer, intent(in) :: degree(2), Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
integer, intent(in) :: holes_array(100,2),particles_array(100,2)
double precision, intent(out) :: phase
integer :: i,ispin,h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end

0
plugins/local/slater_tc/h_tc_bi_ortho_psi.irp.f → plugins/local/slater_tc_no_opt/h_tc_bi_ortho_psi.irp.f
View File

2
plugins/local/slater_tc/slater_tc_3e_slow.irp.f → plugins/local/slater_tc_no_opt/slater_tc_3e_slow.irp.f
View File

@ -1,7 +1,7 @@
! ---
subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
subroutine diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
BEGIN_DOC
! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS

2
plugins/local/slater_tc/slater_tc.irp.f → plugins/local/slater_tc_no_opt/slater_tc_no_opt.irp.f
View File

@ -1,4 +1,4 @@
program slater_tc
program slater_tc_no_opt
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here

73
plugins/local/slater_tc/slater_tc_slow.irp.f → plugins/local/slater_tc_no_opt/slater_tc_slow.irp.f
View File

@ -61,7 +61,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
if(degree.gt.2) return
if(degree == 0) then
call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
call diag_htc_bi_orth_2e_brute(Nint, key_i, hmono, htwoe, htot)
else if (degree == 1) then
call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
else if(degree == 2) then
@ -76,7 +76,7 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then
call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then
call diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
call diag_htc_bi_orth_3e_brute(Nint, key_i, hthree)
endif
endif
@ -95,75 +95,6 @@ end
! ---
subroutine diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
BEGIN_DOC
!
! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
double precision :: get_mo_two_e_integral_tc_int
integer(bit_kind) :: key_i_core(Nint,2)
PROVIDE mo_bi_ortho_tc_two_e
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
do ispin = 1, 2
do i = 1, Ne(ispin)
ii = occ(i,ispin)
hmono += mo_bi_ortho_tc_one_e(ii,ii)
enddo
enddo
! alpha/beta two-body
ispin = 1
jspin = 2
do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
ii = occ(i,ispin)
do j = 1, Ne(jspin) ! electron 2
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
enddo
enddo
! alpha/alpha two-body
do i = 1, Ne(ispin)
ii = occ(i,ispin)
do j = i+1, Ne(ispin)
jj = occ(j,ispin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
! beta/beta two-body
do i = 1, Ne(jspin)
ii = occ(i,jspin)
do j = i+1, Ne(jspin)
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
htot = hmono + htwoe
end
! ---
subroutine double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
BEGIN_DOC

192
src/determinants/slater_rules_general.irp.f Normal file
View File

@ -0,0 +1,192 @@
subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
use bitmasks
BEGIN_DOC
! returns the array, for each spin, of holes/particles between key_i and key_j
!
! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2)
double precision, intent(out) :: phase
integer :: ispin,k,i,pos
integer(bit_kind) :: key_hole, key_particle
integer(bit_kind) :: xorvec(N_int_max,2)
holes_array = -1
particles_array = -1
degree_array = 0
do i = 1, N_int
xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
degree_array(1) += popcnt(xorvec(i,1))
degree_array(2) += popcnt(xorvec(i,2))
enddo
degree_array(1) = shiftr(degree_array(1),1)
degree_array(2) = shiftr(degree_array(2),1)
do ispin = 1, 2
k = 1
!!! GETTING THE HOLES
do i = 1, N_int
key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
do ispin = 1, 2
k = 1
!!! GETTING THE PARTICLES
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_excitation_general '
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
integer :: h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree_array(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end
subroutine get_holes_general(key_i, key_j,Nint, holes_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of holes between key_i and key_j
!
! with the following convention: a_{hole}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: holes_array(100,2)
integer(bit_kind) :: key_hole
integer :: ispin,k,i,pos
holes_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
do while(key_hole .ne.0_bit_kind)
pos = trailz(key_hole)
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_hole = ibclr(key_hole,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_particles_general(key_i, key_j,Nint,particles_array)
use bitmasks
BEGIN_DOC
! returns the array, per spin, of particles between key_i and key_j
!
! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
END_DOC
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
integer, intent(out) :: particles_array(100,2)
integer(bit_kind) :: key_particle
integer :: ispin,k,i,pos
particles_array = -1
do ispin = 1, 2
k = 1
do i = 1, N_int
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
do while(key_particle .ne.0_bit_kind)
pos = trailz(key_particle)
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
key_particle = ibclr(key_particle,pos)
k += 1
if(k .gt.100)then
print*,'WARNING in get_holes_general'
print*,'More than a 100-th excitation for spin ',ispin
print*,'Those are the two determinants'
call debug_det(key_i, N_int)
call debug_det(key_j, N_int)
print*,'stoping ...'
stop
endif
enddo
enddo
enddo
end
subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
implicit none
integer, intent(in) :: degree(2), Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
integer, intent(in) :: holes_array(100,2),particles_array(100,2)
double precision, intent(out) :: phase
integer :: i,ispin,h,p, i_ok
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
integer :: exc(0:2,2,2)
double precision :: phase_tmp
allocate(det_i(Nint,2),det_ip(N_int,2))
det_i = key_i
phase = 1.d0
do ispin = 1, 2
do i = 1, degree(ispin)
h = holes_array(i,ispin)
p = particles_array(i,ispin)
det_ip = det_i
call do_single_excitation(det_ip,h,p,ispin,i_ok)
if(i_ok == -1)then
print*,'excitation was not possible '
stop
endif
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
phase *= phase_tmp
det_i = det_ip
enddo
enddo
end