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

Add the FOBOCI routines

Browse Source
This commit is contained in:
Emmanuel Giner 2016-02-17 17:15:54 +01:00
parent cc71d6f83d
commit bc6c26fb73
34 changed files with 4912 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ocaml/.gitignore vendored
View File

@ -4,6 +4,7 @@ ezfio.ml
Git.ml
Input_auto_generated.ml
Input_determinants.ml
Input_foboci.ml
Input_hartree_fock.ml
Input_integrals_bielec.ml
Input_perturbation.ml

6
ocaml/qp_edit.ml
View File

@ -22,6 +22,7 @@ type keyword =
| Integrals_bielec
| Perturbation
| Properties
| Foboci
| Determinants
;;
@ -37,6 +38,7 @@ let keyword_to_string = function
| Integrals_bielec -> "Integrals_bielec"
| Perturbation -> "Perturbation"
| Properties -> "Properties"
| Foboci -> "Foboci"
| Determinants -> "Determinants"
;;
@ -96,6 +98,8 @@ let get s =
f Perturbation.(read, to_rst)
| Properties ->
f Properties.(read, to_rst)
| Foboci ->
f Foboci.(read, to_rst)
| Determinants ->
f Determinants.(read, to_rst)
end
@ -140,6 +144,7 @@ let set str s =
| Integrals_bielec -> write Integrals_bielec.(of_rst, write) s
| Perturbation -> write Perturbation.(of_rst, write) s
| Properties -> write Properties.(of_rst, write) s
| Foboci -> write Foboci.(of_rst, write) s
| Determinants -> write Determinants.(of_rst, write) s
| Electrons -> write Electrons.(of_rst, write) s
| Determinants_by_hand -> write Determinants_by_hand.(of_rst, write) s
@ -193,6 +198,7 @@ let run check_only ezfio_filename =
Integrals_bielec ;
Perturbation ;
Properties ;
Foboci ;
Determinants ;
Mo_basis;
Determinants_by_hand ;

5
plugins/CAS_SD/cas_sd_selected.irp.f
View File

@ -59,8 +59,9 @@ program full_ci
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
! threshold_selectors = 1.d0
! threshold_generators = 0.999d0
! soft_touch threshold_selectors
call H_apply_CAS_SD_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'

30
plugins/FOBOCI/EZFIO.cfg Normal file
View File

@ -0,0 +1,30 @@
[threshold_singles]
type: double precision
doc: threshold to select the pertinent single excitations at second order
interface: ezfio,provider,ocaml
default: 0.01
[threshold_fobo_dm]
type: double precision
doc: threshold to eliminate small density matrix elements in the fobo procedure
interface: ezfio,provider,ocaml
default: 0.00001
[do_it_perturbative]
type: logical
doc: if true, you do the FOBOCI calculation perturbatively
interface: ezfio,provider,ocaml
default: .False.
[second_order_h]
type: logical
doc: if true, you do the FOBOCI calculation using second order intermediate Hamiltonian
interface: ezfio,provider,ocaml
default: .False.
[do_all_2p]
type: logical
doc: if true, you do all 2p type excitation on the LMCT
interface: ezfio,provider,ocaml
default: .True.

53
plugins/FOBOCI/H_apply.irp.f Normal file
View File

@ -0,0 +1,53 @@
use bitmasks
BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply("just_1h_1p")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.filter_only_1h1p()
print s
s = H_apply("all_but_1h_and_1p")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.filter_1h()
s.filter_1p()
print s
s = H_apply("standard")
s.set_selection_pt2("epstein_nesbet")
s.unset_skip()
print s
s = H_apply("just_mono")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.unset_double_excitations()
print s
s = H_apply("just_mono_no_1h_no_1p")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.unset_double_excitations()
s.filter_1h()
s.filter_1p()
print s
s = H_apply("just_mono_no_1h_no_1p_no_2p")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.unset_double_excitations()
s.filter_1h()
s.filter_1p()
s.filter_2p()
print s
END_SHELL

605
plugins/FOBOCI/H_apply_dressed_autonom.irp.f Normal file
View File

@ -0,0 +1,605 @@
subroutine H_apply_dressed_pert_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref )
use omp_lib
use bitmasks
implicit none
BEGIN_DOC
! Generate all double excitations of key_in using the bit masks of holes and
! particles.
! Assume N_int is already provided.
END_DOC
integer,parameter :: size_max = 3072
integer, intent(in) :: Ndet_generators
double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
integer ,intent(in) :: i_generator
integer(bit_kind),intent(in) :: key_in(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), intent(in) :: hole_2(N_int,2), particl_2(N_int,2)
integer, intent(in) :: iproc_in
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 :: key_union_hole_part(:)
integer :: ii,i,jj,j,k,ispin,l
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
logical :: is_a_two_holes_two_particles
integer, allocatable :: ia_ja_pairs(:,:,:)
integer, allocatable :: ib_jb_pairs(:,:)
double precision :: diag_H_mat_elem
integer :: iproc
integer :: jtest_vvvv
integer(omp_lock_kind), save :: lck, ifirst=0
if (ifirst == 0) then
!$ call omp_init_lock(lck)
ifirst=1
endif
logical :: check_double_excitation
logical :: is_a_1h1p
logical :: b_cycle
check_double_excitation = .True.
iproc = iproc_in
PROVIDE elec_num_tab
! !$OMP PARALLEL DEFAULT(SHARED) &
! !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
! !$OMP occ_particle,occ_hole,j_a,k_a,other_spin, &
! !$OMP hole_save,ispin,jj,l_a,ib_jb_pairs,array_pairs, &
! !$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, &
! !$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
! !$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, &
! !$OMP N_elec_in_key_hole_2,ia_ja_pairs,key_union_hole_part) &
! !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
! !$OMP hole_1, particl_1, hole_2, particl_2, &
! !$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc)
!$ iproc = omp_get_thread_num()
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),key_union_hole_part(N_int))
!!!! First couple hole particle
do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1))
hole(j,2) = iand(hole_1(j,2),key_in(j,2))
particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1))
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole(1,1),occ_hole(1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole(1,2),occ_hole(1,2),N_elec_in_key_hole_1(2),N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2), &
ib_jb_pairs(2,0:(elec_alpha_num)*mo_tot_num))
do ispin=1,2
i=0
do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole
i_a = occ_hole(ii,ispin)
ASSERT (i_a > 0)
ASSERT (i_a <= mo_tot_num)
do jj=1,N_elec_in_key_part_1(ispin) !particle
j_a = occ_particle(jj,ispin)
ASSERT (j_a > 0)
ASSERT (j_a <= mo_tot_num)
i += 1
ia_ja_pairs(1,i,ispin) = i_a
ia_ja_pairs(2,i,ispin) = j_a
enddo
enddo
ia_ja_pairs(1,0,ispin) = i
enddo
key_idx = 0
integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b
integer(bit_kind) :: test(N_int,2)
double precision :: accu
logical, allocatable :: array_pairs(:,:)
allocate(array_pairs(mo_tot_num,mo_tot_num))
accu = 0.d0
do ispin=1,2
other_spin = iand(ispin,1)+1
if (abort_here) then
exit
endif
! !$OMP DO SCHEDULE (static)
do ii=1,ia_ja_pairs(1,0,ispin)
if (abort_here) then
cycle
endif
i_a = ia_ja_pairs(1,ii,ispin)
ASSERT (i_a > 0)
ASSERT (i_a <= mo_tot_num)
j_a = ia_ja_pairs(2,ii,ispin)
ASSERT (j_a > 0)
ASSERT (j_a <= mo_tot_num)
hole = key_in
k = ishft(i_a-1,-bit_kind_shift)+1
j = i_a-ishft(k-1,bit_kind_shift)-1
hole(k,ispin) = ibclr(hole(k,ispin),j)
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
!!!! Second couple hole particle
do j = 1, N_int
hole_tmp(j,1) = iand(hole_2(j,1),hole(j,1))
hole_tmp(j,2) = iand(hole_2(j,2),hole(j,2))
particle_tmp(j,1) = iand(xor(particl_2(j,1),hole(j,1)),particl_2(j,1))
particle_tmp(j,2) = iand(xor(particl_2(j,2),hole(j,2)),particl_2(j,2))
enddo
call bitstring_to_list(particle_tmp(1,1),occ_particle_tmp(1,1),N_elec_in_key_part_2(1),N_int)
call bitstring_to_list(particle_tmp(1,2),occ_particle_tmp(1,2),N_elec_in_key_part_2(2),N_int)
call bitstring_to_list(hole_tmp (1,1),occ_hole_tmp (1,1),N_elec_in_key_hole_2(1),N_int)
call bitstring_to_list(hole_tmp (1,2),occ_hole_tmp (1,2),N_elec_in_key_hole_2(2),N_int)
! hole = a^(+)_j_a(ispin) a_i_a(ispin)|key_in> : mono exc :: orb(i_a,ispin) --> orb(j_a,ispin)
hole_save = hole
! Build array of the non-zero integrals of second excitation
array_pairs = .True.
if (ispin == 1) then
integer :: jjj
i=0
do kk = 1,N_elec_in_key_hole_2(other_spin)
i_b = occ_hole_tmp(kk,other_spin)
ASSERT (i_b > 0)
ASSERT (i_b <= mo_tot_num)
do jjj=1,N_elec_in_key_part_2(other_spin) ! particule
j_b = occ_particle_tmp(jjj,other_spin)
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if (array_pairs(i_b,j_b)) then
i+= 1
ib_jb_pairs(1,i) = i_b
ib_jb_pairs(2,i) = j_b
endif
enddo
enddo
ib_jb_pairs(1,0) = i
do kk = 1,ib_jb_pairs(1,0)
hole = hole_save
i_b = ib_jb_pairs(1,kk)
j_b = ib_jb_pairs(2,kk)
k = ishft(i_b-1,-bit_kind_shift)+1
j = i_b-ishft(k-1,bit_kind_shift)-1
hole(k,other_spin) = ibclr(hole(k,other_spin),j)
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,other_spin) = ibset(key(k,other_spin),l)
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref)
key_idx = 0
endif
if (abort_here) then
exit
endif
enddo
endif
! does all the mono excitations of the same spin
i=0
do kk = 1,N_elec_in_key_hole_2(ispin)
i_b = occ_hole_tmp(kk,ispin)
if (i_b <= i_a.or.i_b == j_a) cycle
ASSERT (i_b > 0)
ASSERT (i_b <= mo_tot_num)
do jjj=1,N_elec_in_key_part_2(ispin) ! particule
j_b = occ_particle_tmp(jjj,ispin)
ASSERT (j_b > 0)
ASSERT (j_b <= mo_tot_num)
if (j_b <= j_a) cycle
if (array_pairs(i_b,j_b)) then
i+= 1
ib_jb_pairs(1,i) = i_b
ib_jb_pairs(2,i) = j_b
endif
enddo
enddo
ib_jb_pairs(1,0) = i
do kk = 1,ib_jb_pairs(1,0)
hole = hole_save
i_b = ib_jb_pairs(1,kk)
j_b = ib_jb_pairs(2,kk)
k = ishft(i_b-1,-bit_kind_shift)+1
j = i_b-ishft(k-1,bit_kind_shift)-1
hole(k,ispin) = ibclr(hole(k,ispin),j)
key = hole
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,ispin) = ibset(key(k,ispin),l)
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref)
key_idx = 0
endif
if (abort_here) then
exit
endif
enddo ! kk
enddo ! ii
! !$OMP ENDDO NOWAIT
enddo ! ispin
call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref)
deallocate (ia_ja_pairs, ib_jb_pairs, &
keys_out, hole_save, &
key,hole, particle, hole_tmp,&
particle_tmp, occ_particle, &
occ_hole, occ_particle_tmp,&
occ_hole_tmp,array_pairs,key_union_hole_part)
! !$OMP END PARALLEL
end
subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref )
use omp_lib
use bitmasks
implicit none
BEGIN_DOC
! Generate all single excitations of key_in using the bit masks of holes and
! particles.
! Assume N_int is already provided.
END_DOC
integer,parameter :: size_max = 3072
integer, intent(in) :: Ndet_generators
double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
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, intent(in) :: iproc_in
integer(bit_kind),allocatable :: keys_out(:,:,:)
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,allocatable :: occ_particle(:,:), occ_hole(:,:)
integer,allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:)
integer,allocatable :: ib_jb_pairs(:,:)
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)
logical :: is_a_two_holes_two_particles
integer(bit_kind), allocatable :: key_union_hole_part(:)
integer, allocatable :: ia_ja_pairs(:,:,:)
logical, allocatable :: array_pairs(:,:)
double precision :: diag_H_mat_elem
integer(omp_lock_kind), save :: lck, ifirst=0
integer :: iproc
logical :: check_double_excitation
logical :: is_a_1h1p
logical :: is_a_1h
logical :: is_a_1p
iproc = iproc_in
check_double_excitation = .True.
check_double_excitation = .False.
if (ifirst == 0) then
ifirst=1
!!$ call omp_init_lock(lck)
endif
PROVIDE elec_num_tab
! !$OMP PARALLEL DEFAULT(SHARED) &
! !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
! !$OMP occ_particle,occ_hole,j_a,k_a,other_spin, &
! !$OMP hole_save,ispin,jj,l_a,ib_jb_pairs,array_pairs, &
! !$OMP accu,i_a,hole_tmp,particle_tmp,occ_particle_tmp, &
! !$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
! !$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, &
! !$OMP N_elec_in_key_hole_2,ia_ja_pairs,key_union_hole_part) &
! !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
! !$OMP hole_1, particl_1, hole_2, particl_2, &
! !$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc)
!!$ iproc = omp_get_thread_num()
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),key_union_hole_part(N_int))
!!!! First couple hole particle
do j = 1, N_int
hole(j,1) = iand(hole_1(j,1),key_in(j,1))
hole(j,2) = iand(hole_1(j,2),key_in(j,2))
particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1))
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
do ispin=1,2
i=0
do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole
i_a = occ_hole(ii,ispin)
do jj=1,N_elec_in_key_part_1(ispin) !particule
j_a = occ_particle(jj,ispin)
i += 1
ia_ja_pairs(1,i,ispin) = i_a
ia_ja_pairs(2,i,ispin) = j_a
enddo
enddo
ia_ja_pairs(1,0,ispin) = i
enddo
key_idx = 0
integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b
integer(bit_kind) :: test(N_int,2)
double precision :: accu
accu = 0.d0
integer :: jjtest,na,nb
do ispin=1,2
other_spin = iand(ispin,1)+1
! !$OMP DO SCHEDULE (static)
do ii=1,ia_ja_pairs(1,0,ispin)
i_a = ia_ja_pairs(1,ii,ispin)
j_a = ia_ja_pairs(2,ii,ispin)
hole = key_in
k = ishft(i_a-1,-bit_kind_shift)+1
j = i_a-ishft(k-1,bit_kind_shift)-1
hole(k,ispin) = ibclr(hole(k,ispin),j)
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
na = 0
nb = 0
! if (is_a_1h(hole)) then
! cycle
! endif
! if (is_a_1p(hole)) then
! cycle
! endif
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2)
enddo
if (key_idx == size_max) then
call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref)
key_idx = 0
endif
enddo ! ii
! !$OMP ENDDO NOWAIT
enddo ! ispin
call standard_dress(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref)
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,key_union_hole_part)
! !$OMP END PARALLEL
end
subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref)
implicit none
use omp_lib
use bitmasks
BEGIN_DOC
! Calls H_apply on the HF determinant and selects all connected single and double
! excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
END_DOC
integer, intent(in) :: Ndet_generators
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators),E_ref
double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
integer :: i_generator, nmax
double precision :: wall_0, wall_1
integer(omp_lock_kind) :: lck
integer(bit_kind), allocatable :: mask(:,:,:)
integer :: ispin, k
integer :: iproc
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map
nmax = mod( Ndet_generators,nproc )
! !$ call omp_init_lock(lck)
call start_progress(Ndet_generators,'Selection (norm)',0.d0)
call wall_time(wall_0)
iproc = 0
allocate( mask(N_int,2,6) )
do i_generator=1,nmax
progress_bar(1) = i_generator
if (abort_here) then
exit
endif
! ! Create bit masks for holes and particles
do ispin=1,2
do k=1,N_int
mask(k,ispin,s_hole) = &
iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), &
psi_det_generators_input(k,ispin,i_generator) )
mask(k,ispin,s_part) = &
iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), &
not(psi_det_generators_input(k,ispin,i_generator)) )
mask(k,ispin,d_hole1) = &
iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), &
psi_det_generators_input(k,ispin,i_generator) )
mask(k,ispin,d_part1) = &
iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), &
not(psi_det_generators_input(k,ispin,i_generator)) )
mask(k,ispin,d_hole2) = &
iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), &
psi_det_generators_input(k,ispin,i_generator) )
mask(k,ispin,d_part2) = &
iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), &
not(psi_det_generators_input(k,ispin,i_generator)) )
enddo
enddo
if(.False.)then
call H_apply_dressed_pert_diexc(psi_det_generators_input(1,1,i_generator), &
mask(1,1,d_hole1), mask(1,1,d_part1), &
mask(1,1,d_hole2), mask(1,1,d_part2), &
i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref)
endif
if(.True.)then
call H_apply_dressed_pert_monoexc(psi_det_generators_input(1,1,i_generator), &
mask(1,1,s_hole ), mask(1,1,s_part ), &
i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref)
endif
call wall_time(wall_1)
if (wall_1 - wall_0 > 2.d0) then
write(output_determinants,*) &
100.*float(i_generator)/float(Ndet_generators), '% in ', wall_1-wall_0, 's'
wall_0 = wall_1
endif
enddo
deallocate( mask )
! !$OMP PARALLEL DEFAULT(SHARED) &
! !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc)
call wall_time(wall_0)
! !$ iproc = omp_get_thread_num()
allocate( mask(N_int,2,6) )
! !$OMP DO SCHEDULE(dynamic,1)
do i_generator=nmax+1,Ndet_generators
if (iproc == 0) then
progress_bar(1) = i_generator
endif
if (abort_here) then
cycle
endif
! Create bit masks for holes and particles
do ispin=1,2
do k=1,N_int
mask(k,ispin,s_hole) = &
iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), &
psi_det_generators_input(k,ispin,i_generator) )
mask(k,ispin,s_part) = &
iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), &
not(psi_det_generators_input(k,ispin,i_generator)) )
mask(k,ispin,d_hole1) = &
iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), &
psi_det_generators_input(k,ispin,i_generator) )
mask(k,ispin,d_part1) = &
iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), &
not(psi_det_generators_input(k,ispin,i_generator)) )
mask(k,ispin,d_hole2) = &
iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), &
psi_det_generators_input(k,ispin,i_generator) )
mask(k,ispin,d_part2) = &
iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), &
not (psi_det_generators_input(k,ispin,i_generator)) )
enddo
enddo
if(.False.)then
call H_apply_dressed_pert_diexc(psi_det_generators_input(1,1,i_generator), &
mask(1,1,d_hole1), mask(1,1,d_part1), &
mask(1,1,d_hole2), mask(1,1,d_part2), &
i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref)
endif
if(.True.)then
call H_apply_dressed_pert_monoexc(psi_det_generators_input(1,1,i_generator), &
mask(1,1,s_hole ), mask(1,1,s_part ), &
i_generator, iproc , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref)
endif
! !$ call omp_set_lock(lck)
call wall_time(wall_1)
if (wall_1 - wall_0 > 2.d0) then
write(output_determinants,*) &
100.*float(i_generator)/float(Ndet_generators), '% in ', wall_1-wall_0, 's'
wall_0 = wall_1
endif
! !$ call omp_unset_lock(lck)
enddo
! !$OMP END DO
deallocate( mask )
! !$OMP END PARALLEL
! !$ call omp_destroy_lock(lck)
abort_here = abort_all
call stop_progress
end

1
plugins/FOBOCI/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Perturbation Generators_restart Selectors_no_sorted

12
plugins/FOBOCI/README.rst Normal file
View File

@ -0,0 +1,12 @@
======
FOBOCI
======
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

362
plugins/FOBOCI/all_singles.irp.f Normal file
View File

@ -0,0 +1,362 @@
subroutine all_single
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision,allocatable :: E_before(:)
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
selection_criterion = 0.d0
soft_touch selection_criterion
threshold_davidson = 1.d-5
soft_touch threshold_davidson davidson_criterion
i = 0
print*,'Doing all the mono excitations !'
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
print*,'N_det_generators = ',N_det_generators
pt2=-1.d0
E_before = ref_bitmask_energy
print*,'Initial Step '
print*,'Inital determinants '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
n_det_max = 100000
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_mono(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
E_before = CI_energy
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
print*,'Final Step '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
! call save_wavefunction
deallocate(pt2,norm_pert,E_before)
end
subroutine all_single_no_1h_or_1p
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision,allocatable :: E_before(:)
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
threshold_davidson = 1.d-5
soft_touch threshold_davidson davidson_criterion
i = 0
print*,'Doing all the mono excitations !'
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
print*,'N_det_generators = ',N_det_generators
pt2=-1.d0
E_before = ref_bitmask_energy
print*,'Initial Step '
print*,'Inital determinants '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
n_det_max = 100000
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_mono_no_1h_no_1p(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
E_before = CI_energy
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
print*,'Final Step '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
! call save_wavefunction
deallocate(pt2,norm_pert,E_before)
end
subroutine all_single_no_1h_or_1p_or_2p
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision,allocatable :: E_before(:)
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
selection_criterion = 0.d0
soft_touch selection_criterion
threshold_davidson = 1.d-5
soft_touch threshold_davidson davidson_criterion
i = 0
print*,'Doing all the mono excitations !'
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
print*,'N_det_generators = ',N_det_generators
pt2=-1.d0
E_before = ref_bitmask_energy
print*,'Initial Step '
print*,'Inital determinants '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
n_det_max = 100000
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_mono_no_1h_no_1p_no_2p(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
E_before = CI_energy
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
print*,'Final Step '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
! call save_wavefunction
deallocate(pt2,norm_pert,E_before)
end
subroutine all_2p
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision,allocatable :: E_before(:)
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
selection_criterion = 0.d0
soft_touch selection_criterion
threshold_davidson = 1.d-5
soft_touch threshold_davidson davidson_criterion
i = 0
print*,''
print*,''
print*,''
print*,''
print*,''
print*,'*****************************'
print*,'Doing all the 2P excitations'
print*,'*****************************'
print*,''
print*,''
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
print*,'N_det_generators = ',N_det_generators
pt2=-1.d0
E_before = ref_bitmask_energy
print*,'Initial Step '
print*,'Inital determinants '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
n_det_max = 100000
i = 0
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_standard(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
E_before = CI_energy
enddo
print*,'Final Step '
print*,'N_det = ',N_det
do i = 1, N_states_diag
print*,''
print*,'i = ',i
print*,'E = ',CI_energy(i)
print*,'S^2 = ',CI_eigenvectors_s2(i)
enddo
! call save_wavefunction
deallocate(pt2,norm_pert,E_before)
end
subroutine all_1h_1p_routine
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision :: E_before
integer :: n_det_before
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
i = 0
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
pt2=-1.d0
E_before = ref_bitmask_energy
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_1h_1p(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
E_before = CI_energy(1)
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
enddo
deallocate(pt2,norm_pert)
end
subroutine all_but_1h_1p_routine
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
double precision :: E_before
integer :: n_det_before
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
i = 0
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
pt2=-1.d0
E_before = ref_bitmask_energy
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_all_but_1h_and_1p(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
E_before = CI_energy(1)
if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0
endif
enddo
deallocate(pt2,norm_pert)
end

243
plugins/FOBOCI/all_singles_split.irp.f Normal file
View File

@ -0,0 +1,243 @@
subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,Ndet_generators_input,dressing_matrix)
implicit none
use bitmasks
integer, intent(in) :: Ndet_generators_input
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators_input)
double precision, intent(inout) :: dressing_matrix(Ndet_generators_input,Ndet_generators_input)
double precision, intent(in) :: psi_coef_generators_input(ndet_generators_input,n_states)
integer :: i,i_hole
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
do i = 1, n_inact_orb
i_hole = list_inact(i)
print*,''
print*,'Doing all the single excitations from the orbital '
print*,i_hole
print*,''
print*,''
threshold_davidson = 1.d-4
soft_touch threshold_davidson davidson_criterion
call modify_bitmasks_for_hole(i_hole)
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_generators_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
call set_psi_det_as_input_psi(ndet_generators_input,psi_det_generators_input,psi_coef_generators_input)
call all_single
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
enddo
n_det_max_jacobi = 1000
soft_touch n_det_max_jacobi
end
subroutine all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
implicit none
use bitmasks
double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators)
integer :: i,i_hole
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
integer :: n_det_1h1p,n_det_2h1p
integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
integer(bit_kind), allocatable :: psi_2h1p(:,:,:)
double precision, allocatable :: psi_ref_coef_out(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_2h1p(:,:)
call all_single_no_1h_or_1p
threshold_davidson = 1.d-12
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p)
allocate(psi_ref_out(N_int,2,N_det_generators))
allocate(psi_1h1p(N_int,2,n_det_1h1p))
allocate(psi_2h1p(N_int,2,n_det_2h1p))
allocate(psi_ref_coef_out(N_det_generators,N_states))
allocate(psi_coef_1h1p(n_det_1h1p,N_states))
allocate(psi_coef_2h1p(n_det_2h1p,N_states))
call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p)
call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_1h1p,psi_coef_1h1p,n_det_1h1p)
call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_2h1p,psi_coef_2h1p,n_det_2h1p)
deallocate(psi_ref_out)
deallocate(psi_1h1p)
deallocate(psi_2h1p)
deallocate(psi_ref_coef_out)
deallocate(psi_coef_1h1p)
deallocate(psi_coef_2h1p)
end
subroutine all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
implicit none
use bitmasks
double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
double precision, intent(inout) :: dressing_matrix_2h1p(N_det_generators,N_det_generators)
integer :: i,i_hole
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
integer :: n_det_1h1p,n_det_2h1p
integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
integer(bit_kind), allocatable :: psi_2h1p(:,:,:)
double precision, allocatable :: psi_ref_coef_out(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_2h1p(:,:)
do i = 1, n_inact_orb
i_hole = list_inact(i)
print*,''
print*,'Doing all the single excitations from the orbital '
print*,i_hole
print*,''
print*,''
threshold_davidson = 1.d-4
soft_touch threshold_davidson davidson_criterion
selection_criterion_factor = 1.d-4
soft_touch selection_criterion_factor selection_criterion selection_criterion_min
call modify_bitmasks_for_hole(i_hole)
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_generators_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators)
call set_psi_det_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators)
call all_single_no_1h_or_1p
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p)
allocate(psi_ref_out(N_int,2,N_det_generators))
allocate(psi_1h1p(N_int,2,n_det_1h1p))
allocate(psi_2h1p(N_int,2,n_det_2h1p))
allocate(psi_ref_coef_out(N_det_generators,N_states))
allocate(psi_coef_1h1p(n_det_1h1p,N_states))
allocate(psi_coef_2h1p(n_det_2h1p,N_states))
call split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p)
call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_1h1p,psi_coef_1h1p,n_det_1h1p)
call provide_matrix_dressing_general(dressing_matrix_2h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_2h1p,psi_coef_2h1p,n_det_2h1p)
deallocate(psi_ref_out)
deallocate(psi_1h1p)
deallocate(psi_2h1p)
deallocate(psi_ref_coef_out)
deallocate(psi_coef_1h1p)
deallocate(psi_coef_2h1p)
enddo
n_det_max_jacobi = 1000
soft_touch n_det_max_jacobi
end
subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
implicit none
use bitmasks
double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators)
integer :: i,i_hole
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
integer :: n_det_1h1p,n_det_1h2p
integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
integer(bit_kind), allocatable :: psi_1h2p(:,:,:)
double precision, allocatable :: psi_ref_coef_out(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_1h2p(:,:)
do i = 1, n_inact_orb
i_hole = list_inact(i)
print*,''
print*,'Doing all the single excitations from the orbital '
print*,i_hole
print*,''
print*,''
threshold_davidson = 1.d-4
soft_touch threshold_davidson davidson_criterion
selection_criterion_factor = 1.d-4
soft_touch selection_criterion_factor selection_criterion selection_criterion_min
call modify_bitmasks_for_hole(i_hole)
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_generators_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators)
call set_psi_det_as_input_psi(n_det_generators,psi_det_generators,psi_coef_generators)
call all_single_no_1h_or_1p
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p)
allocate(psi_ref_out(N_int,2,N_det_generators))
allocate(psi_1h1p(N_int,2,n_det_1h1p))
allocate(psi_1h2p(N_int,2,n_det_1h2p))
allocate(psi_ref_coef_out(N_det_generators,N_states))
allocate(psi_coef_1h1p(n_det_1h1p,N_states))
allocate(psi_coef_1h2p(n_det_1h2p,N_states))
call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p)
call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_1h1p,psi_coef_1h1p,n_det_1h1p)
call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_1h2p,psi_coef_1h2p,n_det_1h2p)
deallocate(psi_ref_out)
deallocate(psi_1h1p)
deallocate(psi_1h2p)
deallocate(psi_ref_coef_out)
deallocate(psi_coef_1h1p)
deallocate(psi_coef_1h2p)
enddo
n_det_max_jacobi = 1000
soft_touch n_det_max_jacobi
end
subroutine all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
implicit none
use bitmasks
double precision, intent(inout) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
double precision, intent(inout) :: dressing_matrix_1h2p(N_det_generators,N_det_generators)
integer :: i,i_hole
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
integer :: n_det_1h1p,n_det_1h2p
integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
integer(bit_kind), allocatable :: psi_1h2p(:,:,:)
double precision, allocatable :: psi_ref_coef_out(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_1h2p(:,:)
call all_single_no_1h_or_1p_or_2p
threshold_davidson = 1.d-12
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p)
allocate(psi_ref_out(N_int,2,N_det_generators))
allocate(psi_1h1p(N_int,2,n_det_1h1p))
allocate(psi_1h2p(N_int,2,n_det_1h2p))
allocate(psi_ref_coef_out(N_det_generators,N_states))
allocate(psi_coef_1h1p(n_det_1h1p,N_states))
allocate(psi_coef_1h2p(n_det_1h2p,N_states))
call split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p)
call provide_matrix_dressing_general(dressing_matrix_1h1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_1h1p,psi_coef_1h1p,n_det_1h1p)
call provide_matrix_dressing_general(dressing_matrix_1h2p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_1h2p,psi_coef_1h2p,n_det_1h2p)
deallocate(psi_ref_out)
deallocate(psi_1h1p)
deallocate(psi_1h2p)
deallocate(psi_ref_coef_out)
deallocate(psi_coef_1h1p)
deallocate(psi_coef_1h2p)
end

218
plugins/FOBOCI/create_1h_or_1p.irp.f Normal file
View File

@ -0,0 +1,218 @@
subroutine create_restart_and_1h(i_hole)
implicit none
use bitmasks
integer, intent(in) :: i_hole
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,i_part_act,ispin,k,l,i_ok
integer :: n_new_det
integer(bit_kind), allocatable :: new_det(:,:,:)
integer(bit_kind), allocatable :: old_psi_det(:,:,:)
allocate (old_psi_det(N_int,2,n_det))
do i = 1, N_det
do j = 1, N_int
old_psi_det(j,1,i) = psi_det(j,1,i)
old_psi_det(j,2,i) = psi_det(j,2,i)
enddo
enddo
n_new_det = 0
do j = 1, n_act_orb
i_part_act = list_act(j) ! index of the particle in the active space
do i = 1, N_det
do ispin = 1,2
do k = 1, N_int
key_tmp(k,1) = psi_det(k,1,i)
key_tmp(k,2) = psi_det(k,2,i)
enddo
call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok)
if(i_ok .ne. 1)cycle
n_new_det +=1
enddo
enddo
enddo
integer :: N_det_old
N_det_old = N_det
N_det += n_new_det
allocate (new_det(N_int,2,n_new_det))
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, N_det_old
do k = 1, N_int
psi_det(k,1,i) = old_psi_det(k,1,i)
psi_det(k,2,i) = old_psi_det(k,2,i)
enddo
enddo
n_new_det = 0
do j = 1, n_act_orb
i_part_act = list_act(j) ! index of the particle in the active space
do i = 1, N_det_old
do ispin = 1,2
do k = 1, N_int
key_tmp(k,1) = psi_det(k,1,i)
key_tmp(k,2) = psi_det(k,2,i)
enddo
call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok)
if(i_ok .ne. 1)cycle
n_new_det +=1
do k = 1, N_int
psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1)
psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2)
enddo
psi_coef(n_det_old+n_new_det,:) = 0.d0
enddo
enddo
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
call remove_duplicates_in_psi_det(found_duplicates)
end
subroutine create_restart_and_1p(i_particle)
implicit none
integer, intent(in) :: i_particle
use bitmasks
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,i_hole_act,ispin,k,l,i_ok
integer :: n_new_det
integer(bit_kind), allocatable :: new_det(:,:,:)
integer(bit_kind), allocatable :: old_psi_det(:,:,:)
allocate (old_psi_det(N_int,2,n_det))
do i = 1, N_det
do j = 1, N_int
old_psi_det(j,1,i) = psi_det(j,1,i)
old_psi_det(j,2,i) = psi_det(j,2,i)
enddo
enddo
n_new_det = 0
do j = 1, n_act_orb
i_hole_act = list_act(j) ! index of the particle in the active space
do i = 1, N_det
do ispin = 1,2
do k = 1, N_int
key_tmp(k,1) = psi_det(k,1,i)
key_tmp(k,2) = psi_det(k,2,i)
enddo
call do_mono_excitation(key_tmp,i_hole_act,i_particle,ispin,i_ok)
if(i_ok .ne. 1)cycle
n_new_det +=1
enddo
enddo
enddo
integer :: N_det_old
N_det_old = N_det
N_det += n_new_det
allocate (new_det(N_int,2,n_new_det))
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, N_det_old
do k = 1, N_int
psi_det(k,1,i) = old_psi_det(k,1,i)
psi_det(k,2,i) = old_psi_det(k,2,i)
enddo
enddo
n_new_det = 0
do j = 1, n_act_orb
i_hole_act = list_act(j) ! index of the particle in the active space
do i = 1, N_det_old
do ispin = 1,2
do k = 1, N_int
key_tmp(k,1) = psi_det(k,1,i)
key_tmp(k,2) = psi_det(k,2,i)
enddo
call do_mono_excitation(key_tmp,i_hole_act,i_particle,ispin,i_ok)
if(i_ok .ne. 1)cycle
n_new_det +=1
do k = 1, N_int
psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1)
psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2)
enddo
psi_coef(n_det_old+n_new_det,:) = 0.d0
enddo
enddo
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
call remove_duplicates_in_psi_det(found_duplicates)
end
subroutine create_restart_1h_1p(i_hole,i_part)
implicit none
use bitmasks
integer, intent(in) :: i_hole
integer, intent(in) :: i_part
integer :: i,j,i_part_act,ispin,k,l,i_ok
integer(bit_kind) :: key_tmp(N_int,2)
integer :: n_new_det
integer(bit_kind), allocatable :: new_det(:,:,:)
integer(bit_kind), allocatable :: old_psi_det(:,:,:)
allocate (old_psi_det(N_int,2,n_det))
do i = 1, N_det
do j = 1, N_int
old_psi_det(j,1,i) = psi_det(j,1,i)
old_psi_det(j,2,i) = psi_det(j,2,i)
enddo
enddo
n_new_det = 0
i_part_act = i_part ! index of the particle in the active space
do i = 1, N_det
do ispin = 1,2
do k = 1, N_int
key_tmp(k,1) = psi_det(k,1,i)
key_tmp(k,2) = psi_det(k,2,i)
enddo
call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok)
if(i_ok .ne. 1)cycle
n_new_det +=1
enddo
enddo
integer :: N_det_old
N_det_old = N_det
N_det += n_new_det
allocate (new_det(N_int,2,n_new_det))
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
do i = 1, N_det_old
do k = 1, N_int
psi_det(k,1,i) = old_psi_det(k,1,i)
psi_det(k,2,i) = old_psi_det(k,2,i)
enddo
enddo
n_new_det = 0
i_part_act = i_part ! index of the particle in the active space
do i = 1, N_det_old
do ispin = 1,2
do k = 1, N_int
key_tmp(k,1) = psi_det(k,1,i)
key_tmp(k,2) = psi_det(k,2,i)
enddo
call do_mono_excitation(key_tmp,i_hole,i_part_act,ispin,i_ok)
if(i_ok .ne. 1)cycle
n_new_det +=1
do k = 1, N_int
psi_det(k,1,n_det_old+n_new_det) = key_tmp(k,1)
psi_det(k,2,n_det_old+n_new_det) = key_tmp(k,2)
enddo
psi_coef(n_det_old+n_new_det,:) = 0.d0
enddo
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
call remove_duplicates_in_psi_det(found_duplicates)
end

133
plugins/FOBOCI/density_matrix.irp.f Normal file
View File

@ -0,0 +1,133 @@
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_generators_restart, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_generators_restart, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, norm_generators_restart]
implicit none
BEGIN_DOC
! Alpha and beta one-body density matrix for the generators restart
END_DOC
integer :: j,k,l,m
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer :: exc(0:2,2,2),n_occ_alpha
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
integer :: degree_respect_to_HF_k
integer :: degree_respect_to_HF_l,index_ref_generators_restart
double precision :: inv_coef_ref_generators_restart
integer :: i
do i = 1, N_det_generators_restart
! Find the reference determinant for intermediate normalization
call get_excitation_degree(ref_generators_restart,psi_det_generators_restart(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
inv_coef_ref_generators_restart = 1.d0/psi_coef_generators_restart(i,1)
exit
endif
enddo
norm_generators_restart = 0.d0
do i = 1, N_det_generators_restart
psi_coef_generators_restart(i,1) = psi_coef_generators_restart(i,1) * inv_coef_ref_generators_restart
norm_generators_restart += psi_coef_generators_restart(i,1)**2
enddo
one_body_dm_mo_alpha_generators_restart = 0.d0
one_body_dm_mo_beta_generators_restart = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
!$OMP tmp_a, tmp_b, n_occ_alpha)&
!$OMP SHARED(psi_det_generators_restart,psi_coef_generators_restart,N_int,elec_alpha_num,&
!$OMP elec_beta_num,one_body_dm_mo_alpha_generators_restart,one_body_dm_mo_beta_generators_restart,N_det_generators_restart,mo_tot_num_align,&
!$OMP mo_tot_num,N_states, state_average_weight)
allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) )
tmp_a = 0.d0
tmp_b = 0.d0
!$OMP DO SCHEDULE(dynamic)
do k=1,N_det_generators_restart
call bitstring_to_list(psi_det_generators_restart(1,1,k), occ(1,1), n_occ_alpha, N_int)
call bitstring_to_list(psi_det_generators_restart(1,2,k), occ(1,2), n_occ_alpha, N_int)
do m=1,N_states
ck = psi_coef_generators_restart(k,m)*psi_coef_generators_restart(k,m) * state_average_weight(m)
do l=1,elec_alpha_num
j = occ(l,1)
tmp_a(j,j) += ck
enddo
do l=1,elec_beta_num
j = occ(l,2)
tmp_b(j,j) += ck
enddo
enddo
do l=1,k-1
call get_excitation_degree(psi_det_generators_restart(1,1,k),psi_det_generators_restart(1,1,l),degree,N_int)
if (degree /= 1) then
cycle
endif
call get_mono_excitation(psi_det_generators_restart(1,1,k),psi_det_generators_restart(1,1,l),exc,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do m=1,N_states
ckl = psi_coef_generators_restart(k,m) * psi_coef_generators_restart(l,m) * phase * state_average_weight(m)
if (s1==1) then
tmp_a(h1,p1) += ckl
tmp_a(p1,h1) += ckl
else
tmp_b(h1,p1) += ckl
tmp_b(p1,h1) += ckl
endif
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
one_body_dm_mo_alpha_generators_restart = one_body_dm_mo_alpha_generators_restart + tmp_a
!$OMP END CRITICAL
!$OMP CRITICAL
one_body_dm_mo_beta_generators_restart = one_body_dm_mo_beta_generators_restart + tmp_b
!$OMP END CRITICAL
deallocate(tmp_a,tmp_b)
!$OMP BARRIER
!$OMP END PARALLEL
do i = 1, mo_tot_num
print*,'DM restat',i,one_body_dm_mo_beta_generators_restart(i,i) + one_body_dm_mo_alpha_generators_restart(i,i)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_mo_generators_restart, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! One-body density matrix for the generators_restart
END_DOC
one_body_dm_mo_generators_restart = one_body_dm_mo_alpha_generators_restart + one_body_dm_mo_beta_generators_restart
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_spin_density_mo_generators_restart, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! rho(alpha) - rho(beta)
END_DOC
one_body_spin_density_mo_generators_restart = one_body_dm_mo_alpha_generators_restart - one_body_dm_mo_beta_generators_restart
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_osoci, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_osoci, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Alpha and beta one-body density matrix that will be used for the OSOCI approach
END_DOC
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha_1h1p, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta_1h1p, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Alpha and beta one-body density matrix that will be used for the 1h1p approach
END_DOC
END_PROVIDER

35
plugins/FOBOCI/diag_fock_inactiv_virt.irp.f Normal file
View File

@ -0,0 +1,35 @@
subroutine diag_inactive_virt_and_update_mos
implicit none
integer :: i,j,i_inact,j_inact,i_virt,j_virt
double precision :: tmp(mo_tot_num_align,mo_tot_num)
character*(64) :: label
tmp = 0.d0
do i = 1, mo_tot_num
tmp(i,i) = Fock_matrix_mo(i,i)
enddo
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = i+1, n_inact_orb
j_inact = list_inact(j)
tmp(i_inact,j_inact) = Fock_matrix_mo(i_inact,j_inact)
tmp(j_inact,i_inact) = Fock_matrix_mo(j_inact,i_inact)
enddo
enddo
do i = 1, n_virt_orb
i_virt = list_virt(i)
do j = i+1, n_virt_orb
j_virt = list_virt(j)
tmp(i_virt,j_virt) = Fock_matrix_mo(i_virt,j_virt)
tmp(j_virt,i_virt) = Fock_matrix_mo(j_virt,i_virt)
enddo
enddo
label = "Canonical"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
soft_touch mo_coef
end

358
plugins/FOBOCI/dress_simple.irp.f Normal file
View File

@ -0,0 +1,358 @@
subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_generator,n_selected,det_buffer,Nint,iproc,psi_det_generators_input,E_ref)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint, iproc
integer, intent(in) :: Ndet_generators,size_buffer
double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
integer(bit_kind), intent(in) :: det_buffer(Nint,2,size_buffer)
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
integer :: i,j,k,m
integer :: new_size
integer :: degree(Ndet_generators)
integer :: idx(0:Ndet_generators)
logical :: good
integer :: c_ref
integer :: connected_to_ref
double precision :: hka, haa
double precision :: haj
double precision :: f
integer :: connected_to_ref_by_mono
logical :: is_in_wavefunction
double precision :: H_array(Ndet_generators)
double precision :: H_matrix_tmp(Ndet_generators+1,Ndet_generators+1)
double precision :: eigenvectors(Ndet_generators+1,Ndet_generators+1), eigenvalues(Ndet_generators+1)
double precision :: contrib,lambda_i,accu
do k = 1, Ndet_generators
call i_h_j(psi_det_generators_input(1,1,k),psi_det_generators_input(1,1,k),Nint,hka)
H_matrix_tmp(k,k) = hka
do j = k+1, Ndet_generators
call i_h_j(psi_det_generators_input(1,1,k),psi_det_generators_input(1,1,j),Nint,hka)
H_matrix_tmp(k,j) = hka
H_matrix_tmp(j,k) = hka
enddo
H_matrix_tmp(k,Ndet_generators+1) = 0.d0
enddo
do i=1,n_selected
c_ref = connected_to_ref_by_mono(det_buffer(1,1,i),psi_det_generators_input,N_int,i_generator,Ndet_generators)
if (c_ref /= 0) then
cycle
endif
if (is_in_wavefunction(det_buffer(1,1,i),Nint)) then
cycle
endif
call get_excitation_degree_vector(psi_det_generators_input,det_buffer(1,1,i),degree,N_int,Ndet_generators,idx)
H_array = 0.d0
do k=1,idx(0)
call i_h_j(det_buffer(1,1,i),psi_det_generators_input(1,1,idx(k)),Nint,hka)
H_array(idx(k)) = hka
enddo
call i_h_j(det_buffer(1,1,i),det_buffer(1,1,i),Nint,haa)
f = 1.d0/(E_ref-haa)
if(second_order_h)then
lambda_i = f
else
! You write the new Hamiltonian matrix
do k = 1, Ndet_generators
H_matrix_tmp(k,Ndet_generators+1) = H_array(k)
H_matrix_tmp(Ndet_generators+1,k) = H_array(k)
enddo
H_matrix_tmp(Ndet_generators+1,Ndet_generators+1) = haa
! Then diagonalize it
call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,Ndet_generators+1,Ndet_generators+1)
! Then you extract the effective denominator
accu = 0.d0
do k = 1, Ndet_generators
accu += eigenvectors(k,1) * H_array(k)
enddo
lambda_i = eigenvectors(Ndet_generators+1,1)/accu
endif
do k=1,idx(0)
contrib = H_array(idx(k)) * H_array(idx(k)) * lambda_i
delta_ij_generators_(idx(k), idx(k)) += contrib
do j=k+1,idx(0)
contrib = H_array(idx(k)) * H_array(idx(j)) * lambda_i
delta_ij_generators_(idx(k), idx(j)) += contrib
delta_ij_generators_(idx(j), idx(k)) += contrib
enddo
enddo
! H_matrix_tmp_bis(idx(k),idx(k)) += contrib
! H_matrix_tmp_bis(idx(k),idx(j)) += contrib
! H_matrix_tmp_bis(idx(j),idx(k)) += contrib
! do k = 1, Ndet_generators
! do j = 1, Ndet_generators
! H_matrix_tmp_bis(k,j) = H_matrix_tmp(k,j)
! enddo
! enddo
! double precision :: H_matrix_tmp_bis(Ndet_generators,Ndet_generators)
! double precision :: eigenvectors_bis(Ndet_generators,Ndet_generators), eigenvalues_bis(Ndet_generators)
! call lapack_diag(eigenvalues_bis,eigenvectors_bis,H_matrix_tmp_bis,Ndet_generators,Ndet_generators)
! print*,'f,lambda_i = ',f,lambda_i
! print*,'eigenvalues_bi(1)',eigenvalues_bis(1)
! print*,'eigenvalues ',eigenvalues(1)
! do k = 1, Ndet_generators
! print*,'coef,coef_dres = ', eigenvectors(k,1), eigenvectors_bis(k,1)
! enddo
! pause
! accu = 0.d0
! do k = 1, Ndet_generators
! do j = 1, Ndet_generators
! accu += eigenvectors(k,1) * eigenvectors(j,1) * (H_matrix_tmp(k,j) + delta_ij_generators_(k,j))
! enddo
! enddo
! print*,'accu,eigv = ',accu,eigenvalues(1)
! pause
enddo
end
subroutine is_a_good_candidate(threshold,is_ok,verbose)
use bitmasks
implicit none
double precision, intent(in) :: threshold
logical, intent(out) :: is_ok
logical, intent(in) :: verbose
integer :: l,k,m
double precision,allocatable :: dressed_H_matrix(:,:)
double precision,allocatable :: psi_coef_diagonalized_tmp(:,:)
integer(bit_kind), allocatable :: psi_det_generators_input(:,:,:)
allocate(psi_det_generators_input(N_int,2,N_det_generators),dressed_H_matrix(N_det_generators,N_det_generators))
allocate(psi_coef_diagonalized_tmp(N_det_generators,N_states))
dressed_H_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_int
psi_det_generators_input(l,1,k) = psi_det_generators(l,1,k)
psi_det_generators_input(l,2,k) = psi_det_generators(l,2,k)
enddo
enddo
!call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input)
call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose)
if(do_it_perturbative)then
if(is_ok)then
N_det = N_det_generators
do m = 1, N_states
do k = 1, N_det_generators
do l = 1, N_int
psi_det(l,1,k) = psi_det_generators_input(l,1,k)
psi_det(l,2,k) = psi_det_generators_input(l,2,k)
enddo
psi_coef(k,m) = psi_coef_diagonalized_tmp(k,m)
enddo
enddo
touch psi_coef psi_det N_det
endif
endif
deallocate(psi_det_generators_input,dressed_H_matrix,psi_coef_diagonalized_tmp)
end
subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
integer, intent(in) :: Ndet_generators
double precision, intent(in) :: threshold
logical, intent(in) :: verbose
logical, intent(out) :: is_ok
double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states)
double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators)
integer :: i,j,degree,index_ref_generators_restart,i_count,k,i_det_no_ref
double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij
double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average
logical :: is_a_ref_det(Ndet_generators)
is_a_ref_det = .False.
do i = 1, N_det_generators
do j = 1, N_det_generators_restart
call get_excitation_degree(psi_det_generators_input(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int)
if(degree == 0)then
is_a_ref_det(i) = .True.
exit
endif
enddo
enddo
do i = 1, Ndet_generators
call get_excitation_degree(ref_generators_restart,psi_det_generators_input(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
endif
do j = 1, Ndet_generators
call i_h_j(psi_det_generators_input(1,1,j),psi_det_generators_input(1,1,i),N_int,hij) ! Fill the zeroth order H matrix
dressed_H_matrix(i,j) = hij
enddo
enddo
i_det_no_ref = 0
diag_h_mat_average = 0.d0
do i = 1, Ndet_generators
if(is_a_ref_det(i))cycle
i_det_no_ref +=1
diag_h_mat_average+=dressed_H_matrix(i,i)
enddo
diag_h_mat_average = diag_h_mat_average/dble(i_det_no_ref)
print*,'diag_h_mat_average = ',diag_h_mat_average
print*,'ref h_mat = ',dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart)
integer :: number_of_particles, number_of_holes
! Filter the the MLCT that are higher than 27.2 eV in energy with respect to the reference determinant
do i = 1, Ndet_generators
if(is_a_ref_det(i))cycle
if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then
if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then
is_ok = .False.
return
endif
endif
! Filter the the LMCT that are higher than 54.4 eV in energy with respect to the reference determinant
if(number_of_holes(psi_det_generators_input(1,1,i)).eq.1 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.0)then
if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then
is_ok = .False.
return
endif
endif
exit
enddo
call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix
double precision :: s2,E_ref(N_states)
integer :: i_state(N_states)
integer :: n_state_good
n_state_good = 0
if(s2_eig)then
do i = 1, Ndet_generators
call get_s2_u0(psi_det_generators_input,eigvectors(1,i),Ndet_generators,Ndet_generators,s2)
print*,'s2 = ',s2
print*,dabs(s2-expected_s2)
if(dabs(s2-expected_s2).le.0.3d0)then
n_state_good +=1
i_state(n_state_good) = i
E_ref(n_state_good) = eigvalues(i)
endif
if(n_state_good==N_states)then
exit
endif
enddo
else
do i = 1, N_states
i_state(i) = i
E_ref(i) = eigvalues(i)
enddo
endif
do i = 1,N_states
print*,'i_state = ',i_state(i)
enddo
do k = 1, N_states
print*,'state ',k
do i = 1, Ndet_generators
psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k))
psi_coef_ref(i,k) = eigvectors(i,i_state(k))
print*,'psi_coef_ref(i) = ',psi_coef_ref(i,k)
enddo
enddo
if(verbose)then
print*,'Zeroth order space :'
do i = 1, Ndet_generators
write(*,'(10(F16.8),X)')dressed_H_matrix(i,:)
enddo
print*,''
print*,'Zeroth order space Diagonalized :'
do k = 1, N_states
print*,'state ',k
do i = 1, Ndet_generators
print*,'coef, <I|H|I> = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i)
enddo
enddo
endif
double precision :: E_ref_average
E_ref_average = 0.d0
do i = 1, N_states
E_ref_average += E_ref(i)
enddo
E_ref_average = E_ref_average / dble(N_states)
call H_apply_dressed_pert(dressed_H_matrix,Ndet_generators,psi_det_generators_input,E_ref_average) ! Calculate the dressing of the H matrix
if(verbose)then
print*,'Zeroth order space Dressed by outer space:'
do i = 1, Ndet_generators
write(*,'(10(F16.8),X)')dressed_H_matrix(i,:)
enddo
endif
call lapack_diagd(eigvalues,eigvectors,dressed_H_matrix,Ndet_generators,Ndet_generators) ! Diagonalize the Dressed_H_matrix
integer :: i_good_state(0:N_states)
i_good_state(0) = 0
do i = 1, Ndet_generators
call get_s2_u0(psi_det_generators_input,eigvectors(1,i),Ndet_generators,Ndet_generators,s2)
! State following
do k = 1, N_states
accu = 0.d0
do j =1, Ndet_generators
accu += eigvectors(j,i) * psi_coef_ref(j,k)
enddo
if(dabs(accu).ge.0.8d0)then
i_good_state(0) +=1
i_good_state(i_good_state(0)) = i
endif
enddo
if(i_good_state(0)==N_states)then
exit
endif
enddo
do i = 1, N_states
i_state(i) = i_good_state(i)
E_ref(i) = eigvalues(i_good_state(i))
enddo
double precision :: accu
accu = 0.d0
do k = 1, N_states
do i = 1, Ndet_generators
psi_coef_diagonalized_tmp(i,k) = eigvectors(i,i_state(k)) / eigvectors(index_ref_generators_restart,i_state(k))
enddo
enddo
if(verbose)then
do k = 1, N_states
print*,'state ',k
do i = 1, Ndet_generators
print*,'coef, <I|H+Delta H|I> = ',psi_coef_diagonalized_tmp(i,k),dressed_H_matrix(i,i)-dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart),is_a_ref_det(i)
enddo
enddo
endif
is_ok = .False.
do i = 1, Ndet_generators
if(is_a_ref_det(i))cycle
do k = 1, N_states
if(dabs(psi_coef_diagonalized_tmp(i,k)) .gt.threshold)then
is_ok = .True.
exit
endif
enddo
if(is_ok)then
exit
endif
enddo
if(verbose)then
print*,'is_ok = ',is_ok
endif
end

5
plugins/FOBOCI/fobo_coupled_ci.irp.f Normal file
View File

@ -0,0 +1,5 @@
program osoci_program
implicit none
call new_approach
! call save_natural_mos
end

18
plugins/FOBOCI/fobo_diff_dm.irp.f Normal file
View File

@ -0,0 +1,18 @@
program osoci_program
call debug_det(ref_bitmask,N_int)
implicit none
call FOBOCI_lmct_mlct_old_thr
call provide_all_the_rest
end
subroutine provide_all_the_rest
implicit none
integer :: i
call update_one_body_dm_mo
call provide_properties
call save_osoci_natural_mos
call save_mos
end

315
plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f Normal file
View File

@ -0,0 +1,315 @@
subroutine FOBOCI_lmct_mlct_old_thr
use bitmasks
implicit none
integer :: i,j,k,l
integer(bit_kind),allocatable :: unpaired_bitmask(:,:)
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha, n_occ_beta
double precision :: norm_tmp(N_states),norm_total(N_states)
logical :: test_sym
double precision :: thr,hij
double precision :: threshold
double precision, allocatable :: dressing_matrix(:,:)
logical :: verbose,is_ok
verbose = .True.
threshold = threshold_singles
print*,'threshold = ',threshold
thr = 1.d-12
allocate(unpaired_bitmask(N_int,2))
allocate (occ(N_int*bit_kind_size,2))
do i = 1, N_int
unpaired_bitmask(i,1) = unpaired_alpha_electrons(i)
unpaired_bitmask(i,2) = unpaired_alpha_electrons(i)
enddo
norm_total = 0.d0
call initialize_density_matrix_osoci
call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int)
print*,''
print*,''
print*,'mulliken spin population analysis'
accu =0.d0
do i = 1, nucl_num
accu += mulliken_spin_densities(i)
print*,i,nucl_charge(i),mulliken_spin_densities(i)
enddo
print*,''
print*,''
print*,'DOING FIRST LMCT !!'
do i = 1, n_inact_orb
integer :: i_hole_osoci
i_hole_osoci = list_inact(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_hole_osoci,thr,test_sym)
if(.not.test_sym)cycle
print*,'i_hole_osoci = ',i_hole_osoci
call create_restart_and_1h(i_hole_osoci)
call set_generators_to_psi_det
print*,'Passed set generators'
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call is_a_good_candidate(threshold,is_ok,verbose)
print*,'is_ok = ',is_ok
if(.not.is_ok)cycle
! so all the mono excitation on the new generators
allocate(dressing_matrix(N_det_generators,N_det_generators))
if(.not.do_it_perturbative)then
! call all_single
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
double precision :: hkl
! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
call debug_det(reunion_of_bitmask,N_int)
call all_single
endif
call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci)
do k = 1, N_states
print*,'norm_tmp = ',norm_tmp(k)
norm_total(k) += norm_tmp(k)
enddo
call update_density_matrix_osoci
deallocate(dressing_matrix)
enddo
if(.True.)then
print*,''
print*,'DOING THEN THE MLCT !!'
do i = 1, n_virt_orb
integer :: i_particl_osoci
i_particl_osoci = list_virt(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_particl_osoci,thr,test_sym)
if(.not.test_sym)cycle
print*,'i_particl_osoci= ',i_particl_osoci
! Initialize the bitmask to the restart ones
call initialize_bitmask_to_restart_ones
! Impose that only the hole i_hole_osoci can be done
call modify_bitmasks_for_particl(i_particl_osoci)
call print_generators_bitmasks_holes
! Impose that only the active part can be reached
call set_bitmask_hole_as_input(unpaired_bitmask)
!! call all_single_h_core
call create_restart_and_1p(i_particl_osoci)
!! ! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
!! ! so all the mono excitation on the new generators
call is_a_good_candidate(threshold,is_ok,verbose)
print*,'is_ok = ',is_ok
if(.not.is_ok)cycle
allocate(dressing_matrix(N_det_generators,N_det_generators))
if(.not.do_it_perturbative)then
dressing_matrix = 0.d0
do k = 1, N_det_generators
do l = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,k),psi_det_generators(1,1,l),N_int,hkl)
dressing_matrix(k,l) = hkl
enddo
enddo
! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
call all_single
endif
call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci)
do k = 1, N_states
print*,'norm_tmp = ',norm_tmp(k)
norm_total(k) += norm_tmp(k)
enddo
call update_density_matrix_osoci
deallocate(dressing_matrix)
enddo
endif
if(.False.)then
print*,'LAST loop for all the 1h-1p'
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call initialize_bitmask_to_restart_ones
! Impose that only the hole i_hole_osoci can be done
call set_bitmask_particl_as_input(inact_virt_bitmask)
call set_bitmask_hole_as_input(inact_virt_bitmask)
! call set_bitmask_particl_as_input(reunion_of_bitmask)
! call set_bitmask_hole_as_input(reunion_of_bitmask)
call all_single
call set_intermediate_normalization_1h1p(norm_tmp)
norm_total += norm_tmp
call update_density_matrix_osoci
endif
print*,'norm_total = ',norm_total
norm_total = norm_generators_restart
norm_total = 1.d0/norm_total
! call rescale_density_matrix_osoci(norm_total)
double precision :: accu
accu = 0.d0
do i = 1, mo_tot_num
accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i)
enddo
print*,'accu = ',accu
end
subroutine FOBOCI_mlct_old
use bitmasks
implicit none
integer :: i,j,k,l
integer(bit_kind),allocatable :: unpaired_bitmask(:,:)
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha, n_occ_beta
double precision :: norm_tmp,norm_total
logical :: test_sym
double precision :: thr
double precision :: threshold
logical :: verbose,is_ok
verbose = .False.
threshold = 1.d-2
thr = 1.d-12
allocate(unpaired_bitmask(N_int,2))
allocate (occ(N_int*bit_kind_size,2))
do i = 1, N_int
unpaired_bitmask(i,1) = unpaired_alpha_electrons(i)
unpaired_bitmask(i,2) = unpaired_alpha_electrons(i)
enddo
norm_total = 0.d0
call initialize_density_matrix_osoci
call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int)
print*,''
print*,''
print*,''
print*,'DOING FIRST MLCT !!'
do i = 1, n_virt_orb
integer :: i_particl_osoci
i_particl_osoci = list_virt(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_particl_osoci,thr,test_sym)
if(.not.test_sym)cycle
print*,'i_particl_osoci= ',i_particl_osoci
! Initialize the bitmask to the restart ones
call initialize_bitmask_to_restart_ones
! Impose that only the hole i_hole_osoci can be done
call modify_bitmasks_for_particl(i_particl_osoci)
call print_generators_bitmasks_holes
! Impose that only the active part can be reached
call set_bitmask_hole_as_input(unpaired_bitmask)
! call all_single_h_core
call create_restart_and_1p(i_particl_osoci)
! ! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
! ! so all the mono excitation on the new generators
call is_a_good_candidate(threshold,is_ok,verbose)
print*,'is_ok = ',is_ok
is_ok =.True.
if(.not.is_ok)cycle
call all_single
call set_intermediate_normalization_mlct_old(norm_tmp,i_particl_osoci)
print*,'norm_tmp = ',norm_tmp
norm_total += norm_tmp
call update_density_matrix_osoci
enddo
print*,'norm_total = ',norm_total
norm_total += 1.d0
norm_total = 1.d0/norm_total
call rescale_density_matrix_osoci(norm_total)
double precision :: accu
accu = 0.d0
do i = 1, mo_tot_num
accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i)
enddo
print*,'accu = ',accu
end
subroutine FOBOCI_lmct_old
use bitmasks
implicit none
integer :: i,j,k,l
integer(bit_kind),allocatable :: unpaired_bitmask(:,:)
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha, n_occ_beta
double precision :: norm_tmp,norm_total
logical :: test_sym
double precision :: thr
double precision :: threshold
logical :: verbose,is_ok
verbose = .False.
threshold = 1.d-2
thr = 1.d-12
allocate(unpaired_bitmask(N_int,2))
allocate (occ(N_int*bit_kind_size,2))
do i = 1, N_int
unpaired_bitmask(i,1) = unpaired_alpha_electrons(i)
unpaired_bitmask(i,2) = unpaired_alpha_electrons(i)
enddo
norm_total = 0.d0
call initialize_density_matrix_osoci
call bitstring_to_list(inact_bitmask(1,1), occ(1,1), n_occ_beta, N_int)
print*,''
print*,''
print*,'DOING FIRST LMCT !!'
do i = 1, n_inact_orb
integer :: i_hole_osoci
i_hole_osoci = list_inact(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_hole_osoci,thr,test_sym)
if(.not.test_sym)cycle
print*,'i_hole_osoci = ',i_hole_osoci
! Initialize the bitmask to the restart ones
call initialize_bitmask_to_restart_ones
! Impose that only the hole i_hole_osoci can be done
call modify_bitmasks_for_hole(i_hole_osoci)
call print_generators_bitmasks_holes
! Impose that only the active part can be reached
call set_bitmask_particl_as_input(unpaired_bitmask)
! call all_single_h_core
call create_restart_and_1h(i_hole_osoci)
! ! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call is_a_good_candidate(threshold,is_ok,verbose)
print*,'is_ok = ',is_ok
if(.not.is_ok)cycle
! ! so all the mono excitation on the new generators
call all_single
! call set_intermediate_normalization_lmct_bis(norm_tmp,i_hole_osoci)
call set_intermediate_normalization_lmct_old(norm_tmp,i_hole_osoci)
print*,'norm_tmp = ',norm_tmp
norm_total += norm_tmp
call update_density_matrix_osoci
enddo
print*,'norm_total = ',norm_total
norm_total += 1.d0
norm_total = 1.d0/norm_total
call rescale_density_matrix_osoci(norm_total)
double precision :: accu
accu = 0.d0
do i = 1, mo_tot_num
accu += one_body_dm_mo_alpha_osoci(i,i) + one_body_dm_mo_beta_osoci(i,i)
enddo
print*,'accu = ',accu
end

126
plugins/FOBOCI/generators_restart_save.irp.f Normal file
View File

@ -0,0 +1,126 @@
use bitmasks
BEGIN_PROVIDER [ integer, N_det_generators_restart ]
implicit none
BEGIN_DOC
! Number of determinants in the wave function
END_DOC
logical :: exists
character*64 :: label
integer, save :: ifirst = 0
!if(ifirst == 0)then
PROVIDE ezfio_filename
call ezfio_has_determinants_n_det(exists)
print*,'exists = ',exists
if(.not.exists)then
print*,'The OSOCI needs a restart WF'
print*,'There are none in the EZFIO file ...'
print*,'Stopping ...'
stop
endif
print*,'passed N_det_generators_restart'
call ezfio_get_determinants_n_det(N_det_generators_restart)
ASSERT (N_det_generators_restart > 0)
ifirst = 1
!endif
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ]
implicit none
BEGIN_DOC
! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
! is empty
END_DOC
integer :: i
logical :: exists
character*64 :: label
integer, save :: ifirst = 0
!if(ifirst == 0)then
provide N_det_generators_restart
if(.True.)then
call ezfio_has_determinants_N_int(exists)
if (exists) then
call ezfio_has_determinants_bit_kind(exists)
if (exists) then
call ezfio_has_determinants_N_det(exists)
if (exists) then
call ezfio_has_determinants_N_states(exists)
if (exists) then
call ezfio_has_determinants_psi_det(exists)
endif
endif
endif
endif
if(.not.exists)then
print*,'The OSOCI needs a restart WF'
print*,'There are none in the EZFIO file ...'
print*,'Stopping ...'
stop
endif
print*,'passed psi_det_generators_restart'
call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart)
do i = 1, N_int
ref_generators_restart(i,1) = psi_det_generators_restart(i,1,1)
ref_generators_restart(i,2) = psi_det_generators_restart(i,2,1)
enddo
endif
ifirst = 1
!endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (psi_det_size,N_states_diag) ]
implicit none
BEGIN_DOC
! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
! is empty
END_DOC
integer :: i,k, N_int2
logical :: exists
double precision, allocatable :: psi_coef_read(:,:)
character*(64) :: label
integer, save :: ifirst = 0
!if(ifirst == 0)then
psi_coef_generators_restart = 0.d0
do i=1,N_states_diag
psi_coef_generators_restart(i,i) = 1.d0
enddo
call ezfio_has_determinants_psi_coef(exists)
if(.not.exists)then
print*,'The OSOCI needs a restart WF'
print*,'There are none in the EZFIO file ...'
print*,'Stopping ...'
stop
endif
print*,'passed psi_coef_generators_restart'
if (exists) then
allocate (psi_coef_read(N_det_generators_restart,N_states))
call ezfio_get_determinants_psi_coef(psi_coef_read)
do k=1,N_states
do i=1,N_det_generators_restart
psi_coef_generators_restart(i,k) = psi_coef_read(i,k)
enddo
enddo
deallocate(psi_coef_read)
endif
ifirst = 1
!endif
END_PROVIDER

157
plugins/FOBOCI/modify_generators.irp.f Normal file
View File

@ -0,0 +1,157 @@
subroutine set_generators_to_psi_det
implicit none
BEGIN_DOC
! subroutines that sets psi_det_generators to
! the current psi_det
END_DOC
N_det_generators = N_det
integer :: i,k
do i=1,N_det_generators
do k=1,N_int
psi_det_generators(k,1,i) = psi_det(k,1,i)
psi_det_generators(k,2,i) = psi_det(k,2,i)
enddo
do k = 1, N_states
psi_coef_generators(i,k) = psi_coef(i,k)
enddo
enddo
touch N_det_generators psi_coef_generators psi_det_generators
end
subroutine set_generators_as_input_psi(ndet_input,psi_det_input,psi_coef_input)
implicit none
integer, intent(in) :: ndet_input
integer(bit_kind), intent(in) :: psi_det_input(N_int,2,ndet_input)
double precision, intent(in) :: psi_coef_input(ndet_input,N_states)
BEGIN_DOC
! subroutines that sets psi_det_generators to
! the current psi_det
END_DOC
N_det_generators = ndet_input
integer :: i,k
do i=1,N_det_generators
do k=1,N_int
psi_det_generators(k,1,i) = psi_det_input(k,1,i)
psi_det_generators(k,2,i) = psi_det_input(k,2,i)
enddo
do k = 1, N_states
psi_coef_generators(i,k) = psi_coef_input(i,k)
enddo
enddo
touch N_det_generators psi_coef_generators psi_det_generators
end
subroutine set_psi_det_as_input_psi(ndet_input,psi_det_input,psi_coef_input)
implicit none
integer, intent(in) :: ndet_input
integer(bit_kind), intent(in) :: psi_det_input(N_int,2,ndet_input)
double precision, intent(in) :: psi_coef_input(ndet_input,N_states)
BEGIN_DOC
! subroutines that sets psi_det_generators to
! the current psi_det
END_DOC
N_det= ndet_input
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
integer :: i,k
do i=1,N_det
do k=1,N_int
psi_det(k,1,i) = psi_det_input(k,1,i)
psi_det(k,2,i) = psi_det_input(k,2,i)
enddo
do k = 1, N_states
psi_coef(i,k) = psi_coef_input(i,k)
enddo
enddo
soft_touch N_det psi_coef psi_det
end
subroutine set_psi_det_to_generators
implicit none
BEGIN_DOC
! subroutines that sets psi_det_generators to
! the current psi_det
END_DOC
N_det= N_det_generators
integer :: i,k
do i = 1, psi_det_size
do k=1,N_int
psi_det(k,1,i) = 0_bit_kind
psi_det(k,2,i) = 0_bit_kind
enddo
do k = 1, N_states
psi_coef(i,k) = 0.d0
enddo
enddo
do i=1,N_det_generators
do k=1,N_int
psi_det(k,1,i) = psi_det_generators(k,1,i)
psi_det(k,2,i) = psi_det_generators(k,2,i)
enddo
do k = 1, N_states
psi_coef(i,k) = psi_coef_generators(i,k)
enddo
enddo
touch N_det psi_coef psi_det
end
subroutine set_generators_to_generators_restart
implicit none
BEGIN_DOC
! subroutines that sets psi_det_generators to
! the current psi_det
END_DOC
N_det_generators = N_det_generators_restart
integer :: i,k
do i=1,N_det_generators
do k=1,N_int
psi_det_generators(k,1,i) = psi_det_generators_restart(k,1,i)
psi_det_generators(k,2,i) = psi_det_generators_restart(k,2,i)
enddo
do k = 1, N_states
psi_coef_generators(i,k) = psi_coef_generators_restart(i,k)
enddo
enddo
touch N_det_generators psi_coef_generators psi_det_generators
end
subroutine set_psi_det_to_generators_restart
implicit none
BEGIN_DOC
! subroutines that sets psi_det_generators to
! the current psi_det
END_DOC
N_det = N_det_generators_restart
integer :: i,k
do i=1,N_det_generators
do k=1,N_int
psi_det(k,1,i) = psi_det_generators_restart(k,1,i)
psi_det(k,2,i) = psi_det_generators_restart(k,2,i)
enddo
do k = 1, N_states
psi_coef(i,k) = psi_coef_generators_restart(i,k)
enddo
enddo
touch N_det psi_coef psi_det
end

414
plugins/FOBOCI/new_approach.irp.f Normal file
View File

@ -0,0 +1,414 @@
subroutine new_approach
use bitmasks
implicit none
integer :: n_max_good_det
n_max_good_det = n_inact_orb * n_act_orb *n_det_generators_restart + n_virt_orb * n_act_orb * n_det_generators_restart
integer :: n_good_det,n_good_hole, n_good_particl
n_good_det = 0
n_good_hole = 0
n_good_particl = 0
integer(bit_kind), allocatable :: psi_good_det(:,:,:)
double precision, allocatable :: dressing_restart_good_det(:,:)
double precision, allocatable :: dressing_matrix_restart_1h1p(:,:)
double precision, allocatable :: dressing_matrix_restart_2h1p(:,:)
double precision, allocatable :: dressing_matrix_restart_1h2p(:,:)
double precision, allocatable :: dressing_diag_good_det(:)
double precision :: hjk
integer :: i,j,k,l,i_hole_foboci
logical :: test_sym
double precision :: thr,hij
double precision :: threshold,accu
double precision, allocatable :: dressing_matrix_1h1p(:,:)
double precision, allocatable :: dressing_matrix_2h1p(:,:)
double precision, allocatable :: dressing_matrix_1h2p(:,:)
double precision, allocatable :: H_matrix_tmp(:,:)
logical :: verbose,is_ok
double precision,allocatable :: eigenvectors(:,:), eigenvalues(:)
allocate(psi_good_det(N_int,2,n_max_good_det))
allocate(dressing_restart_good_det(n_max_good_det,n_det_generators_restart))
allocate(dressing_matrix_restart_1h1p(N_det_generators_restart, N_det_generators_restart))
allocate(dressing_matrix_restart_2h1p(N_det_generators_restart, N_det_generators_restart))
allocate(dressing_matrix_restart_1h2p(N_det_generators_restart, N_det_generators_restart))
allocate(dressing_diag_good_det(n_max_good_det))
dressing_restart_good_det = 0.d0
dressing_matrix_restart_1h1p = 0.d0
dressing_matrix_restart_2h1p = 0.d0
dressing_matrix_restart_1h2p = 0.d0
dressing_diag_good_det = 0.d0
verbose = .True.
threshold = threshold_singles
print*,'threshold = ',threshold
thr = 1.d-12
print*,''
print*,''
print*,'mulliken spin population analysis'
accu =0.d0
do i = 1, nucl_num
accu += mulliken_spin_densities(i)
print*,i,nucl_charge(i),mulliken_spin_densities(i)
enddo
print*,''
print*,''
print*,'DOING FIRST LMCT !!'
integer :: i_particl_osoci
do i = 1, n_inact_orb
i_hole_foboci = list_inact(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_hole_foboci,thr,test_sym)
if(.not.test_sym)cycle
print*,'i_hole_foboci = ',i_hole_foboci
call create_restart_and_1h(i_hole_foboci)
! ! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call is_a_good_candidate(threshold,is_ok,verbose)
print*,'is_ok = ',is_ok
if(.not.is_ok)cycle
! so all the mono excitation on the new generators
allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators))
allocate(dressing_matrix_2h1p(N_det_generators,N_det_generators))
dressing_matrix_1h1p = 0.d0
dressing_matrix_2h1p = 0.d0
if(.not.do_it_perturbative)then
n_good_hole +=1
! call all_single_split_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
call all_single_for_1h(dressing_matrix_1h1p,dressing_matrix_2h1p)
allocate(H_matrix_tmp(N_det_generators,N_det_generators))
do j = 1,N_det_generators
do k = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,j),psi_det_generators(1,1,k),N_int,hjk)
H_matrix_tmp(j,k) = hjk
enddo
enddo
do j = 1, N_det_generators
do k = 1, N_det_generators
H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_2h1p(j,k)
enddo
enddo
hjk = H_matrix_tmp(1,1)
do j = 1, N_det_generators
H_matrix_tmp(j,j) -= hjk
enddo
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'Dressed matrix :'
do j = 1, N_det_generators
write(*,'(100(X,F8.5))') H_matrix_tmp(j,:)
enddo
allocate(eigenvectors(N_det_generators,N_det_generators), eigenvalues(N_det_generators))
call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,N_det_generators,N_det_generators)
print*,'Eigenvector of the dressed matrix :'
do j = 1, N_det_generators
print*,'coef = ',eigenvectors(j,1)
enddo
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
deallocate(eigenvectors, eigenvalues)
deallocate(H_matrix_tmp)
call update_dressing_matrix(dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_restart_good_det,dressing_matrix_restart_1h1p, &
dressing_matrix_restart_2h1p,dressing_diag_good_det,psi_good_det,n_good_det,n_max_good_det)
endif
deallocate(dressing_matrix_1h1p)
deallocate(dressing_matrix_2h1p)
enddo
print*,''
print*,''
print*,'DOING THEN THE MLCT !!'
do i = 1, n_virt_orb
i_particl_osoci = list_virt(i)
print*,'--------------------------'
! First set the current generators to the one of restart
call set_generators_to_generators_restart
call set_psi_det_to_generators
call check_symetry(i_particl_osoci,thr,test_sym)
if(.not.test_sym)cycle
print*,'i_part_foboci = ',i_particl_osoci
call create_restart_and_1p(i_particl_osoci)
! Update the generators
call set_generators_to_psi_det
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call is_a_good_candidate(threshold,is_ok,verbose)
print*,'is_ok = ',is_ok
if(.not.is_ok)cycle
! so all the mono excitation on the new generators
allocate(dressing_matrix_1h1p(N_det_generators,N_det_generators))
allocate(dressing_matrix_1h2p(N_det_generators,N_det_generators))
dressing_matrix_1h1p = 0.d0
dressing_matrix_1h2p = 0.d0
if(.not.do_it_perturbative)then
n_good_hole +=1
! call all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
call all_single_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
allocate(H_matrix_tmp(N_det_generators,N_det_generators))
do j = 1,N_det_generators
do k = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,j),psi_det_generators(1,1,k),N_int,hjk)
H_matrix_tmp(j,k) = hjk
enddo
enddo
do j = 1, N_det_generators
do k = 1, N_det_generators
H_matrix_tmp(j,k) += dressing_matrix_1h1p(j,k) + dressing_matrix_1h2p(j,k)
enddo
enddo
hjk = H_matrix_tmp(1,1)
do j = 1, N_det_generators
H_matrix_tmp(j,j) -= hjk
enddo
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'Dressed matrix :'
do j = 1, N_det_generators
write(*,'(100(F8.5))') H_matrix_tmp(j,:)
enddo
allocate(eigenvectors(N_det_generators,N_det_generators), eigenvalues(N_det_generators))
call lapack_diag(eigenvalues,eigenvectors,H_matrix_tmp,N_det_generators,N_det_generators)
print*,'Eigenvector of the dressed matrix :'
do j = 1, N_det_generators
print*,'coef = ',eigenvectors(j,1)
enddo
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
print*,'-----------------------'
deallocate(eigenvectors, eigenvalues)
deallocate(H_matrix_tmp)
call update_dressing_matrix(dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_restart_good_det,dressing_matrix_restart_1h1p, &
dressing_matrix_restart_1h2p,dressing_diag_good_det,psi_good_det,n_good_det,n_max_good_det)
endif
deallocate(dressing_matrix_1h1p)
deallocate(dressing_matrix_1h2p)
enddo
double precision, allocatable :: H_matrix_total(:,:)
integer :: n_det_total
n_det_total = N_det_generators_restart + n_good_det
allocate(H_matrix_total(n_det_total, n_det_total))
! Building of the effective Hamiltonian
! We assume that the first determinants are the n_det_generators_restart ones
! and then come the n_good_det determinants in psi_good_det
H_matrix_total = 0.d0
do i = 1, N_det_generators_restart
do j = 1, N_det_generators_restart
call i_H_j(psi_det_generators_restart(1,1,i),psi_det_generators_restart(1,1,j),N_int,hij)
H_matrix_total(i,j) = hij
!!! Adding the averaged dressing coming from the 1h1p that are redundant for each of the "n_good_hole" 1h
H_matrix_total(i,j) += dressing_matrix_restart_1h1p(i,j)/dble(n_good_hole+n_good_particl)
!!! Adding the dressing coming from the 2h1p that are not redundant for the any of CI calculations
H_matrix_total(i,j) += dressing_matrix_restart_2h1p(i,j)
enddo
enddo
do i = 1, n_good_det
call i_H_j(psi_good_det(1,1,i),psi_good_det(1,1,i),N_int,hij)
!!! Adding the diagonal dressing coming from the singles
H_matrix_total(n_det_generators_restart+i,n_det_generators_restart+i) = hij + dressing_diag_good_det(i)
do j = 1, N_det_generators_restart
!!! Adding the extra diagonal dressing between the references and the singles
print*,' dressing_restart_good_det = ',dressing_restart_good_det(i,j)
call i_H_j(psi_good_det(1,1,i),psi_det_generators_restart(1,1,j),N_int,hij)
H_matrix_total(n_det_generators_restart+i,j) += hij
H_matrix_total(j,n_det_generators_restart+i) += hij
H_matrix_total(j,n_det_generators_restart+i) += dressing_restart_good_det(i,j)
H_matrix_total(n_det_generators_restart+i,j) += dressing_restart_good_det(i,j)
enddo
do j = i+1, n_good_det
!!! Building the naked Hamiltonian matrix between the singles
call i_H_j(psi_good_det(1,1,i),psi_good_det(1,1,j),N_int,hij)
H_matrix_total(n_det_generators_restart+i,n_det_generators_restart+j) = hij
H_matrix_total(n_det_generators_restart+j,n_det_generators_restart+i) = hij
enddo
enddo
print*,'H matrix to diagonalize'
double precision :: href
href = H_matrix_total(1,1)
do i = 1, n_det_total
H_matrix_total(i,i) -= href
enddo
do i = 1, n_det_total
write(*,'(100(X,F16.8))')H_matrix_total(i,:)
enddo
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
allocate(eigvalues(n_det_total),eigvectors(n_det_total,n_det_total))
call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total)
print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion + href
do i = 1, n_det_total
print*,'coef = ',eigvectors(i,1)
enddo
integer(bit_kind), allocatable :: psi_det_final(:,:,:)
double precision, allocatable :: psi_coef_final(:,:)
double precision :: norm
allocate(psi_coef_final(n_det_total, N_states))
allocate(psi_det_final(N_int,2,n_det_total))
do i = 1, N_det_generators_restart
do j = 1,N_int
psi_det_final(j,1,i) = psi_det_generators_restart(j,1,i)
psi_det_final(j,2,i) = psi_det_generators_restart(j,2,i)
enddo
enddo
do i = 1, n_good_det
do j = 1,N_int
psi_det_final(j,1,n_det_generators_restart+i) = psi_good_det(j,1,i)
psi_det_final(j,2,n_det_generators_restart+i) = psi_good_det(j,2,i)
enddo
enddo
norm = 0.d0
do i = 1, n_det_total
do j = 1, N_states
psi_coef_final(i,j) = eigvectors(i,j)
enddo
norm += psi_coef_final(i,1)**2
! call debug_det(psi_det_final(1, 1, i), N_int)
enddo
print*,'norm = ',norm
call set_psi_det_as_input_psi(n_det_total,psi_det_final,psi_coef_final)
print*,''
!do i = 1, N_det
! call debug_det(psi_det(1,1,i),N_int)
! print*,'coef = ',psi_coef(i,1)
!enddo
provide one_body_dm_mo
integer :: i_core,iorb,jorb,i_inact,j_inact,i_virt,j_virt,j_core
do i = 1, n_core_orb
i_core = list_core(i)
one_body_dm_mo(i_core,i_core) = 10.d0
do j = i+1, n_core_orb
j_core = list_core(j)
one_body_dm_mo(i_core,j_core) = 0.d0
one_body_dm_mo(j_core,i_core) = 0.d0
enddo
do j = 1, n_inact_orb
iorb = list_inact(j)
one_body_dm_mo(i_core,iorb) = 0.d0
one_body_dm_mo(iorb,i_core) = 0.d0
enddo
do j = 1, n_act_orb
iorb = list_act(j)
one_body_dm_mo(i_core,iorb) = 0.d0
one_body_dm_mo(iorb,i_core) = 0.d0
enddo
do j = 1, n_virt_orb
iorb = list_virt(j)
one_body_dm_mo(i_core,iorb) = 0.d0
one_body_dm_mo(iorb,i_core) = 0.d0
enddo
enddo
! Set to Zero the inact-inact part to avoid arbitrary rotations
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = i+1, n_inact_orb
j_inact = list_inact(j)
one_body_dm_mo(i_inact,j_inact) = 0.d0
one_body_dm_mo(j_inact,i_inact) = 0.d0
enddo
enddo
! Set to Zero the inact-virt part to avoid arbitrary rotations
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = 1, n_virt_orb
j_virt = list_virt(j)
one_body_dm_mo(i_inact,j_virt) = 0.d0
one_body_dm_mo(j_virt,i_inact) = 0.d0
enddo
enddo
! Set to Zero the virt-virt part to avoid arbitrary rotations
do i = 1, n_virt_orb
i_virt = list_virt(i)
do j = i+1, n_virt_orb
j_virt = list_virt(j)
one_body_dm_mo(i_virt,j_virt) = 0.d0
one_body_dm_mo(j_virt,i_virt) = 0.d0
enddo
enddo
print*,''
print*,'Inactive-active Part of the One body DM'
print*,''
do i = 1,n_act_orb
iorb = list_act(i)
print*,''
print*,'ACTIVE ORBITAL ',iorb
do j = 1, n_inact_orb
jorb = list_inact(j)
if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then
print*,'INACTIVE '
print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb))
endif
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
if(dabs(one_body_dm_mo(iorb,jorb)).gt.threshold_singles)then
print*,'VIRT '
print*,'DM ',iorb,jorb,dabs(one_body_dm_mo(iorb,jorb))
endif
enddo
enddo
do i = 1, mo_tot_num
do j = i+1, mo_tot_num
if(dabs(one_body_dm_mo(i,j)).le.threshold_fobo_dm)then
one_body_dm_mo(i,j) = 0.d0
one_body_dm_mo(j,i) = 0.d0
endif
enddo
enddo
label = "Natural"
character*(64) :: label
integer :: sign
sign = -1
call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign)
soft_touch mo_coef
call save_mos
deallocate(eigvalues,eigvectors,psi_det_final,psi_coef_final)
deallocate(H_matrix_total)
deallocate(psi_good_det)
deallocate(dressing_restart_good_det)
deallocate(dressing_matrix_restart_1h1p)
deallocate(dressing_matrix_restart_2h1p)
deallocate(dressing_diag_good_det)
print*,'YOUPI C EST FINI !!'
end

56
plugins/FOBOCI/routine_new_approach.irp.f Normal file
View File

@ -0,0 +1,56 @@
subroutine update_dressing_matrix(dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_restart_good_det,dressing_matrix_restart_1h1p, &
dressing_matrix_restart_2h1p,dressing_diag_good_det,psi_good_det,n_good_det,n_max_good_det)
implicit none
integer, intent(in) :: n_max_good_det
integer, intent(inout) :: n_good_det
integer(bit_kind), intent(inout) :: psi_good_det(N_int,2,n_max_good_det)
double precision, intent(in) :: dressing_matrix_1h1p(N_det_generators,N_det_generators)
double precision, intent(in) :: dressing_matrix_2h1p(N_det_generators,N_det_generators)
double precision, intent(inout) :: dressing_matrix_restart_1h1p(N_det_generators_restart, N_det_generators_restart)
double precision, intent(inout) :: dressing_matrix_restart_2h1p(N_det_generators_restart, N_det_generators_restart)
double precision, intent(inout) :: dressing_restart_good_det(n_max_good_det,N_det_generators_restart)
double precision, intent(inout) :: dressing_diag_good_det(n_max_good_det)
use bitmasks
integer :: k,l,degree
logical, allocatable :: is_a_ref_det(:)
integer, allocatable :: index_restart_generators(:)
allocate(is_a_ref_det(N_det_generators),index_restart_generators(N_det_generators))
is_a_ref_det = .False.
do k = 1, N_det_generators
do l = 1, N_det_generators_restart
call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators_restart(1,1,l), degree, N_int)
if(degree==0)then
is_a_ref_det(k) = .True.
index_restart_generators(k) = l
endif
enddo
enddo
do k = 1, N_det_generators
if(is_a_ref_det(k))then
do l = 1, N_det_generators
if(.not.is_a_ref_det(l))cycle
!!!! Dressing of the reference space in the order of the restart determinants
dressing_matrix_restart_1h1p(index_restart_generators(l),index_restart_generators(k)) += dressing_matrix_1h1p(k,l)
print*,' dressing_matrix_1h1p(k,l) = ',dressing_matrix_1h1p(k,l)
dressing_matrix_restart_2h1p(index_restart_generators(l),index_restart_generators(k)) += dressing_matrix_2h1p(k,l)
enddo
else
!!!! Incrementing the counting of the good determinants
n_good_det +=1
!!!! Adding the good determinant to the global_list (psi_good_det)
do l = 1, N_int
psi_good_det(l,1,n_good_det) = psi_det_generators(l,1,k)
psi_good_det(l,2,n_good_det) = psi_det_generators(l,2,k)
enddo
!!! Storing the diagonal dressing of the good det
dressing_diag_good_det(n_good_det) = dressing_matrix_1h1p(k,k) + dressing_matrix_2h1p(k,k)
do l = 1, N_det_generators
if(.not.is_a_ref_det(l))cycle
!!! Storing the extra diagonal dressing of the good det with the restart determinants
dressing_restart_good_det(n_good_det,index_restart_generators(l)) = dressing_matrix_1h1p(k,l) + dressing_matrix_2h1p(k,l)
enddo
endif
enddo
deallocate(is_a_ref_det,index_restart_generators)
end

457
plugins/FOBOCI/routines_dressing.irp.f Normal file
View File

@ -0,0 +1,457 @@
subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
use bitmasks
implicit none
integer, intent(in) :: ndet_generators_input
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input)
double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
double precision :: H_array(N_det),hka
logical :: is_a_ref_det(N_det)
integer :: i,j,n_det_ref_tmp
integer :: connected_to_ref_by_mono,degree
double precision :: coef_ref(Ndet_generators_input)
double precision :: accu,lambda_i
integer :: k
integer :: index_ref_tmp(N_det)
is_a_ref_det = .False.
n_det_ref_tmp = 0
do i = 1, N_det
do j = 1, Ndet_generators_input
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_input(1,1,j),degree,N_int)
if(degree == 0)then
is_a_ref_det(i) = .True.
n_det_ref_tmp +=1
index_ref_tmp(n_det_ref_tmp) = i
coef_ref(n_det_ref_tmp) = psi_coef(i,1)
exit
endif
enddo
enddo
if( ndet_generators_input .ne. n_det_ref_tmp)then
print*,'Problem !!!! '
print*,' ndet_generators .ne. n_det_ref_tmp !!!'
print*,'ndet_generators,n_det_ref_tmp'
print*,ndet_generators_input,n_det_ref_tmp
stop
endif
call i_h_j(psi_det_generators_input(1,1,1),psi_det_generators_input(1,1,1),N_int,href)
integer :: i_pert, i_pert_count
i_pert_count = 0
do i = 1, N_det
if(is_a_ref_det(i))cycle
call i_h_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hka)
double precision :: f,href
f = 1.d0/(href - hka)
H_array = 0.d0
accu = 0.d0
do j=1,ndet_generators_input
call i_h_j(psi_det(1,1,i),psi_det_generators_input(1,1,j),N_int,hka)
H_array(j) = hka
accu += coef_ref(j) * hka
enddo
lambda_i = psi_coef(i,1)/accu
i_pert = 1
if(accu * f / psi_coef(i,1) .gt. 0.5d0 .and. accu * f/psi_coef(i,1).gt.0.d0)then
i_pert = 0
endif
do j = 1, ndet_generators_input
if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
i_pert = 1
exit
endif
enddo
! print*,''
! print*,'lambda_i,f = ',lambda_i,f
! print*,'i_pert = ',i_pert
! print*,''
if(i_pert==1)then
lambda_i = f
i_pert_count +=1
endif
do k=1,ndet_generators_input
double precision :: contrib
contrib = H_array(k) * H_array(k) * lambda_i
dressing_matrix(k, k) += contrib
do j=k+1,ndet_generators_input
contrib = H_array(k) * H_array(j) * lambda_i
dressing_matrix(k, j) += contrib
dressing_matrix(j, k) += contrib
enddo
enddo
enddo
!print*,'i_pert_count = ',i_pert_count
end
subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi_coef_ref_input,n_det_ref_input, &
psi_det_outer_input,psi_coef_outer_input,n_det_outer_input)
use bitmasks
implicit none
integer, intent(in) :: n_det_ref_input
integer(bit_kind), intent(in) :: psi_det_ref_input(N_int,2,n_det_ref_input)
double precision, intent(in) :: psi_coef_ref_input(n_det_ref_input,N_states)
integer, intent(in) :: n_det_outer_input
integer(bit_kind), intent(in) :: psi_det_outer_input(N_int,2,n_det_outer_input)
double precision, intent(in) :: psi_coef_outer_input(n_det_outer_input,N_states)
double precision, intent(inout) :: dressing_matrix(n_det_ref_input,n_det_ref_input)
integer :: i_pert, i_pert_count,i,j,k
double precision :: f,href,hka,lambda_i
double precision :: H_array(n_det_ref_input),accu
call i_h_j(psi_det_ref_input(1,1,1),psi_det_ref_input(1,1,1),N_int,href)
i_pert_count = 0
do i = 1, n_det_outer_input
call i_h_j(psi_det_outer_input(1,1,i),psi_det_outer_input(1,1,i),N_int,hka)
f = 1.d0/(href - hka)
H_array = 0.d0
accu = 0.d0
do j=1,n_det_ref_input
call i_h_j(psi_det_outer_input(1,1,i),psi_det_ref_input(1,1,j),N_int,hka)
H_array(j) = hka
accu += psi_coef_ref_input(j,1) * hka
enddo
lambda_i = psi_coef_outer_input(i,1)/accu
i_pert = 1
if(accu * f / psi_coef_outer_input(i,1) .gt. 0.5d0 .and. accu * f/psi_coef_outer_input(i,1).gt.0.d0)then
i_pert = 0
endif
do j = 1, n_det_ref_input
if(dabs(H_array(j)*lambda_i).gt.0.3d0)then
i_pert = 1
exit
endif
enddo
if(i_pert==1)then
lambda_i = f
i_pert_count +=1
endif
do k=1,n_det_ref_input
double precision :: contrib
contrib = H_array(k) * H_array(k) * lambda_i
dressing_matrix(k, k) += contrib
do j=k+1,n_det_ref_input
contrib = H_array(k) * H_array(j) * lambda_i
dressing_matrix(k, j) += contrib
dressing_matrix(j, k) += contrib
enddo
enddo
enddo
end
subroutine diag_dressed_matrix_and_set_to_psi_det(psi_det_generators_input,Ndet_generators_input,dressing_matrix)
use bitmasks
implicit none
integer, intent(in) :: ndet_generators_input
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,ndet_generators_input)
double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
integer :: i,j
double precision :: eigenvectors(Ndet_generators_input,Ndet_generators_input), eigenvalues(Ndet_generators_input)
call lapack_diag(eigenvalues,eigenvectors,dressing_matrix,Ndet_generators_input,Ndet_generators_input)
print*,'Dressed eigenvalue, to be compared with the CI one'
print*,'E = ',eigenvalues(1)+nuclear_repulsion
print*,'Dressed matrix, to be compared to the intermediate Hamiltonian one'
do i = 1, Ndet_generators_input
write(*,'(100(F12.5,X))')dressing_matrix(i,:)
enddo
n_det = Ndet_generators_input
do i = 1, Ndet_generators_input
psi_coef(i,1) = eigenvectors(i,1)
do j = 1, N_int
psi_det(j,1,i) = psi_det_generators_input(j,1,i)
psi_det(j,2,i) = psi_det_generators_input(j,2,i)
enddo
enddo
touch N_det psi_coef psi_det
end
subroutine give_n_1h1p_and_n_2h1p_in_psi_det(n_det_1h1p,n_det_2h1p)
use bitmasks
implicit none
integer, intent(out) :: n_det_1h1p, n_det_2h1p
integer :: i
integer :: n_det_ref_restart_tmp,n_det_1h
integer :: number_of_holes,n_h, number_of_particles,n_p
n_det_ref_restart_tmp = 0
n_det_1h = 0
n_det_1h1p = 0
n_det_2h1p = 0
do i = 1, N_det
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if(n_h == 0 .and. n_p == 0)then
n_det_ref_restart_tmp +=1
else if (n_h ==1 .and. n_p==0)then
n_det_1h +=1
else if (n_h ==1 .and. n_p==1)then
n_det_1h1p +=1
else if (n_h ==2 .and. n_p==1)then
n_det_2h1p +=1
else
print*,'PB !!!!'
print*,'You have something else than a 1h, 1h1p or 2h1p'
call debug_det(psi_det(1,1,i),N_int)
stop
endif
enddo
! if(n_det_1h.ne.1)then
! print*,'PB !! You have more than one 1h'
! stop
! endif
if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then
print*,'PB !!!!'
print*,'You have forgotten something in your generators ... '
stop
endif
end
subroutine give_n_1h1p_and_n_1h2p_in_psi_det(n_det_1h1p,n_det_1h2p)
use bitmasks
implicit none
integer, intent(out) :: n_det_1h1p, n_det_1h2p
integer :: i
integer :: n_det_ref_restart_tmp,n_det_1h
integer :: number_of_holes,n_h, number_of_particles,n_p
n_det_ref_restart_tmp = 0
n_det_1h = 0
n_det_1h1p = 0
n_det_1h2p = 0
do i = 1, N_det
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if(n_h == 0 .and. n_p == 0)then
n_det_ref_restart_tmp +=1
else if (n_h ==0 .and. n_p==1)then
n_det_1h +=1
else if (n_h ==1 .and. n_p==1)then
n_det_1h1p +=1
else if (n_h ==1 .and. n_p==2)then
n_det_1h2p +=1
else
print*,'PB !!!!'
print*,'You have something else than a 1p, 1h1p or 1h2p'
call debug_det(psi_det(1,1,i),N_int)
stop
endif
enddo
if(n_det_ref_restart_tmp + n_det_1h .ne. n_det_generators)then
print*,'PB !!!!'
print*,'You have forgotten something in your generators ... '
stop
endif
end
subroutine split_wf_generators_and_1h1p_and_2h1p(n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p)
use bitmasks
implicit none
integer, intent(in) :: n_det_1h1p,n_det_2h1p
integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators)
integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p)
integer(bit_kind), intent(out) :: psi_2h1p(N_int,2,n_det_2h1p)
double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states)
double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states)
double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p, N_states)
integer :: i,j
integer :: degree
integer :: number_of_holes,n_h, number_of_particles,n_p
integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_2h1p_tmp
integer, allocatable :: index_generator(:)
integer, allocatable :: index_1h1p(:)
integer, allocatable :: index_2h1p(:)
allocate(index_1h1p(n_det))
allocate(index_2h1p(n_det))
allocate(index_generator(N_det))
n_det_generators_tmp = 0
n_det_1h1p_tmp = 0
n_det_2h1p_tmp = 0
do i = 1, n_det
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if (n_h ==1 .and. n_p==1)then
n_det_1h1p_tmp +=1
index_1h1p(n_det_1h1p_tmp) = i
else if (n_h ==2 .and. n_p==1)then
n_det_2h1p_tmp +=1
index_2h1p(n_det_2h1p_tmp) = i
endif
do j = 1, N_det_generators
call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int)
if(degree == 0)then
n_det_generators_tmp +=1
index_generator(n_det_generators_tmp) = i
endif
enddo
enddo
if(n_det_1h1p_tmp.ne.n_det_1h1p)then
print*,'PB !!!'
print*,'n_det_1h1p_tmp.ne.n_det_1h1p)'
stop
endif
if(n_det_2h1p_tmp.ne.n_det_2h1p)then
print*,'PB !!!'
print*,'n_det_2h1p_tmp.ne.n_det_2h1p)'
stop
endif
if(N_det_generators_tmp.ne.n_det_generators)then
print*,'PB !!!'
print*,'N_det_generators_tmp.ne.n_det_generators'
stop
endif
do i = 1,N_det_generators
do j = 1, N_int
psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i))
psi_ref_out(j,2,i) = psi_det(j,2,index_generator(i))
enddo
do j = 1, N_states
psi_ref_coef_out(i,j) = psi_coef(index_generator(i),j)
enddo
enddo
do i = 1, n_det_1h1p
do j = 1, N_int
psi_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
psi_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
enddo
do j = 1, N_states
psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
enddo
enddo
do i = 1, n_det_2h1p
do j = 1, N_int
psi_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i))
psi_2h1p(j,2,i) = psi_det(j,2,index_2h1p(i))
enddo
do j = 1, N_states
psi_coef_2h1p(i,j) = psi_coef(index_2h1p(i),j)
enddo
enddo
deallocate(index_generator)
deallocate(index_1h1p)
deallocate(index_2h1p)
end
subroutine split_wf_generators_and_1h1p_and_1h2p(n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p)
use bitmasks
implicit none
integer, intent(in) :: n_det_1h1p,n_det_1h2p
integer(bit_kind), intent(out) :: psi_ref_out(N_int,2,N_det_generators)
integer(bit_kind), intent(out) :: psi_1h1p(N_int,2,n_det_1h1p)
integer(bit_kind), intent(out) :: psi_1h2p(N_int,2,n_det_1h2p)
double precision, intent(out) :: psi_ref_coef_out(N_det_generators,N_states)
double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p, N_states)
double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p, N_states)
integer :: i,j
integer :: degree
integer :: number_of_holes,n_h, number_of_particles,n_p
integer :: n_det_generators_tmp,n_det_1h1p_tmp,n_det_1h2p_tmp
integer, allocatable :: index_generator(:)
integer, allocatable :: index_1h1p(:)
integer, allocatable :: index_1h2p(:)
allocate(index_1h1p(n_det))
allocate(index_1h2p(n_det))
allocate(index_generator(N_det))
n_det_generators_tmp = 0
n_det_1h1p_tmp = 0
n_det_1h2p_tmp = 0
do i = 1, n_det
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if (n_h ==1 .and. n_p==1)then
n_det_1h1p_tmp +=1
index_1h1p(n_det_1h1p_tmp) = i
else if (n_h ==1 .and. n_p==2)then
n_det_1h2p_tmp +=1
index_1h2p(n_det_1h2p_tmp) = i
endif
do j = 1, N_det_generators
call get_excitation_degree(psi_det_generators(1,1,j),psi_det(1,1,i), degree, N_int)
if(degree == 0)then
n_det_generators_tmp +=1
index_generator(n_det_generators_tmp) = i
endif
enddo
enddo
if(n_det_1h1p_tmp.ne.n_det_1h1p)then
print*,'PB !!!'
print*,'n_det_1h1p_tmp.ne.n_det_1h1p)'
stop
endif
if(n_det_1h2p_tmp.ne.n_det_1h2p)then
print*,'PB !!!'
print*,'n_det_1h2p_tmp.ne.n_det_1h2p)'
stop
endif
if(N_det_generators_tmp.ne.n_det_generators)then
print*,'PB !!!'
print*,'N_det_generators_tmp.ne.n_det_generators'
stop
endif
do i = 1,N_det_generators
do j = 1, N_int
psi_ref_out(j,1,i) = psi_det(j,1,index_generator(i))
psi_ref_out(j,2,i) = psi_det(j,2,index_generator(i))
enddo
do j = 1, N_states
psi_ref_coef_out(i,j) = psi_coef(index_generator(i),j)
enddo
enddo
do i = 1, n_det_1h1p
do j = 1, N_int
psi_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
psi_1h1p(j,2,i) = psi_det(j,2,index_1h1p(i))
enddo
do j = 1, N_states
psi_coef_1h1p(i,j) = psi_coef(index_1h1p(i),j)
enddo
enddo
do i = 1, n_det_1h2p
do j = 1, N_int
psi_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i))
psi_1h2p(j,2,i) = psi_det(j,2,index_1h2p(i))
enddo
do j = 1, N_states
psi_coef_1h2p(i,j) = psi_coef(index_1h2p(i),j)
enddo
enddo
deallocate(index_generator)
deallocate(index_1h1p)
deallocate(index_1h2p)
end

616
plugins/FOBOCI/routines_foboci.irp.f Normal file
View File

@ -0,0 +1,616 @@
subroutine set_intermediate_normalization_lmct_old(norm,i_hole)
implicit none
integer, intent(in) :: i_hole
double precision, intent(out) :: norm(N_states)
integer :: i,j,degree,index_ref_generators_restart,k
integer:: number_of_holes,n_h, number_of_particles,n_p
integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:)
integer, allocatable :: index_one_p(:)
integer :: n_one_hole,n_one_hole_one_p,n_two_hole_one_p,n_two_hole,n_one_p
logical :: is_the_hole_in_det
double precision :: inv_coef_ref_generators_restart(N_states),hij,hii,accu
integer :: index_good_hole(1000)
integer :: n_good_hole
logical,allocatable :: is_a_ref_det(:)
allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det))
n_one_hole = 0
n_one_hole_one_p = 0
n_two_hole_one_p = 0
n_two_hole = 0
n_one_p = 0
n_good_hole = 0
! Find the one holes and one hole one particle
is_a_ref_det = .False.
do i = 1, N_det
! Find the reference determinant for intermediate normalization
call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
do k = 1, N_states
inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k)
enddo
! cycle
endif
! Find all the determinants present in the reference wave function
do j = 1, N_det_generators_restart
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int)
if(degree == 0)then
is_a_ref_det(i) = .True.
exit
endif
enddo
if(is_a_ref_det(i))cycle
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if(n_h == 1 .and. n_p == 0)then
if(is_the_hole_in_det(psi_det(1,1,i),1,i_hole).or.is_the_hole_in_det(psi_det(1,1,i),2,i_hole))then
n_good_hole +=1
index_good_hole(n_good_hole) = i
else
do k = 1, N_states
psi_coef(i,k) = 0.d0
enddo
endif
else
do k = 1, N_states
psi_coef(i,k) = 0.d0
enddo
endif
enddo
!do k = 1, N_det
! call debug_det(psi_det(1,1,k),N_int)
! print*,'k,coef = ',k,psi_coef(k,1)/psi_coef(index_ref_generators_restart,1)
!enddo
print*,''
print*,'n_good_hole = ',n_good_hole
do k = 1,N_states
print*,'state ',k
do i = 1, n_good_hole
print*,'psi_coef(index_good_hole) = ',psi_coef(index_good_hole(i),k)/psi_coef(index_ref_generators_restart,k)
enddo
print*,''
enddo
norm = 0.d0
! Set the wave function to the intermediate normalization
do k = 1, N_states
do i = 1, N_det
psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k)
enddo
enddo
do k = 1,N_states
print*,'state ',k
do i = 1, N_det
!! print*,'psi_coef(i_ref) = ',psi_coef(i,1)
if (is_a_ref_det(i))then
print*,'i,psi_coef_ref = ',psi_coef(i,k)
cycle
endif
norm(k) += psi_coef(i,k) * psi_coef(i,k)
enddo
print*,'norm = ',norm(k)
enddo
deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det)
soft_touch psi_coef
end
subroutine set_intermediate_normalization_mlct_old(norm,i_particl)
implicit none
integer, intent(in) :: i_particl
double precision, intent(out) :: norm(N_states)
integer :: i,j,degree,index_ref_generators_restart,k
integer:: number_of_holes,n_h, number_of_particles,n_p
integer, allocatable :: index_one_hole(:),index_one_hole_one_p(:),index_two_hole_one_p(:),index_two_hole(:)
integer, allocatable :: index_one_p(:),index_one_hole_two_p(:)
integer :: n_one_hole,n_one_hole_one_p,n_two_hole_one_p,n_two_hole,n_one_p,n_one_hole_two_p
logical :: is_the_particl_in_det
double precision :: inv_coef_ref_generators_restart(N_states)
integer :: exc(0:2,2,2)
double precision :: phase,hij,hii,accu
integer :: h1,p1,h2,p2,s1,s2
integer :: index_good_particl(1000)
integer :: n_good_particl
logical,allocatable :: is_a_ref_det(:)
integer :: i_count
allocate(index_one_hole(n_det),index_one_hole_one_p(n_det),index_two_hole_one_p(N_det),index_two_hole(N_det),index_one_p(N_det),is_a_ref_det(N_det))
allocate(index_one_hole_two_p(n_det))
n_one_hole = 0
n_one_hole_one_p = 0
n_two_hole_one_p = 0
n_two_hole = 0
n_one_p = 0
n_one_hole_two_p = 0
n_good_particl = 0
! Find the one holes and one hole one particle
i_count = 0
is_a_ref_det = .False.
do i = 1, N_det
call get_excitation_degree(ref_generators_restart,psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_ref_generators_restart = i
do k = 1, N_states
inv_coef_ref_generators_restart(k) = 1.d0/psi_coef(i,k)
enddo
! cycle
endif
! Find all the determinants present in the reference wave function
do j = 1, N_det_generators_restart
call get_excitation_degree(psi_det(1,1,i),psi_det_generators_restart(1,1,j),degree,N_int)
if(degree == 0)then
is_a_ref_det(i) = .True.
exit
endif
enddo
if(is_a_ref_det(i))cycle
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if(n_h == 0 .and. n_p == 1)then ! 1p
if(is_the_particl_in_det(psi_det(1,1,i),1,i_particl).or.is_the_particl_in_det(psi_det(1,1,i),2,i_particl))then
n_good_particl += 1
index_good_particl(n_good_particl) = i
else
do k = 1, N_states
psi_coef(i,k) = 0.d0
enddo
endif
else
do k = 1, N_states
psi_coef(i,k) = 0.d0
enddo
endif
enddo
norm = 0.d0
print*,''
print*,'n_good_particl = ',n_good_particl
do k = 1, N_states
print*,'state ',k
do i = 1, n_good_particl
print*,'psi_coef(index_good_particl,1) = ',psi_coef(index_good_particl(i),k)/psi_coef(index_ref_generators_restart,k)
enddo
print*,''
enddo
! Set the wave function to the intermediate normalization
do k = 1, N_states
do i = 1, N_det
psi_coef(i,k) = psi_coef(i,k) * inv_coef_ref_generators_restart(k)
enddo
enddo
do k = 1, N_states
print*,'state ',k
do i = 1, N_det
!! print*,'i = ',i, psi_coef(i,1)
if (is_a_ref_det(i))then
print*,'i,psi_coef_ref = ',psi_coef(i,k)
cycle
endif
norm(k) += psi_coef(i,k) * psi_coef(i,k)
enddo
print*,'norm = ',norm
enddo
soft_touch psi_coef
deallocate(index_one_hole,index_one_hole_one_p,index_two_hole_one_p,index_two_hole,index_one_p,is_a_ref_det)
end
subroutine update_density_matrix_osoci
implicit none
BEGIN_DOC
! one_body_dm_mo_alpha_osoci += Delta rho alpha
! one_body_dm_mo_beta_osoci += Delta rho beta
END_DOC
integer :: i,j
integer :: iorb,jorb
do i = 1, mo_tot_num
do j = 1, mo_tot_num
one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha(i,j) - one_body_dm_mo_alpha_generators_restart(i,j))
one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta(i,j) - one_body_dm_mo_beta_generators_restart(i,j))
enddo
enddo
end
subroutine initialize_density_matrix_osoci
implicit none
one_body_dm_mo_alpha_osoci = one_body_dm_mo_alpha_generators_restart
one_body_dm_mo_beta_osoci = one_body_dm_mo_beta_generators_restart
end
subroutine rescale_density_matrix_osoci(norm)
implicit none
double precision, intent(in) :: norm(N_states)
integer :: i,j
double precision :: norm_tmp
norm_tmp = 0.d0
do i = 1, N_states
norm_tmp += norm(i)
enddo
print*,'norm = ',norm_tmp
do i = 1, mo_tot_num
do j = 1,mo_tot_num
one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) * norm_tmp
one_body_dm_mo_beta_osoci(j,i) = one_body_dm_mo_beta_osoci(j,i) * norm_tmp
enddo
enddo
end
subroutine save_osoci_natural_mos
implicit none
BEGIN_DOC
! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
END_DOC
character*(64) :: label
double precision, allocatable :: tmp(:,:),tmp_bis(:,:)
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha,i,i_core,j_core,iorb,jorb,j,i_inact,j_inact,i_virt,j_virt
allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
allocate(tmp_bis(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
allocate (occ(N_int*bit_kind_size,2))
! Negation to have the occupied MOs first after the diagonalization
tmp_bis = -one_body_dm_mo_alpha_osoci - one_body_dm_mo_beta_osoci
! Set to Zero the core-inact-act-virt part
do i = 1, n_core_orb
i_core = list_core(i)
tmp_bis(i_core,i_core) = -10.d0
do j = i+1, n_core_orb
j_core = list_core(j)
tmp_bis(i_core,j_core) = 0.d0
tmp_bis(j_core,i_core) = 0.d0
enddo
do j = 1, n_inact_orb
iorb = list_inact(j)
tmp_bis(i_core,iorb) = 0.d0
tmp_bis(iorb,i_core) = 0.d0
enddo
do j = 1, n_act_orb
iorb = list_act(j)
tmp_bis(i_core,iorb) = 0.d0
tmp_bis(iorb,i_core) = 0.d0
enddo
do j = 1, n_virt_orb
iorb = list_virt(j)
tmp_bis(i_core,iorb) = 0.d0
tmp_bis(iorb,i_core) = 0.d0
enddo
enddo
do i = 1, n_core_orb
print*,'dm core = ',list_core(i),tmp_bis(list_core(i),list_core(i))
enddo
! Set to Zero the inact-inact part to avoid arbitrary rotations
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = i+1, n_inact_orb
j_inact = list_inact(j)
tmp_bis(i_inact,j_inact) = 0.d0
tmp_bis(j_inact,i_inact) = 0.d0
enddo
enddo
! Set to Zero the inact-virt part to avoid arbitrary rotations
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = 1, n_virt_orb
j_virt = list_virt(j)
tmp_bis(i_inact,j_virt) = 0.d0
tmp_bis(j_virt,i_inact) = 0.d0
enddo
enddo
! Set to Zero the virt-virt part to avoid arbitrary rotations
do i = 1, n_virt_orb
i_virt = list_virt(i)
do j = i+1, n_virt_orb
j_virt = list_virt(j)
tmp_bis(i_virt,j_virt) = 0.d0
tmp_bis(j_virt,i_virt) = 0.d0
enddo
enddo
double precision :: accu
! Set to Zero the act-act part to avoid arbitrary rotations
do i = 1,n_act_orb
iorb = list_act(i)
do j = i+1,n_act_orb
jorb = list_act(j)
tmp_bis(iorb,jorb) = 0.d0
tmp_bis(jorb,iorb) = 0.d0
enddo
enddo
tmp = tmp_bis
!! Symetrization act-virt
do j = 1, n_virt_orb
j_virt= list_virt(j)
accu = 0.d0
do i = 1, n_act_orb
jorb = list_act(i)
accu += dabs(tmp_bis(j_virt,jorb))
enddo
do i = 1, n_act_orb
iorb = list_act(i)
tmp(j_virt,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb))
tmp(iorb,j_virt) = dsign(accu/dble(n_act_orb),tmp_bis(j_virt,iorb))
enddo
enddo
!! Symetrization act-inact
!do j = 1, n_inact_orb
! j_inact = list_inact(j)
! accu = 0.d0
! do i = 1, n_act_orb
! jorb = list_act(i)
! accu += dabs(tmp_bis(j_inact,jorb))
! enddo
! do i = 1, n_act_orb
! iorb = list_act(i)
! tmp(j_inact,iorb) = dsign(accu/dble(n_act_orb),tmp_bis(j_inact,iorb))
! tmp(iorb,j_inact) = dsign(accu/dble(n_act_orb),tmp_bis(j_inact,iorb))
! enddo
!enddo
!!! Symetrization act-act
!!accu = 0.d0
!!do i = 1, n_act_orb
!! iorb = list_act(i)
!! accu += tmp_bis(iorb,iorb)
!!enddo
!!do i = 1, n_act_orb
!! iorb = list_act(i)
!! tmp(iorb,iorb) = accu/dble(n_act_orb)
!!enddo
call bitstring_to_list(reunion_of_bitmask(1,1), occ(1,1), n_occ_alpha, N_int)
double precision :: maxvaldm,imax,jmax
maxvaldm = 0.d0
imax = 1
jmax = 1
print*,''
print*,'Inactive-active Part of the One body DM'
print*,''
do i = 1,n_act_orb
iorb = list_act(i)
print*,''
print*,'ACTIVE ORBITAL ',iorb
do j = 1, n_inact_orb
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
print*,'INACTIVE '
print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb))
endif
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
print*,'VIRT '
print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb))
endif
enddo
enddo
do i = 1, mo_tot_num
do j = i+1, mo_tot_num
if(dabs(tmp(i,j)).le.threshold_fobo_dm)then
tmp(i,j) = 0.d0
tmp(j,i) = 0.d0
endif
enddo
enddo
label = "Natural"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
soft_touch mo_coef
deallocate(tmp,occ)
end
subroutine set_osoci_natural_mos
implicit none
BEGIN_DOC
! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
END_DOC
character*(64) :: label
double precision, allocatable :: tmp(:,:),tmp_bis(:,:)
integer, allocatable :: occ(:,:)
integer :: n_occ_alpha,i,i_core,j_core,iorb,jorb,j,i_inact,j_inact,i_virt,j_virt
allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
allocate(tmp_bis(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
allocate (occ(N_int*bit_kind_size,2))
! Negation to have the occupied MOs first after the diagonalization
tmp_bis = -one_body_dm_mo_alpha_osoci - one_body_dm_mo_beta_osoci
! Set to Zero the core-inact-act-virt part
do i = 1, n_core_orb
i_core = list_core(i)
tmp_bis(i_core,i_core) = -10.d0
do j = i+1, n_core_orb
j_core = list_core(j)
tmp_bis(i_core,j_core) = 0.d0
tmp_bis(j_core,i_core) = 0.d0
enddo
do j = 1, n_inact_orb
iorb = list_inact(j)
tmp_bis(i_core,iorb) = 0.d0
tmp_bis(iorb,i_core) = 0.d0
enddo
do j = 1, n_act_orb
iorb = list_act(j)
tmp_bis(i_core,iorb) = 0.d0
tmp_bis(iorb,i_core) = 0.d0
enddo
do j = 1, n_virt_orb
iorb = list_virt(j)
tmp_bis(i_core,iorb) = 0.d0
tmp_bis(iorb,i_core) = 0.d0
enddo
enddo
do i = 1, n_core_orb
print*,'dm core = ',list_core(i),tmp_bis(list_core(i),list_core(i))
enddo
! Set to Zero the inact-inact part to avoid arbitrary rotations
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = i+1, n_inact_orb
j_inact = list_inact(j)
tmp_bis(i_inact,j_inact) = 0.d0
tmp_bis(j_inact,i_inact) = 0.d0
enddo
enddo
! Set to Zero the inact-virt part to avoid arbitrary rotations
do i = 1, n_inact_orb
i_inact = list_inact(i)
do j = 1, n_virt_orb
j_virt = list_virt(j)
tmp_bis(i_inact,j_virt) = 0.d0
tmp_bis(j_virt,i_inact) = 0.d0
enddo
enddo
! Set to Zero the virt-virt part to avoid arbitrary rotations
do i = 1, n_virt_orb
i_virt = list_virt(i)
do j = i+1, n_virt_orb
j_virt = list_virt(j)
tmp_bis(i_virt,j_virt) = 0.d0
tmp_bis(j_virt,i_virt) = 0.d0
enddo
enddo
double precision :: accu
! Set to Zero the act-act part to avoid arbitrary rotations
do i = 1,n_act_orb
iorb = list_act(i)
do j = i+1,n_act_orb
jorb = list_act(j)
tmp_bis(iorb,jorb) = 0.d0
tmp_bis(jorb,iorb) = 0.d0
enddo
enddo
tmp = tmp_bis
call bitstring_to_list(reunion_of_bitmask(1,1), occ(1,1), n_occ_alpha, N_int)
double precision :: maxvaldm,imax,jmax
maxvaldm = 0.d0
imax = 1
jmax = 1
print*,''
print*,'Inactive-active Part of the One body DM'
print*,''
do i = 1,n_act_orb
iorb = list_act(i)
print*,''
print*,'ACTIVE ORBITAL ',iorb
do j = 1, n_inact_orb
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
print*,'INACTIVE '
print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb))
endif
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
print*,'VIRT '
print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb))
endif
enddo
enddo
do i = 1, mo_tot_num
do j = i+1, mo_tot_num
if(dabs(tmp(i,j)).le.threshold_fobo_dm)then
tmp(i,j) = 0.d0
tmp(j,i) = 0.d0
endif
enddo
enddo
label = "Natural"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
soft_touch mo_coef
deallocate(tmp,occ)
end
subroutine check_symetry(i_hole,thr,test)
implicit none
integer, intent(in) :: i_hole
double precision, intent(in) :: thr
logical, intent(out) :: test
integer :: i,j,k,l
double precision :: accu
accu = 0.d0
do i = 1, n_act_orb
accu += dabs(mo_mono_elec_integral(i_hole,list_act(i)))
enddo
if(accu.gt.thr)then
test = .True.
else
test = .false.
endif
end
subroutine check_symetry_1h1p(i_hole,i_part,thr,test)
implicit none
integer, intent(in) :: i_hole,i_part
double precision, intent(in) :: thr
logical, intent(out) :: test
integer :: i,j,k,l
double precision :: accu
accu = dabs(mo_mono_elec_integral(i_hole,i_part))
if(accu.gt.thr)then
test = .True.
else
test = .false.
endif
end
subroutine update_one_body_dm_mo
implicit none
integer :: i
double precision :: accu_tot,accu_sd
print*,'touched the one_body_dm_mo_beta'
one_body_dm_mo_alpha = one_body_dm_mo_alpha_osoci
one_body_dm_mo_beta = one_body_dm_mo_beta_osoci
touch one_body_dm_mo_alpha one_body_dm_mo_beta
accu_tot = 0.d0
accu_sd = 0.d0
do i = 1, mo_tot_num
accu_tot += one_body_dm_mo_alpha(i,i) + one_body_dm_mo_beta(i,i)
accu_sd += one_body_dm_mo_alpha(i,i) - one_body_dm_mo_beta(i,i)
enddo
print*,'accu_tot = ',accu_tot
print*,'accu_sdt = ',accu_sd
end
subroutine provide_properties
implicit none
integer :: i
double precision :: accu
if(.True.)then
accu= 0.d0
do i = 1, nucl_num
accu += mulliken_spin_densities(i)
print*,i,nucl_charge(i),mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
endif
end

7
plugins/FOBOCI/save_fock_diag_inactiv_virt.irp.f Normal file
View File

@ -0,0 +1,7 @@
program save_fock_inactiv_virt_mos
implicit none
call diag_inactive_virt_and_update_mos
call save_mos
end

6
plugins/Hartree_Fock/EZFIO.cfg
View File

@ -33,3 +33,9 @@ doc: If true, skip the (inactive+core) --> (active) and the (active) --> (virtua
interface: ezfio,provider,ocaml
default: False
[take_input_guess]
type: logical
doc: If true, the MOs in the directory are taken as input
interface: ezfio,provider,ocaml
default: False

2
plugins/loc_cele/loc.f
View File

@ -17,7 +17,7 @@ C
data small/1.d-6/
zprt=.true.
niter=100
niter=500
conv=1.d-8
write (6,5) n,m,conv

299
plugins/loc_cele/loc_cele.irp.f
View File

@ -9,7 +9,7 @@
! id1=max is the number of MO in a given symmetry.
END_DOC
integer id1
integer id1,i_atom,shift,shift_h
parameter (id1=300)
@ -92,7 +92,7 @@
nrot(1) = 6 ! number of orbitals to be localized
nrot(1) = 6 ! number of orbitals to be localized
integer :: index_rot(1000,1)
@ -101,12 +101,30 @@
cmoref = 0.d0
! Definition of the index of the MO to be rotated
irot(1,1) = 20 ! the first mo to be rotated is the 19 th MO
irot(2,1) = 21 ! the first mo to be rotated is the 20 th MO
irot(3,1) = 22 ! etc....
irot(4,1) = 23 !
irot(5,1) = 24 !
irot(6,1) = 25 !
! irot(2,1) = 21 ! the first mo to be rotated is the 21 th MO
! irot(3,1) = 22 ! etc....
! irot(4,1) = 23 !
! irot(5,1) = 24 !
! irot(6,1) = 25 !
! do i = 1,12
! irot(i,1) = i+6
! enddo
irot(1,1) = 5
irot(2,1) = 6
irot(3,1) = 7
irot(4,1) = 8
irot(5,1) = 9
irot(6,1) = 10
do i = 1, nrot(1)
print*,'irot(i,1) = ',irot(i,1)
enddo
pause
cmoref(4,1,1) = 1.d0 ! 2S function
cmoref(5,2,1) = 1.d0 ! 2S function
cmoref(6,3,1) = 1.d0 ! 2S function
cmoref(19,4,1) = 1.d0 ! 2S function
cmoref(20,5,1) = 1.d0 ! 2S function
cmoref(21,6,1) = 1.d0 ! 2S function
! you define the guess vectors that you want
! the new MO to be close to
@ -120,23 +138,221 @@
! own guess vectors for the MOs
! The new MOs are provided in output
! in the same order than the guess MOs
cmoref(3,1,1) = 1.d0 !
cmoref(12,1,1) = 1.d0 !
! C-C bonds
! 1-2
! i_atom = 1
! shift = (i_atom -1) * 15
! cmoref(1+shift,1,1) = -0.012d0 ! 2S function
! cmoref(2+shift,1,1) = 0.18d0 !
! cmoref(3+shift,1,1) = 0.1d0 !
cmoref(21,2,1) = 1.d0 !
cmoref(30,2,1) = 1.d0 !
! cmoref(5+shift,1,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,1,1) = -0.1d0 ! 2pZ function
cmoref(39,3,1) = 1.d0 !
cmoref(48,3,1) = 1.d0 !
! i_atom = 2
! shift = (i_atom -1) * 15
! cmoref(1+shift,1,1) = -0.012d0 ! 2S function
! cmoref(2+shift,1,1) = 0.18d0 !
! cmoref(3+shift,1,1) = 0.1d0 !
cmoref(3,4,1) = 1.d0 !
cmoref(12,4,1) =-1.d0 !
! cmoref(5+shift,1,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,1,1) = 0.1d0 ! 2pZ function
cmoref(21,5,1) = 1.d0 !
cmoref(30,5,1) =-1.d0 !
cmoref(39,6,1) = 1.d0 !
cmoref(48,6,1) =-1.d0 !
! ! 1-3
! i_atom = 1
! shift = (i_atom -1) * 15
! cmoref(1+shift,2,1) = -0.012d0 ! 2S function
! cmoref(2+shift,2,1) = 0.18d0 !
! cmoref(3+shift,2,1) = 0.1d0 !
! cmoref(5+shift,2,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,2,1) = -0.1d0 ! 2pZ function
! i_atom = 3
! shift = (i_atom -1) * 15
! cmoref(1+shift,2,1) = -0.012d0 ! 2S function
! cmoref(2+shift,2,1) = 0.18d0 !
! cmoref(3+shift,2,1) = 0.1d0 !
! cmoref(5+shift,2,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,2,1) = 0.1d0 ! 2pZ function
! ! 4-6
! i_atom = 4
! shift = (i_atom -1) * 15
! cmoref(1+shift,3,1) = -0.012d0 ! 2S function
! cmoref(2+shift,3,1) = 0.18d0 !
! cmoref(3+shift,3,1) = 0.1d0 !
! cmoref(5+shift,3,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,3,1) = -0.1d0 ! 2pZ function
! i_atom = 6
! shift = (i_atom -1) * 15
! cmoref(1+shift,3,1) = -0.012d0 ! 2S function
! cmoref(2+shift,3,1) = 0.18d0 !
! cmoref(3+shift,3,1) = 0.1d0 !
! cmoref(5+shift,3,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,3,1) = 0.1d0 ! 2pZ function
! ! 6-5
! i_atom = 6
! shift = (i_atom -1) * 15
! cmoref(1+shift,4,1) = -0.012d0 ! 2S function
! cmoref(2+shift,4,1) = 0.18d0 !
! cmoref(3+shift,4,1) = 0.1d0 !
! cmoref(5+shift,4,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,4,1) = 0.1d0 ! 2pZ function
! i_atom = 5
! shift = (i_atom -1) * 15
! cmoref(1+shift,4,1) = -0.012d0 ! 2S function
! cmoref(2+shift,4,1) = 0.18d0 !
! cmoref(3+shift,4,1) = 0.1d0 !
! cmoref(5+shift,4,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,4,1) = -0.1d0 ! 2pZ function
! ! 2-4
! i_atom = 2
! shift = (i_atom -1) * 15
! cmoref(1+shift,5,1) = -0.012d0 ! 2S function
! cmoref(2+shift,5,1) = 0.18d0 !
! cmoref(3+shift,5,1) = 0.1d0 !
! cmoref(6+shift,5,1) = 0.1d0 ! 2pZ function
! i_atom = 4
! shift = (i_atom -1) * 15
! cmoref(1+shift,5,1) = -0.012d0 ! 2S function
! cmoref(2+shift,5,1) = 0.18d0 !
! cmoref(3+shift,5,1) = 0.1d0 !
! cmoref(6+shift,5,1) = -0.1d0 ! 2pZ function
! ! 3-5
! i_atom = 3
! shift = (i_atom -1) * 15
! cmoref(1+shift,6,1) = -0.012d0 ! 2S function
! cmoref(2+shift,6,1) = 0.18d0 !
! cmoref(3+shift,6,1) = 0.1d0 !
! cmoref(6+shift,6,1) = 0.1d0 ! 2pZ function
! i_atom = 5
! shift = (i_atom -1) * 15
! cmoref(1+shift,6,1) = -0.012d0 ! 2S function
! cmoref(2+shift,6,1) = 0.18d0 !
! cmoref(3+shift,6,1) = 0.1d0 !
! cmoref(6+shift,6,1) = -0.1d0 ! 2pZ function
! ! C-H bonds
! ! 2-7
! i_atom = 2
! shift = (i_atom -1) * 15
! cmoref(1+shift,7,1) = -0.012d0 ! 2S function
! cmoref(2+shift,7,1) = 0.18d0 !
! cmoref(3+shift,7,1) = 0.1d0 !
! cmoref(5+shift,7,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,7,1) = 0.1d0 ! 2pZ function
!
! i_atom = 7
! shift_h = (6-1) * 15 + (i_atom - 6)*5
! cmoref(1+shift_h,7,1) = 0.12d0 ! 1S function
! ! 4-10
! i_atom = 4
! shift = (i_atom -1) * 15
! cmoref(1+shift,8,1) = -0.012d0 ! 2S function
! cmoref(2+shift,8,1) = 0.18d0 !
! cmoref(3+shift,8,1) = 0.1d0 !
! cmoref(5+shift,8,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,8,1) = -0.1d0 ! 2pZ function
!
! i_atom = 10
! shift_h = (6-1) * 15 + (i_atom - 6)*5
! cmoref(1+shift_h,8,1) = 0.12d0 ! 1S function
! ! 5-11
! i_atom = 5
! shift = (i_atom -1) * 15
! cmoref(1+shift,9,1) = -0.012d0 ! 2S function
! cmoref(2+shift,9,1) = 0.18d0 !
! cmoref(3+shift,9,1) = 0.1d0 !
! cmoref(5+shift,9,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,9,1) = -0.1d0 ! 2pZ function
!
! i_atom = 11
! shift_h = (6-1) * 15 + (i_atom - 6)*5
! cmoref(1+shift_h,9,1) = 0.12d0 ! 1S function
! ! 3-8
! i_atom = 3
! shift = (i_atom -1) * 15
! cmoref(1+shift,10,1) = -0.012d0 ! 2S function
! cmoref(2+shift,10,1) = 0.18d0 !
! cmoref(3+shift,10,1) = 0.1d0 !
!
! cmoref(5+shift,10,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,10,1) = 0.1d0 ! 2pZ function
!
! i_atom = 8
! shift_h = (6-1) * 15 + (i_atom - 6)*5
! cmoref(1+shift_h,10,1) = 0.12d0 ! 1S function
! ! 1-9
! i_atom = 1
! shift = (i_atom -1) * 15
! cmoref(1+shift,11,1) = -0.012d0 ! 2S function
! cmoref(2+shift,11,1) = 0.18d0 !
! cmoref(3+shift,11,1) = 0.1d0 !
!
! cmoref(6+shift,11,1) = 0.1d0 ! 2pZ function
! i_atom = 9
! shift_h = (6-1) * 15 + (i_atom - 6)*5
! cmoref(1+shift_h,11,1) = 0.12d0 ! 1S function
!
! ! 6-12
! i_atom = 6
! shift = (i_atom -1) * 15
! cmoref(1+shift,12,1) = -0.012d0 ! 2S function
! cmoref(2+shift,12,1) = 0.18d0 !
! cmoref(3+shift,12,1) = 0.1d0 !
!
! cmoref(6+shift,12,1) = -0.1d0 ! 2pZ function
! i_atom = 12
! shift_h = (6-1) * 15 + (i_atom - 6)*5
! cmoref(1+shift_h,12,1) = 0.12d0 ! 1S function
! cmoref(12,1,1) = 1.d0 !
! cmoref(21,2,1) = 1.d0 !
! cmoref(30,2,1) = 1.d0 !
! cmoref(39,3,1) = 1.d0 !
! cmoref(48,3,1) = 1.d0 !
! cmoref(3,4,1) = 1.d0 !
! cmoref(12,4,1) =-1.d0 !
! cmoref(21,5,1) = 1.d0 !
! cmoref(30,5,1) =-1.d0 !
! cmoref(39,6,1) = 1.d0 !
! cmoref(48,6,1) =-1.d0 !
@ -146,48 +362,11 @@
do isym=1,nsym
if (nrot(isym).eq.0) cycle
do i=1,ao_num
s(i,i,isym)=1.d0
do j=1,ao_num
if (i.ne.j) s(i,j,isym)=0.d0
ddum(i,j)=0.d0
do k=1,nmo(isym)
ddum(i,j)=ddum(i,j)+cmo(i,k,isym)*cmo(j,k,isym)
enddo
do i = 1, ao_num
do j = 1, ao_num
s(i,j,1) = ao_overlap(i,j)
enddo
enddo
call dgesv(ao_num,ao_num,ddum,id1,ipiv,s(1,1,isym),id1,info)
if (info.ne.0) then
write (6,*) 'Something wrong in dgsev',isym
stop
endif
enddo
!Now big loop over symmetry

44
scripts/generate_h_apply.py
View File

@ -23,7 +23,12 @@ deinit_thread
skip
init_main
filter_integrals
filter2p
filter2h2p
filter1h
filter1p
only_2p_single
only_2p_double
filter_only_1h1p_single
filter_only_1h1p_double
filterhole
@ -152,6 +157,32 @@ class H_apply(object):
self["filterparticle"] = """
if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind )cycle
"""
def filter_1h(self):
self["filter1h"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h(hole)) cycle
"""
def filter_2p(self):
self["filter2p"] = """
! ! DIR$ FORCEINLINE
if (is_a_2p(hole)) cycle
"""
def filter_1p(self):
self["filter0p"] = """
! ! DIR$ FORCEINLINE
if (is_a_1p(hole)) cycle
"""
def filter_only_2p(self):
self["only_2p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_2p(hole).eq..False.) cycle
"""
self["only_2p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_2p(key).eq..False.) cycle
"""
def filter_only_1h1p(self):
self["filter_only_1h1p_single"] = """
@ -217,9 +248,16 @@ class H_apply(object):
PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
"""
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,N_int,key_mask,fock_diag_tmp)
"""%(pert,)
if(check_double_excitation)then
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,N_int,key_mask,fock_diag_tmp)
else
call perturb_buffer_by_mono_%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,N_int,key_mask,fock_diag_tmp)
endif
"""%(pert,pert)
self.data["finalization"] = """
"""
self.data["copy_buffer"] = ""

55
src/Bitmask/bitmasks.irp.f
View File

@ -399,18 +399,6 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)]
implicit none
BEGIN_DOC
! Bitmask of the core orbitals that are never excited in post CAS method
END_DOC
integer :: i,j
do i = 1, N_int
core_bitmask(i,1) = iand(ref_bitmask(i,1),reunion_of_bitmask(i,1))
core_bitmask(i,2) = iand(ref_bitmask(i,2),reunion_of_bitmask(i,2))
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, list_core, (n_core_orb)]
BEGIN_DOC
! List of the core orbitals that are never excited in post CAS method
@ -426,20 +414,21 @@ END_PROVIDER
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, n_core_orb ]
BEGIN_PROVIDER [ integer(bit_kind), core_bitmask, (N_int,2)]
&BEGIN_PROVIDER [ integer, n_core_orb]
implicit none
BEGIN_DOC
! Number of core orbitals that are never excited in post CAS method
! Core orbitals bitmask
END_DOC
logical :: exists
integer :: j,i
integer :: i_hole,i_part,i_gen
integer :: i,j
n_core_orb = 0
do j = 1, N_int
n_core_orb += popcnt(core_bitmask(j,1))
do i = 1, N_int
core_bitmask(i,1) = xor(closed_shell_ref_bitmask(i,1),reunion_of_cas_inact_bitmask(i,1))
core_bitmask(i,2) = xor(closed_shell_ref_bitmask(i,2),reunion_of_cas_inact_bitmask(i,2))
n_core_orb += popcnt(core_bitmask(i,1))
enddo
END_PROVIDER
print*,'n_core_orb = ',n_core_orb
END_PROVIDER
BEGIN_PROVIDER [ integer, i_bitmask_gen ]
@ -490,3 +479,27 @@ BEGIN_PROVIDER [integer, list_act, (n_act_orb)]
enddo
END_PROVIDER
BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)]
implicit none
integer :: i,j
do i = 1, N_int
closed_shell_ref_bitmask(i,1) = ior(ref_bitmask(i,1),cas_bitmask(i,1,1))
closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),cas_bitmask(i,2,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), reunion_of_cas_inact_bitmask, (N_int,2)]
implicit none
BEGIN_DOC
! Reunion of the inactive, active and virtual bitmasks
END_DOC
integer :: i,j
do i = 1, N_int
reunion_of_cas_inact_bitmask(i,1) = ior(cas_bitmask(i,1,1),inact_bitmask(i,1))
reunion_of_cas_inact_bitmask(i,2) = ior(cas_bitmask(i,2,1),inact_bitmask(i,2))
enddo
END_PROVIDER

55
src/Bitmask/find_hole.irp.f Normal file
View File

@ -0,0 +1,55 @@
logical function is_the_hole_in_det(key_in,ispin,i_hole)
use bitmasks
! returns true if the electron ispin is absent from i_hole
implicit none
integer, intent(in) :: i_hole,ispin
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer(bit_kind) :: key_tmp(N_int)
integer(bit_kind) :: itest(N_int)
integer :: i,j,k
do i = 1, N_int
itest(i) = 0_bit_kind
enddo
k = ishft(i_hole-1,-bit_kind_shift)+1
j = i_hole-ishft(k-1,bit_kind_shift)-1
itest(k) = ibset(itest(k),j)
j = 0
do i = 1, N_int
key_tmp(i) = iand(itest(i),key_in(i,ispin))
j += popcnt(key_tmp(i))
enddo
if(j==0)then
is_the_hole_in_det = .True.
else
is_the_hole_in_det = .False.
endif
end
logical function is_the_particl_in_det(key_in,ispin,i_particl)
use bitmasks
! returns true if the electron ispin is absent from i_particl
implicit none
integer, intent(in) :: i_particl,ispin
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer(bit_kind) :: key_tmp(N_int)
integer(bit_kind) :: itest(N_int)
integer :: i,j,k
do i = 1, N_int
itest(i) = 0_bit_kind
enddo
k = ishft(i_particl-1,-bit_kind_shift)+1
j = i_particl-ishft(k-1,bit_kind_shift)-1
itest(k) = ibset(itest(k),j)
j = 0
do i = 1, N_int
key_tmp(i) = iand(itest(i),key_in(i,ispin))
j += popcnt(key_tmp(i))
enddo
if(j==0)then
is_the_particl_in_det = .False.
else
is_the_particl_in_det = .True.
endif
end

280
src/Bitmask/modify_bitmasks.irp.f Normal file
View File

@ -0,0 +1,280 @@
use bitmasks
subroutine initialize_bitmask_to_restart_ones
implicit none
integer :: i,j,k,l,m
integer :: ispin
BEGIN_DOC
! Initialization of the generators_bitmask to the restart bitmask
END_DOC
do i = 1, N_int
do k=1,N_generators_bitmask
do ispin=1,2
generators_bitmask(i,ispin,s_hole ,k) = generators_bitmask_restart(i,ispin,s_hole ,k)
generators_bitmask(i,ispin,s_part ,k) = generators_bitmask_restart(i,ispin,s_part ,k)
generators_bitmask(i,ispin,d_hole1,k) = generators_bitmask_restart(i,ispin,d_hole1,k)
generators_bitmask(i,ispin,d_part1,k) = generators_bitmask_restart(i,ispin,d_part1,k)
generators_bitmask(i,ispin,d_hole2,k) = generators_bitmask_restart(i,ispin,d_hole2,k)
generators_bitmask(i,ispin,d_part2,k) = generators_bitmask_restart(i,ispin,d_part2,k)
enddo
enddo
enddo
end
subroutine modify_bitmasks_for_hole(i_hole)
implicit none
integer, intent(in) :: i_hole
integer :: i,j,k,l,m
integer :: ispin
BEGIN_DOC
! modify the generators_bitmask in order that one can only excite
! the electrons occupying i_hole
END_DOC
! Set to Zero the holes
do k=1,N_generators_bitmask
do l = 1, 3
i = index_holes_bitmask(l)
do ispin=1,2
do j = 1, N_int
generators_bitmask(j,ispin,i,k) = 0_bit_kind
enddo
enddo
enddo
enddo
k = ishft(i_hole-1,-bit_kind_shift)+1
j = i_hole-ishft(k-1,bit_kind_shift)-1
do m = 1, N_generators_bitmask
do l = 1, 3
i = index_holes_bitmask(l)
do ispin=1,2
generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j)
enddo
enddo
enddo
end
subroutine modify_bitmasks_for_hole_in_out(i_hole)
implicit none
integer, intent(in) :: i_hole
integer :: i,j,k,l,m
integer :: ispin
BEGIN_DOC
! modify the generators_bitmask in order that one can only excite
! the electrons occupying i_hole
END_DOC
k = ishft(i_hole-1,-bit_kind_shift)+1
j = i_hole-ishft(k-1,bit_kind_shift)-1
do m = 1, N_generators_bitmask
do l = 1, 3
i = index_holes_bitmask(l)
do ispin=1,2
generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j)
enddo
enddo
enddo
end
subroutine modify_bitmasks_for_particl(i_part)
implicit none
integer, intent(in) :: i_part
integer :: i,j,k,l,m
integer :: ispin
BEGIN_DOC
! modify the generators_bitmask in order that one can only excite
! the electrons to the orbital i_part
END_DOC
! Set to Zero the particles
do k=1,N_generators_bitmask
do l = 1, 3
i = index_particl_bitmask(l)
do ispin=1,2
do j = 1, N_int
generators_bitmask(j,ispin,i,k) = 0_bit_kind
enddo
enddo
enddo
enddo
k = ishft(i_part-1,-bit_kind_shift)+1
j = i_part-ishft(k-1,bit_kind_shift)-1
do m = 1, N_generators_bitmask
do l = 1, 3
i = index_particl_bitmask(l)
do ispin=1,2
generators_bitmask(k,ispin,i,m) = ibset(generators_bitmask(k,ispin,i,m),j)
enddo
enddo
enddo
end
subroutine set_bitmask_particl_as_input(input_bimask)
implicit none
integer(bit_kind), intent(in) :: input_bimask(N_int,2)
integer :: i,j,k,l,m
integer :: ispin
BEGIN_DOC
! set the generators_bitmask for the particles
! as the input_bimask
END_DOC
do k=1,N_generators_bitmask
do l = 1, 3
i = index_particl_bitmask(l)
do ispin=1,2
do j = 1, N_int
generators_bitmask(j,ispin,i,k) = input_bimask(j,ispin)
enddo
enddo
enddo
enddo
touch generators_bitmask
end
subroutine set_bitmask_hole_as_input(input_bimask)
implicit none
integer(bit_kind), intent(in) :: input_bimask(N_int,2)
integer :: i,j,k,l,m
integer :: ispin
BEGIN_DOC
! set the generators_bitmask for the holes
! as the input_bimask
END_DOC
do k=1,N_generators_bitmask
do l = 1, 3
i = index_holes_bitmask(l)
do ispin=1,2
do j = 1, N_int
generators_bitmask(j,ispin,i,k) = input_bimask(j,ispin)
enddo
enddo
enddo
enddo
touch generators_bitmask
end
subroutine print_generators_bitmasks_holes
implicit none
integer :: i,j,k,l
integer(bit_kind),allocatable :: key_tmp(:,:)
allocate(key_tmp(N_int,2))
do l = 1, 3
k = 1
i = index_holes_bitmask(l)
do j = 1, N_int
key_tmp(j,1) = generators_bitmask(j,1,i,k)
key_tmp(j,2) = generators_bitmask(j,2,i,k)
enddo
print*,''
print*,'index hole = ',i
call print_det(key_tmp,N_int)
print*,''
enddo
deallocate(key_tmp)
end
subroutine print_generators_bitmasks_particles
implicit none
integer :: i,j,k,l
integer(bit_kind),allocatable :: key_tmp(:,:)
allocate(key_tmp(N_int,2))
do l = 1, 3
k = 1
i = index_particl_bitmask(l)
do j = 1, N_int
key_tmp(j,1) = generators_bitmask(j,1,i,k)
key_tmp(j,2) = generators_bitmask(j,2,i,k)
enddo
print*,''
print*,'index particl ',i
call print_det(key_tmp,N_int)
print*,''
enddo
deallocate(key_tmp)
end
subroutine print_generators_bitmasks_holes_for_one_generator(i_gen)
implicit none
integer, intent(in) :: i_gen
integer :: i,j,k,l
integer(bit_kind),allocatable :: key_tmp(:,:)
allocate(key_tmp(N_int,2))
do l = 1, 3
k = i_gen
i = index_holes_bitmask(l)
do j = 1, N_int
key_tmp(j,1) = generators_bitmask(j,1,i,k)
key_tmp(j,2) = generators_bitmask(j,2,i,k)
enddo
print*,''
print*,'index hole = ',i
call print_det(key_tmp,N_int)
print*,''
enddo
deallocate(key_tmp)
end
subroutine print_generators_bitmasks_particles_for_one_generator(i_gen)
implicit none
integer, intent(in) :: i_gen
integer :: i,j,k,l
integer(bit_kind),allocatable :: key_tmp(:,:)
allocate(key_tmp(N_int,2))
do l = 1, 3
k = i_gen
i = index_particl_bitmask(l)
do j = 1, N_int
key_tmp(j,1) = generators_bitmask(j,1,i,k)
key_tmp(j,2) = generators_bitmask(j,2,i,k)
enddo
print*,''
print*,'index particl ',i
call print_det(key_tmp,N_int)
print*,''
enddo
deallocate(key_tmp)
end
BEGIN_PROVIDER [integer, index_holes_bitmask, (3)]
implicit none
BEGIN_DOC
! Index of the holes in the generators_bitmasks
END_DOC
index_holes_bitmask(1) = d_hole1
index_holes_bitmask(2) = d_hole2
index_holes_bitmask(3) = s_hole
END_PROVIDER
BEGIN_PROVIDER [integer, index_particl_bitmask, (3)]
implicit none
BEGIN_DOC
! Index of the holes in the generators_bitmasks
END_DOC
index_particl_bitmask(1) = d_part1
index_particl_bitmask(2) = d_part2
index_particl_bitmask(3) = s_part
END_PROVIDER

12
src/Determinants/H_apply.template.f
View File

@ -166,6 +166,9 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
logical :: check_double_excitation
logical :: is_a_1h1p
logical :: is_a_1h
logical :: is_a_1p
logical :: is_a_2p
logical :: b_cycle
check_double_excitation = .True.
iproc = iproc_in
@ -300,6 +303,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
key(k,other_spin) = ibset(key(k,other_spin),l)
$filter2h2p
$filter_only_1h1p_double
$only_2p_double
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
@ -351,6 +355,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
key(k,ispin) = ibset(key(k,ispin),l)
$filter2h2p
$filter_only_1h1p_double
$only_2p_double
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
@ -422,6 +427,9 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
logical :: check_double_excitation
logical :: is_a_1h1p
logical :: is_a_1h
logical :: is_a_1p
logical :: is_a_2p
key_mask(:,:) = 0_bit_kind
@ -493,6 +501,10 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
$filterparticle
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
$only_2p_single
$filter1h
$filter1p
$filter2p
$filter2h2p
$filter_only_1h1p_single
key_idx += 1

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

@ -642,6 +642,14 @@ subroutine read_dets(det,Nint,Ndet)
end
subroutine save_ref_determinant
implicit none
use bitmasks
call save_wavefunction_general(1,1,ref_bitmask,1,1.d0)
end
subroutine save_wavefunction
implicit none

5
src/Determinants/save_HF_determinant.irp.f Normal file
View File

@ -0,0 +1,5 @@
program save_HF
implicit none
call save_ref_determinant
end