LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-01-12 04:58:08 +01:00
Code Releases Activity

Merge branch 'dev' into features_spack

Browse Source
This commit is contained in:
Anthony Scemama 2020-02-27 10:19:40 +01:00
commit a712441069
5 changed files with 148 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
configure vendored
View File

@ -1,4 +1,4 @@
#!/bin/bash #!/bin/bash
# #
# Quantum Package configuration script # Quantum Package configuration script
# #
@ -76,7 +76,7 @@ Usage:
Options: Options:
-c, --config=<file> Define a COMPILATION configuration file, -c, --config=<file> Define a COMPILATION configuration file,
in "${QP_ROOT}/config/". in "${QP_ROOT}/config/".
-h, --help Print the HELP message -h, --help Print the HELP message
Example: Example:
@ -103,7 +103,7 @@ function execute () {
while read -r line; do while read -r line; do
echo " " $line echo " " $line
_command+="${line} ;" _command+="${line} ;"
done done
sleep 1 sleep 1
echo "" echo ""
printf "\e[0;94m" printf "\e[0;94m"
@ -314,14 +314,12 @@ printf "\e[m\n"
if [[ -f ${QP_ROOT}/build.ninja ]] ; then if [[ -f ${QP_ROOT}/build.ninja ]] ; then
[[ -z ${TRAVIS} ]] && echo "You can now run ${QP_ROOT}/bin/qpsh to enter in the QP shell mode :)" [[ -z ${TRAVIS} ]] && echo "You can now run ${QP_ROOT}/bin/qpsh to enter in the QP shell mode :)"
else
echo "" echo ""
echo "${QP_ROOT}/build.ninja does not exist," echo "${QP_ROOT}/build.ninja does not exist,"
echo "you need to specify the COMPILATION configuration file." echo "you need to specify the COMPILATION configuration file."
echo "See ./configure --help for more detail." echo "See ./configure --help for more detail."
echo "" echo ""
fi fi
exit 0 exit 0

7
scripts/generate_h_apply.py
View File

@ -11,7 +11,6 @@ declarations
decls_main decls_main
deinit_thread deinit_thread
init_main init_main
filter_integrals
filter2p filter2p
filter2h2p_double filter2h2p_double
filter2h2p_single filter2h2p_single
@ -106,12 +105,6 @@ class H_apply(object):
s["do_double_excitations"] = d[do_double_exc] s["do_double_excitations"] = d[do_double_exc]
s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)" s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)"
s["filter_integrals"] = "array_pairs = .True."
if SingleRef:
s["filter_integrals"] = """
call get_mo_bielec_integrals_existing_ik(i_a,j_a,mo_num,array_pairs,mo_integrals_map)
"""
s["generate_psi_guess"] = """ s["generate_psi_guess"] = """
! Sort H_jj to find the N_states lowest states ! Sort H_jj to find the N_states lowest states
integer :: i integer :: i

160
src/cipsi/selection.irp.f
View File

@ -13,7 +13,7 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Weights adjusted along the selection to make the variances ! Weights adjusted along the selection to make the variances
! of each state coincide. ! of each state coincide.
END_DOC END_DOC
variance_match_weight(:) = 1.d0 variance_match_weight(:) = 1.d0
@ -46,10 +46,10 @@ subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
i_iter = 1 i_iter = 1
endif endif
dt = 4.d0 dt = 4.d0
do k=1,N_st do k=1,N_st
rpt2(k) = pt2(k)/(1.d0 + norm(k)) rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo enddo
avg = sum(rpt2(1:N_st)) / dble(N_st) avg = sum(rpt2(1:N_st)) / dble(N_st)
@ -90,7 +90,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
case (1) case (1)
print *, 'Using 1/c_max^2 weight in selection' print *, 'Using 1/c_max^2 weight in selection'
selection_weight(1:N_states) = c0_weight(1:N_states) selection_weight(1:N_states) = c0_weight(1:N_states)
case (2) case (2)
print *, 'Using pt2-matching weight in selection' print *, 'Using pt2-matching weight in selection'
@ -114,7 +114,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
print *, '# var weight ', real(variance_match_weight(:),4) print *, '# var weight ', real(variance_match_weight(:),4)
case (6) case (6)
print *, 'Using CI coefficient-based selection' print *, 'Using CI coefficient-based selection'
selection_weight(1:N_states) = c0_weight(1:N_states) selection_weight(1:N_states) = c0_weight(1:N_states)
case (7) case (7)
@ -289,34 +289,34 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
monoAdo = .true. monoAdo = .true.
monoBdo = .true. monoBdo = .true.
do k=1,N_int do k=1,N_int
hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1)) hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1))
hole (k,2) = iand(psi_det_generators(k,2,i_generator), hole_mask(k,2)) hole (k,2) = iand(psi_det_generators(k,2,i_generator), hole_mask(k,2))
particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1)) particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), particle_mask(k,1))
particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2)) particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2))
enddo enddo
integer :: N_holes(2), N_particles(2) integer :: N_holes(2), N_particles(2)
integer :: hole_list(N_int*bit_kind_size,2) integer :: hole_list(N_int*bit_kind_size,2)
integer :: particle_list(N_int*bit_kind_size,2) integer :: particle_list(N_int*bit_kind_size,2)
call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) call bitstring_to_list_ab(hole , hole_list , N_holes , N_int)
call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) call bitstring_to_list_ab(particle, particle_list, N_particles, N_int)
integer :: l_a, nmax, idx integer :: l_a, nmax, idx
integer, allocatable :: indices(:), exc_degree(:), iorder(:) integer, allocatable :: indices(:), exc_degree(:), iorder(:)
allocate (indices(N_det), & allocate (indices(N_det), &
exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
k=1 k=1
do i=1,N_det_alpha_unique do i=1,N_det_alpha_unique
call get_excitation_degree_spin(psi_det_alpha_unique(1,i), & call get_excitation_degree_spin(psi_det_alpha_unique(1,i), &
psi_det_generators(1,1,i_generator), exc_degree(i), N_int) psi_det_generators(1,1,i_generator), exc_degree(i), N_int)
enddo enddo
do j=1,N_det_beta_unique do j=1,N_det_beta_unique
call get_excitation_degree_spin(psi_det_beta_unique(1,j), & call get_excitation_degree_spin(psi_det_beta_unique(1,j), &
psi_det_generators(1,2,i_generator), nt, N_int) psi_det_generators(1,2,i_generator), nt, N_int)
@ -332,12 +332,12 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
endif endif
enddo enddo
enddo enddo
do i=1,N_det_beta_unique do i=1,N_det_beta_unique
call get_excitation_degree_spin(psi_det_beta_unique(1,i), & call get_excitation_degree_spin(psi_det_beta_unique(1,i), &
psi_det_generators(1,2,i_generator), exc_degree(i), N_int) psi_det_generators(1,2,i_generator), exc_degree(i), N_int)
enddo enddo
do j=1,N_det_alpha_unique do j=1,N_det_alpha_unique
call get_excitation_degree_spin(psi_det_alpha_unique(1,j), & call get_excitation_degree_spin(psi_det_alpha_unique(1,j), &
psi_det_generators(1,1,i_generator), nt, N_int) psi_det_generators(1,1,i_generator), nt, N_int)
@ -356,7 +356,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
endif endif
enddo enddo
enddo enddo
deallocate(exc_degree) deallocate(exc_degree)
nmax=k-1 nmax=k-1
@ -366,18 +366,18 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
enddo enddo
call isort(indices,iorder,nmax) call isort(indices,iorder,nmax)
deallocate(iorder) deallocate(iorder)
! Start with 32 elements. Size will double along with the filtering. ! Start with 32 elements. Size will double along with the filtering.
allocate(preinteresting(0:32), prefullinteresting(0:32), & allocate(preinteresting(0:32), prefullinteresting(0:32), &
interesting(0:32), fullinteresting(0:32)) interesting(0:32), fullinteresting(0:32))
preinteresting(:) = 0 preinteresting(:) = 0
prefullinteresting(:) = 0 prefullinteresting(:) = 0
do i=1,N_int do i=1,N_int
negMask(i,1) = not(psi_det_generators(i,1,i_generator)) negMask(i,1) = not(psi_det_generators(i,1,i_generator))
negMask(i,2) = not(psi_det_generators(i,2,i_generator)) negMask(i,2) = not(psi_det_generators(i,2,i_generator))
end do end do
do k=1,nmax do k=1,nmax
i = indices(k) i = indices(k)
mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i)) mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
@ -388,10 +388,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i)) mobMask(j,2) = iand(negMask(j,2), psi_det_sorted(j,2,i))
nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) nt = nt + popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2))
end do end do
if(nt <= 4) then if(nt <= 4) then
if(i <= N_det_selectors) then if(i <= N_det_selectors) then
sze = preinteresting(0) sze = preinteresting(0)
if (sze+1 == size(preinteresting)) then if (sze+1 == size(preinteresting)) then
allocate (tmp_array(0:sze)) allocate (tmp_array(0:sze))
tmp_array(0:sze) = preinteresting(0:sze) tmp_array(0:sze) = preinteresting(0:sze)
@ -403,7 +403,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
preinteresting(0) = sze+1 preinteresting(0) = sze+1
preinteresting(sze+1) = i preinteresting(sze+1) = i
else if(nt <= 2) then else if(nt <= 2) then
sze = prefullinteresting(0) sze = prefullinteresting(0)
if (sze+1 == size(prefullinteresting)) then if (sze+1 == size(prefullinteresting)) then
allocate (tmp_array(0:sze)) allocate (tmp_array(0:sze))
tmp_array(0:sze) = prefullinteresting(0:sze) tmp_array(0:sze) = prefullinteresting(0:sze)
@ -422,17 +422,17 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
! !$OMP CRITICAL ! !$OMP CRITICAL
! print *, 'Step1: ', i_generator, preinteresting(0) ! print *, 'Step1: ', i_generator, preinteresting(0)
! !$OMP END CRITICAL ! !$OMP END CRITICAL
allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2)) allocate(banned(mo_num, mo_num,2), bannedOrb(mo_num, 2))
allocate (mat(N_states, mo_num, mo_num)) allocate (mat(N_states, mo_num, mo_num))
maskInd = -1 maskInd = -1
integer :: nb_count, maskInd_save integer :: nb_count, maskInd_save
logical :: monoBdo_save logical :: monoBdo_save
logical :: found logical :: found
do s1=1,2 do s1=1,2
do i1=N_holes(s1),1,-1 ! Generate low excitations first do i1=N_holes(s1),1,-1 ! Generate low excitations first
found = .False. found = .False.
monoBdo_save = monoBdo monoBdo_save = monoBdo
maskInd_save = maskInd maskInd_save = maskInd
@ -447,21 +447,21 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
enddo enddo
if(s1 /= s2) monoBdo = .false. if(s1 /= s2) monoBdo = .false.
enddo enddo
if (.not.found) cycle if (.not.found) cycle
monoBdo = monoBdo_save monoBdo = monoBdo_save
maskInd = maskInd_save maskInd = maskInd_save
h1 = hole_list(i1,s1) h1 = hole_list(i1,s1)
call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int) call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int)
negMask = not(pmask) negMask = not(pmask)
interesting(0) = 0 interesting(0) = 0
fullinteresting(0) = 0 fullinteresting(0) = 0
do ii=1,preinteresting(0) do ii=1,preinteresting(0)
i = preinteresting(ii) i = preinteresting(ii)
select case (N_int) select case (N_int)
case (1) case (1)
mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i)) mobMask(1,1) = iand(negMask(1,1), psi_det_sorted(1,1,i))
@ -515,9 +515,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
endif endif
end do end do
end select end select
if(nt <= 4) then if(nt <= 4) then
sze = interesting(0) sze = interesting(0)
if (sze+1 == size(interesting)) then if (sze+1 == size(interesting)) then
allocate (tmp_array(0:sze)) allocate (tmp_array(0:sze))
tmp_array(0:sze) = interesting(0:sze) tmp_array(0:sze) = interesting(0:sze)
@ -529,7 +529,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
interesting(0) = sze+1 interesting(0) = sze+1
interesting(sze+1) = i interesting(sze+1) = i
if(nt <= 2) then if(nt <= 2) then
sze = fullinteresting(0) sze = fullinteresting(0)
if (sze+1 == size(fullinteresting)) then if (sze+1 == size(fullinteresting)) then
allocate (tmp_array(0:sze)) allocate (tmp_array(0:sze))
tmp_array(0:sze) = fullinteresting(0:sze) tmp_array(0:sze) = fullinteresting(0:sze)
@ -542,9 +542,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
fullinteresting(sze+1) = i fullinteresting(sze+1) = i
end if end if
end if end if
end do end do
do ii=1,prefullinteresting(0) do ii=1,prefullinteresting(0)
i = prefullinteresting(ii) i = prefullinteresting(ii)
nt = 0 nt = 0
@ -560,7 +560,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
end do end do
if(nt <= 2) then if(nt <= 2) then
sze = fullinteresting(0) sze = fullinteresting(0)
if (sze+1 == size(fullinteresting)) then if (sze+1 == size(fullinteresting)) then
allocate (tmp_array(0:sze)) allocate (tmp_array(0:sze))
tmp_array(0:sze) = fullinteresting(0:sze) tmp_array(0:sze) = fullinteresting(0:sze)
@ -577,7 +577,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
allocate (fullminilist (N_int, 2, fullinteresting(0)), & allocate (fullminilist (N_int, 2, fullinteresting(0)), &
minilist (N_int, 2, interesting(0)) ) minilist (N_int, 2, interesting(0)) )
if(pert_2rdm)then if(pert_2rdm)then
allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
do i=1,fullinteresting(0) do i=1,fullinteresting(0)
do j = 1, N_states do j = 1, N_states
coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
@ -588,7 +588,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
do i=1,fullinteresting(0) do i=1,fullinteresting(0)
fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i)) fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i))
enddo enddo
do i=1,interesting(0) do i=1,interesting(0)
minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i)) minilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,interesting(i))
enddo enddo
@ -624,21 +624,21 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
monoAdo = .false. monoAdo = .false.
end if end if
end if end if
maskInd = maskInd + 1 maskInd = maskInd + 1
if(mod(maskInd, csubset) == (subset-1)) then if(mod(maskInd, csubset) == (subset-1)) then
call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
if(fullMatch) cycle if(fullMatch) cycle
! !$OMP CRITICAL ! !$OMP CRITICAL
! print *, 'Step3: ', i_generator, h1, interesting(0) ! print *, 'Step3: ', i_generator, h1, interesting(0)
! !$OMP END CRITICAL ! !$OMP END CRITICAL
call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
if(.not.pert_2rdm)then if(.not.pert_2rdm)then
call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf) call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf)
else else
call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
endif endif
end if end if
@ -682,7 +682,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
double precision, allocatable :: values(:) double precision, allocatable :: values(:)
integer, allocatable :: keys(:,:) integer, allocatable :: keys(:,:)
integer :: nkeys integer :: nkeys
if(sp == 3) then if(sp == 3) then
s1 = 1 s1 = 1
@ -712,7 +712,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! to a determinant of the future. In that case, the determinant will be ! to a determinant of the future. In that case, the determinant will be
! detected as already generated when generating in the future with a ! detected as already generated when generating in the future with a
! double excitation. ! double excitation.
! !
! if (.not.do_singles) then ! if (.not.do_singles) then
! if ((h1 == p1) .or. (h2 == p2)) then ! if ((h1 == p1) .or. (h2 == p2)) then
! cycle ! cycle
@ -783,6 +783,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
norm(istate) = norm(istate) + coef * coef norm(istate) = norm(istate) + coef * coef
!!!DEBUG !!!DEBUG
! pt2(istate) = pt2(istate) - e_pert + alpha_h_psi**2/delta_E
!
! integer :: k ! integer :: k
! double precision :: alpha_h_psi_2,hij ! double precision :: alpha_h_psi_2,hij
! alpha_h_psi_2 = 0.d0 ! alpha_h_psi_2 = 0.d0
@ -792,7 +794,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! enddo ! enddo
! if(dabs(alpha_h_psi_2 - alpha_h_psi).gt.1.d-12)then ! if(dabs(alpha_h_psi_2 - alpha_h_psi).gt.1.d-12)then
! call debug_det(psi_det_generators(1,1,i_generator),N_int) ! call debug_det(psi_det_generators(1,1,i_generator),N_int)
! call debug_det(det,N_int) ! call debug_det(det,N_int)
! print*,'alpha_h_psi,alpha_h_psi_2 = ',alpha_h_psi,alpha_h_psi_2 ! print*,'alpha_h_psi,alpha_h_psi_2 = ',alpha_h_psi,alpha_h_psi_2
! stop ! stop
! endif ! endif
@ -816,7 +818,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
if(pseudo_sym)then if(pseudo_sym)then
if(dabs(mat(1, p1, p2)).lt.thresh_sym)then if(dabs(mat(1, p1, p2)).lt.thresh_sym)then
w = 0.d0 w = 0.d0
endif endif
endif endif
@ -857,7 +859,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
negMask(i,2) = not(mask(i,2)) negMask(i,2) = not(mask(i,2))
end do end do
do i=1, N_sel do i=1, N_sel
if (interesting(i) < 0) then if (interesting(i) < 0) then
stop 'prefetch interesting(i) and det(i)' stop 'prefetch interesting(i) and det(i)'
endif endif
@ -1210,7 +1212,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
endif endif
end if end if
! enddo ! enddo
! !
putj = p2 putj = p2
! do puti=1,mo_num !HOT ! do puti=1,mo_num !HOT
if(.not. banned(putj,puti,bant)) then if(.not. banned(putj,puti,bant)) then
@ -1573,15 +1575,15 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
double precision, intent(in) :: coefs(N_states) double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num) double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp integer, intent(in) :: h(0:2,2), p(0:4,2), sp
integer :: i, j, s, h1, h2, p1, p2, puti, putj integer :: i, j, s, h1, h2, p1, p2, puti, putj
double precision :: hij, phase double precision :: hij, phase
double precision, external :: get_phase_bi, mo_two_e_integral double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok logical :: ok
integer :: bant integer :: bant
bant = 1 bant = 1
if(sp == 3) then ! AB if(sp == 3) then ! AB
h1 = p(1,1) h1 = p(1,1)
@ -1619,7 +1621,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
end do end do
end do end do
end if end if
end end
subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks use bitmasks
@ -1639,27 +1641,27 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
logical, allocatable :: lbanned(:,:) logical, allocatable :: lbanned(:,:)
integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j
integer :: hfix, pfix, h1, h2, p1, p2, ib integer :: hfix, pfix, h1, h2, p1, p2, ib
integer, parameter :: turn2(2) = (/2,1/) integer, parameter :: turn2(2) = (/2,1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant integer :: bant
allocate (lbanned(mo_num, 2)) allocate (lbanned(mo_num, 2))
lbanned = bannedOrb lbanned = bannedOrb
do i=1, p(0,1) do i=1, p(0,1)
lbanned(p(i,1), 1) = .true. lbanned(p(i,1), 1) = .true.
end do end do
do i=1, p(0,2) do i=1, p(0,2)
lbanned(p(i,2), 2) = .true. lbanned(p(i,2), 2) = .true.
end do end do
ma = 1 ma = 1
if(p(0,2) >= 2) ma = 2 if(p(0,2) >= 2) ma = 2
mi = turn2(ma) mi = turn2(ma)
bant = 1 bant = 1
if(sp == 3) then if(sp == 3) then
@ -1682,7 +1684,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
tmp_row(1:N_states,putj) += hij * coefs(1:N_states) tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
end do end do
if(ma == 1) then if(ma == 1) then
mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num) mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num)
else else
mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num) mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num)
@ -1701,14 +1703,14 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
tmp_row(:,puti) += hij * coefs(:) tmp_row(:,puti) += hij * coefs(:)
end if end if
putj = p2 putj = p2
if(.not. banned(putj,puti,bant)) then if(.not. banned(putj,puti,bant)) then
hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
tmp_row2(:,puti) += hij * coefs(:) tmp_row2(:,puti) += hij * coefs(:)
end if end if
end do end do
if(mi == 1) then if(mi == 1) then
mat(:,:,p1) += tmp_row(:,:) mat(:,:,p1) += tmp_row(:,:)
mat(:,:,p2) += tmp_row2(:,:) mat(:,:,p2) += tmp_row2(:,:)
@ -1752,7 +1754,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
tmp_row(:,puti) += hij * coefs(:) tmp_row(:,puti) += hij * coefs(:)
end if end if
putj = p1 putj = p1
if(.not. banned(puti,putj,1)) then if(.not. banned(puti,putj,1)) then
hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
@ -1788,7 +1790,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
mat(:, p1, p2) += coefs(:) * hij mat(:, p1, p2) += coefs(:) * hij
end do end do
end do end do
end end
subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs) subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks use bitmasks
@ -1800,30 +1802,30 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
double precision, intent(in) :: coefs(N_states) double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num) double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi, mo_two_e_integral double precision, external :: get_phase_bi, mo_two_e_integral
integer :: i, j, tip, ma, mi, puti, putj integer :: i, j, tip, ma, mi, puti, putj
integer :: h1, h2, p1, p2, i1, i2 integer :: h1, h2, p1, p2, i1, i2
double precision :: hij, phase double precision :: hij, phase
integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/)) integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/))
integer, parameter :: turn2(2) = (/2, 1/) integer, parameter :: turn2(2) = (/2, 1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant integer :: bant
bant = 1 bant = 1
tip = p(0,1) * p(0,2) tip = p(0,1) * p(0,2)
ma = sp ma = sp
if(p(0,1) > p(0,2)) ma = 1 if(p(0,1) > p(0,2)) ma = 1
if(p(0,1) < p(0,2)) ma = 2 if(p(0,1) < p(0,2)) ma = 2
mi = mod(ma, 2) + 1 mi = mod(ma, 2) + 1
if(sp == 3) then if(sp == 3) then
if(ma == 2) bant = 2 if(ma == 2) bant = 2
if(tip == 3) then if(tip == 3) then
puti = p(1, mi) puti = p(1, mi)
do i = 1, 3 do i = 1, 3
@ -1835,7 +1837,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p2 = p(i2, ma) p2 = p(i2, ma)
h1 = h(1, ma) h1 = h(1, ma)
h2 = h(2, ma) h2 = h(2, ma)
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
if(ma == 1) then if(ma == 1) then
mat(:, putj, puti) += coefs(:) * hij mat(:, putj, puti) += coefs(:) * hij
@ -1851,10 +1853,10 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p2 = p(turn2(j), 2) p2 = p(turn2(j), 2)
do i = 1,2 do i = 1,2
puti = p(i, 1) puti = p(i, 1)
if(banned(puti,putj,bant)) cycle if(banned(puti,putj,bant)) cycle
p1 = p(turn2(i), 1) p1 = p(turn2(i), 1)
hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij mat(:, puti, putj) += coefs(:) * hij
end do end do
@ -1870,7 +1872,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
do j=i+1,4 do j=i+1,4
putj = p(j, ma) putj = p(j, ma)
if(banned(puti,putj,1)) cycle if(banned(puti,putj,1)) cycle
i1 = turn2d(1, i, j) i1 = turn2d(1, i, j)
i2 = turn2d(2, i, j) i2 = turn2d(2, i, j)
p1 = p(i1, ma) p1 = p(i1, ma)
@ -1888,7 +1890,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
putj = p(turn3(2,i), ma) putj = p(turn3(2,i), ma)
if(banned(puti,putj,1)) cycle if(banned(puti,putj,1)) cycle
p2 = p(i, ma) p2 = p(i, ma)
hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int)
mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij
end do end do
@ -1905,6 +1907,6 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
end if end if
end if end if
end if end if
end end

34
src/tools/print_energy.irp.f Normal file
View File

@ -0,0 +1,34 @@
program print_energy
implicit none
BEGIN_DOC
! Prints the energy of the wave function stored in the |EZFIO| directory.
END_DOC
! this has to be done in order to be sure that N_det, psi_det and
! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
read_wf = .True.
touch read_wf
call run
end
subroutine run
implicit none
integer :: i
double precision :: i_H_psi_array(N_states)
double precision :: E(N_states)
double precision :: norm(N_states)
E(:) = nuclear_repulsion
norm(:) = 0.d0
do i=1,N_det
call i_H_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
size(psi_coef,1), N_states, i_H_psi_array)
norm(:) += psi_coef(i,:)**2
E(:) += i_H_psi_array(:) * psi_coef(i,:)
enddo
print *, 'Energy:'
do i=1,N_states
print *, E(i)/norm(i)
enddo
end

29
src/tools/print_hamiltonian.irp.f Normal file
View File

@ -0,0 +1,29 @@
program print_hamiltonian
implicit none
BEGIN_DOC
! Prints the Hamiltonian matrix defined in the space of determinants
! present in the |EZFIO| directory.
END_DOC
! this has to be done in order to be sure that N_det, psi_det and
! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
read_wf = .True.
touch read_wf
call run
end
subroutine run
implicit none
integer :: i, j
double precision :: hij
do j=1,N_det
do i=1,N_det
call i_H_j(psi_det(1,1,i), psi_det(1,1,j), N_int, hij)
if (dabs(hij) > 1.d-20) then
print *, i, j, hij
endif
enddo
enddo
end