LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00
Code Issues Releases Wiki Activity

Extracting davidson from slater rules

Browse Source
This commit is contained in:
Anthony Scemama 2016-09-20 17:29:02 +02:00
parent c99b0cb299
commit 98059e87f5
26 changed files with 291 additions and 382 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -32,7 +32,7 @@ OPENMP : 1 ; Append OpenMP flags
#
[OPT]
FC : -traceback
FCFLAGS : -xHost -O2 -ip -ftz -g
FCFLAGS : -xHost -O0 -ip -ftz -g
# Profiling flags
#################

2
plugins/All_singles/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Generators_restart Perturbation Properties Selectors_no_sorted Utils
Generators_restart Perturbation Properties Selectors_no_sorted Utils Davidson

2
plugins/CAS_SD/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_CAS
Perturbation Selectors_full Generators_CAS Davidson

2
plugins/CID/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Selectors_full SingleRefMethod
Selectors_full SingleRefMethod Davidson

2
plugins/CIS/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Selectors_full SingleRefMethod
Selectors_full SingleRefMethod Davidson

2
plugins/CISD/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Selectors_full SingleRefMethod
Selectors_full SingleRefMethod Davidson

2
plugins/Casino/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants
Determinants Davidson

2
plugins/DDCI_selected/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_CAS
Perturbation Selectors_full Generators_CAS Davidson

2
plugins/FCIdump/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants
Determinants Davidson

2
plugins/FOBOCI/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_no_sorted Hartree_Fock
Perturbation Selectors_no_sorted Hartree_Fock Davidson

22
plugins/Full_CI/H_apply.irp.f
View File

@ -30,28 +30,6 @@ s.unset_openmp()
print s
s = H_apply("select_mono_delta_rho")
s.unset_double_excitations()
s.set_selection_pt2("delta_rho_one_point")
s.unset_openmp()
print s
s = H_apply("pt2_mono_delta_rho")
s.unset_double_excitations()
s.set_perturbation("delta_rho_one_point")
s.unset_openmp()
print s
s = H_apply("select_mono_di_delta_rho")
s.set_selection_pt2("delta_rho_one_point")
s.unset_openmp()
print s
s = H_apply("pt2_mono_di_delta_rho")
s.set_perturbation("delta_rho_one_point")
s.unset_openmp()
print s
END_SHELL

2
plugins/Full_CI/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_full
Perturbation Selectors_full Generators_full Davidson

2
plugins/Perturbation/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Properties Hartree_Fock
Properties Hartree_Fock Davidson

77
plugins/Perturbation/delta_rho_perturbation.irp.f
View File

@ -1,77 +0,0 @@
subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st,minilist,idx_minilist,N_minilist)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
double precision :: i_O1_psi_array(N_st)
double precision :: i_H_psi_array(N_st)
integer, intent(in) :: N_minilist
integer, intent(in) :: idx_minilist(0:N_det_selectors)
integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors)
BEGIN_DOC
! compute the perturbatibe contribution to the Integrated Spin density at z = z_one point of one determinant
!
! for the various n_st states, at various level of theory.
!
! c_pert(i) = <psi(i)|H|det_pert>/(<psi(i)|H|psi(i)> - <det_pert|H|det_pert>)
!
! e_2_pert(i) = c_pert(i) * <det_pert|O|psi(i)>
!
! H_pert_diag(i) = c_pert(i)^2 * <det_pert|O|det_pert>
!
! To get the contribution of the first order :
!
! <O_1> = sum(over i) e_2_pert(i)
!
! To get the contribution of the diagonal elements of the second order :
!
! [ <O_0> + <O_1> + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2]
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,diag_o1_mat_elem_alpha_beta
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e,h,oii,diag_o1_mat_elem
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
! call get_excitation_degree(HF_bitmask,det_pert,degree,N_int)
! if(degree.gt.degree_max_generators+1)then
! H_pert_diag = 0.d0
! e_2_pert = 0.d0
! c_pert = 0.d0
! return
! endif
call i_O1_psi_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int)
do i =1,N_st
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
e_2_pert(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) ) + c_pert(i) * c_pert(i) * oii
H_pert_diag(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) )
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = c_pert(i) * i_O1_psi_array(i)
endif
enddo
end

2
plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Psiref_Utils
Psiref_Utils Davidson

2
plugins/QmcChem/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Determinants
Determinants Davidson

3
src/Davidson/davidson.irp.f Normal file
View File

@ -0,0 +1,3 @@
program davidson
stop 1
end

65
src/Determinants/davidson.irp.f → src/Davidson/diagonalization.irp.f
View File

@ -1,20 +1,3 @@
BEGIN_PROVIDER [ integer, davidson_iter_max ]
implicit none
BEGIN_DOC
! Max number of Davidson iterations
END_DOC
davidson_iter_max = 100
END_PROVIDER
BEGIN_PROVIDER [ integer, davidson_sze_max ]
implicit none
BEGIN_DOC
! Max number of Davidson sizes
END_DOC
ASSERT (davidson_sze_max <= davidson_iter_max)
davidson_sze_max = max(8,2*N_states_diag)
END_PROVIDER
subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)
use bitmasks
implicit none
@ -69,6 +52,9 @@ end
logical function det_inf(key1, key2, Nint)
use bitmasks
implicit none
BEGIN_DOC
! Ordering function for determinants
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: key1(Nint, 2), key2(Nint, 2)
integer :: i,j
@ -91,7 +77,6 @@ end function
subroutine tamiser(key, idx, no, n, Nint, N_key)
use bitmasks
implicit none
BEGIN_DOC
! Uncodumented : TODO
END_DOC
@ -619,47 +604,3 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
)
end
BEGIN_PROVIDER [ character(64), davidson_criterion ]
implicit none
BEGIN_DOC
! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
END_DOC
davidson_criterion = 'residual'
END_PROVIDER
subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)
implicit none
BEGIN_DOC
! True if the Davidson algorithm is converged
END_DOC
integer, intent(in) :: N_st, iterations
logical, intent(out) :: converged
double precision, intent(in) :: energy(N_st), residual(N_st)
double precision, intent(in) :: wall, cpu
double precision :: E(N_st), time
double precision, allocatable, save :: energy_old(:)
if (.not.allocated(energy_old)) then
allocate(energy_old(N_st))
energy_old = 0.d0
endif
E = energy - energy_old
energy_old = energy
if (davidson_criterion == 'energy') then
converged = dabs(maxval(E(1:N_st))) < threshold_davidson
else if (davidson_criterion == 'residual') then
converged = dabs(maxval(residual(1:N_st))) < threshold_davidson
else if (davidson_criterion == 'both') then
converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) &
< threshold_davidson
else if (davidson_criterion == 'wall_time') then
call wall_time(time)
converged = time - wall > threshold_davidson
else if (davidson_criterion == 'cpu_time') then
call cpu_time(time)
converged = time - cpu > threshold_davidson
else if (davidson_criterion == 'iterations') then
converged = iterations >= int(threshold_davidson)
endif
end

16
src/Determinants/diagonalize_CI.irp.f → src/Davidson/diagonalize_CI.irp.f
View File

@ -1,19 +1,3 @@
BEGIN_PROVIDER [ character*(64), diag_algorithm ]
implicit none
BEGIN_DOC
! Diagonalization algorithm (Davidson or Lapack)
END_DOC
if (N_det > N_det_max_jacobi) then
diag_algorithm = "Davidson"
else
diag_algorithm = "Lapack"
endif
if (N_det < N_states_diag) then
diag_algorithm = "Lapack"
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
implicit none

0
src/Determinants/diagonalize_CI_SC2.irp.f → src/Davidson/diagonalize_CI_SC2.irp.f
View File

0
src/Determinants/diagonalize_CI_mono.irp.f → src/Davidson/diagonalize_CI_mono.irp.f
View File

62
src/Davidson/parameters.irp.f Normal file
View File

@ -0,0 +1,62 @@
BEGIN_PROVIDER [ integer, davidson_iter_max ]
implicit none
BEGIN_DOC
! Max number of Davidson iterations
END_DOC
davidson_iter_max = 100
END_PROVIDER
BEGIN_PROVIDER [ integer, davidson_sze_max ]
implicit none
BEGIN_DOC
! Max number of Davidson sizes
END_DOC
ASSERT (davidson_sze_max <= davidson_iter_max)
davidson_sze_max = max(8,2*N_states_diag)
END_PROVIDER
BEGIN_PROVIDER [ character(64), davidson_criterion ]
implicit none
BEGIN_DOC
! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
END_DOC
davidson_criterion = 'residual'
END_PROVIDER
subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)
implicit none
BEGIN_DOC
! True if the Davidson algorithm is converged
END_DOC
integer, intent(in) :: N_st, iterations
logical, intent(out) :: converged
double precision, intent(in) :: energy(N_st), residual(N_st)
double precision, intent(in) :: wall, cpu
double precision :: E(N_st), time
double precision, allocatable, save :: energy_old(:)
if (.not.allocated(energy_old)) then
allocate(energy_old(N_st))
energy_old = 0.d0
endif
E = energy - energy_old
energy_old = energy
if (davidson_criterion == 'energy') then
converged = dabs(maxval(E(1:N_st))) < threshold_davidson
else if (davidson_criterion == 'residual') then
converged = dabs(maxval(residual(1:N_st))) < threshold_davidson
else if (davidson_criterion == 'both') then
converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) &
< threshold_davidson
else if (davidson_criterion == 'wall_time') then
call wall_time(time)
converged = time - wall > threshold_davidson
else if (davidson_criterion == 'cpu_time') then
call cpu_time(time)
converged = time - cpu > threshold_davidson
else if (davidson_criterion == 'iterations') then
converged = iterations >= int(threshold_davidson)
endif
end

190
src/Davidson/u0Hu0.irp.f Normal file
View File

@ -0,0 +1,190 @@
subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint
double precision, intent(out) :: e_0
double precision, intent(in) :: u_0(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
call u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,1,n)
end
subroutine u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint, N_st, sze_8
double precision, intent(out) :: e_0(N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision, allocatable :: H_jj(:), v_0(:,:)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
allocate (H_jj(n), v_0(sze_8,N_st))
do i = 1, n
H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint)
enddo
call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
do i=1,N_st
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
enddo
end
subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
END_DOC
integer, intent(in) :: n,Nint
double precision, intent(out) :: v_0(n)
double precision, intent(in) :: u_0(n)
double precision, intent(in) :: H_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,1,n)
end
subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8
double precision, intent(out) :: v_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij
double precision, allocatable :: vt(:,:)
integer :: i,j,k,l, jj,ii
integer :: i0, j0
integer, allocatable :: shortcut(:,:), sort_idx(:,:)
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
v_0 = 0.d0
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)&
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze_8)
allocate(vt(sze_8,N_st))
Vt = 0.d0
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,1)
do sh2=sh,shortcut(0,1)
exa = 0
do ni=1,Nint
exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1)))
end do
if(exa > 2) then
cycle
end if
do i=shortcut(sh,1),shortcut(sh+1,1)-1
org_i = sort_idx(i,1)
if(sh==sh2) then
endi = i-1
else
endi = shortcut(sh2+1,1)-1
end if
do ni=1,Nint
sorted_i(ni) = sorted(ni,i,1)
enddo
do j=shortcut(sh2,1),endi
org_j = sort_idx(j,1)
ext = exa
do ni=1,Nint
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
end do
if(ext <= 4) then
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
do istate=1,N_st
vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate)
vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate)
enddo
endif
enddo
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
do j=shortcut(sh,2),i-1
org_j = sort_idx(j,2)
ext = 0
do ni=1,Nint
ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2)))
end do
if(ext == 4) then
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
do istate=1,N_st
vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate)
vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate)
enddo
end if
end do
end do
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
do istate=1,N_st
do i=n,1,-1
v_0(i,istate) = v_0(i,istate) + vt(i,istate)
enddo
enddo
!$OMP END CRITICAL
deallocate(vt)
!$OMP END PARALLEL
do istate=1,N_st
do i=1,n
v_0(i,istate) += H_jj(i) * u_0(i,istate)
enddo
enddo
deallocate (shortcut, sort_idx, sorted, version)
end

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

@ -1,5 +1,22 @@
use bitmasks
BEGIN_PROVIDER [ character*(64), diag_algorithm ]
implicit none
BEGIN_DOC
! Diagonalization algorithm (Davidson or Lapack)
END_DOC
if (N_det > N_det_max_jacobi) then
diag_algorithm = "Davidson"
else
diag_algorithm = "Lapack"
endif
if (N_det < N_states_diag) then
diag_algorithm = "Lapack"
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det ]
implicit none
BEGIN_DOC

2
src/Determinants/s2.irp.f
View File

@ -162,7 +162,7 @@ subroutine S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy davidson_criterion
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
v_0 = 0.d0

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

@ -1634,195 +1634,6 @@ subroutine get_occ_from_key(key,occ,Nint)
end
subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint
double precision, intent(out) :: e_0
double precision, intent(in) :: u_0(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
call u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,1,n)
end
subroutine u_0_H_u_0_nstates(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes e_0 = <u_0|H|u_0>/<u_0|u_0>
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint, N_st, sze_8
double precision, intent(out) :: e_0(N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision, allocatable :: H_jj(:), v_0(:,:)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
allocate (H_jj(n), v_0(sze_8,N_st))
do i = 1, n
H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint)
enddo
call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
do i=1,N_st
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
enddo
end
subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
END_DOC
integer, intent(in) :: n,Nint
double precision, intent(out) :: v_0(n)
double precision, intent(in) :: u_0(n)
double precision, intent(in) :: H_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
call H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,1,n)
end
subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
use bitmasks
implicit none
BEGIN_DOC
! Computes v_0 = H|u_0>
!
! n : number of determinants
!
! H_jj : array of <j|H|j>
END_DOC
integer, intent(in) :: N_st,n,Nint, sze_8
double precision, intent(out) :: v_0(sze_8,N_st)
double precision, intent(in) :: u_0(sze_8,N_st)
double precision, intent(in) :: H_jj(n)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision :: hij
double precision, allocatable :: vt(:,:)
integer :: i,j,k,l, jj,ii
integer :: i0, j0
integer, allocatable :: shortcut(:,:), sort_idx(:,:)
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy davidson_criterion
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
v_0 = 0.d0
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)&
!$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze_8)
allocate(vt(sze_8,N_st))
Vt = 0.d0
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,1)
do sh2=sh,shortcut(0,1)
exa = 0
do ni=1,Nint
exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1)))
end do
if(exa > 2) then
cycle
end if
do i=shortcut(sh,1),shortcut(sh+1,1)-1
org_i = sort_idx(i,1)
if(sh==sh2) then
endi = i-1
else
endi = shortcut(sh2+1,1)-1
end if
do ni=1,Nint
sorted_i(ni) = sorted(ni,i,1)
enddo
do j=shortcut(sh2,1),endi
org_j = sort_idx(j,1)
ext = exa
do ni=1,Nint
ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
end do
if(ext <= 4) then
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
do istate=1,N_st
vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate)
vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate)
enddo
endif
enddo
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,2)
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
do j=shortcut(sh,2),i-1
org_j = sort_idx(j,2)
ext = 0
do ni=1,Nint
ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2)))
end do
if(ext == 4) then
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
do istate=1,N_st
vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate)
vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate)
enddo
end if
end do
end do
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
do istate=1,N_st
do i=n,1,-1
v_0(i,istate) = v_0(i,istate) + vt(i,istate)
enddo
enddo
!$OMP END CRITICAL
deallocate(vt)
!$OMP END PARALLEL
do istate=1,N_st
do i=1,n
v_0(i,istate) += H_jj(i) * u_0(i,istate)
enddo
enddo
deallocate (shortcut, sort_idx, sorted, version)
end
subroutine get_double_excitation_phase(det1,det2,exc,phase,Nint)
use bitmasks