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

Merge pull request #147 from eginer/master

Browse Source 
Add new modules for OVB analysis, Dressed Hamiltonian and FOBOCI programs. Also, I corrected the CAS+S program
This commit is contained in:
Anthony Scemama 2016-02-18 18:06:39 +01:00
commit e5c668deef
48 changed files with 5550 additions and 88 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 ;

4
plugins/CAS_SD/H_apply.irp.f
View File

@ -20,18 +20,22 @@ print s
s = H_apply("CAS_S",do_double_exc=False)
s.unset_double_excitations()
print s
s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
s.unset_double_excitations()
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
print s
s = H_apply("CAS_S_selected",do_double_exc=False)
s.unset_double_excitations()
s.set_selection_pt2("epstein_nesbet_2x2")
print s
s = H_apply("CAS_S_PT2",do_double_exc=False)
s.unset_double_excitations()
s.set_perturbation("epstein_nesbet_2x2")
print s

37
plugins/Dressed_Ref_Hamiltonian/Dressed_Ref_Hamiltonian.main.irp.f Normal file
View File

@ -0,0 +1,37 @@
program Dressed_Ref_Hamiltonian implicit none
BEGIN_DOC
! TODO
END_DOC
print *, ' _/ '
print *, ' -:\_?, _Jm####La '
print *, 'J"(:" > _]#AZ#Z#UUZ##, '
print *, '_,::./ %(|i%12XmX1*1XL _?, '
print *, ' \..\ _\(vmWQwodY+ia%lnL _",/ ( '
print *, ' .:< ]J=mQD?WXn<uQWmmvd, -.-:=!'
print *, ' "{Z jC]QW|=3Zv)Bi3BmXv3 = _7'
print *, ' ]h[Z6)WQ;)jZs]C;|$BZv+, : ./ '
print *, ' -#sJX%$Wmm#ev]hinW#Xi:` c ; '
print *, ' #X#X23###1}vI$WWmX1>|,)nr" '
print *, ' 4XZ#Xov1v}=)vnXAX1nnv;1n" '
print *, ' ]XX#ZXoovvvivnnnlvvo2*i7 '
print *, ' "23Z#1S2oo2XXSnnnoSo2>v" '
print *, ' miX#L -~`""!!1}oSoe|i7 '
print *, ' 4cn#m, v221=|v[ '
print *, ' ]hI3Zma,;..__wXSe=+vo '
print *, ' ]Zov*XSUXXZXZXSe||vo2 '
print *, ' ]Z#><iiii|i||||==vn2( '
print *, ' ]Z#i<ii||+|=||=:{no2[ '
print *, ' ]ZUsiiiiivi|=||=vo22[ '
print *, ' ]XZvlliiIi|i=|+|vooo '
print *, ' =v1llli||||=|||||lii( '
print *, ' ]iillii||||||||=>=|< '
print *, ' -ziiiii||||||+||==+> '
print *, ' -%|+++||=|=+|=|==/ '
print *, ' -a>====+|====-:- '
print *, ' "~,- -- /- '
print *, ' -. )> '
print *, ' .~ +- '
print *, ' . .... : . '
print *, ' -------~ '
print *, ''
end

1
plugins/Dressed_Ref_Hamiltonian/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
MRCC_Utils

16
plugins/Dressed_Ref_Hamiltonian/README.rst Normal file
View File

@ -0,0 +1,16 @@
=======================
Dressed_Ref_Hamiltonian
=======================
The following modules proposes to build an effective Hamiltonian
spanned on the reference determinants supposed to be the CAS ones.
The effective matrix Hamiltonian are built using the multi parentage
proposal used in the MR-CCSD formalism of Giner et. al. (JCP, 144, 064101 (2016); doi: 10.1063/1.4940781)
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.

41
plugins/Dressed_Ref_Hamiltonian/dressed_eigenvectors.irp.f Normal file
View File

@ -0,0 +1,41 @@
BEGIN_PROVIDER [double precision, psi_ref_coef_dressed, (n_det_ref,N_states) ]
&BEGIN_PROVIDER [double precision, energies_ref_dressed, (N_states) ]
implicit none
integer :: i,j,k,l,istate,igoodstate
double precision, allocatable :: H_matrix_tmp(:,:)
double precision, allocatable :: eigvalues(:),eigvectors(:,:),psi_coef_ref_tmp(:)
double precision :: accu, accu1
allocate(H_matrix_tmp(n_det_ref,n_det_ref))
allocate(eigvalues(n_det_ref))
allocate(eigvectors(n_det_ref,n_det_ref))
allocate(psi_coef_ref_tmp(n_det_ref))
do istate = 1, N_states
accu1 = 0.d0
do j = 1, n_det_ref
accu1 += psi_ref_coef(j,istate)**2 ! norm of the "istate" eigenvector in the projected in the reference space
do k = 1, n_det_ref
H_matrix_tmp(j,k) = hamiltonian_total_dressed(j,k,istate)
enddo
enddo
accu1 = 1.d0/dsqrt(accu1)
do j = 1, n_det_ref
psi_coef_ref_tmp(j) = psi_ref_coef(j,istate) * accu1
enddo
call lapack_diagd(eigvalues,eigvectors,H_matrix_tmp,n_det_ref,n_det_ref)
do j = 1, n_det_ref
accu = 0.d0
do k = 1, n_det_ref
accu += eigvectors(k,j) * psi_coef_ref_tmp(k)
enddo
if(dabs(accu).gt.0.9d0)then
igoodstate = j
exit
endif
enddo
energies_ref_dressed(istate) = eigvalues(igoodstate)
do j = 1,n_det_ref
psi_ref_coef_dressed(j,istate) = eigvectors(j,igoodstate)
enddo
enddo
END_PROVIDER

46
plugins/Dressed_Ref_Hamiltonian/dressed_hamiltonian.irp.f Normal file
View File

@ -0,0 +1,46 @@
BEGIN_PROVIDER [double precision, dressing_ref_hamiltonian, (n_det_ref,n_det_ref,N_states)]
implicit none
integer :: i,j,k,l
integer :: ii,jj,istate
double precision :: hij,sec_order,H_ref(N_det_ref),hik,hkl
integer :: idx(0:N_det_ref)
double precision :: accu_negative,accu_positive,phase
integer :: degree_exc_ionic,degree_exc_neutral,exc(0:2,2,2)
dressing_ref_hamiltonian = 0.d0
accu_negative = 0.d0
accu_positive = 0.d0
integer :: h1,p1,h2,p2,s1,s2
do istate = 1, N_states
do i = 1, N_det_non_ref
call filter_connected_i_H_psi0(psi_ref,psi_non_ref(1,1,i),N_int,N_det_ref,idx)
H_ref = 0.d0
do ii=1,idx(0)
k = idx(ii)
!DEC$ FORCEINLINE
call i_H_j(psi_ref(1,1,k),psi_non_ref(1,1,i),N_int,hij)
H_ref(k) = hij
enddo
do ii= 1, idx(0)
k = idx(ii)
hik = H_ref(k) * lambda_mrcc(istate,i)
do jj = 1, idx(0)
l = idx(jj)
dressing_ref_hamiltonian(k,l,istate) += hik * H_ref(l)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, hamiltonian_total_dressed, (n_det_ref,n_det_ref,N_states)]
implicit none
integer :: i,j,k
do k = 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
hamiltonian_total_dressed(j,i,k) = dressing_ref_hamiltonian(j,i,k) + ref_hamiltonian_matrix(j,i)
enddo
enddo
enddo
END_PROVIDER

108
plugins/Dressed_Ref_Hamiltonian/print_CAS_effective_Hamiltonian.irp.f Normal file
View File

@ -0,0 +1,108 @@
program print
read_wf = .True.
touch read_wf
call provide_all_stuffs
end
subroutine provide_all_stuffs
implicit none
provide ref_hamiltonian_matrix dressing_ref_hamiltonian
integer :: i,j,istate
double precision, allocatable :: psi_restart_ref_normalized(:),psi_ref_zeroth_order(:),psi_ref_dressed(:)
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
double precision, allocatable :: H_dressed(:,:)
double precision, allocatable :: H_print(:,:)
double precision :: accu_norm
allocate (H_dressed(N_det_ref,N_det_ref))
allocate (H_print(N_det_ref,N_det_ref))
allocate (psi_restart_ref_normalized(N_det_ref))
allocate (psi_ref_zeroth_order(N_det_ref))
print*,'# nuclear_repulsion = ',nuclear_repulsion
allocate (psi_ref_dressed(N_det_ref))
allocate (eigvalues(N_det_ref))
allocate (eigvectors(N_det_ref,N_det_ref))
do istate= 1, N_states
do i = 1, N_det_ref
do j = 1, N_det_ref
H_print(i,j) = ref_hamiltonian_matrix(j,i)
enddo
enddo
do i = 1, N_det_ref
H_print(i,i) -= ref_hamiltonian_matrix(1,1)
enddo
print*,'Ref Hamiltonian matrix emelent = ',ref_hamiltonian_matrix(1,1)
print*,'ISTATE = ',istate
accu_norm = 0.d0
do i = 1, N_det_ref
accu_norm += psi_ref_coef(i,1) * psi_ref_coef(i,1)
enddo
print*,'accu_norm = ',accu_norm
accu_norm = 1.d0/dsqrt(accu_norm)
do i = 1, N_det_ref
psi_restart_ref_normalized(i) = psi_ref_coef(i,istate)* accu_norm
enddo
print*,'-------------------'
print*,'-------------------'
print*,'CAS MATRIX '
print*,''
do i = 1, N_det_ref
write(*,'(10(F8.5 ,4X))') H_print(i,:)
enddo
print*,''
print*,'-------------------'
print*,'-------------------'
print*,'CAS MATRIX DRESSING'
print*,''
do i = 1, N_det_ref
write(*,'(10(F8.5 ,4X))') dressing_ref_hamiltonian(i,:,istate)
enddo
print*,''
print*,'-------------------'
print*,'-------------------'
do i = 1, N_det_ref
do j = 1, N_det_ref
H_dressed(j,i) = ref_hamiltonian_matrix(j,i) + dressing_ref_hamiltonian(j,i,istate)
H_print(i,j) += dressing_ref_hamiltonian(j,i,istate)
enddo
enddo
print*,''
print*,'-------------------'
print*,'-------------------'
print*,'TOTAL DRESSED H MATRIX '
print*,''
do i = 1, N_det_ref
write(*,'(10(F8.5 ,4X))') H_print(i,:)
enddo
print*,''
print*,''
print*,''
call lapack_diagd(eigvalues,eigvectors,ref_hamiltonian_matrix,n_det_ref,n_det_ref)
do i = 1, N_det_ref
psi_ref_zeroth_order(i) = eigvectors(i,istate)
enddo
call lapack_diagd(eigvalues,eigvectors,H_dressed,n_det_ref,n_det_ref)
do i = 1, N_det_ref
psi_ref_dressed(i) = eigvectors(i,istate)
enddo
print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion
do i = 1, N_det_ref
write(*,'(10(F10.7 ,4X))') psi_ref_coef(i,istate)/psi_ref_coef(1,istate), psi_ref_dressed(i)/psi_ref_dressed(1),psi_ref_zeroth_order(i)/psi_ref_zeroth_order(1)
enddo
enddo
deallocate (H_dressed)
deallocate (H_print)
deallocate (psi_restart_ref_normalized)
deallocate (psi_ref_zeroth_order)
deallocate (psi_ref_dressed)
deallocate (eigvalues)
deallocate (eigvectors)
end

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.

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

@ -0,0 +1,50 @@
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",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
print s
s = H_apply("just_mono_no_1h_no_1p",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.filter_1h()
s.filter_1p()
print s
s = H_apply("just_mono_no_1h_no_1p_no_2p",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
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 = 1.d-8
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

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

@ -0,0 +1,413 @@
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)
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

2
plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_full Generators_full Psiref_Utils
Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS

4
plugins/MRCC_Utils/mrcc_dummy.irp.f Normal file
View File

@ -0,0 +1,4 @@
program pouet
end

1
plugins/OVB_effective_Hamiltonian/NEEDED_CHILDREN_MODULES Normal file
View File

@ -0,0 +1 @@
Dressed_Ref_Hamiltonian OVB

59
plugins/OVB_effective_Hamiltonian/OVB_effective_H.irp.f Normal file
View File

@ -0,0 +1,59 @@
BEGIN_PROVIDER [double precision, H_OVB_dressing, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)]
BEGIN_DOC
! Hamiltonian matrix expressed in the basis of all the
END_DOC
implicit none
integer :: i,j,istate,k,l
double precision :: accu,hij
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
accu = 0.d0
do istate = 1, N_states
do k = 1, ionic_index(i,0)
do l = 1, ionic_index(j,0)
hij = dressing_ref_hamiltonian(ionic_index(i,k),ionic_index(j,l),istate)
! accu += psi_ref_coef(ionic_index(i,k),istate) * normalization_factor_ionic(i,istate) * &
! psi_ref_coef(ionic_index(j,l),istate) * normalization_factor_ionic(j,istate) * hij
accu += psi_ref_coef_dressed(ionic_index(i,k),istate) * normalization_factor_ionic_dressed(i,istate) * &
psi_ref_coef_dressed(ionic_index(j,l),istate) * normalization_factor_ionic_dressed(j,istate) * hij
enddo
enddo
H_OVB_dressing(i,j,istate) = accu
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, H_OVB_total_dressed, (min_number_ionic:max_number_ionic, min_number_ionic:max_number_ionic, n_states)]
BEGIN_DOC
! Hamiltonian matrix expressed in the basis of all the
END_DOC
implicit none
integer :: i,j,istate
double precision :: accu,hij
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
do istate = 1, N_states
H_OVB_total_dressed(i,j,istate) = H_OVB_dressing(i,j,istate) + H_OVB_naked(i,j,istate)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, normalization_factor_ionic_dressed, (min_number_ionic:max_number_ionic, N_states) ]
implicit none
integer :: i,j,istate,k
double precision :: accu
do j = min_number_ionic, max_number_ionic
do istate = 1, N_states
accu = 0.d0
do k = 1, ionic_index(j,0)
accu += psi_ref_coef_dressed(ionic_index(j,k),istate) **2
enddo
normalization_factor_ionic_dressed(j,istate) = 1.d0/dsqrt(accu)
enddo
enddo
END_PROVIDER

13
plugins/OVB_effective_Hamiltonian/README.rst Normal file
View File

@ -0,0 +1,13 @@
=========================
OVB_effective_Hamiltonian
=========================
Dressing of the OVB matrix by use of the Dressed_Ref_Hamiltonian dressing
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.

101
plugins/OVB_effective_Hamiltonian/print_OVB_effective_H_diagonalized.irp.f Normal file
View File

@ -0,0 +1,101 @@
program print
read_wf = .True.
touch read_wf
call provide_all_stuffs
end
subroutine provide_all_stuffs
implicit none
provide ref_hamiltonian_matrix dressing_ref_hamiltonian
integer :: i,j,istate
double precision, allocatable :: psi_restart_ref_normalized(:),psi_ref_zeroth_order(:),psi_ref_dressed(:)
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
double precision, allocatable :: H_naked(:,:)
double precision, allocatable :: H_dressed(:,:)
double precision, allocatable :: H_print(:,:)
double precision :: accu_norm
allocate (H_dressed(max_number_ionic+1,max_number_ionic+1))
allocate (H_print(min_number_ionic:max_number_ionic,min_number_ionic:max_number_ionic))
allocate (H_naked(max_number_ionic+1,max_number_ionic+1))
allocate (psi_restart_ref_normalized(min_number_ionic:max_number_ionic))
allocate (psi_ref_zeroth_order(min_number_ionic:max_number_ionic))
print*,'# nuclear_repulsion = ',nuclear_repulsion
allocate (psi_ref_dressed(min_number_ionic:max_number_ionic))
allocate (eigvalues(max_number_ionic+1))
allocate (eigvectors(max_number_ionic+1,max_number_ionic+1))
do istate= 1, N_states
print*,'ISTATE = ',istate
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
H_print(i,j) = H_OVB_naked(j,i,istate)
enddo
enddo
do i = min_number_ionic,max_number_ionic
H_print(i,i) -= H_OVB_naked(min_number_ionic,min_number_ionic,istate)
enddo
print*,'Ref Hamiltonian matrix emelent = ',H_OVB_naked(min_number_ionic,min_number_ionic,istate)
print*,'-------------------'
print*,'-------------------'
print*,'CAS MATRIX '
print*,''
do i = min_number_ionic,max_number_ionic
write(*,'(I4,X,10(F8.5 ,4X))')i, H_print(i,:)
enddo
print*,'CAS MATRIX DRESSING'
print*,''
do i = min_number_ionic,max_number_ionic
write(*,'(I4,X,10(F8.5 ,4X))')i, H_OVB_dressing(i,:,istate)
enddo
print*,''
print*,'-------------------'
print*,'-------------------'
print*,'CAS MATRIX DRESSED '
print*,''
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
H_print(i,j) = H_OVB_total_dressed(j,i,istate)
enddo
enddo
do i = min_number_ionic,max_number_ionic
H_print(i,i) -= H_OVB_total_dressed(min_number_ionic,min_number_ionic,istate)
enddo
do i = min_number_ionic,max_number_ionic
write(*,'(I4,X,10(F8.5 ,4X))')i, H_print(i,:)
enddo
print*,''
do i = min_number_ionic,max_number_ionic
do j = min_number_ionic,max_number_ionic
H_dressed(j+1,i+1) = H_OVB_total_dressed(i,j,istate)
H_naked(j+1,i+1) = H_OVB_naked(i,j,istate)
enddo
enddo
call lapack_diagd(eigvalues,eigvectors,H_naked,max_number_ionic+1,max_number_ionic+1)
print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion
do i = min_number_ionic,max_number_ionic
psi_ref_zeroth_order(i) = eigvectors(i+1,istate)
enddo
call lapack_diagd(eigvalues,eigvectors,H_dressed,max_number_ionic+1,max_number_ionic+1)
do i = min_number_ionic,max_number_ionic
psi_ref_dressed(i) = eigvectors(i+1,istate)
enddo
print*,'E+PT2 = ',eigvalues(istate) + nuclear_repulsion
do i = min_number_ionic,max_number_ionic
write(*,'(10(F10.7 ,4X))') psi_ref_dressed(i)/psi_ref_dressed(min_number_ionic) ,psi_ref_zeroth_order(i)/psi_ref_zeroth_order(min_number_ionic)
enddo
enddo
deallocate (H_dressed)
deallocate (H_naked)
deallocate (psi_restart_ref_normalized)
deallocate (psi_ref_zeroth_order)
deallocate (psi_ref_dressed)
deallocate (eigvalues)
deallocate (eigvectors)
end

90
plugins/OVB_effective_Hamiltonian/save_wf_only_ionic_and_1p_amplitudes.irp.f Normal file
View File

@ -0,0 +1,90 @@
program save_wf
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
use bitmasks
integer :: i,j,k,l
integer(bit_kind), allocatable :: psi_save_final(:,:,:)
double precision, allocatable :: psi_coef_save_final(:,:)
integer :: index_ref_determinants_save(psi_det_size)
integer :: n_det_ref_determinants_save
integer :: index_non_ref_determinants_save(psi_det_size)
integer :: n_det_non_ref_determinants_save
integer :: n_det_save_final
integer :: number_of_particles
n_det_ref_determinants_save = 0
integer :: ionic_level
ionic_level = 1
do i = 1, ionic_index(ionic_level,0) ! number of determinants in the ref wf that are neutrals
n_det_ref_determinants_save +=1
index_ref_determinants_save(n_det_ref_determinants_save) = ionic_index(ionic_level,i)
enddo
! save all the 1p determinants in order to have the single excitations
! on the top of the neutral structures
n_det_non_ref_determinants_save = 0
do i = 1, N_det_non_ref
if(number_of_particles(psi_non_ref(1,1,i))==1)then
n_det_non_ref_determinants_save +=1
index_non_ref_determinants_save(n_det_non_ref_determinants_save) = i
endif
enddo
print*,'n_det_ref_determinants_save = ',n_det_ref_determinants_save
print*,'n_det_non_ref_determinants_save = ',n_det_non_ref_determinants_save
n_det_save_final = n_det_ref_determinants_save + n_det_non_ref_determinants_save
allocate (psi_save_final(N_int,2,n_det_save_final))
allocate (psi_coef_save_final(n_det_save_final,1))
integer :: n_det_tmp
n_det_tmp = 0
do i = 1, n_det_ref_determinants_save ! set the CAS determinants to psi_save_final
n_det_tmp +=1
do j = 1, N_int
psi_save_final(j,1,n_det_tmp) = psi_ref(j,1,index_ref_determinants_save(i))
psi_save_final(j,2,n_det_tmp) = psi_ref(j,2,index_ref_determinants_save(i))
enddo
psi_coef_save_final(n_det_tmp,1) = psi_ref_coef(index_ref_determinants_save(i),1)
enddo
pause
do i = 1, n_det_non_ref_determinants_save ! set the non ref determinants to psi_save_final
n_det_tmp +=1
do j = 1, N_int
psi_save_final(j,1,n_det_tmp) = psi_non_ref(j,1,index_non_ref_determinants_save(i))
psi_save_final(j,2,n_det_tmp) = psi_non_ref(j,2,index_non_ref_determinants_save(i))
enddo
accu = 0.d0
double precision :: t_ik,hij
do j = 1, n_det_ref_determinants_save
call i_H_j(psi_non_ref(1,1,index_non_ref_determinants_save(i)),psi_ref(1,1,index_ref_determinants_save(j)),N_int,hij)
t_ik = hij * lambda_mrcc(1,index_non_ref_determinants_save(i))
accu += psi_ref_coef(index_ref_determinants_save(j),1) * t_ik
enddo
psi_coef_save_final(n_det_tmp,1) = accu
enddo
double precision :: accu
accu = 0.d0
do i = 1, n_det_save_final
accu += psi_coef_save_final(i,1) * psi_coef_save_final(i,1)
enddo
accu = 1.d0/dsqrt(accu)
do i = 1, n_det_save_final
psi_coef_save_final(i,1) = accu * psi_coef_save_final(i,1)
enddo
do i = 1, n_det_save_final
print*,''
print*,'Det'
call debug_det(psi_save_final(1,1,i),N_int)
print*,'coef = ',psi_coef_save_final(i,1)
enddo
call save_wavefunction_general(n_det_save_final,1,psi_save_final,n_det_save_final,psi_coef_save_final)
deallocate (psi_save_final)
deallocate (psi_coef_save_final)
end

88
plugins/OVB_effective_Hamiltonian/save_wf_only_neutral_and_1p_amplitudes.irp.f Normal file
View File

@ -0,0 +1,88 @@
program save_wf
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
use bitmasks
integer :: i,j,k,l
integer(bit_kind), allocatable :: psi_save_final(:,:,:)
double precision, allocatable :: psi_coef_save_final(:,:)
integer :: index_ref_determinants_save(psi_det_size)
integer :: n_det_ref_determinants_save
integer :: index_non_ref_determinants_save(psi_det_size)
integer :: n_det_non_ref_determinants_save
integer :: n_det_save_final
integer :: number_of_particles
n_det_ref_determinants_save = 0
do i = 1, ionic_index(0,0) ! number of determinants in the ref wf that are neutrals
n_det_ref_determinants_save +=1
index_ref_determinants_save(n_det_ref_determinants_save) = ionic_index(0,i)
enddo
! save all the 1p determinants in order to have the single excitations
! on the top of the neutral structures
n_det_non_ref_determinants_save = 0
do i = 1, N_det_non_ref
if(number_of_particles(psi_non_ref(1,1,i))==1)then
n_det_non_ref_determinants_save +=1
index_non_ref_determinants_save(n_det_non_ref_determinants_save) = i
endif
enddo
print*,'n_det_ref_determinants_save = ',n_det_ref_determinants_save
print*,'n_det_non_ref_determinants_save = ',n_det_non_ref_determinants_save
n_det_save_final = n_det_ref_determinants_save + n_det_non_ref_determinants_save
allocate (psi_save_final(N_int,2,n_det_save_final))
allocate (psi_coef_save_final(n_det_save_final,1))
integer :: n_det_tmp
n_det_tmp = 0
do i = 1, n_det_ref_determinants_save ! set the CAS determinants to psi_save_final
n_det_tmp +=1
do j = 1, N_int
psi_save_final(j,1,n_det_tmp) = psi_ref(j,1,index_ref_determinants_save(i))
psi_save_final(j,2,n_det_tmp) = psi_ref(j,2,index_ref_determinants_save(i))
enddo
psi_coef_save_final(n_det_tmp,1) = psi_ref_coef(index_ref_determinants_save(i),1)
enddo
pause
do i = 1, n_det_non_ref_determinants_save ! set the non ref determinants to psi_save_final
n_det_tmp +=1
do j = 1, N_int
psi_save_final(j,1,n_det_tmp) = psi_non_ref(j,1,index_non_ref_determinants_save(i))
psi_save_final(j,2,n_det_tmp) = psi_non_ref(j,2,index_non_ref_determinants_save(i))
enddo
accu = 0.d0
double precision :: t_ik,hij
do j = 1, n_det_ref_determinants_save
call i_H_j(psi_non_ref(1,1,index_non_ref_determinants_save(i)),psi_ref(1,1,index_ref_determinants_save(j)),N_int,hij)
t_ik = hij * lambda_mrcc(1,index_non_ref_determinants_save(i))
accu += psi_ref_coef(index_ref_determinants_save(j),1) * t_ik
enddo
psi_coef_save_final(n_det_tmp,1) = accu
enddo
double precision :: accu
accu = 0.d0
do i = 1, n_det_save_final
accu += psi_coef_save_final(i,1) * psi_coef_save_final(i,1)
enddo
accu = 1.d0/dsqrt(accu)
do i = 1, n_det_save_final
psi_coef_save_final(i,1) = accu * psi_coef_save_final(i,1)
enddo
do i = 1, n_det_save_final
print*,''
print*,'Det'
call debug_det(psi_save_final(1,1,i),N_int)
print*,'coef = ',psi_coef_save_final(i,1)
enddo
call save_wavefunction_general(n_det_save_final,1,psi_save_final,n_det_save_final,psi_coef_save_final)
deallocate (psi_save_final)
deallocate (psi_coef_save_final)
end

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

297
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)
@ -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 !
cmoref(21,2,1) = 1.d0 !
cmoref(30,2,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(39,3,1) = 1.d0 !
cmoref(48,3,1) = 1.d0 !
! cmoref(5+shift,1,1) = -0.1d0 ! 2pX function
! cmoref(6+shift,1,1) = -0.1d0 ! 2pZ function
cmoref(3,4,1) = 1.d0 !
cmoref(12,4,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(21,5,1) = 1.d0 !
cmoref(30,5,1) =-1.d0 !
! cmoref(5+shift,1,1) = 0.1d0 ! 2pX function
! cmoref(6+shift,1,1) = 0.1d0 ! 2pZ function
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)
do i = 1, ao_num
do j = 1, ao_num
s(i,j,1) = ao_overlap(i,j)
enddo
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

48
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
@ -44,7 +49,7 @@ class H_apply(object):
self.selection_pt2 = None
self.perturbation = None
self.do_double_exc = do_double_exc
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
s["omp_parallel"] = """ PROVIDE elec_num_tab
!$OMP PARALLEL DEFAULT(SHARED) &
@ -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"] = """
@ -216,10 +247,23 @@ class H_apply(object):
self.data["initialization"] = """
PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
"""
if self.do_double_exc == True:
self.data["keys_work"] = """
! 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)
"""%(pert,)
! 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)
else:
self.data["keys_work"] = """
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)
"""%(pert)
self.data["finalization"] = """
"""
self.data["copy_buffer"] = ""

53
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,19 +414,20 @@ 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
print*,'n_core_orb = ',n_core_orb
END_PROVIDER
@ -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

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

@ -189,6 +189,39 @@ logical function is_connected_to(key,keys,Nint,Ndet)
enddo
end
logical function is_connected_to_by_mono(key,keys,Nint,Ndet)
use bitmasks
implicit none
integer, intent(in) :: Nint, Ndet
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
integer :: i, l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
is_connected_to_by_mono = .false.
do i=1,Ndet
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
popcnt(xor( key(l,2), keys(l,2,i)))
enddo
if (degree_x2 > 2) then
cycle
else
is_connected_to_by_mono = .true.
return
endif
enddo
end
integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
use bitmasks

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