2016-03-04 16:52:46 +01:00
|
|
|
use bitmasks
|
|
|
|
|
|
|
|
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-01 12:00:03 +02:00
|
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ]
|
2016-03-11 19:35:57 +01:00
|
|
|
use bitmasks
|
|
|
|
implicit none
|
2016-04-01 12:00:03 +02:00
|
|
|
integer :: i, j, i_state
|
|
|
|
|
|
|
|
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-01 12:00:03 +02:00
|
|
|
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
|
2016-04-29 16:23:20 +02:00
|
|
|
! do i = 1, N_det_ref
|
|
|
|
! delta_ii(i_state,i)= delta_ii_old(i,i_state)
|
|
|
|
! do j = 1, N_det_non_ref
|
|
|
|
! delta_ij(i_state,j,i) = delta_ij_old(i,j,i_state)
|
|
|
|
! end do
|
|
|
|
! end do
|
2016-04-01 12:00:03 +02:00
|
|
|
do i = 1, N_det_ref
|
2016-04-29 16:23:20 +02:00
|
|
|
delta_ii(i_state,i)= delta_ii_old(i_state,i)
|
2016-04-01 12:00:03 +02:00
|
|
|
do j = 1, N_det_non_ref
|
2016-04-29 16:23:20 +02:00
|
|
|
delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i)
|
2016-04-01 12:00:03 +02:00
|
|
|
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)
|
|
|
|
do j = 1, N_det_non_ref
|
|
|
|
delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state)
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
else
|
|
|
|
stop "invalid mrmode"
|
|
|
|
end if
|
|
|
|
end do
|
2016-04-29 16:23:20 +02:00
|
|
|
! do i=1,N_det_ref
|
|
|
|
! print *, delta_ii(1,i)
|
|
|
|
! end do
|
|
|
|
! do i=1,min(N_det_non_ref,100)
|
|
|
|
! print *, delta_ij(1,i,:)
|
|
|
|
! end do
|
2016-04-26 17:27:21 +02:00
|
|
|
! stop
|
|
|
|
|
|
|
|
|
2016-03-11 19:35:57 +01:00
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
|
|
|
|
2016-04-01 12:00:03 +02:00
|
|
|
|
2016-03-04 16:52:46 +01:00
|
|
|
BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ]
|
|
|
|
&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ]
|
2016-04-04 15:51:32 +02:00
|
|
|
&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_int,2,N_det_non_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_int,2,N_det_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ integer(bit_kind), active_sorb, (N_int,2) ]
|
2016-04-15 15:16:46 +02:00
|
|
|
&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0, (N_int,2,N_det_non_ref) ]
|
2016-04-29 16:23:20 +02:00
|
|
|
&BEGIN_PROVIDER [ integer, nlink, (N_det_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ integer, linked, (N_det_non_ref,N_det_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ integer, blokMwen, (N_det_non_ref,N_det_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ double precision, searchance, (N_det_ref) ]
|
|
|
|
&BEGIN_PROVIDER [ integer, child_num, (N_det_non_ref,N_det_ref) ]
|
|
|
|
|
2016-03-04 16:52:46 +01:00
|
|
|
use bitmasks
|
|
|
|
implicit none
|
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), nonactive_sorb(N_int,2), det(N_int, 2)
|
2016-04-29 16:23:20 +02:00
|
|
|
integer i, II, j, k, n, ni, idx(N_det_non_ref), shortcut(0:N_det_non_ref+1), blok, degree
|
2016-03-04 16:52:46 +01:00
|
|
|
logical, external :: detEq
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
active_sorb(:,:) = 0_8
|
|
|
|
nonactive_sorb(:,:) = not(0_8)
|
2016-03-11 19:35:57 +01:00
|
|
|
|
|
|
|
if(N_det_ref > 1) then
|
|
|
|
do i=1, N_det_ref
|
2016-04-04 15:51:32 +02:00
|
|
|
do k=1, N_int
|
|
|
|
active_sorb(k,1) = ior(psi_ref(k,1,i), active_sorb(k,1))
|
|
|
|
active_sorb(k,2) = ior(psi_ref(k,2,i), active_sorb(k,2))
|
|
|
|
nonactive_sorb(k,1) = iand(psi_ref(k,1,i), nonactive_sorb(k,1))
|
|
|
|
nonactive_sorb(k,2) = iand(psi_ref(k,2,i), nonactive_sorb(k,2))
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
do k=1, N_int
|
|
|
|
active_sorb(k,1) = iand(active_sorb(k,1), not(nonactive_sorb(k,1)))
|
|
|
|
active_sorb(k,2) = iand(active_sorb(k,2), not(nonactive_sorb(k,2)))
|
2016-03-11 19:35:57 +01:00
|
|
|
end do
|
|
|
|
end if
|
|
|
|
|
2016-03-04 16:52:46 +01:00
|
|
|
do i=1, N_det_non_ref
|
2016-04-04 15:51:32 +02:00
|
|
|
do k=1, N_int
|
|
|
|
det_noactive(k,1,i) = iand(psi_non_ref(k,1,i), not(active_sorb(k,1)))
|
|
|
|
det_noactive(k,2,i) = iand(psi_non_ref(k,2,i), not(active_sorb(k,2)))
|
|
|
|
end do
|
2016-03-04 16:52:46 +01:00
|
|
|
end do
|
|
|
|
|
|
|
|
call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int)
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
do i=1,N_det_non_ref
|
|
|
|
det_cepa0(:,:,i) = psi_non_ref(:,:,det_cepa0_idx(i))
|
2016-03-04 16:52:46 +01:00
|
|
|
end do
|
2016-04-26 17:27:21 +02:00
|
|
|
|
2016-03-04 16:52:46 +01:00
|
|
|
cepa0_shortcut(0) = 1
|
|
|
|
cepa0_shortcut(1) = 1
|
|
|
|
do i=2,N_det_non_ref
|
|
|
|
if(.not. detEq(det_noactive(1,1,i), det_noactive(1,1,i-1), N_int)) then
|
|
|
|
cepa0_shortcut(0) += 1
|
|
|
|
cepa0_shortcut(cepa0_shortcut(0)) = i
|
|
|
|
end if
|
|
|
|
end do
|
2016-03-11 19:35:57 +01:00
|
|
|
cepa0_shortcut(cepa0_shortcut(0)+1) = N_det_non_ref+1
|
2016-04-15 15:16:46 +02:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
if(.true.) then
|
|
|
|
do i=1,cepa0_shortcut(0)
|
|
|
|
n = cepa0_shortcut(i+1) - cepa0_shortcut(i)
|
|
|
|
call sort_dets_ab(det_cepa0(1,1,cepa0_shortcut(i)), idx, shortcut, n, N_int)
|
|
|
|
do k=1,n
|
|
|
|
idx(k) = det_cepa0_idx(cepa0_shortcut(i)-1+idx(k))
|
|
|
|
end do
|
|
|
|
det_cepa0_idx(cepa0_shortcut(i):cepa0_shortcut(i)+n-1) = idx(:n)
|
|
|
|
end do
|
|
|
|
end if
|
|
|
|
|
|
|
|
|
|
|
|
do i=1,N_det_ref
|
|
|
|
do k=1, N_int
|
|
|
|
det_ref_active(k,1,i) = iand(psi_ref(k,1,i), active_sorb(k,1))
|
|
|
|
det_ref_active(k,2,i) = iand(psi_ref(k,2,i), active_sorb(k,2))
|
|
|
|
end do
|
2016-04-15 15:16:46 +02:00
|
|
|
end do
|
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
do i=1,N_det_non_ref
|
|
|
|
do k=1, N_int
|
|
|
|
det_cepa0_active(k,1,i) = iand(psi_non_ref(k,1,det_cepa0_idx(i)), active_sorb(k,1))
|
|
|
|
det_cepa0_active(k,2,i) = iand(psi_non_ref(k,2,det_cepa0_idx(i)), active_sorb(k,2))
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
do i=1,N_det_non_ref
|
|
|
|
if(.not. detEq(psi_non_ref(1,1,det_cepa0_idx(i)), det_cepa0(1,1,i),N_int)) stop "STOOOP"
|
|
|
|
end do
|
2016-04-29 16:23:20 +02:00
|
|
|
|
|
|
|
searchance = 0d0
|
|
|
|
child_num = 0
|
|
|
|
do J = 1, N_det_ref
|
|
|
|
nlink(J) = 0
|
|
|
|
do blok=1,cepa0_shortcut(0)
|
|
|
|
do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
|
|
|
|
call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
|
|
|
|
if(degree <= 2) then
|
|
|
|
nlink(J) += 1
|
|
|
|
linked(nlink(J),J) = k
|
|
|
|
child_num(k, J) = nlink(J)
|
|
|
|
blokMwen(nlink(J),J) = blok
|
|
|
|
searchance(J) += 1d0 + log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
print *, "pre done"
|
2016-03-04 16:52:46 +01:00
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
BEGIN_PROVIDER [ double precision, delta_cas_old, (N_det_ref, N_det_ref, N_states) ]
|
2016-03-04 16:52:46 +01:00
|
|
|
use bitmasks
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k
|
2016-04-01 12:00:03 +02:00
|
|
|
double precision :: Hjk, Hki, Hij
|
|
|
|
integer i_state, degree
|
2016-03-04 16:52:46 +01:00
|
|
|
|
2016-03-11 19:35:57 +01:00
|
|
|
provide lambda_mrcc
|
2016-04-04 15:51:32 +02:00
|
|
|
do i_state = 1, N_states
|
2016-04-08 13:25:55 +02:00
|
|
|
!$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Hjk,Hki,degree) shared(no_mono_dressing,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref)
|
2016-04-04 15:51:32 +02:00
|
|
|
do i=1,N_det_ref
|
|
|
|
do j=1,i
|
2016-04-15 09:16:31 +02:00
|
|
|
!call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int)
|
2016-04-04 15:51:32 +02:00
|
|
|
delta_cas(i,j,i_state) = 0d0
|
2016-04-11 17:42:15 +02:00
|
|
|
!if(no_mono_dressing .and. degree == 1) cycle
|
2016-04-04 15:51:32 +02:00
|
|
|
do k=1,N_det_non_ref
|
|
|
|
|
|
|
|
call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
|
|
|
|
call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki)
|
|
|
|
|
|
|
|
delta_cas(i,j,i_state) += Hjk * Hki * lambda_mrcc(i_state, k)
|
|
|
|
end do
|
2016-03-04 16:52:46 +01:00
|
|
|
end do
|
|
|
|
end do
|
2016-04-04 15:51:32 +02:00
|
|
|
!$OMP END PARALLEL DO
|
2016-04-10 12:30:22 +02:00
|
|
|
do i=1,N_det_ref
|
|
|
|
do j=1,i
|
|
|
|
delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
|
|
|
|
end do
|
|
|
|
end do
|
2016-03-04 16:52:46 +01:00
|
|
|
end do
|
|
|
|
END_PROVIDER
|
|
|
|
|
2016-04-01 12:00:03 +02:00
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ]
|
|
|
|
use bitmasks
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k
|
|
|
|
double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall
|
2016-04-15 15:16:46 +02:00
|
|
|
integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref)
|
2016-04-15 09:16:31 +02:00
|
|
|
|
|
|
|
provide lambda_mrcc
|
2016-04-15 15:16:46 +02:00
|
|
|
npres = 0
|
2016-04-15 09:16:31 +02:00
|
|
|
delta_cas = 0d0
|
|
|
|
call wall_time(wall)
|
2016-04-15 15:16:46 +02:00
|
|
|
print *, "dcas ", wall
|
2016-04-15 09:16:31 +02:00
|
|
|
do i_state = 1, N_states
|
2016-04-29 16:23:20 +02:00
|
|
|
!!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref)
|
2016-04-15 09:16:31 +02:00
|
|
|
do k=1,N_det_non_ref
|
|
|
|
if(lambda_mrcc(i_state, k) == 0d0) cycle
|
|
|
|
npre = 0
|
|
|
|
do i=1,N_det_ref
|
|
|
|
call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki)
|
|
|
|
if(Hki /= 0d0) then
|
2016-04-29 16:23:20 +02:00
|
|
|
!!$OMP ATOMIC
|
2016-04-15 15:16:46 +02:00
|
|
|
npres(i) += 1
|
2016-04-15 09:16:31 +02:00
|
|
|
npre += 1
|
|
|
|
ipre(npre) = i
|
|
|
|
pre(npre) = Hki
|
|
|
|
end if
|
|
|
|
end do
|
2016-04-15 15:16:46 +02:00
|
|
|
|
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
do i=1,npre
|
|
|
|
do j=1,i
|
2016-04-29 16:23:20 +02:00
|
|
|
!!$OMP ATOMIC
|
2016-04-15 09:16:31 +02:00
|
|
|
delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k)
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
2016-04-29 16:23:20 +02:00
|
|
|
!!$OMP END PARALLEL DO
|
2016-04-15 15:16:46 +02:00
|
|
|
print *, npres
|
|
|
|
npre=0
|
|
|
|
do i=1,N_det_ref
|
|
|
|
npre += npres(i)
|
|
|
|
end do
|
|
|
|
print *, npre
|
2016-04-26 17:27:21 +02:00
|
|
|
!stop
|
2016-04-15 09:16:31 +02:00
|
|
|
do i=1,N_det_ref
|
|
|
|
do j=1,i
|
|
|
|
delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
call wall_time(wall)
|
2016-04-15 15:16:46 +02:00
|
|
|
print *, "dcas", wall
|
2016-04-15 09:16:31 +02:00
|
|
|
! stop
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
2016-03-04 16:52:46 +01:00
|
|
|
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
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
logical function isInCassd(a,Nint)
|
|
|
|
use bitmasks
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: Nint
|
|
|
|
integer(bit_kind), intent(in) :: a(Nint,2)
|
|
|
|
integer(bit_kind) :: inac, virt
|
|
|
|
integer :: ni, i, deg
|
|
|
|
|
|
|
|
|
|
|
|
isInCassd = .false.
|
|
|
|
|
|
|
|
|
|
|
|
deg = 0
|
|
|
|
do i=1,2
|
|
|
|
do ni=1,Nint
|
|
|
|
virt = iand(not(HF_bitmask(ni,i)), not(active_sorb(ni,i)))
|
|
|
|
deg += popcnt(iand(virt, a(ni,i)))
|
|
|
|
if(deg > 2) return
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
deg = 0
|
|
|
|
do i=1,2
|
|
|
|
do ni=1,Nint
|
|
|
|
inac = iand(HF_bitmask(ni,i), not(active_sorb(ni,i)))
|
|
|
|
deg += popcnt(xor(iand(inac,a(ni,i)), inac))
|
|
|
|
if(deg > 2) return
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
|
|
|
|
isInCassd = .true.
|
|
|
|
end function
|
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
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
|
2016-04-26 17:27:21 +02:00
|
|
|
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
|
|
|
|
|
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
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
|
2016-03-04 16:52:46 +01:00
|
|
|
|
|
|
|
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-01 12:00:03 +02:00
|
|
|
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) ]
|
|
|
|
use bitmasks
|
|
|
|
implicit none
|
|
|
|
|
|
|
|
integer :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni
|
2016-04-15 09:16:31 +02:00
|
|
|
integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref)
|
2016-04-01 12:00:03 +02:00
|
|
|
logical :: ok
|
|
|
|
double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
|
2016-04-15 15:16:46 +02:00
|
|
|
double precision :: contrib, HIIi, HJk, wall
|
2016-04-01 12:00:03 +02:00
|
|
|
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
|
2016-04-15 09:16:31 +02:00
|
|
|
integer(bit_kind) :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2), sortRef(N_int,2,N_det_ref)
|
2016-04-01 12:00:03 +02:00
|
|
|
integer, allocatable :: idx_sorted_bit(:)
|
2016-04-15 09:16:31 +02:00
|
|
|
integer, external :: get_index_in_psi_det_sorted_bit, searchDet
|
|
|
|
logical, external :: is_in_wavefunction, detEq
|
2016-04-01 12:00:03 +02:00
|
|
|
|
|
|
|
integer :: II, blok
|
2016-04-26 17:27:21 +02:00
|
|
|
integer*8, save :: notf = 0
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
call wall_time(wall)
|
|
|
|
print *, "cepa0", wall
|
2016-04-01 12:00:03 +02:00
|
|
|
provide det_cepa0_active delta_cas lambda_mrcc
|
|
|
|
provide mo_bielec_integrals_in_map
|
|
|
|
allocate(idx_sorted_bit(N_det))
|
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
sortRef(:,:,:) = det_ref_active(:,:,:)
|
2016-04-15 09:16:31 +02:00
|
|
|
call sort_det(sortRef, sortRefIdx, N_det_ref, N_int)
|
|
|
|
|
2016-04-01 12:00:03 +02:00
|
|
|
idx_sorted_bit(:) = -1
|
|
|
|
do i=1,N_det_non_ref
|
|
|
|
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
|
|
|
|
enddo
|
2016-04-04 15:51:32 +02:00
|
|
|
|
|
|
|
|
|
|
|
do i_state = 1, N_states
|
|
|
|
delta_mrcepa0_ii(:,:) = 0d0
|
|
|
|
delta_mrcepa0_ij(:,:,:) = 0d0
|
|
|
|
|
|
|
|
!$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii) &
|
2016-04-15 09:16:31 +02:00
|
|
|
!$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib) &
|
2016-04-04 15:51:32 +02:00
|
|
|
!$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) &
|
|
|
|
!$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas) &
|
2016-04-26 17:27:21 +02:00
|
|
|
!$OMP shared(notf,i_state, sortRef, sortRefIdx)
|
2016-04-04 15:51:32 +02:00
|
|
|
do blok=1,cepa0_shortcut(0)
|
|
|
|
do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
|
|
|
|
do II=1,N_det_ref
|
|
|
|
call get_excitation_degree(psi_ref(1,1,II),psi_non_ref(1,1,det_cepa0_idx(i)),degree,N_int)
|
|
|
|
if (degree > 2 ) cycle
|
|
|
|
|
|
|
|
do ni=1,N_int
|
|
|
|
made_hole(ni,1) = iand(det_ref_active(ni,1,II), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II)))
|
|
|
|
made_hole(ni,2) = iand(det_ref_active(ni,2,II), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II)))
|
2016-04-15 09:16:31 +02:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
made_particle(ni,1) = iand(det_cepa0_active(ni,1,i), xor(det_cepa0_active(ni,1,i), det_ref_active(ni,1,II)))
|
|
|
|
made_particle(ni,2) = iand(det_cepa0_active(ni,2,i), xor(det_cepa0_active(ni,2,i), det_ref_active(ni,2,II)))
|
|
|
|
end do
|
2016-04-01 12:00:03 +02:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
kloop: do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 !i
|
2016-04-04 15:51:32 +02:00
|
|
|
if(lambda_mrcc(i_state, det_cepa0_idx(k)) == 0d0) cycle
|
2016-04-01 12:00:03 +02:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
do ni=1,N_int
|
|
|
|
if(iand(made_hole(ni,1), det_cepa0_active(ni,1,k)) /= 0) cycle kloop
|
|
|
|
if(iand(made_particle(ni,1), det_cepa0_active(ni,1,k)) /= made_particle(ni,1)) cycle kloop
|
|
|
|
if(iand(made_hole(ni,2), det_cepa0_active(ni,2,k)) /= 0) cycle kloop
|
|
|
|
if(iand(made_particle(ni,2), det_cepa0_active(ni,2,k)) /= made_particle(ni,2)) cycle kloop
|
|
|
|
end do
|
|
|
|
do ni=1,N_int
|
|
|
|
myActive(ni,1) = xor(det_cepa0_active(ni,1,k), made_hole(ni,1))
|
|
|
|
myActive(ni,1) = xor(myActive(ni,1), made_particle(ni,1))
|
|
|
|
myActive(ni,2) = xor(det_cepa0_active(ni,2,k), made_hole(ni,2))
|
|
|
|
myActive(ni,2) = xor(myActive(ni,2), made_particle(ni,2))
|
|
|
|
end do
|
2016-04-01 12:00:03 +02:00
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
j = searchDet(sortRef, myActive, N_det_ref, N_int)
|
2016-04-26 17:27:21 +02:00
|
|
|
if(j == -1) then
|
|
|
|
cycle
|
|
|
|
end if
|
2016-04-15 09:16:31 +02:00
|
|
|
j = sortRefIdx(j)
|
2016-04-26 17:27:21 +02:00
|
|
|
!$OMP ATOMIC
|
|
|
|
notf = notf+1
|
|
|
|
!if(i/=k .and. dabs(psi_non_ref_coef(det_cepa0_idx(i),i_state)) < dabs(psi_non_ref_coef(det_cepa0_idx(k),i_state))) cycle
|
|
|
|
! if(dabs(lambda_mrcc(i_state,det_cepa0_idx(i))) > dabs(lambda_mrcc(i_state,det_cepa0_idx(k)))) cycle
|
|
|
|
! if(dabs(lambda_mrcc(i_state,det_cepa0_idx(i))) == dabs(lambda_mrcc(i_state,det_cepa0_idx(k))) .and. i < k) cycle
|
|
|
|
!if(.not. j==II .and. dabs(psi_ref_coef(II,i_state)) < dabs(psi_ref_coef(j,i_state))) cycle
|
|
|
|
|
2016-04-15 09:16:31 +02:00
|
|
|
call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk)
|
|
|
|
contrib = delta_cas(II, J, i_state) * HJk * lambda_mrcc(i_state, det_cepa0_idx(k))
|
|
|
|
!$OMP ATOMIC
|
|
|
|
delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib
|
|
|
|
|
|
|
|
if(dabs(psi_ref_coef(J,i_state)).ge.5.d-5) then
|
2016-04-04 15:51:32 +02:00
|
|
|
!$OMP ATOMIC
|
2016-04-26 17:27:21 +02:00
|
|
|
delta_mrcepa0_ii(J,i_state) -= contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state)
|
2016-04-15 09:16:31 +02:00
|
|
|
end if
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
end do kloop
|
2016-04-01 12:00:03 +02:00
|
|
|
end do
|
|
|
|
end do
|
2016-04-04 15:51:32 +02:00
|
|
|
end do
|
|
|
|
!$OMP END PARALLEL DO
|
2016-04-01 12:00:03 +02:00
|
|
|
end do
|
|
|
|
deallocate(idx_sorted_bit)
|
2016-04-15 15:16:46 +02:00
|
|
|
call wall_time(wall)
|
2016-04-26 17:27:21 +02:00
|
|
|
print *, "cepa0", wall, notf
|
|
|
|
!stop
|
2016-04-01 12:00:03 +02:00
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
2016-03-11 19:35:57 +01:00
|
|
|
BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ]
|
|
|
|
&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ]
|
2016-03-04 16:52:46 +01:00
|
|
|
use bitmasks
|
|
|
|
implicit none
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
integer :: i_state, i, i_I, J, k, degree, degree2, l, deg, ni
|
2016-03-11 19:35:57 +01:00
|
|
|
integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_
|
|
|
|
logical :: ok
|
|
|
|
double precision :: phase_Ji, phase_Ik, phase_Ii
|
|
|
|
double precision :: contrib, delta_IJk, HJk, HIk, HIl
|
|
|
|
integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii
|
|
|
|
integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2)
|
|
|
|
integer, allocatable :: idx_sorted_bit(:)
|
|
|
|
integer, external :: get_index_in_psi_det_sorted_bit
|
|
|
|
|
|
|
|
integer :: II, blok
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
provide delta_cas lambda_mrcc
|
2016-03-11 19:35:57 +01:00
|
|
|
allocate(idx_sorted_bit(N_det))
|
|
|
|
idx_sorted_bit(:) = -1
|
|
|
|
do i=1,N_det_non_ref
|
|
|
|
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
|
|
|
|
enddo
|
2016-04-04 15:51:32 +02:00
|
|
|
|
|
|
|
do i_state = 1, N_states
|
|
|
|
delta_sub_ij(:,:,:) = 0d0
|
|
|
|
delta_sub_ii(:,:) = 0d0
|
|
|
|
|
|
|
|
provide mo_bielec_integrals_in_map
|
|
|
|
|
|
|
|
|
|
|
|
!$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii) &
|
|
|
|
!$OMP private(i, J, k, degree, degree2, l, deg, ni) &
|
|
|
|
!$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_) &
|
|
|
|
!$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) &
|
|
|
|
!$OMP private(det_tmp, det_tmp2, II, blok) &
|
|
|
|
!$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) &
|
|
|
|
!$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb)
|
|
|
|
do i=1,N_det_non_ref
|
|
|
|
if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref
|
|
|
|
do J=1,N_det_ref
|
|
|
|
call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_Ji,degree,phase_Ji,N_int)
|
|
|
|
if(degree == -1) cycle
|
2016-03-11 19:35:57 +01:00
|
|
|
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
do II=1,N_det_ref
|
|
|
|
call apply_excitation(psi_ref(1,1,II),exc_Ji,det_tmp,ok,N_int)
|
2016-04-15 09:16:31 +02:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
if(.not. ok) cycle
|
|
|
|
l = get_index_in_psi_det_sorted_bit(det_tmp, N_int)
|
|
|
|
if(l == 0) cycle
|
|
|
|
l = idx_sorted_bit(l)
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
call i_h_j(psi_ref(1,1,II), det_tmp, N_int, HIl)
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
do k=1,N_det_non_ref
|
|
|
|
if(lambda_mrcc(i_state, k) == 0d0) cycle
|
|
|
|
call get_excitation(psi_ref(1,1,II),psi_non_ref(1,1,k),exc_Ik,degree2,phase_Ik,N_int)
|
|
|
|
|
|
|
|
det_tmp(:,:) = 0_bit_kind
|
|
|
|
det_tmp2(:,:) = 0_bit_kind
|
|
|
|
|
|
|
|
ok = .true.
|
|
|
|
do ni=1,N_int
|
|
|
|
det_tmp(ni,1) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,k)), not(active_sorb(ni,1)))
|
|
|
|
det_tmp(ni,2) = iand(xor(HF_bitmask(ni,1), psi_non_ref(ni,1,i)), not(active_sorb(ni,1)))
|
|
|
|
ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2))))
|
|
|
|
|
|
|
|
det_tmp(ni,1) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,k)), not(active_sorb(ni,2)))
|
|
|
|
det_tmp(ni,2) = iand(xor(HF_bitmask(ni,2), psi_non_ref(ni,2,i)), not(active_sorb(ni,2)))
|
|
|
|
ok = ok .and. (popcnt(det_tmp(ni,1)) + popcnt(det_tmp(ni,2)) == popcnt(xor(det_tmp(ni,1), det_tmp(ni,2))))
|
|
|
|
end do
|
|
|
|
|
|
|
|
if(ok) cycle
|
|
|
|
|
|
|
|
|
|
|
|
call i_h_j(psi_ref(1,1,J), psi_non_ref(1,1,k), N_int, HJk)
|
|
|
|
call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,k), N_int, HIk)
|
|
|
|
if(HJk == 0) cycle
|
|
|
|
!assert HIk == 0
|
|
|
|
delta_IJk = HJk * HIk * lambda_mrcc(i_state, k)
|
|
|
|
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
|
|
|
|
if(ok) cycle
|
|
|
|
contrib = delta_IJk * HIl * lambda_mrcc(i_state,l)
|
|
|
|
!$OMP ATOMIC
|
|
|
|
delta_sub_ij(II, i, i_state) += contrib
|
|
|
|
if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
|
|
|
|
!$OMP ATOMIC
|
|
|
|
delta_sub_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state)
|
|
|
|
endif
|
|
|
|
end do
|
2016-03-11 19:35:57 +01:00
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
2016-04-04 15:51:32 +02:00
|
|
|
!$OMP END PARALLEL DO
|
2016-03-11 19:35:57 +01:00
|
|
|
end do
|
|
|
|
deallocate(idx_sorted_bit)
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
|
|
|
subroutine set_det_bit(det, p, s)
|
|
|
|
implicit none
|
|
|
|
integer(bit_kind),intent(inout) :: det(N_int, 2)
|
|
|
|
integer, intent(in) :: p, s
|
|
|
|
integer :: ni, pos
|
|
|
|
|
|
|
|
ni = (p-1)/bit_kind_size + 1
|
|
|
|
pos = mod(p-1, bit_kind_size)
|
|
|
|
det(ni,s) = ibset(det(ni,s), pos)
|
2016-03-04 16:52:46 +01:00
|
|
|
end subroutine
|
2016-03-11 19:35:57 +01:00
|
|
|
|
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ]
|
|
|
|
integer :: i,j
|
|
|
|
do i=1,N_det_ref
|
|
|
|
do j=1,N_det_non_ref
|
|
|
|
call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,j), N_int, h_(i,j))
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
END_PROVIDER
|
|
|
|
|
|
|
|
|
2016-04-28 09:55:07 +02:00
|
|
|
|
|
|
|
BEGIN_PROVIDER [ double precision, delta_ij_older, (N_det_ref,N_det_non_ref,N_states) ]
|
|
|
|
&BEGIN_PROVIDER [ double precision, delta_ii_older, (N_det_ref,N_states) ]
|
2016-04-15 15:16:46 +02:00
|
|
|
implicit none
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
integer :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni
|
2016-04-26 17:27:21 +02:00
|
|
|
integer :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s
|
2016-04-29 16:23:20 +02:00
|
|
|
! integer, allocatable :: linked(:,:), blokMwen(:, :), nlink(:)
|
2016-03-11 19:35:57 +01:00
|
|
|
logical :: ok
|
2016-04-26 17:27:21 +02:00
|
|
|
double precision :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states)
|
2016-04-29 16:23:20 +02:00
|
|
|
double precision :: contrib, wall, iwall ! , searchance(N_det_ref)
|
2016-04-28 09:55:07 +02:00
|
|
|
double precision, allocatable :: deltaMwen(:,:,:), deltaIImwen(:,:)
|
2016-03-11 19:35:57 +01:00
|
|
|
integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
|
2016-04-15 09:16:31 +02:00
|
|
|
integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt
|
2016-04-26 17:27:21 +02:00
|
|
|
integer, allocatable :: idx_sorted_bit(:)
|
|
|
|
integer, external :: get_index_in_psi_det_sorted_bit, searchDet, detCmp
|
|
|
|
logical, external :: is_in_wavefunction, isInCassd, detEq
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-28 09:55:07 +02:00
|
|
|
! -459.6346665282306
|
|
|
|
! -459.6346665282306
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
call wall_time(iwall)
|
2016-03-11 19:35:57 +01:00
|
|
|
allocate(idx_sorted_bit(N_det))
|
2016-04-29 16:23:20 +02:00
|
|
|
! allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref))
|
2016-03-11 19:35:57 +01:00
|
|
|
|
|
|
|
idx_sorted_bit(:) = -1
|
|
|
|
do i=1,N_det_non_ref
|
|
|
|
idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
|
|
|
|
enddo
|
2016-03-04 16:52:46 +01:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
|
|
|
|
do i_state = 1, N_states
|
2016-03-11 19:35:57 +01:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
delta_ii_old(:,:) = 0d0
|
|
|
|
delta_ij_old(:,:,:) = 0d0
|
2016-04-28 09:55:07 +02:00
|
|
|
searchance = 0d0
|
|
|
|
!$OMP PARALLEL DO default(none) schedule(dynamic) private(blok,k,degree) shared(searchance,nlink,linked,blokMwen,psi_ref, det_cepa0,cepa0_shortcut, N_det_ref, N_int)
|
2016-04-15 15:16:46 +02:00
|
|
|
do J = 1, N_det_ref
|
|
|
|
nlink(J) = 0
|
|
|
|
do blok=1,cepa0_shortcut(0)
|
|
|
|
do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
|
|
|
|
call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
|
|
|
|
if(degree <= 2) then
|
|
|
|
nlink(J) += 1
|
|
|
|
linked(nlink(J),J) = k
|
|
|
|
blokMwen(nlink(J),J) = blok
|
2016-04-28 09:55:07 +02:00
|
|
|
searchance(J) += log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
|
2016-04-15 15:16:46 +02:00
|
|
|
end if
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
end do
|
|
|
|
!$OMP END PARALLEL DO
|
|
|
|
|
|
|
|
|
2016-04-28 09:55:07 +02:00
|
|
|
! do i=1,cepa0_shortcut(0)
|
|
|
|
! print *, cepa0_shortcut(i+1) - cepa0_shortcut(i)
|
|
|
|
! end do
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
!$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_ij_old, delta_ii_old) &
|
2016-04-26 17:27:21 +02:00
|
|
|
!$OMP private(m,kk, i_I, i, J, k, degree, degree2, l, deg, ni, inac, virt) &
|
|
|
|
!$OMP private(ok,p1,p2,h1,h2,s1,s2, blok, wall, I_s, J_s) &
|
|
|
|
!$OMP private(phase_iI, phase_Ik, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv, dia_hla) &
|
2016-04-04 15:51:32 +02:00
|
|
|
!$OMP private(contrib, exc_iI, exc_Ik, exc_IJ, det_tmp, det_tmp2) &
|
|
|
|
!$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) &
|
2016-04-15 15:16:46 +02:00
|
|
|
!$OMP shared(i_state, lambda_mrcc, hf_bitmask, active_sorb,cepa0_shortcut,det_cepa0) &
|
2016-04-28 09:55:07 +02:00
|
|
|
!$OMP shared(h_,det_cepa0_idx, linked, blokMwen, nlink, iwall, searchance)
|
2016-04-04 15:51:32 +02:00
|
|
|
do i = 1 , N_det_non_ref
|
2016-04-15 15:16:46 +02:00
|
|
|
if(mod(i,100) == 0) then
|
|
|
|
call wall_time(wall)
|
|
|
|
wall = wall-iwall
|
|
|
|
print *, i, "/", N_det_non_ref, wall * (dfloat(N_det_non_ref) / dfloat(i)), wall, wall * (dfloat(N_det_non_ref) / dfloat(i))-wall
|
|
|
|
end if
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
if(lambda_mrcc(i_state, i) == 0d0) cycle
|
2016-04-15 15:16:46 +02:00
|
|
|
|
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
do I_s = 1, N_det_ref
|
2016-04-28 09:55:07 +02:00
|
|
|
do J_s = 1, I_s
|
2016-04-15 15:16:46 +02:00
|
|
|
|
2016-04-28 09:55:07 +02:00
|
|
|
if(nlink(I_s) < nlink(J_s)) then
|
|
|
|
!if(searchance(I_s) < searchance(J_s)) then
|
2016-04-26 17:27:21 +02:00
|
|
|
i_I = I_s
|
|
|
|
J = J_s
|
|
|
|
else
|
|
|
|
i_I = J_s
|
|
|
|
J = I_s
|
|
|
|
end if
|
|
|
|
|
|
|
|
!call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,J),N_int,hJi)
|
|
|
|
!!!!
|
2016-04-28 09:55:07 +02:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
!!!!!!
|
|
|
|
hJi = h_(J,i)
|
2016-04-15 09:16:31 +02:00
|
|
|
if(hJi == 0) cycle
|
2016-04-26 17:27:21 +02:00
|
|
|
hIi = h_(i_I,i)
|
|
|
|
if(hIi == 0) cycle
|
|
|
|
|
|
|
|
diI = hIi * lambda_mrcc(i_state, i)
|
2016-04-04 15:51:32 +02:00
|
|
|
delta_JI = hJi * diI
|
2016-04-26 17:27:21 +02:00
|
|
|
|
|
|
|
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,i),exc_iI,degree2,phase_iI,N_int)
|
|
|
|
if(degree2 == -1) cycle
|
|
|
|
!call i_h_j(psi_non_ref(1,1,i), psi_ref(1,1,i_I),N_int,hIi)
|
2016-04-15 15:16:46 +02:00
|
|
|
|
|
|
|
do kk = 1 , nlink(i_I)
|
|
|
|
k = linked(kk,i_I)
|
|
|
|
blok = blokMwen(kk,i_I)
|
2016-04-04 15:51:32 +02:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
!if(lambda_mrcc(i_state, det_cepa0_idx(k)) == 0d0) cycle
|
2016-04-04 15:51:32 +02:00
|
|
|
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
call get_excitation(psi_ref(1,1,i_I),det_cepa0(1,1,k),exc_Ik,degree,phase_Ik,N_int)
|
|
|
|
!if(degree == -1) cycle
|
|
|
|
if(degree == -1) stop "STOP; ( linked )"
|
2016-04-15 09:16:31 +02:00
|
|
|
|
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
|
2016-04-15 09:16:31 +02:00
|
|
|
if(.not. ok) cycle
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int)
|
2016-04-04 15:51:32 +02:00
|
|
|
if(.not. ok) cycle
|
2016-04-15 15:16:46 +02:00
|
|
|
|
|
|
|
if(isInCassd(det_tmp, N_int)) cycle
|
|
|
|
|
2016-04-04 15:51:32 +02:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2, cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int)
|
|
|
|
if(l == -1) cycle
|
|
|
|
l = det_cepa0_idx(cepa0_shortcut(blok)-1+l)
|
|
|
|
!call i_h_j(det_cepa0(1,1,k), det_tmp, N_int, HiI)
|
|
|
|
!call i_h_j(psi_non_ref(1,1,l), det_tmp, N_int, HJi)
|
2016-04-04 15:51:32 +02:00
|
|
|
|
2016-04-26 17:27:21 +02:00
|
|
|
|
|
|
|
|
2016-04-28 09:55:07 +02:00
|
|
|
call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int)
|
|
|
|
call get_excitation(det_tmp,psi_non_ref(1,1,l),exc_IJ,degree2,phase_al,N_int)
|
|
|
|
delta_JI = hJi * hIi * lambda_mrcc(i_state, i) * phase_al * phase_Ji
|
2016-04-26 17:27:21 +02:00
|
|
|
ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
|
|
|
|
HkI = h_(i_I,det_cepa0_idx(k))
|
|
|
|
dkI(i_state) = HkI * lambda_mrcc(i_state, det_cepa0_idx(k))
|
|
|
|
contrib = dkI(i_state) * delta_JI
|
2016-04-28 09:55:07 +02:00
|
|
|
!$OMP ATOMIC
|
2016-04-26 17:27:21 +02:00
|
|
|
delta_ij_old(i_I,l,i_state) += contrib
|
|
|
|
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
|
2016-04-28 09:55:07 +02:00
|
|
|
!$OMP ATOMIC
|
2016-04-26 17:27:21 +02:00
|
|
|
delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state)
|
|
|
|
endif
|
|
|
|
!
|
2016-04-28 09:55:07 +02:00
|
|
|
if(l == det_cepa0_idx(k)) cycle
|
|
|
|
call get_excitation(psi_ref(1,1,I_i),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int)
|
|
|
|
call get_excitation(det_tmp,det_cepa0(1,1,k),exc_IJ,degree2,phase_al,N_int)
|
|
|
|
delta_JI = hJi * hIi * lambda_mrcc(i_state, i) * phase_al * phase_Ji
|
|
|
|
|
|
|
|
ci_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state)
|
|
|
|
HkI = h_(J,l)
|
|
|
|
dkI(i_state) = HkI * lambda_mrcc(i_state, l)
|
|
|
|
contrib = dkI(i_state) * delta_JI
|
|
|
|
!$OMP ATOMIC
|
|
|
|
delta_ij_old(J,det_cepa0_idx(k),i_state) += contrib
|
|
|
|
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
|
|
|
|
!$OMP ATOMIC
|
|
|
|
delta_ii_old(J,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
|
|
|
|
end if
|
2016-04-04 15:51:32 +02:00
|
|
|
|
|
|
|
enddo
|
2016-03-11 19:35:57 +01:00
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
enddo
|
2016-04-26 17:27:21 +02:00
|
|
|
!$OMP END PARALLEL DO
|
|
|
|
|
|
|
|
! do i=1,N_det_non_ref
|
|
|
|
! print *, delta_ij_old(:,i,i_state)
|
|
|
|
! end do
|
|
|
|
! stop
|
2016-04-04 15:51:32 +02:00
|
|
|
end do
|
2016-04-01 12:00:03 +02:00
|
|
|
deallocate(idx_sorted_bit)
|
2016-04-26 17:27:21 +02:00
|
|
|
|
|
|
|
|
2016-04-15 15:16:46 +02:00
|
|
|
! call wall_time(wall)
|
|
|
|
! print *, "old ", wall
|
2016-04-01 12:00:03 +02:00
|
|
|
END_PROVIDER
|
2016-03-04 16:52:46 +01:00
|
|
|
|
2016-04-28 09:55:07 +02:00
|
|
|
!
|
|
|
|
! BEGIN_PROVIDER [ double precision, delta_ij_old, (N_det_ref,N_det_non_ref,N_states) ]
|
|
|
|
! &BEGIN_PROVIDER [ double precision, delta_ii_old, (N_det_ref,N_states) ]
|
|
|
|
! implicit none
|
|
|
|
!
|
|
|
|
! integer :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2
|
|
|
|
! integer :: p1,p2,h1,h2,s1,s2, blok, I_s, J_s
|
2016-04-29 16:23:20 +02:00
|
|
|
! ! integer, allocatable :: linked(:,:), blokMwen(:, :), nlink(:)
|
2016-04-28 09:55:07 +02:00
|
|
|
! logical :: ok
|
|
|
|
! double precision :: phase_iI, phase_Ik, phase_Jl, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI(N_states), HkI, ci_inv(N_states), dia_hla(N_states)
|
2016-04-29 16:23:20 +02:00
|
|
|
! double precision :: contrib, wall, iwall !, searchance(N_det_ref)
|
2016-04-28 09:55:07 +02:00
|
|
|
! double precision, allocatable :: deltaMwen(:,:,:), deltaIImwen(:,:)
|
|
|
|
! integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ
|
|
|
|
! integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2), inac, virt
|
|
|
|
! integer, external :: get_index_in_psi_det_sorted_bit, searchDet, detCmp
|
|
|
|
! logical, external :: is_in_wavefunction, isInCassd, detEq
|
|
|
|
!
|
|
|
|
! ! -459.6346665282306
|
|
|
|
! ! -459.6346665282306
|
|
|
|
!
|
|
|
|
! call wall_time(iwall)
|
2016-04-29 16:23:20 +02:00
|
|
|
! !allocate(linked(N_det_non_ref, N_det_ref), blokMwen(N_det_non_ref, N_det_ref), nlink(N_det_ref))
|
2016-04-28 09:55:07 +02:00
|
|
|
!
|
|
|
|
!
|
|
|
|
! delta_ii_old(:,:) = 0d0
|
|
|
|
! delta_ij_old(:,:,:) = 0d0
|
|
|
|
!
|
2016-04-29 16:23:20 +02:00
|
|
|
! ! searchance = 0d0
|
|
|
|
! ! do J = 1, N_det_ref
|
|
|
|
! ! nlink(J) = 0
|
|
|
|
! ! do blok=1,cepa0_shortcut(0)
|
|
|
|
! ! do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
|
|
|
|
! ! call get_excitation_degree(psi_ref(1,1,J),det_cepa0(1,1,k),degree,N_int)
|
|
|
|
! ! if(degree <= 2) then
|
|
|
|
! ! nlink(J) += 1
|
|
|
|
! ! linked(nlink(J),J) = k
|
|
|
|
! ! blokMwen(nlink(J),J) = blok
|
|
|
|
! ! searchance(J) += log(dfloat(cepa0_shortcut(blok+1) - cepa0_shortcut(blok)))
|
|
|
|
! ! end if
|
|
|
|
! ! end do
|
|
|
|
! ! end do
|
|
|
|
! ! end do
|
2016-04-28 09:55:07 +02:00
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! do I_s = 1, N_det_ref
|
|
|
|
! if(mod(I_s,1) == 0) then
|
|
|
|
! call wall_time(wall)
|
|
|
|
! wall = wall-iwall
|
|
|
|
! print *, I_s, "/", N_det_ref, wall * (dfloat(N_det_ref) / dfloat(I_s)), wall, wall * (dfloat(N_det_ref) / dfloat(I_s))-wall
|
|
|
|
! end if
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! do J_s = 1, I_s
|
|
|
|
!
|
2016-04-29 16:23:20 +02:00
|
|
|
! call get_excitation_degree(psi_ref(1,1,J_s), psi_ref(1,1,I_s), degree, N_int)
|
|
|
|
! if(degree > 3) cycle
|
|
|
|
!
|
2016-04-28 09:55:07 +02:00
|
|
|
! if(searchance(I_s) < searchance(J_s)) then
|
|
|
|
! i_I = I_s
|
|
|
|
! J = J_s
|
|
|
|
! else
|
|
|
|
! i_I = J_s
|
|
|
|
! J = I_s
|
|
|
|
! end if
|
2016-04-29 16:23:20 +02:00
|
|
|
!
|
|
|
|
! !$OMP PARALLEL DO default(none) schedule(dynamic,1) shared(delta_ij_old, delta_ii_old) &
|
|
|
|
! !$OMP private(m,m2,kk, i, k, degree, degree2, l, deg, ni, inac, virt) &
|
|
|
|
! !$OMP private(ok,p1,p2,h1,h2,s1,s2, blok, wall, I_s, J_s) &
|
|
|
|
! !$OMP private(phase_iI, phase_Ik, phase_Ji, phase_al, diI, hIi, hJi, delta_JI, dkI, HkI, ci_inv, dia_hla) &
|
|
|
|
! !$OMP private(i_state, contrib, exc_iI, exc_Ik, exc_IJ, det_tmp, det_tmp2) &
|
|
|
|
! !$OMP shared(N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) &
|
|
|
|
! !$OMP shared(lambda_mrcc, hf_bitmask, active_sorb,cepa0_shortcut,det_cepa0,N_states) &
|
|
|
|
! !$OMP shared(i_I, J, h_,det_cepa0_idx, linked, blokMwen, nlink, iwall, searchance)
|
2016-04-28 09:55:07 +02:00
|
|
|
! do kk = 1 , nlink(i_I)
|
|
|
|
! k = linked(kk,i_I)
|
|
|
|
! blok = blokMwen(kk,i_I)
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! call get_excitation(psi_ref(1,1,i_I),det_cepa0(1,1,k),exc_Ik,degree,phase_Ik,N_int)
|
|
|
|
!
|
|
|
|
! call apply_excitation(psi_ref(1,1,J),exc_Ik,det_tmp2,ok,N_int)
|
|
|
|
! if(.not. ok) cycle
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! l = searchDet(det_cepa0(1,1,cepa0_shortcut(blok)), det_tmp2, cepa0_shortcut(blok+1)-cepa0_shortcut(blok), N_int)
|
|
|
|
! if(l == -1) cycle
|
|
|
|
! l = det_cepa0_idx(cepa0_shortcut(blok)-1+l)
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! m = 1
|
|
|
|
! m2 = 1
|
|
|
|
! do while(m <= nlink(i_I) .and. m2 <= nlink(J))
|
|
|
|
! if(linked(m, i_I) < linked(m2, J)) then
|
|
|
|
! m += 1
|
|
|
|
! cycle
|
|
|
|
! else if(linked(m, i_I) > linked(m2, J)) then
|
|
|
|
! m2 += 1
|
|
|
|
! cycle
|
|
|
|
! end if
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! i = det_cepa0_idx(linked(m,i_I))
|
|
|
|
! m += 1
|
|
|
|
! m2 += 1
|
|
|
|
!
|
|
|
|
! do i_state = 1, N_states
|
|
|
|
! if(lambda_mrcc(i_state, i) == 0d0) cycle
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! hJi = h_(J,i)
|
|
|
|
! if(hJi == 0) cycle
|
|
|
|
! hIi = h_(i_I,i)
|
|
|
|
! if(hIi == 0) cycle
|
|
|
|
!
|
|
|
|
! call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int)
|
|
|
|
! if(.not. ok) cycle
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! if(isInCassd(det_tmp, N_int)) cycle
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! call get_excitation(psi_ref(1,1,J),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int)
|
|
|
|
! call get_excitation(det_tmp,psi_non_ref(1,1,l),exc_IJ,degree2,phase_al,N_int)
|
|
|
|
! delta_JI = hJi * hIi * lambda_mrcc(i_state, i) * phase_al * phase_Ji
|
|
|
|
! ci_inv(i_state) = 1.d0 / psi_ref_coef(i_I,i_state)
|
|
|
|
! HkI = h_(i_I,det_cepa0_idx(k))
|
|
|
|
! dkI(i_state) = HkI * lambda_mrcc(i_state, det_cepa0_idx(k))
|
|
|
|
! contrib = dkI(i_state) * delta_JI
|
2016-04-29 16:23:20 +02:00
|
|
|
! !$OMP ATOMIC
|
2016-04-28 09:55:07 +02:00
|
|
|
! delta_ij_old(i_I,l,i_state) += contrib
|
|
|
|
! if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
|
2016-04-29 16:23:20 +02:00
|
|
|
! !$OMP ATOMIC
|
2016-04-28 09:55:07 +02:00
|
|
|
! delta_ii_old(i_I,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(l,i_state)
|
|
|
|
! endif
|
|
|
|
! !
|
|
|
|
! if(l == det_cepa0_idx(k)) cycle
|
|
|
|
! call get_excitation(psi_ref(1,1,I_i),psi_non_ref(1,1,i),exc_IJ,degree2,phase_Ji,N_int)
|
|
|
|
! call get_excitation(det_tmp,det_cepa0(1,1,k),exc_IJ,degree2,phase_al,N_int)
|
|
|
|
! delta_JI = hJi * hIi * lambda_mrcc(i_state, i) * phase_al * phase_Ji
|
|
|
|
!
|
|
|
|
! ci_inv(i_state) = 1.d0 / psi_ref_coef(J,i_state)
|
|
|
|
! HkI = h_(J,l)
|
|
|
|
! dkI(i_state) = HkI * lambda_mrcc(i_state, l)
|
|
|
|
! contrib = dkI(i_state) * delta_JI
|
2016-04-29 16:23:20 +02:00
|
|
|
! !$OMP ATOMIC
|
2016-04-28 09:55:07 +02:00
|
|
|
! delta_ij_old(J,det_cepa0_idx(k),i_state) += contrib
|
|
|
|
! if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5) then
|
2016-04-29 16:23:20 +02:00
|
|
|
! !$OMP ATOMIC
|
2016-04-28 09:55:07 +02:00
|
|
|
! delta_ii_old(J,i_state) -= contrib * ci_inv(i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
|
|
|
|
! end if
|
|
|
|
! enddo !i_state
|
|
|
|
! end do ! while
|
|
|
|
! enddo !kk
|
|
|
|
! enddo !J
|
|
|
|
!
|
|
|
|
! enddo !I
|
|
|
|
!
|
|
|
|
! END_PROVIDER
|
2016-04-29 16:23:20 +02:00
|
|
|
! !
|
|
|
|
!
|
2016-03-04 16:52:46 +01:00
|
|
|
|