subroutine fill_cc_amplitudes(n_excits)
! design decisions: since the singles are not so many in general
! and because i could need them to correct the doubles amplitudes
! store all of the possible ones! and encode them specifically through
! (i) and (a)
integer, intent(in) :: n_excits(cc_order)
character(*), parameter :: this_routine = "fill_cc_amplitudes"
integer :: idet, ic, ex(2, cc_order), j, i
integer :: ia, ib, ja, jb, ind
integer :: a, b, elec_i, elec_j, orb_a, orb_b
logical :: t_par, t_store_full_doubles = .true.
real(dp) :: amp
real(dp), allocatable :: temp_amps(:)
integer, allocatable :: temp_ops(:, :, :)
integer :: ab(2), ac(2), bc(2), ij(2), ik(2), jk(2)
integer :: c, k, ij_ac, ij_bc, ij_ab, ik_bc, ik_ac, ik_ab
integer :: jk_bc, jk_ac, jk_ab, elec_k, orb_c, jc, ka, kb, kc, n
integer(n_int) :: temp_ilut(0:nifd)
integer :: dummy_ind, dummy_hash, temp_nI(nel)
logical :: tSuccess
integer :: i_a(2, 1), j_b(2, 1)
integer :: n_disc_1, n_disc_2
integer :: quad_ind(8), hash_ind
logical :: t_found
integer(n_int) :: temp_int(0:nifd)
real(dp) :: sign_tmp(lenof_sign)
! for this it is helpful to have an upper limit of the number of
! possible amplitudes, but just do it for the singles for now..
allocate(cc_ops(cc_order))
if (t_store_full_doubles .and. cc_order > 2 .and. t_store_hash_quadrupels) then
! also store the intermediate T1^2 and T2^2 in operator form to
! make my life easier
allocate(disc_cc_ops(2))
disc_cc_ops(1)%order = 2
disc_cc_ops(2)%order = 4
end if
! and do a nice initialization depending on the order
do i = 1, cc_order
cc_ops(i)%order = i
end do
allocate(cc_ops(1)%amplitudes(n_excits(1)))
allocate(cc_ops(1)%operators(n_excits(1), 2, 1))
allocate(cc_ops(1)%set_flag(n_excits(1)))
cc_ops(1)%operators = 0
cc_ops(1)%amplitudes = 0.0_dp
cc_ops(1)%set_flag = 0
cc_ops(1)%n_ops = n_excits(1)
! first figure out the number of double and fill the singles
! amplitudes!
n_doubles = 0
do idet = 1, int(totwalkers)
! for now also count the exact number of triples and quadrupels
ic = FindBitExcitLevel(ilutRef(:, 1), CurrentDets(:, idet))
call extract_sign(CurrentDets(:, idet), sign_tmp)
! for some strange reason there are zero elements stored..
if (abs(sign_tmp(1)) > EPS) then
select case (ic)
case (1)
n_singles = n_singles + 1
call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)
ind = cc_ops(1)%get_ind(ex(1, 1), ex(2, 1))
! for now only do it for single runs.. do i need the normalising?
amp = sign_tmp(1)
cc_ops(1)%amplitudes(ind) = amp
cc_ops(1)%operators(ind, 1, 1) = ex(1, 1)
cc_ops(1)%operators(ind, 2, 1) = ex(2, 1)
cc_ops(1)%set_flag(ind) = 1
cc_amp_norm(0, 1) = cc_amp_norm(0, 1) + 1
cc_amp_norm(1, 1) = cc_amp_norm(1, 1) + abs(amp)
cc_amp_norm(2, 1) = cc_amp_norm(2, 1) + amp**2
case (2)
n_doubles = n_doubles + 1
case (3)
n_triples = n_triples + 1
case (4)
n_quads = n_quads + 1
end select
end if
end do
cc_amp_norm(2, 1) = sqrt(cc_amp_norm(2, 1))
print *, "direct sampled doubles: ", n_doubles, "on proc: ", iProcIndex
if (allocated(cc_ops(2)%operators)) deallocate(cc_ops(2)%operators)
if (allocated(cc_ops(2)%amplitudes)) deallocate(cc_ops(2)%amplitudes)
! here i have a choice to only store the non-zero contributions
! or encode them in the full-list to better access them later on.
! especially when we want to calculate the coupled-cluster
! triples and quadrupels..
if (t_store_full_doubles) then
allocate(cc_ops(2)%operators(n_excits(2), 2, 2))
allocate(cc_ops(2)%amplitudes(n_excits(2)))
allocate(cc_ops(2)%set_flag(n_excits(2)))
cc_ops(2)%n_ops = n_excits(2)
! then also store the intermediate
allocate(disc_cc_ops(1)%operators(n_excits(2), 2, 2))
allocate(disc_cc_ops(1)%amplitudes(n_excits(2)))
allocate(disc_cc_ops(1)%set_flag(n_excits(2)))
disc_cc_ops(1)%n_ops = n_excits(2)
disc_cc_ops(1)%operators = 0
disc_cc_ops(1)%amplitudes = 0
disc_cc_ops(1)%set_flag = 0
else
allocate(cc_ops(2)%operators(n_doubles, 2, 2))
allocate(cc_ops(2)%amplitudes(n_doubles))
allocate(cc_ops(2)%set_flag(n_doubles))
cc_ops(2)%n_ops = n_doubles
cc_ops(2)%operators = 0
cc_ops(2)%amplitudes = 0.0_dp
cc_ops(2)%set_flag = 0
end if
if (t_store_disc_ops) then
n_disc_1 = 0
n_disc_2 = 0
do i = 1, cc_ops(1)%n_ops
if (cc_ops(1)%set_flag(i) == 1) then
do j = i + 1, cc_ops(1)%n_ops
if (cc_ops(1)%set_flag(j) == 1) then
i_a = cc_ops(1)%get_ex(i)
j_b = cc_ops(1)%get_ex(j)
! ensure ordering, they should be ordered by
! electrons but electron index can still be the
! same ..
if (i_a(1, 1) > j_b(1, 1)) then
call stop_all(this_routine, &
"i did smth wrong with the ordering")
end if
if (i_a(1, 1) /= j_b(1, 1)) then
amp = cc_ops(1)%get_amp(i) * cc_ops(1)%get_amp(j)
if (i_a(2, 1) < j_b(2, 1)) then
ind = cc_ops(2)%get_ind([i_a(1, 1), j_b(1, 1)], &
[i_a(2, 1), j_b(2, 1)])
amp = -amp
else if (i_a(2, 1) > j_b(2, 1)) then
ind = cc_ops(2)%get_ind([i_a(1, 1), j_b(1, 1)], &
[j_b(2, 1), i_a(2, 1)])
else
! orbitals are the same..
end if
! did i already find that:
if (disc_cc_ops(1)%set_flag(ind) == 1) then
disc_cc_ops(1)%amplitudes(ind) = &
disc_cc_ops(1)%amplitudes(ind) + amp
! if it cancelled: remove the entry:
if (abs(disc_cc_ops(1)%amplitudes(ind)) < EPS) then
disc_cc_ops(1)%set_flag(ind) = 0
disc_cc_ops(1)%operators(ind, :, :) = 0
disc_cc_ops(1)%amplitudes(ind) = 0.0_dp
! and lower ofc..
n_disc_1 = n_disc_1 - 1
end if
else
! otherwise set it..
disc_cc_ops(1)%amplitudes(ind) = amp
disc_cc_ops(1)%set_flag(ind) = 1
disc_cc_ops(1)%operators(ind, 1, :) = [i_a(1, 1), j_b(1, 1)]
disc_cc_ops(1)%operators(ind, 2, :) = &
[min(i_a(2, 1), j_b(2, 1)), max(i_a(2, 1), j_b(2, 1))]
n_disc_1 = n_disc_1 + 1
end if
end if
end if
end do
end if
end do
end if
print *, "number of disconnected T1^2: ", n_disc_1, "on proc: ", iProcIndex
cc_ops(2)%operators = 0
cc_ops(2)%amplitudes = 0.0_dp
cc_ops(2)%set_flag = 0
do idet = 1, int(totwalkers)
ic = FindBitExcitLevel(ilutRef(:, 1), CurrentDets(:, idet))
if (ic == 2) then
call extract_sign(CurrentDets(:, idet), sign_tmp)
! i hope this routine is not buggy..
call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)
ia = cc_ops(1)%get_ind(ex(1, 1), ex(2, 1))
ib = cc_ops(1)%get_ind(ex(1, 1), ex(2, 2))
ja = cc_ops(1)%get_ind(ex(1, 2), ex(2, 1))
jb = cc_ops(1)%get_ind(ex(1, 2), ex(2, 2))
! check first if the amp gets 0 due to the singles
amp = sign_tmp(1) + &
cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(ib) - &
cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb)
if (abs(amp) > EPS) then
if (t_store_full_doubles) then
ind = cc_ops(2)%get_ind(ex(1, :), ex(2, :))
else
ind = j
j = j + 1
end if
cc_ops(2)%operators(ind, 1, :) = ex(1, 1:2)
cc_ops(2)%operators(ind, 2, :) = ex(2, 1:2)
cc_ops(2)%amplitudes(ind) = amp
cc_ops(2)%set_flag(ind) = 1
! i think i can already sum up the norms here..
! since we have all the necessary contribs
cc_amp_norm(0, 2) = cc_amp_norm(0, 2) + 1
cc_amp_norm(1, 2) = cc_amp_norm(1, 2) + abs(amp)
cc_amp_norm(2, 2) = cc_amp_norm(2, 2) + amp**2
end if
end if
end do
! i just realise that i have to run over all the possible double
! excitations, since t_ij^ab = c_ij^ab - t_i^a... could be
! non-zero, just from the single contribution..
if (t_store_full_doubles) then
do i = 1, nel
elec_i = projedet(i, 1)
do j = i + 1, nel
elec_j = projedet(j, 1)
t_par = same_spin(elec_i, elec_j)
if (t_par) then
do a = 1, nvirt
orb_a = mask_virt_ni(a, 1)
do b = a + 1, nvirt
orb_b = mask_virt_ni(b, 1)
if (same_spin(orb_a, orb_b) .and. &
same_spin(orb_a, elec_i)) then
ind = cc_ops(2)%get_ind([elec_i, elec_j], &
[orb_a, orb_b])
ASSERT(elec_i < elec_j)
ASSERT(orb_a < orb_b)
! do the disc op storage here!
ia = cc_ops(1)%get_ind([elec_i], [orb_a])
ib = cc_ops(1)%get_ind([elec_i], [orb_b])
ja = cc_ops(1)%get_ind([elec_j], [orb_a])
jb = cc_ops(1)%get_ind([elec_j], [orb_b])
amp = cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(ib) - &
cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb)
if (abs(amp) > EPS) then
if (.not. cc_ops(2)%set_flag(ind) == 1) then
print *, "para: ", elec_i, elec_j, orb_a, orb_b
cc_ops(2)%operators(ind, 1, :) = [elec_i, elec_j]
cc_ops(2)%operators(ind, 2, :) = [orb_a, orb_b]
cc_ops(2)%amplitudes(ind) = amp
cc_ops(2)%set_flag(ind) = 1
n_doubles = n_doubles + 1
cc_amp_norm(0, 2) = cc_amp_norm(0, 2) + 1
cc_amp_norm(1, 2) = cc_amp_norm(1, 2) + abs(amp)
cc_amp_norm(2, 2) = cc_amp_norm(2, 2) + amp**2
end if
if (t_store_disc_ops) then
disc_cc_ops(1)%operators(ind, 1, :) = [elec_i, elec_j]
disc_cc_ops(1)%operators(ind, 2, :) = [orb_a, orb_b]
disc_cc_ops(1)%amplitudes(ind) = amp
disc_cc_ops(1)%set_flag(ind) = 1
n_disc_2 = n_disc_2 + 1
end if
end if
end if
end do
end do
else
do a = 1, nvirt
orb_a = mask_virt_ni(a, 1)
do b = a + 1, nvirt
orb_b = mask_virt_ni(b, 1)
if (.not. same_spin(orb_a, orb_b)) then
ind = cc_ops(2)%get_ind([elec_i, elec_j], &
[orb_a, orb_b])
ASSERT(elec_i < elec_j)
ASSERT(orb_a < orb_b)
ia = cc_ops(1)%get_ind([elec_i], [orb_a])
ib = cc_ops(1)%get_ind([elec_i], [orb_b])
ja = cc_ops(1)%get_ind([elec_j], [orb_a])
jb = cc_ops(1)%get_ind([elec_j], [orb_b])
amp = cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(ib) - &
cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb)
if (abs(amp) > EPS) then
if (.not. cc_ops(2)%set_flag(ind) == 1) then
print *, "anti: ", elec_i, elec_j, orb_a, orb_b
cc_ops(2)%operators(ind, 1, :) = [elec_i, elec_j]
cc_ops(2)%operators(ind, 2, :) = [orb_a, orb_b]
cc_ops(2)%amplitudes(ind) = amp
cc_ops(2)%set_flag(ind) = 1
n_doubles = n_doubles + 1
cc_amp_norm(0, 2) = cc_amp_norm(0, 2) + 1
cc_amp_norm(1, 2) = cc_amp_norm(1, 2) + abs(amp)
cc_amp_norm(2, 2) = cc_amp_norm(2, 2) + amp**2
end if
if (t_store_disc_ops) then
disc_cc_ops(1)%operators(ind, 1, :) = [elec_i, elec_j]
disc_cc_ops(1)%operators(ind, 2, :) = [orb_a, orb_b]
disc_cc_ops(1)%amplitudes(ind) = amp
disc_cc_ops(1)%set_flag(ind) = 1
n_disc_2 = n_disc_2 + 1
end if
end if
end if
end do
end do
end if
end do
end do
! end if
end if
print *, "disconnecte T1^2 calc. differently: ", n_disc_2, "on proc: ", &
iProcIndex
cc_amp_norm(2, 2) = sqrt(cc_amp_norm(2, 2))
print *, "doubles after singles contribution: ", n_doubles, "on proc:", &
iProcIndex
! and should i also do the triples.. just to check maybe?
! but here i definetly only want to store the non-zero
! contributions.. but.. i actually would need to run over
! all the possible ones also since the contribution could be
! approximated by T1*T2 or T1^3..
! i could loop.. but i do not want to encode all of them!
! but essentially i have to first allocate a vector of all possible
! ones to be sure..
if (cc_order >= 3) then
allocate(temp_amps(n_excits(3)))
allocate(temp_ops(n_excits(3), 2, 3))
n = 1
! first analyze the wavefunction:
do idet = 1, int(totwalkers)
ic = FindBitExcitLevel(iLutRef(:, 1), CurrentDets(:, idet))
if (ic == 3) then
call extract_sign(CurrentDets(:, idet), sign_tmp)
call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)
ia = cc_ops(1)%get_ind(ex(1, 1), ex(2, 1))
ib = cc_ops(1)%get_ind(ex(1, 1), ex(2, 2))
ic = cc_ops(1)%get_ind(ex(1, 1), ex(2, 3))
ja = cc_ops(1)%get_ind(ex(1, 2), ex(2, 1))
jb = cc_ops(1)%get_ind(ex(1, 2), ex(2, 2))
jc = cc_ops(1)%get_ind(ex(1, 2), ex(2, 3))
ka = cc_ops(1)%get_ind(ex(1, 3), ex(2, 1))
kb = cc_ops(1)%get_ind(ex(1, 3), ex(2, 2))
kc = cc_ops(1)%get_ind(ex(1, 3), ex(2, 3))
ij = [ex(1, 1), ex(1, 2)]
ik = [ex(1, 1), ex(1, 3)]
jk = [ex(1, 2), ex(1, 3)]
ab = [ex(2, 1), ex(2, 2)]
ac = [ex(2, 1), ex(2, 3)]
bc = [ex(2, 2), ex(2, 3)]
jk_bc = cc_ops(2)%get_ind(jk, bc)
jk_ac = cc_ops(2)%get_ind(jk, ac)
jk_ab = cc_ops(2)%get_ind(jk, ab)
ik_bc = cc_ops(2)%get_ind(ik, bc)
ik_ac = cc_ops(2)%get_ind(ik, ac)
ik_ab = cc_ops(2)%get_ind(ik, ab)
ij_bc = cc_ops(2)%get_ind(ij, bc)
ij_ac = cc_ops(2)%get_ind(ij, ac)
ij_ab = cc_ops(2)%get_ind(ij, ab)
amp = sign_tmp(1) &
- cc_ops(1)%get_amp(ia) * cc_ops(2)%get_amp(jk_bc) &
+ cc_ops(1)%get_amp(ib) * cc_ops(2)%get_amp(jk_ac) &
- cc_ops(1)%get_amp(ic) * cc_ops(2)%get_amp(jk_ab) &
+ cc_ops(1)%get_amp(ja) * cc_ops(2)%get_amp(ik_bc) &
- cc_ops(1)%get_amp(jb) * cc_ops(2)%get_amp(ik_ac) &
+ cc_ops(1)%get_amp(jc) * cc_ops(2)%get_amp(ik_ab) &
- cc_ops(1)%get_amp(ka) * cc_ops(2)%get_amp(ij_bc) &
+ cc_ops(1)%get_amp(kb) * cc_ops(2)%get_amp(ij_ac) &
- cc_ops(1)%get_amp(kc) * cc_ops(2)%get_amp(ij_ab) &
- cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(kc) &
+ cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(kb) &
+ cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kc) &
- cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(ka) &
- cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kb) &
+ cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(ka)
if (abs(amp) > EPS) then
temp_amps(n) = amp
temp_ops(n, :, :) = ex
n = n + 1
! also here i can add the norms up already..
cc_amp_norm(0, 3) = cc_amp_norm(0, 3) + 1
cc_amp_norm(1, 3) = cc_amp_norm(1, 3) + abs(amp)
cc_amp_norm(2, 3) = cc_amp_norm(2, 3) + amp**2
end if
end if
end do
print *, "direct sampled triples: ", n - 1, "on Proc: ", iProcIndex
do i = 1, nel
elec_i = projedet(i, 1)
do j = i + 1, nel
elec_j = projedet(j, 1)
ij = [elec_i, elec_j]
do k = j + 1, nel
elec_k = projedet(k, 1)
ik = [elec_i, elec_k]
jk = [elec_j, elec_k]
do a = 1, nvirt
orb_a = mask_virt_ni(a, 1)
ia = cc_ops(1)%get_ind([elec_i], [orb_a])
ja = cc_ops(1)%get_ind([elec_j], [orb_a])
ka = cc_ops(1)%get_ind([elec_k], [orb_a])
do b = a + 1, nvirt
orb_b = mask_virt_ni(b, 1)
ab = [orb_a, orb_b]
ib = cc_ops(1)%get_ind([elec_i], [orb_b])
jb = cc_ops(1)%get_ind([elec_j], [orb_b])
kb = cc_ops(1)%get_ind([elec_k], [orb_b])
ij_ab = cc_ops(2)%get_ind(ij, ab)
ik_ab = cc_ops(2)%get_ind(ik, ab)
jk_ab = cc_ops(2)%get_ind(jk, ab)
do c = b + 1, nvirt
orb_c = mask_virt_ni(c, 1)
ac = [orb_a, orb_c]
bc = [orb_b, orb_c]
ic = cc_ops(1)%get_ind([elec_i], [orb_c])
jc = cc_ops(1)%get_ind([elec_j], [orb_c])
kc = cc_ops(1)%get_ind([elec_k], [orb_c])
ij_ac = cc_ops(2)%get_ind(ij, ac)
ik_ac = cc_ops(2)%get_ind(ik, ac)
jk_ac = cc_ops(2)%get_ind(jk, ac)
ij_bc = cc_ops(2)%get_ind(ij, bc)
ik_bc = cc_ops(2)%get_ind(ik, bc)
jk_bc = cc_ops(2)%get_ind(jk, bc)
! and now i have to check if the
! triple excitation (ijk|abc) is in the
! wallker list! use the hash-table!
! todo
temp_ilut = apply_excit_ops(iLutRef, &
[elec_i, elec_j, elec_k], [orb_a, orb_b, orb_c])
call decode_bit_det(temp_nI, temp_ilut)
call hash_table_lookup(temp_nI, temp_ilut, &
nifd, HashIndex, CurrentDets, dummy_ind, &
dummy_hash, tSuccess)
if (.not. tSuccess) then
amp = -cc_ops(1)%get_amp(ia) * cc_ops(2)%get_amp(jk_bc) &
+ cc_ops(1)%get_amp(ib) * cc_ops(2)%get_amp(jk_ac) &
- cc_ops(1)%get_amp(ic) * cc_ops(2)%get_amp(jk_ab) &
+ cc_ops(1)%get_amp(ja) * cc_ops(2)%get_amp(ik_bc) &
- cc_ops(1)%get_amp(jb) * cc_ops(2)%get_amp(ik_ac) &
+ cc_ops(1)%get_amp(jc) * cc_ops(2)%get_amp(ik_ab) &
- cc_ops(1)%get_amp(ka) * cc_ops(2)%get_amp(ij_bc) &
+ cc_ops(1)%get_amp(kb) * cc_ops(2)%get_amp(ij_ac) &
- cc_ops(1)%get_amp(kc) * cc_ops(2)%get_amp(ij_ab) &
- cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(kc) &
+ cc_ops(1)%get_amp(ia) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(kb) &
+ cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kc) &
- cc_ops(1)%get_amp(ib) * cc_ops(1)%get_amp(jc) * cc_ops(1)%get_amp(ka) &
- cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(ja) * cc_ops(1)%get_amp(kb) &
+ cc_ops(1)%get_amp(ic) * cc_ops(1)%get_amp(jb) * cc_ops(1)%get_amp(ka)
if (abs(amp) > EPS) then
! then add this to the list!
temp_amps(n) = amp
temp_ops(n, 1, :) = [elec_i, elec_j, elec_k]
temp_ops(n, 2, :) = [orb_a, orb_b, orb_c]
n = n + 1
cc_amp_norm(0, 3) = cc_amp_norm(0, 3) + 1
cc_amp_norm(1, 3) = cc_amp_norm(1, 3) + abs(amp)
cc_amp_norm(2, 3) = cc_amp_norm(2, 3) + amp**2
end if
end if
end do
end do
end do
end do
end do
end do
! end if
cc_amp_norm(2, 3) = sqrt(cc_amp_norm(2, 3))
print *, "triples after T1*T2 and T1^3: ", n - 1, "on proc: ", &
iProcIndex
! and then allocate the actual array:
allocate(cc_ops(3)%amplitudes(n - 1))
allocate(cc_ops(3)%operators(n - 1, 2, 3))
cc_ops(3)%amplitudes = temp_amps(1:n - 1)
cc_ops(3)%operators = temp_ops(1:n - 1, :, :)
cc_ops(3)%n_ops = n - 1
end if
if (cc_order >= 4) then
! to compare the fciqmc amplitudes and the CC amplitudes I
! also need to store the quadrupels of the wavefunction
! NOTE: for now only do that if there are no singles!
! as this allows us to only consider the T_2^2 influence!
if (n_singles /= 0) then
call stop_all(this_routine, &
"T_4 only implemented for zero singles!")
end if
if (t_store_hash_quadrupels) then
quad_hash_size_wf = n_doubles * (n_doubles - 1) / 2
call init_cc_hash(quad_hash_wf, quad_hash_size_wf)
end if
do idet = 1, int(totwalkers)
ic = FindBitExcitLevel(iLutRef(:, 1), CurrentDets(:, idet))
call extract_sign(CurrentDets(:, idet), sign_tmp)
if (abs(sign_tmp(1)) > EPS) then
if (ic == 4) then
call get_bit_excitmat(iLutRef(:, 1), CurrentDets(:, idet), ex, ic)
! convert it to the unique quad-index
quad_ind = [ex(1, 1:2), ex(2, 1:2), ex(1, 3:4), ex(2, 3:4)]
! for some strange reason i decided to have a different
! ordering for this than above.. change that! TODO
temp_int = apply_excit_ops(iLutRef(:, 1), &
[ex(1, :)], [ex(2, :)])
! this should never find..
call cc_hash_look_up(quad_ind, temp_int, quad_hash_wf, &
hash_ind, t_found)
if (t_found) then
call stop_all(this_routine, &
"repeated quad not possible from wavefunction!")
end if
call cc_hash_add(quad_hash_wf, hash_ind, temp_int, sign_tmp(1))
cc_amp_norm(0, 4) = cc_amp_norm(0, 4) + 1
end if
end if
end do
print *, "cc-4 L0 norm of wavefunction: ", n_quads, "on proc: ", iProcIndex
call calc_cc_quad_norm(quad_hash_wf, quad_hash_size_wf)
if (iProcIndex == root) then
print *, "---test quads from wf: "
print *, "direct sampled quads: ", n_quads
print *, "L0 Norm of quadrupels in FCIMC-wavefunction: ", &
cc_amp_norm(0, 4)
print *, "L1 Norm of quadrupels in FCIMC-wavefunction: ", &
cc_amp_norm(1, 4)
print *, "L2 Norm of quadrupels in FCIMC-wavefunction: ", &
cc_amp_norm(2, 4)
end if
end if
end subroutine fill_cc_amplitudes