LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Blanks at end of line

Browse Source
This commit is contained in:
Anthony Scemama 2020-02-27 10:16:27 +01:00
parent eac877dd44
commit 4741f2f544
1 changed files with 81 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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