LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

working mrcc

Browse Source
This commit is contained in:
Yann Garniron 2016-05-20 09:44:22 +02:00
parent 55138005a0
commit 33bd506328
14 changed files with 733 additions and 433 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/gfortran_debug.cfg
View File

@ -51,7 +51,7 @@ FCFLAGS : -Ofast
# -g : Extra debugging information
#
[DEBUG]
FCFLAGS : -g -pedantic -msse4.2
FCFLAGS : -g -msse4.2
# OpenMP flags
#################

2
plugins/MRCC_Utils/davidson.irp.f
View File

@ -269,7 +269,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
to_print(2,k) = residual_norm(k)
enddo
write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st)
write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))') iter, to_print(:,1:N_st)
call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
if (converged) then
exit

72
plugins/MRCC_Utils/mrcc_dress.irp.f
View File

@ -51,9 +51,9 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:)
integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:)
integer :: mobiles(2), smallerlist
logical, external :: is_generable
print *, i_generator
leng = max(N_det_generators, N_det_non_ref)
allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref))
@ -69,7 +69,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
allocate( microlist(Nint,2,N_minilist*4), &
idx_microlist(N_minilist*4))
if(key_mask(1,1) /= 0) then
if(key_mask(1,1) /= 0_8) then
call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
else
@ -87,6 +87,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
! |alpha>
if(N_tq > 0) then
call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint)
if(N_minilist == 0) return
@ -117,8 +118,18 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
do i_alpha=1,N_tq
! ok = .false.
! do i=N_det_generators, 1, -1
! if(is_generable(psi_det_generators(1,1,i), tq(1,1,i_alpha), Nint)) then
! ok = .true.
! exit
! end if
! end do
! if(.not. ok) then
! cycle
! end if
if(key_mask(1,1) /= 0) then
call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint)
@ -138,37 +149,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
do j=1,idx_alpha(0)
idx_alpha(j) = idx_microlist_zero(idx_alpha(j))
end do
! i = 1
! j = 2
! do j = 2, idx_alpha_tmp(0)
! if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit
! end do
!
! m = j
!
! idx_alpha(0) = idx_alpha_tmp(0)
!
! do l = 1, idx_alpha(0)
! if(j > idx_alpha_tmp(0)) then
! k = i
! i += 1
! else if(i >= m) then
! k = j
! j += 1
! else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then
! k = i
! i += 1
! else
! k = j
! j += 1
! end if
! ! k=l
! idx_alpha(l) = idx_alpha_tmp(k)
! degree_alpha(l) = degree_alpha_tmp(k)
! end do
!
else
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
do j=1,idx_alpha(0)
@ -177,12 +157,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
end if
! call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
! do j=1,idx_alpha(0)
! idx_alpha(j) = idx_miniList(idx_alpha(j))
! end do
!print *, idx_alpha(:idx_alpha(0))
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
@ -285,6 +259,8 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
enddo
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do i_state=1,Nstates
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
do l_sd=1,idx_alpha(0)
@ -293,7 +269,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
enddo
else
delta_ii_(i_state,i_I) = 0.d0
!delta_ii_(i_state,i_I) = 0.d0
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
@ -303,15 +279,11 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
!deallocate (dIa_hla,hij_cache)
!deallocate(miniList, idx_miniList)
deallocate (dIa_hla,hij_cache)
deallocate(miniList, idx_miniList)
end
subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
use bitmasks
@ -360,7 +332,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
endif
enddo
if (good) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
N_tq += 1
do k=1,N_int
tq(k,1,N_tq) = det_buffer(k,1,i)
@ -437,7 +409,7 @@ subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,N
endif
enddo
if (good) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
N_tq += 1
do k=1,N_int
tq(k,1,N_tq) = det_buffer(k,1,i)

455
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -112,20 +112,12 @@ END_PROVIDER
lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i
endif
endif
! j = int(lambda_mrcc(k,i) * 100)
! if(j < -200) j = -200
! if(j > 200) j = 200
! histo(j) += 1
enddo
enddo
lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2
end if
! do i=-200,200
! print *, i, histo(i)
! end do
print*,'N_det_non_ref = ',N_det_non_ref
!print*,'Number of ignored determinants = ',i_pert_count
print*,'Number of ignored determinants = ',i_pert_count
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
print*,'lambda max = ',maxval(dabs(lambda_mrcc))
@ -164,6 +156,7 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
implicit none
BEGIN_DOC
@ -288,3 +281,447 @@ subroutine diagonalize_CI_dressed(lambda)
SOFT_TOUCH psi_coef
end
logical function is_generable(det1, det2, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
integer :: degree, f, exc(0:2, 2, 2), h1, h2, p1, p2, s1, s2, t
integer, external :: searchExc
logical, external :: excEq
double precision :: phase
is_generable = .false.
call get_excitation(det1, det2, exc, degree, phase, Nint)
if(degree == -1) return
if(degree == 0) then
is_generable = .true.
return
end if
if(degree > 2) stop "?22??"
!!!!!
! call dec_exc(exc, h1, h2, p1, p2)
! f = searchExc(toutmoun(1,1), (/h1, h2, p1, p2/), hh_shortcut(hh_shortcut(0)+1)-1)
! !print *, toutmoun(:,1), hh_shortcut(hh_shortcut(0)+1)-1, (/h1, h2, p1, p2/)
! if(f /= -1) then
! is_generable = .true.
! if(.not. excEq(toutmoun(1,f), (/h1, h2, p1, p2/))) stop "????"
! end if
! ! print *, f
! return
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
if(degree == 1) then
h2 = h1
p2 = p1
s2 = s1
h1 = 0
p1 = 0
s1 = 0
end if
if(h1 + s1*mo_tot_num < h2 + s2*mo_tot_num) then
f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0))
else
f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0))
end if
if(f == -1) return
if(p1 + s1*mo_tot_num < p2 + s2*mo_tot_num) then
f = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f))
else
f = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f))
end if
if(f /= -1) is_generable = .true.
end function
integer function searchDet(dets, det, n, Nint)
implicit none
use bitmasks
integer(bit_kind),intent(in) :: dets(Nint,2,n), det(Nint,2)
integer, intent(in) :: nint, n
integer :: l, h, c
integer, external :: detCmp
logical, external :: detEq
l = 1
h = n
do while(.true.)
searchDet = (l+h)/2
c = detCmp(dets(1,1,searchDet), det(:,:), Nint)
if(c == 0) return
if(c == 1) then
h = searchDet-1
else
l = searchDet+1
end if
if(l > h) then
searchDet = -1
return
end if
end do
end function
integer function searchExc(excs, exc, n)
implicit none
use bitmasks
integer, intent(in) :: n
integer,intent(in) :: excs(4,n), exc(4)
integer :: l, h, c
integer, external :: excCmp
logical, external :: excEq
l = 1
h = n
do
searchExc = (l+h)/2
c = excCmp(excs(1,searchExc), exc(:))
if(c == 0) return
if(c == 1) then
h = searchExc-1
else
l = searchExc+1
end if
if(l > h) then
searchExc = -1
return
end if
end do
end function
subroutine sort_det(key, idx, N_key, Nint)
implicit none
integer, intent(in) :: Nint, N_key
integer(8),intent(inout) :: key(Nint,2,N_key)
integer,intent(out) :: idx(N_key)
integer(8) :: tmp(Nint, 2)
integer :: tmpidx,i,ni
do i=1,N_key
idx(i) = i
end do
do i=N_key/2,1,-1
call tamiser(key, idx, i, N_key, Nint, N_key)
end do
do i=N_key,2,-1
do ni=1,Nint
tmp(ni,1) = key(ni,1,i)
tmp(ni,2) = key(ni,2,i)
key(ni,1,i) = key(ni,1,1)
key(ni,2,i) = key(ni,2,1)
key(ni,1,1) = tmp(ni,1)
key(ni,2,1) = tmp(ni,2)
enddo
tmpidx = idx(i)
idx(i) = idx(1)
idx(1) = tmpidx
call tamiser(key, idx, 1, i-1, Nint, N_key)
end do
end subroutine
subroutine sort_exc(key, N_key)
implicit none
integer, intent(in) :: N_key
integer,intent(inout) :: key(4,N_key)
integer :: tmp(4)
integer :: i,ni
do i=N_key/2,1,-1
call tamise_exc(key, i, N_key, N_key)
end do
do i=N_key,2,-1
do ni=1,4
tmp(ni) = key(ni,i)
key(ni,i) = key(ni,1)
key(ni,1) = tmp(ni)
enddo
call tamise_exc(key, 1, i-1, N_key)
end do
end subroutine
logical function exc_inf(exc1, exc2)
implicit none
integer,intent(in) :: exc1(4), exc2(4)
integer :: i
exc_inf = .false.
do i=1,4
if(exc1(i) < exc2(i)) then
exc_inf = .true.
return
else if(exc1(i) > exc2(i)) then
return
end if
end do
end function
subroutine tamise_exc(key, no, n, N_key)
use bitmasks
implicit none
BEGIN_DOC
! Uncodumented : TODO
END_DOC
integer,intent(in) :: no, n, N_key
integer,intent(inout) :: key(4, N_key)
integer :: k,j
integer :: tmp(4)
logical :: exc_inf
integer :: ni
k = no
j = 2*k
do while(j <= n)
if(j < n) then
if (exc_inf(key(1,j), key(1,j+1))) then
j = j+1
endif
endif
if(exc_inf(key(1,k), key(1,j))) then
do ni=1,4
tmp(ni) = key(ni,k)
key(ni,k) = key(ni,j)
key(ni,j) = tmp(ni)
enddo
k = j
j = k+k
else
return
endif
enddo
end subroutine
subroutine dec_exc(exc, h1, h2, p1, p2)
implicit none
integer :: exc(0:2,2,2), s1, s2, degree
integer, intent(out) :: h1, h2, p1, p2
degree = exc(0,1,1) + exc(0,1,2)
h1 = 0
h2 = 0
p1 = 0
p2 = 0
if(degree == 0) return
call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
h1 += mo_tot_num * (s1-1)
p1 += mo_tot_num * (s1-1)
if(degree == 2) then
h2 += mo_tot_num * (s2-1)
p2 += mo_tot_num * (s2-1)
if(h1 > h2) then
s1 = h1
h1 = h2
h2 = s1
end if
if(p1 > p2) then
s1 = p1
p1 = p2
p2 = s1
end if
else
h2 = h1
p2 = p1
p1 = 0
h1 = 0
end if
end subroutine
BEGIN_PROVIDER [ integer, hh_exists, (4, N_det_ref * N_det_non_ref) ]
&BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_det_ref * N_det_non_ref + 1) ]
&BEGIN_PROVIDER [ integer, pp_exists, (4, N_det_ref * N_det_non_ref) ]
implicit none
integer,allocatable :: num(:,:)
integer :: exc(0:2, 2, 2), degree, n, on, s, h1, h2, p1, p2, l, i
double precision :: phase
logical, external :: excEq
allocate(num(4, N_det_ref * N_det_non_ref))
hh_shortcut = 0
hh_exists = 0
pp_exists = 0
num = 0
n = 0
do i=1, N_det_ref
do l=1, N_det_non_ref
call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int)
if(degree == -1) cycle
call dec_exc(exc, h1, h2, p1, p2)
n += 1
num(:, n) = (/h1, h2, p1, p2/)
end do
end do
call sort_exc(num, n)
hh_shortcut(0) = 1
hh_shortcut(1) = 1
hh_exists(:,1) = (/1, num(1,1), 1, num(2,1)/)
pp_exists(:,1) = (/1, num(3,1), 1, num(4,1)/)
s = 1
do i=2,n
if(.not. excEq(num(1,i), num(1,s))) then
s += 1
num(:, s) = num(:, i)
pp_exists(:,s) = (/1, num(3,s), 1, num(4,s)/)
if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. &
hh_exists(4, hh_shortcut(0)) /= num(2,s)) then
hh_shortcut(0) += 1
hh_shortcut(hh_shortcut(0)) = s
hh_exists(:,hh_shortcut(0)) = (/1, num(1,s), 1, num(2,s)/)
end if
end if
end do
hh_shortcut(hh_shortcut(0)+1) = s+1
do s=2,4,2
do i=1,hh_shortcut(0)
if(hh_exists(s, i) == 0) then
hh_exists(s-1, i) = 0
else if(hh_exists(s, i) > mo_tot_num) then
hh_exists(s, i) -= mo_tot_num
hh_exists(s-1, i) = 2
end if
end do
do i=1,hh_shortcut(hh_shortcut(0)+1)-1
if(pp_exists(s, i) == 0) then
pp_exists(s-1, i) = 0
else if(pp_exists(s, i) > mo_tot_num) then
pp_exists(s, i) -= mo_tot_num
pp_exists(s-1, i) = 2
end if
end do
end do
END_PROVIDER
logical function excEq(exc1, exc2)
implicit none
integer, intent(in) :: exc1(4), exc2(4)
integer :: i
excEq = .false.
do i=1, 4
if(exc1(i) /= exc2(i)) return
end do
excEq = .true.
end function
integer function excCmp(exc1, exc2)
implicit none
integer, intent(in) :: exc1(4), exc2(4)
integer :: i
excCmp = 0
do i=1, 4
if(exc1(i) > exc2(i)) then
excCmp = 1
return
else if(exc1(i) < exc2(i)) then
excCmp = -1
return
end if
end do
end function
subroutine apply_hole(det, exc, res, ok, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: exc(4)
integer :: s1, s2, h1, h2
integer(bit_kind),intent(in) :: det(Nint, 2)
integer(bit_kind),intent(out) :: res(Nint, 2)
logical, intent(out) :: ok
integer :: ii, pos
ok = .false.
s1 = exc(1)
h1 = exc(2)
s2 = exc(3)
h2 = exc(4)
res = det
if(h1 /= 0) then
ii = (h1-1)/bit_kind_size + 1
pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s1) = ibclr(res(ii, s1), pos)
end if
ii = (h2-1)/bit_kind_size + 1
pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s2) = ibclr(res(ii, s2), pos)
ok = .true.
end subroutine
subroutine apply_particle(det, exc, res, ok, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: exc(4)
integer :: s1, s2, p1, p2
integer(bit_kind),intent(in) :: det(Nint, 2)
integer(bit_kind),intent(out) :: res(Nint, 2)
logical, intent(out) :: ok
integer :: ii, pos
ok = .false.
s1 = exc(1)
p1 = exc(2)
s2 = exc(3)
p2 = exc(4)
res = det
if(p1 /= 0) then
ii = (p1-1)/bit_kind_size + 1
pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s1) = ibset(res(ii, s1), pos)
end if
ii = (p2-1)/bit_kind_size + 1
pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s2) = ibset(res(ii, s2), pos)
ok = .true.
end subroutine

497
plugins/mrcepa0/dressing.irp.f
View File

@ -1,227 +1,32 @@
use bitmasks
subroutine dec_exc(exc, h1, h2, p1, p2)
implicit none
integer :: exc(0:2,2,2), s1, s2, degree
integer, intent(out) :: h1, h2, p1, p2
degree = exc(0,1,1) + exc(0,1,2)
h1 = 0
h2 = 0
p1 = 0
p2 = 0
if(degree == 0) return
call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2)
h1 += mo_tot_num * (s1-1)
p1 += mo_tot_num * (s1-1)
if(degree == 2) then
h2 += mo_tot_num * (s2-1)
p2 += mo_tot_num * (s2-1)
if(h1 > h2) then
s1 = h1
h1 = h2
h2 = s1
end if
if(p1 > p2) then
s1 = p1
p1 = p2
p2 = s1
end if
else
h2 = h1
p2 = p1
p1 = 0
h1 = 0
end if
end subroutine
BEGIN_PROVIDER [ integer, hh_exists, (4, N_det_ref * N_det_non_ref) ]
&BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_det_ref * N_det_non_ref + 1) ]
&BEGIN_PROVIDER [ integer, pp_exists, (4, N_det_ref * N_det_non_ref) ]
implicit none
integer :: num(0:mo_tot_num*2, 0:mo_tot_num*2)
integer :: exc(0:2, 2, 2), degree, n, on, s, h1, h2, p1, p2, l, i
double precision :: phase
hh_shortcut = 0
hh_exists = 0
pp_exists = 0
num = 0
do i=1, N_det_ref
do l=1, N_det_non_ref
call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int)
if(degree == -1) cycle
call dec_exc(exc, h1, h2, p1, p2)
num(h1, h2) += 1
end do
end do
n = 1
do l=0,mo_tot_num*2
do i=0,l
on = num(i,l)
if(on /= 0) then
hh_shortcut(0) += 1
hh_shortcut(hh_shortcut(0)) = n
hh_exists(:, hh_shortcut(0)) = (/1, i, 1, l/)
end if
num(i,l) = n
n += on
end do
end do
hh_shortcut(hh_shortcut(0)+1) = n
do i=1, N_det_ref
do l=1, N_det_non_ref
call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int)
if(degree == -1) cycle
call dec_exc(exc, h1, h2, p1, p2)
pp_exists(:, num(h1, h2)) = (/1,p1,1,p2/)
num(h1, h2) += 1
end do
end do
do s=2,4,2
do i=1,hh_shortcut(0)
if(hh_exists(s, i) == 0) then
hh_exists(s-1, i) = 0
else if(hh_exists(s, i) > mo_tot_num) then
hh_exists(s, i) -= mo_tot_num
hh_exists(s-1, i) = 2
end if
end do
do i=1,hh_shortcut(hh_shortcut(0)+1)-1
if(pp_exists(s, i) == 0) then
pp_exists(s-1, i) = 0
else if(pp_exists(s, i) > mo_tot_num) then
pp_exists(s, i) -= mo_tot_num
pp_exists(s-1, i) = 2
end if
end do
end do
END_PROVIDER
subroutine apply_hole(det, exc, res, ok, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: exc(4)
integer :: s1, s2, h1, h2
integer(bit_kind),intent(in) :: det(Nint, 2)
integer(bit_kind),intent(out) :: res(Nint, 2)
logical, intent(out) :: ok
integer :: ii, pos
ok = .false.
s1 = exc(1)
h1 = exc(2)
s2 = exc(3)
h2 = exc(4)
res = det
if(h1 /= 0) then
ii = (h1-1)/bit_kind_size + 1
pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s1) = ibclr(res(ii, s1), pos)
end if
ii = (h2-1)/bit_kind_size + 1
pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s2) = ibclr(res(ii, s2), pos)
ok = .true.
end subroutine
subroutine apply_particle(det, exc, res, ok, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: exc(4)
integer :: s1, s2, p1, p2
integer(bit_kind),intent(in) :: det(Nint, 2)
integer(bit_kind),intent(out) :: res(Nint, 2)
logical, intent(out) :: ok
integer :: ii, pos
ok = .false.
s1 = exc(1)
p1 = exc(2)
s2 = exc(3)
p2 = exc(4)
res = det
if(p1 /= 0) then
ii = (p1-1)/bit_kind_size + 1
pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s1) = ibset(res(ii, s1), pos)
end if
ii = (p2-1)/bit_kind_size + 1
pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s2) = ibset(res(ii, s2), pos)
ok = .true.
end subroutine
BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ]
use bitmasks
implicit none
integer :: gen, h, p, i_state, n, t
integer :: gen, h, p, i_state, n, t, i, h1, h2, p1, p2, s1, s2
integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2), buf(N_int, 2, N_det_non_ref)
logical :: ok
logical, external :: detEq
delta_ij_mrcc = 0d0
delta_ii_mrcc = 0d0
i_state = 1
do gen= 1, N_det_generators
!print *, gen, "/", N_det_generators
print *, gen, "/", N_det_generators
do h=1, hh_shortcut(0)
call apply_hole(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int)
if(.not. ok) cycle
omask = 0
omask = 0_bit_kind
if(hh_exists(1, h) /= 0) omask = mask
!-459.6378590456251
!-199.0659502581943
n = 1
ploop : do p=hh_shortcut(h), hh_shortcut(h+1)-1
do t=hh_shortcut(h), p-1
if(pp_exists(1, p) == pp_exists(1,t) .and. &
pp_exists(2, p) == pp_exists(2,t) .and. &
pp_exists(3, p) == pp_exists(3,t) .and. &
pp_exists(4, p) == pp_exists(4,t)) cycle ploop
end do
do p=hh_shortcut(h), hh_shortcut(h+1)-1
call apply_particle(mask, pp_exists(1, p), buf(1,1,n), ok, N_int)
!-459.6379081607463
!-199.0659982685706
if(ok) n = n + 1
end do ploop
end do
n = n - 1
if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask)
end do
@ -229,7 +34,6 @@ end subroutine
END_PROVIDER
subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffer,Nint,key_mask)
use bitmasks
implicit none
@ -258,7 +62,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
integer :: i_state, k_sd, l_sd, i_I, i_alpha
integer(bit_kind),allocatable :: miniList(:,:,:)
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
integer(bit_kind) :: key_mask(Nint, 2)
integer,allocatable :: idx_miniList(:)
integer :: N_miniList, ni, leng
double precision, allocatable :: hij_cache(:)
@ -266,7 +70,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:)
integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:)
integer :: mobiles(2), smallerlist
logical, external :: detEq, is_generable
leng = max(N_det_generators, N_det_non_ref)
@ -275,9 +79,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
if(fullMatch) then
return
end if
! if(fullMatch) then
! return
! end if
allocate(ptr_microlist(0:mo_tot_num*2+1), &
N_microlist(0:mo_tot_num*2) )
@ -286,9 +90,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
if(key_mask(1,1) /= 0) then
call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
call filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
else
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
call filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
end if
@ -463,6 +267,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
enddo
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do i_state=1,N_states
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
do l_sd=1,idx_alpha(0)
@ -481,8 +286,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
!deallocate (dIa_hla,hij_cache)
!deallocate(miniList, idx_miniList)
deallocate (dIa_hla,hij_cache)
deallocate(miniList, idx_miniList)
end
@ -494,29 +299,30 @@ end
implicit none
integer :: i, j, i_state
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
do i_state = 1, N_states
if(mrmode == 3) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
end do
end do
else if(mrmode == 2) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_ii_mrcc(i_state,i)
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i)
end do
end do
!
! do i = 1, N_det_ref
! delta_ii(i_state,i)= delta_ii_old(i_state,i)
! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
! do j = 1, N_det_non_ref
! delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i)
! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
! end do
! end do
else if(mrmode == 2) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_ii_old(i_state,i)
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i)
end do
end do
else if(mrmode == 1) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state)
@ -646,7 +452,6 @@ END_PROVIDER
end do
end do
end do
print *, "pre done"
END_PROVIDER
@ -708,21 +513,6 @@ END_PROVIDER
END_PROVIDER
logical function detEq(a,b,Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2)
integer :: ni, i
detEq = .false.
do i=1,2
do ni=1,Nint
if(a(ni,i) /= b(ni,i)) return
end do
end do
detEq = .true.
end function
logical function isInCassd(a,Nint)
@ -793,106 +583,6 @@ subroutine getHP(a,h,p,Nint)
!isInCassd = .true.
end function
integer function detCmp(a,b,Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2)
integer :: ni, i
detCmp = 0
do i=1,2
do ni=Nint,1,-1
if(a(ni,i) < b(ni,i)) then
detCmp = -1
return
else if(a(ni,i) > b(ni,i)) then
detCmp = 1
return
end if
end do
end do
end function
integer function searchDet(dets, det, n, Nint)
implicit none
use bitmasks
integer(bit_kind),intent(in) :: dets(Nint,2,n), det(Nint,2)
integer, intent(in) :: nint, n
integer :: l, h, c
integer, external :: detCmp
logical, external :: detEq
!do l=1,n
! if(detEq(det(1,1), dets(1,1,l),Nint)) then
! searchDet = l
! return
! end if
!end do
!searchDet = -1
!return
l = 1
h = n
do while(.true.)
searchDet = (l+h)/2
c = detCmp(dets(1,1,searchDet), det(:,:), Nint)
if(c == 0) return
if(c == 1) then
h = searchDet-1
else
l = searchDet+1
end if
if(l > h) then
searchDet = -1
return
end if
end do
end function
subroutine sort_det(key, idx, N_key, Nint)
implicit none
integer, intent(in) :: Nint, N_key
integer(8),intent(inout) :: key(Nint,2,N_key)
integer,intent(out) :: idx(N_key)
integer(8) :: tmp(Nint, 2)
integer :: tmpidx,i,ni
do i=1,N_key
idx(i) = i
end do
do i=N_key/2,1,-1
call tamiser(key, idx, i, N_key, Nint, N_key)
end do
do i=N_key,2,-1
do ni=1,Nint
tmp(ni,1) = key(ni,1,i)
tmp(ni,2) = key(ni,2,i)
key(ni,1,i) = key(ni,1,1)
key(ni,2,i) = key(ni,2,1)
key(ni,1,1) = tmp(ni,1)
key(ni,2,1) = tmp(ni,2)
enddo
tmpidx = idx(i)
idx(i) = idx(1)
idx(1) = tmpidx
call tamiser(key, idx, 1, i-1, Nint, N_key)
end do
end subroutine
BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ]
@ -1116,6 +806,7 @@ end subroutine
BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ]
implicit none
integer :: i,j
do i=1,N_det_ref
do j=1,N_det_non_ref
@ -1125,3 +816,135 @@ BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ]
END_PROVIDER
subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,m
logical :: is_in_wavefunction
integer :: degree(psi_det_size)
integer :: idx(0:psi_det_size)
logical :: good
integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
integer, intent(out) :: N_tq
integer :: nt,ni
logical, external :: is_connected_to, is_generable
integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_generators)
integer,intent(in) :: N_miniList
N_tq = 0
i_loop : do i=1,N_selected
do k=1, N_minilist
if(is_generable(miniList(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop
end do
! Select determinants that are triple or quadruple excitations
! from the ref
good = .True.
call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
!good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
do k=1,idx(0)
if (degree(k) < 3) then
good = .False.
exit
endif
enddo
if (good) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
N_tq += 1
do k=1,N_int
tq(k,1,N_tq) = det_buffer(k,1,i)
tq(k,2,N_tq) = det_buffer(k,2,i)
enddo
endif
endif
enddo i_loop
end
subroutine filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,m
logical :: is_in_wavefunction
integer :: degree(psi_det_size)
integer :: idx(0:psi_det_size)
logical :: good
integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
integer, intent(out) :: N_tq
integer :: nt,ni
logical, external :: is_connected_to, is_generable
integer(bit_kind),intent(in) :: microlist(Nint,2,*)
integer,intent(in) :: ptr_microlist(0:*)
integer,intent(in) :: N_microlist(0:*)
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
integer :: mobiles(2), smallerlist
N_tq = 0
i_loop : do i=1,N_selected
call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint)
if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
smallerlist = mobiles(1)
else
smallerlist = mobiles(2)
end if
if(N_microlist(smallerlist) > 0) then
do k=ptr_microlist(smallerlist), ptr_microlist(smallerlist)+N_microlist(smallerlist)-1
if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop
end do
end if
if(N_microlist(0) > 0) then
do k=1, N_microlist(0)
if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop
end do
end if
! Select determinants that are triple or quadruple excitations
! from the ref
good = .True.
call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
!good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
do k=1,idx(0)
if (degree(k) < 3) then
good = .False.
exit
endif
enddo
if (good) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
N_tq += 1
do k=1,N_int
tq(k,1,N_tq) = det_buffer(k,1,i)
tq(k,2,N_tq) = det_buffer(k,2,i)
enddo
endif
endif
enddo i_loop
end

2
plugins/mrcepa0/mrsc2sub.irp.f → plugins/mrcepa0/mrcc.irp.f
View File

@ -3,7 +3,7 @@ program mrsc2sub
double precision, allocatable :: energy(:)
allocate (energy(N_states))
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
mrmode = 3
read_wf = .True.

2
plugins/mrcepa0/mrcepa0.irp.f
View File

@ -3,7 +3,7 @@ program mrcepa0
double precision, allocatable :: energy(:)
allocate (energy(N_states))
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
mrmode = 1
read_wf = .True.

1
plugins/mrcepa0/mrcepa0_general.irp.f
View File

@ -15,6 +15,7 @@ subroutine run(N_st,energy)
integer :: n_it_mrcc_max
double precision :: thresh_mrcc
thresh_mrcc = 1d-7
n_it_mrcc_max = 10

3
plugins/mrcepa0/mrsc2.irp.f
View File

@ -3,9 +3,8 @@ program mrsc2
double precision, allocatable :: energy(:)
allocate (energy(N_states))
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
mrmode = 2
read_wf = .True.
SOFT_TOUCH read_wf
call print_cas_coefs

25
src/Determinants/connected_to_ref.irp.f
View File

@ -165,7 +165,7 @@ logical function is_connected_to(key,keys,Nint,Ndet)
integer :: i, l
integer :: degree_x2
logical, external :: is_generable_cassd
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
@ -183,12 +183,35 @@ logical function is_connected_to(key,keys,Nint,Ndet)
if (degree_x2 > 4) then
cycle
else
! if(.not. is_generable_cassd(keys(1,1,i), key(1,1), Nint)) cycle !!!Nint==1 !!!!!
is_connected_to = .true.
return
endif
enddo
end
logical function is_generable_cassd(det1, det2, Nint) !!! TEST Cl HARD !!!!!
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2)
integer :: degree, f, exc(0:2, 2, 2), h1, h2, p1, p2, s1, s2, t
double precision :: phase
is_generable_cassd = .false.
call get_excitation(det1, det2, exc, degree, phase, Nint)
if(degree == -1) return
if(degree == 0) then
is_generable_cassd = .true.
return
end if
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
if(degree == 1 .and. h1 <= 11) is_generable_cassd = .true.
if(degree == 2 .and. h1 <= 11 .and. h2 <= 11) is_generable_cassd = .true.
end function
logical function is_connected_to_by_mono(key,keys,Nint,Ndet)
use bitmasks
implicit none

41
src/Determinants/determinants.irp.f
View File

@ -664,3 +664,44 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde
end
logical function detEq(a,b,Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2)
integer :: ni, i
detEq = .false.
do i=1,2
do ni=1,Nint
if(a(ni,i) /= b(ni,i)) return
end do
end do
detEq = .true.
end function
integer function detCmp(a,b,Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2)
integer :: ni, i
detCmp = 0
do i=1,2
do ni=Nint,1,-1
if(a(ni,i) < b(ni,i)) then
detCmp = -1
return
else if(a(ni,i) > b(ni,i)) then
detCmp = 1
return
end if
end do
end do
end function

10
src/Determinants/slater_rules.irp.f
View File

@ -914,7 +914,6 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis
fullMatch = .false.
N_miniList = 0
N_subList = 0
l = popcnt(key_mask(1,1)) + popcnt(key_mask(1,2))
do ni = 2,Nint
l = l + popcnt(key_mask(ni,1)) + popcnt(key_mask(ni,2))
@ -947,8 +946,13 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis
miniList(ni,2,N_minilist) = fullList(ni,2,i)
enddo
else if(k == 0) then
fullMatch = .true.
return
N_minilist += 1
do ni=1,Nint
miniList(ni,1,N_minilist) = fullList(ni,1,i)
miniList(ni,2,N_minilist) = fullList(ni,2,i)
enddo
! fullMatch = .true.
! return
end if
end do
end if

20
src/Determinants/spindeterminants.irp.f
View File

@ -10,7 +10,7 @@ integer*8 function spin_det_search_key(det,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Return an integer*8 corresponding to a determinant index for searching
! Return an integer(8) corresponding to a determinant index for searching
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: det(Nint)
@ -64,9 +64,9 @@ BEGIN_TEMPLATE
integer :: i,j,k
integer, allocatable :: iorder(:)
integer*8, allocatable :: bit_tmp(:)
integer*8 :: last_key
integer*8, external :: spin_det_search_key
integer(8), allocatable :: bit_tmp(:)
integer(8) :: last_key
integer(8), external :: spin_det_search_key
logical,allocatable :: duplicate(:)
allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) )
@ -149,8 +149,8 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint)
integer(bit_kind), intent(in) :: key(Nint)
integer :: i, ibegin, iend, istep, l
integer*8 :: det_ref, det_search
integer*8, external :: spin_det_search_key
integer(8) :: det_ref, det_search
integer(8), external :: spin_det_search_key
logical :: in_wavefunction
in_wavefunction = .False.
@ -231,8 +231,8 @@ integer function get_index_in_psi_det_beta_unique(key,Nint)
integer(bit_kind), intent(in) :: key(Nint)
integer :: i, ibegin, iend, istep, l
integer*8 :: det_ref, det_search
integer*8, external :: spin_det_search_key
integer(8) :: det_ref, det_search
integer(8), external :: spin_det_search_key
logical :: in_wavefunction
in_wavefunction = .False.
@ -305,10 +305,10 @@ end
subroutine write_spindeterminants
use bitmasks
implicit none
integer*8, allocatable :: tmpdet(:,:)
integer(8), allocatable :: tmpdet(:,:)
integer :: N_int2
integer :: i,j,k
integer*8 :: det_8(100)
integer(8) :: det_8(100)
integer(bit_kind) :: det_bk((100*8)/bit_kind)
equivalence (det_8, det_bk)