LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-22 05:02:15 +02:00
Code Issues Releases Wiki Activity

Merge branch 'master' of github.com:LCPQ/quantum_package

Browse Source
This commit is contained in:
Manu 2014-05-29 10:54:21 +02:00
commit cbc70ec1db
29 changed files with 459 additions and 159 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
scripts/generate_h_apply.py
View File

@ -16,6 +16,8 @@ keys_work
copy_buffer
finalization
generate_psi_guess
init_thread
deinit_thread
""".split()
class H_apply(object):
@ -42,7 +44,7 @@ class H_apply(object):
!$OMP N_elec_in_key_hole_2,ia_ja_pairs,iproc) &
!$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
!$OMP hole_1, particl_1, hole_2, particl_2, &
!$OMP thresh,elec_alpha_num,i_generator)"""
!$OMP elec_alpha_num,i_generator)"""
s["omp_end_parallel"] = "!$OMP END PARALLEL"
s["omp_master"] = "!$OMP MASTER"
s["omp_end_master"] = "!$OMP END MASTER"
@ -54,7 +56,7 @@ class H_apply(object):
s["generate_psi_guess"] = """
! Sort H_jj to find the N_states lowest states
integer :: i,k
integer :: i
integer, allocatable :: iorder(:)
double precision, allocatable :: H_jj(:)
double precision, external :: diag_h_mat_elem
@ -81,9 +83,6 @@ class H_apply(object):
for k in s:
s[k] = ""
s["size_max"] = str(1024*128)
s["set_i_H_j_threshold"] = """
thresh = H_apply_threshold
"""
s["copy_buffer"] = "call copy_h_apply_buffer_to_wf"
self.data = s
@ -109,28 +108,41 @@ class H_apply(object):
integer, intent(in) :: N_st,Nint
double precision, intent(inout) :: sum_e_2_pert_in(N_st)
double precision, intent(inout) :: sum_norm_pert_in(N_st)
double precision, intent(inout) :: sum_H_pert_diag_in
double precision, intent(inout) :: sum_H_pert_diag_in(N_st)
double precision :: sum_e_2_pert(N_st)
double precision :: sum_norm_pert(N_st)
double precision :: sum_H_pert_diag
double precision :: e_2_pert_buffer(N_st,size_max)
double precision :: coef_pert_buffer(N_st,size_max)
double precision :: sum_H_pert_diag(N_st)
double precision, allocatable :: e_2_pert_buffer(:,:)
double precision, allocatable :: coef_pert_buffer(:,:)
ASSERT (Nint == N_int)
"""
self.data["init_thread"] = """
allocate (e_2_pert_buffer(N_st,size_max), coef_pert_buffer(N_st,size_max))
do k=1,N_st
sum_e_2_pert(k) = 0.d0
sum_norm_pert(k) = 0.d0
sum_H_pert_diag(k) = 0.d0
enddo
"""
self.data["deinit_thread"] = """
!$OMP CRITICAL
do k=1,N_st
sum_e_2_pert_in(k) = sum_e_2_pert_in(k) + sum_e_2_pert(k)
sum_norm_pert_in(k) = sum_norm_pert_in(k) + sum_norm_pert(k)
sum_H_pert_diag_in(k) = sum_H_pert_diag_in(k) + sum_H_pert_diag(k)
enddo
!$OMP END CRITICAL
deallocate (e_2_pert_buffer, coef_pert_buffer)
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
sum_e_2_pert = sum_e_2_pert_in
sum_norm_pert = sum_norm_pert_in
sum_H_pert_diag = sum_H_pert_diag_in
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors
"""
self.data["keys_work"] = """
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,Nint)
sum_norm_pert,sum_H_pert_diag,N_st,N_int)
"""%(pert,)
self.data["finalization"] = """
sum_e_2_pert_in = sum_e_2_pert
sum_norm_pert_in = sum_norm_pert
sum_H_pert_diag_in = sum_H_pert_diag
"""
self.data["copy_buffer"] = ""
self.data["generate_psi_guess"] = ""
@ -140,17 +152,19 @@ class H_apply(object):
self.data["decls_main"] = """ integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag
double precision, intent(inout):: H_pert_diag(N_st)
PROVIDE CI_electronic_energy N_det_generators key_pattern_not_in_ref
pt2 = 0.d0
norm_pert = 0.d0
H_pert_diag = 0.d0
do k=1,N_st
pt2(k) = 0.d0
norm_pert(k) = 0.d0
H_pert_diag(k) = 0.d0
enddo
"""
if self.openmp:
self.data["omp_parallel"] += """&
!$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
!$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
!$OMP SHARED(N_st) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
!$OMP PRIVATE(sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
def set_selection_pt2(self,pert):
if self.selection_pt2 is not None:
@ -163,7 +177,7 @@ class H_apply(object):
call copy_h_apply_buffer_to_wf
selection_criterion_min = selection_criterion_min*0.1d0
selection_criterion = selection_criterion_min
call remove_small_contributions
!call remove_small_contributions
"""
self.data["keys_work"] = """
e_2_pert_buffer = 0.d0

2
setup_environment.sh
View File

@ -6,8 +6,6 @@
QPACKAGE_ROOT=${PWD}
IRPF90=$(which irpf90)
if [[ -z ${IRPF90} ]] ;
then
make irpf90

3
src/CISD/SC2.irp.f
View File

@ -143,7 +143,8 @@ end
subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2),Nint
integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2)
integer :: Nint
integer,intent(out) :: i_ok
integer :: ispin,i_hole,k_hole,j_hole,i_particl,k_particl,j_particl,i_trou,degree,exc(0:2,2,2)
double precision :: phase

3
src/CISD_selected/H_apply.irp.f
View File

@ -14,8 +14,7 @@ subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st
implicit none
character*(64), intent(in) :: perturbation
integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st), norm_pert(N_st)
double precision,intent(inout) :: H_pert_diag
double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st)
BEGIN_SHELL [ /usr/bin/env python ]
from perturbation import perturbations

6
src/CISD_selected/README.rst
View File

@ -8,6 +8,12 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/H_apply.irp.f#L13>`_
Undocumented
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/cisd_selection.irp.f#L1>`_
Undocumented
Needed Modules

12
src/CISD_selected/cisd_selection.irp.f
View File

@ -3,27 +3,25 @@ program cisd
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:)
double precision :: H_pert_diag, E_old
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, iter
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st))
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
pt2 = 1.d0
perturbation = "epstein_nesbet"
do while (maxval(abs(pt2(1:N_st))) > 1.d-6)
E_old = CI_energy(1)
call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', E_old
print *, 'E+PT2 = ', E_old+pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
if (abort_all) then
exit
endif
enddo
deallocate(pt2,norm_pert)
deallocate(pt2,norm_pert,H_pert_diag)
end

2
src/Dets/H_apply.irp.f
View File

@ -19,7 +19,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
! Uninitialized. Filled by H_apply subroutines.
END_DOC
integer :: iproc, sze
sze = 10
sze = 100
if (.not.associated(H_apply_buffer)) then
allocate(H_apply_buffer(0:nproc-1))
iproc = 0

80
src/Dets/H_apply_template.f
View File

@ -14,31 +14,33 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
integer(bit_kind),allocatable :: keys_out(:,:,:)
integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2)
integer(bit_kind) :: hole_save(N_int,2)
integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2)
integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2)
integer(bit_kind), allocatable :: hole_save(:,:)
integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:)
integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:)
integer :: ii,i,jj,j,k,ispin,l
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle_tmp(N_int*bit_kind_size,2)
integer :: occ_hole_tmp(N_int*bit_kind_size,2)
integer, allocatable :: occ_particle(:,:), occ_hole(:,:)
integer, allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:)
integer :: kk,pp,other_spin,key_idx
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
double precision :: mo_bielec_integral, thresh
double precision :: mo_bielec_integral
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem
integer :: iproc
$set_i_H_j_threshold
PROVIDE H_apply_threshold
$initialization
iproc = 0
$omp_parallel
!$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max))
allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), &
key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),&
particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), &
occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),&
occ_hole_tmp(N_int*bit_kind_size,2))
$init_thread
!print*,'key_in !!'
!call print_key(key_in)
@ -142,7 +144,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
j_b = occ_particle_tmp(jjj,other_spin)
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if(dabs( mo_bielec_integral(j_a,j_b,i_a,i_b))<thresh)cycle
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
@ -179,7 +180,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if (j_b <= j_a) cycle
if(dabs( mo_bielec_integral(j_a,j_b,i_b,i_a))<thresh)cycle
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
@ -205,7 +205,13 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
$omp_enddo
enddo ! ispin
$keys_work
deallocate (keys_out,ia_ja_pairs)
$deinit_thread
deallocate (ia_ja_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp,&
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp,&
occ_hole_tmp)
$omp_end_parallel
$finalization
abort_here = abort_all
@ -224,34 +230,35 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator $parameters
$declarations
integer ,intent(in) :: i_generator
integer(bit_kind),intent(in) :: key_in(N_int,2)
integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind),allocatable :: keys_out(:,:,:)
integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2)
integer(bit_kind) :: hole_2(N_int,2), particl_2(N_int,2)
integer(bit_kind) :: hole_save(N_int,2)
integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2)
integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2)
integer(bit_kind),allocatable :: hole_save(:,:)
integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:)
integer(bit_kind),allocatable :: hole_tmp(:,:), particle_tmp(:,:)
integer(bit_kind),allocatable :: hole_2(:,:), particl_2(:,:)
integer :: ii,i,jj,j,k,ispin,l
integer :: occ_particle(N_int*bit_kind_size,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_particle_tmp(N_int*bit_kind_size,2)
integer :: occ_hole_tmp(N_int*bit_kind_size,2)
integer,allocatable :: occ_particle(:,:), occ_hole(:,:)
integer,allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:)
integer :: kk,pp,other_spin,key_idx
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
double precision :: thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem
integer :: iproc
$set_i_H_j_threshold
PROVIDE H_apply_threshold
$initialization
iproc = 0
$omp_parallel
!$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max))
allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), &
key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),&
particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), &
occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),&
occ_hole_tmp(N_int*bit_kind_size,2))
$init_thread
!!!! First couple hole particle
do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1))
@ -312,7 +319,13 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator $parameters
$omp_enddo
enddo ! ispin
$keys_work
deallocate (keys_out,ia_ja_pairs)
$deinit_thread
deallocate (ia_ja_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp,&
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp,&
occ_hole_tmp)
$omp_end_parallel
$finalization
@ -330,20 +343,19 @@ subroutine $subroutine($params_main)
$decls_main
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_reference psi_generators
integer :: i_generator
integer :: i_generator, k
print *, irp_here
do i_generator=1,N_det_generators
call $subroutine_monoexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,s_hole ,i_bitmask_gen), &
generators_bitmask(1,1,s_part ,i_bitmask_gen), &
i_generator $params_post)
call $subroutine_diexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,d_hole1,i_bitmask_gen), &
generators_bitmask(1,1,d_part1,i_bitmask_gen), &
generators_bitmask(1,1,d_hole2,i_bitmask_gen), &
generators_bitmask(1,1,d_part2,i_bitmask_gen), &
i_generator $params_post)
call $subroutine_monoexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,s_hole ,i_bitmask_gen), &
generators_bitmask(1,1,s_part ,i_bitmask_gen), &
i_generator $params_post)
if (abort_here) then
exit
endif

12
src/Dets/connected_to_ref.irp.f
View File

@ -32,6 +32,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo
@ -64,6 +67,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo
@ -101,6 +107,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo
@ -140,6 +149,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
endif
if (degree_x2 > 5) then
cycle
else
connected_to_ref = i
return
endif
enddo

4
src/Dets/davidson.irp.f
View File

@ -374,8 +374,8 @@ end
BEGIN_DOC
! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
END_DOC
davidson_criterion = 'both'
davidson_threshold = 1.d-8
davidson_criterion = 'residual'
davidson_threshold = 1.d-6
END_PROVIDER
subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)

71
src/Dets/determinants.irp.f
View File

@ -57,18 +57,18 @@ END_PROVIDER
ifirst = 1
psi_det = 0_bit_kind
psi_coef = 0.d0
integer :: i
do i=1,N_int
psi_det(i,1,1) = HF_bitmask(i,1)
psi_det(i,2,1) = HF_bitmask(i,2)
enddo
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
endif
integer :: i
do i=1,N_int
psi_det(i,1,1) = HF_bitmask(i,1)
psi_det(i,2,1) = HF_bitmask(i,2)
enddo
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
END_PROVIDER
@ -80,3 +80,54 @@ BEGIN_PROVIDER [ integer, N_det_reference ]
N_det_reference = N_det
ASSERT (N_det_reference > 0)
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (N_det) ]
implicit none
BEGIN_DOC
! Contribution of determinants to the state-averaged density
END_DOC
integer :: i,j,k
double precision :: f
f = 1.d0/dble(N_states)
do i=1,N_det
psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f
enddo
do k=2,N_states
do i=1,N_det
psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + &
psi_coef(i,k)*psi_coef(i,k)*f
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted, (N_int,2,N_det) ]
&BEGIN_PROVIDER [ double precision, psi_coef_sorted, (N_det,N_states) ]
&BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted, (N_det) ]
implicit none
BEGIN_DOC
! Wave function sorted by determinants (state-averaged)
END_DOC
integer :: i,j,k
integer, allocatable :: iorder(:)
allocate ( iorder(N_det) )
do i=1,N_det
psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib(i)
iorder(i) = i
enddo
call dsort(psi_average_norm_contrib_sorted,iorder,N_det)
!DIR$ IVDEP
do i=1,N_det
do j=1,N_int
psi_det_sorted(j,1,i) = psi_det(j,1,iorder(i))
psi_det_sorted(j,2,i) = psi_det(j,2,iorder(i))
enddo
do k=1,N_states
psi_coef_sorted(i,k) = psi_coef(iorder(i),k)
enddo
psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib_sorted(i)
enddo
deallocate(iorder)
END_PROVIDER

4
src/Dets/diagonalize_CI.irp.f
View File

@ -22,8 +22,12 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
END_DOC
integer :: j
character*(8) :: st
call write_time(output_Dets)
do j=1,N_states
CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st))
enddo
END_PROVIDER

0
src/Full_CI/ASSUMPTIONS.rst Normal file
View File

10
src/Full_CI/H_apply.irp.f Normal file
View File

@ -0,0 +1,10 @@
use bitmasks
BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply("FCI",openmp=True)
s.set_selection_pt2("epstein_nesbet_2x2")
print s
END_SHELL

8
src/Full_CI/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Full_CI/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MonoInts MOs Nuclei Output Selectors_full Utils Perturbation

39
src/Full_CI/README.rst Normal file
View File

@ -0,0 +1,39 @@
==============
Full_CI Module
==============
Performs a perturbatively selected Full-CI.
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Full_CI/full_ci.irp.f#L1>`_
Undocumented
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_

28
src/Full_CI/full_ci.irp.f Normal file
View File

@ -0,0 +1,28 @@
program cisd
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
pt2 = 1.d0
diag_algorithm = "Lapack"
do while (maxval(abs(pt2(1:N_st))) > 1.d-4)
call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
if (abort_all) then
exit
endif
enddo
deallocate(pt2,norm_pert)
end

0
src/Generators_full/ASSUMPTIONS.rst Normal file
View File

8
src/Generators_full/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Generators_full/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils

7
src/Generators_full/README.rst Normal file
View File

@ -0,0 +1,7 @@
======================
Generators_full Module
======================
All the determinants of the wave function are generators. In this way, the Full CI
space is explored.

44
src/Generators_full/generators.irp.f Normal file
View File

@ -0,0 +1,44 @@
use bitmasks
BEGIN_PROVIDER [ double precision, generators_threshold ]
implicit none
BEGIN_DOC
! Percentage of the norm of the state-averaged wave function to
! consider for the generators
END_DOC
generators_threshold = 0.97d0
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_generators ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the number of generators is 1 : the
! Hartree-Fock determinant
END_DOC
integer :: i
double precision :: norm
call write_time(output_dets)
norm = 0.d0
N_det_generators = N_det
do i=1,N_det
norm = norm + psi_average_norm_contrib_sorted(i)
if (norm > generators_threshold) then
N_det_generators = i-1
exit
endif
enddo
N_det_generators = max(N_det_generators,1)
call write_int(output_dets,N_det_generators,'Number of generators')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
psi_generators = psi_det_sorted
END_PROVIDER

33
src/Generators_full/tests/Makefile Normal file
View File

@ -0,0 +1,33 @@
OPENMP =1
PROFILE =0
DEBUG = 0
IRPF90+= -I tests
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
.PHONY: clean executables serial_tests parallel_tests
all: clean executables serial_tests parallel_tests
parallel_tests: $(REF_FILES)
@echo ; echo " ---- Running parallel tests ----" ; echo
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
serial_tests: $(REF_FILES)
@echo ; echo " ---- Running serial tests ----" ; echo
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
executables: $(wildcard *.irp.f) veryclean
$(MAKE) -C ..
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
clean:
$(MAKE) -C .. clean
veryclean:
$(MAKE) -C .. veryclean

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 Generators_full Full_CI

7
src/Perturbation/Moller_plesset.irp.f
View File

@ -3,7 +3,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
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
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
@ -21,7 +21,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
double precision :: diag_H_mat_elem
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase,delta_e
double precision :: phase,delta_e,h
integer :: h1,h2,p1,p2,s1,s2
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
@ -32,8 +32,9 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
delta_e = 1.d0/delta_e
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
H_pert_diag(i) = h
c_pert(i) = i_H_psi_array(i) *delta_e
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo

101
src/Perturbation/epstein_nesbet.irp.f
View File

@ -1,15 +1,15 @@
subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
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
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states.
! for the various N_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
@ -18,14 +18,15 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem
double precision :: diag_H_mat_elem, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
H_pert_diag(i) = h
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_H_psi_array(i)
else
c_pert(i) = 0.d0
@ -35,18 +36,18 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
end
subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
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
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various n_st states.
! for the various N_st states.
!
! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
!
@ -55,13 +56,14 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,delta_e
double precision :: diag_H_mat_elem,delta_e, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
delta_e = H_pert_diag - CI_electronic_energy(i)
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
H_pert_diag(i) = h
delta_e = h - CI_electronic_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
@ -72,19 +74,19 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
end
subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
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
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states,
! for the various N_st states,
!
! but with the correction in the denominator
!
@ -105,40 +107,41 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij
double precision :: diag_H_mat_elem,accu_e_corr,hij,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array,idx_repeat)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,n_st
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,N_st
H_pert_diag(i) = h
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag)
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else
c_pert(i) = 0.d0
endif
enddo
end
subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
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
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various n_st states.
! for the various N_st states.
!
! but with the correction in the denominator
!
@ -159,19 +162,20 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e
double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array,idx_repeat)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
do i = 1, idx_repeat(0)
call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij)
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,n_st
delta_e = H_pert_diag - CI_electronic_energy(i)
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
do i =1,N_st
H_pert_diag(i) = h
delta_e = h - CI_electronic_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
@ -182,19 +186,19 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
end
subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
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
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states,
! for the various N_st states,
!
! but with the correction in the denominator
!
@ -219,10 +223,10 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j
double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array,idx_repeat)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
call i_H_j(ref_bitmask,det_pert,Nint,h0j)
do i = 1, idx_repeat(0)
@ -230,13 +234,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1)
enddo
accu_e_corr = accu_e_corr / psi_selectors_coef(1,1)
H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - H_pert_diag)
c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - h)
e_2_pert(1) = c_pert(i) * h0j
do i =2,n_st
if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else
c_pert(i) = 0.d0
endif

4
src/Perturbation/perturbation_template.f
View File

@ -13,7 +13,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
integer :: i,k, c_ref
integer :: connected_to_ref
@ -38,7 +38,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
coef_pert_buffer(k,i) = c_pert(k)
sum_norm_pert(k) += c_pert(k) * c_pert(k)
sum_e_2_pert(k) += e_2_pert(k)
sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag
sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag(k)
enddo
enddo

66
src/Perturbation/selection.irp.f
View File

@ -31,13 +31,12 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
l=H_apply_buffer(iproc)%N_det
do i=1,n_selected
s = 0.d0
is_selected = .False.
do j=1,N_st
s -= e_2_pert_buffer(j,i)
s = -e_2_pert_buffer(j,i)
is_selected = s > selection_criterion*selection_criterion_factor .or. is_selected
enddo
ASSERT (s>=-1.d-8)
is_selected = s > selection_criterion * selection_criterion_factor
if (is_selected) then
l = l+1
@ -60,8 +59,10 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c
ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
selection_criterion = smax
selection_criterion_min = smin
!$OMP CRITICAL
selection_criterion = max(selection_criterion,smax)
selection_criterion_min = min(selection_criterion_min,smin)
!$OMP END CRITICAL
end
BEGIN_PROVIDER [ double precision, selection_criterion ]
@ -71,7 +72,7 @@ end
BEGIN_DOC
! Threshold to select determinants. Set by selection routines.
END_DOC
selection_criterion = .1d0
selection_criterion = 1.d0
selection_criterion_factor = 0.01d0
selection_criterion_min = selection_criterion
@ -80,33 +81,52 @@ END_PROVIDER
subroutine remove_small_contributions
implicit none
BEGIN_DOC
! Remove determinants with small contributions
! Remove determinants with small contributions. N_states is assumed to be
! provided.
END_DOC
integer :: i,j,k, N_removed
logical keep
N_removed = 0
do i=N_det,1,-1
keep = .False.
logical, allocatable :: keep(:)
double precision :: i_H_psi_array(N_states)
allocate (keep(N_det))
call diagonalize_CI
do i=1,N_det
keep(i) = .True.
enddo
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,i_H_psi_array) &
!$OMP SHARED(k,psi_det_sorted,psi_coef_sorted,N_int,N_det,psi_det_size,N_states, &
!$OMP selection_criterion_min,keep,N_det_generators) &
!$OMP REDUCTION(+:N_removed)
!$OMP DO
do i=2*N_det_generators+1, N_det
call i_H_psi(psi_det_sorted(1,1,i),psi_det_sorted,psi_coef_sorted,N_int,min(N_det,2*N_det_generators),psi_det_size,N_states,i_H_psi_array)
keep(i) = .False.
do j=1,N_states
keep = keep .or. (dabs(psi_coef(i,j)) > selection_criterion_min)
keep(i) = keep(i) .or. (-(psi_coef_sorted(i,j)*i_H_psi_array(j)) > selection_criterion_min)
enddo
if (.not.keep) then
do k=i+1,N_det
do j=1,N_int
psi_det(j,1,k-1) = psi_det(j,1,k)
psi_det(j,2,k-1) = psi_det(j,2,k)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
N_removed = 0
k = 0
do i=1, N_det
if (keep(i)) then
k += 1
do j=1,N_int
psi_det(j,1,k) = psi_det_sorted(j,1,i)
psi_det(j,2,k) = psi_det_sorted(j,2,i)
enddo
do j=1,N_states
do k=i+1,N_det
psi_coef(k-1,j) = psi_coef(k,j)
enddo
psi_coef(k,j) = psi_coef_sorted(i,j)
enddo
else
N_removed += 1
endif
enddo
deallocate(keep)
if (N_removed > 0) then
N_det -= N_removed
N_det = N_det - N_removed
SOFT_TOUCH N_det psi_det psi_coef
call write_int(output_dets,N_removed, 'Removed determinants')
endif
end