LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-04-08 17:51:43 +02:00
Code Issues Releases Wiki Activity

FOBO-SCF executable works

Browse Source
This commit is contained in:
Emmanuel Giner 2016-03-11 23:27:39 +01:00
parent 518520d619
commit ac8e530372
43 changed files with 3082 additions and 499 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -51,7 +51,7 @@ FCFLAGS : -xSSE4.2 -O2 -ip -ftz
#
[DEBUG]
FC : -g -traceback
FCFLAGS : -xSSE2 -C -fpe0
FCFLAGS : -xSSE2 -C
IRPF90_FLAGS : --openmp
# OpenMP flags

3
plugins/All_singles/H_apply.irp.f
View File

@ -8,10 +8,9 @@ s.unset_skip()
s.filter_only_1h1p()
print s
s = H_apply("just_mono")
s = H_apply("just_mono",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.unset_double_excitations()
print s
END_SHELL

2
plugins/All_singles/all_singles.irp.f
View File

@ -15,7 +15,7 @@ subroutine routine
integer :: N_st, degree
double precision,allocatable :: E_before(:)
integer :: n_det_before
N_st = N_states
N_st = N_states_diag
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
i = 0
print*,'N_det = ',N_det

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

@ -20,22 +20,18 @@ 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

23
plugins/DDCI_selected/ddci.irp.f
View File

@ -3,10 +3,10 @@ program ddci
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:)
integer :: N_st, degree
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
N_st = N_states_diag
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
character*(64) :: perturbation
pt2 = 1.d0
@ -27,6 +27,8 @@ program ddci
print *, 'E+PT2 = ', CI_energy+pt2
print *, '-----'
endif
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
call H_apply_DDCI_selection(pt2, norm_pert, H_pert_diag, N_st)
@ -47,8 +49,21 @@ program ddci
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', CI_energy
print *, 'E+PT2 = ', CI_energy+pt2
print *, 'E+PT2 = ', E_before+pt2
print *, '-----'
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
call ezfio_set_ddci_selected_energy(CI_energy)
if (abort_all) then
exit

16
plugins/FOBOCI/H_apply.irp.f
View File

@ -18,8 +18,22 @@ print s
s = H_apply("standard")
s = H_apply("only_1h2p")
s.set_selection_pt2("epstein_nesbet")
s.filter_only_1h2p()
s.unset_skip()
print s
s = H_apply("only_2h2p")
s.set_selection_pt2("epstein_nesbet")
s.filter_only_2h2p()
s.unset_skip()
print s
s = H_apply("only_2p")
s.set_selection_pt2("epstein_nesbet")
s.filter_only_2p()
s.unset_skip()
print s

4
plugins/FOBOCI/H_apply_dressed_autonom.irp.f
View File

@ -449,8 +449,8 @@ subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_g
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(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
integer :: i_generator, nmax

385
plugins/FOBOCI/H_apply_dressed_autonom_bis.irp.f Normal file
View File

@ -0,0 +1,385 @@
subroutine H_apply_dressed_pert_monoexc_bis(key_in, hole_1,particl_1,i_generator,iproc_in , delta_ij_generators_, Ndet_generators,psi_det_generators_input,E_ref,n_det_input,psi_det_input )
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,n_det_input
double precision, intent(inout) :: delta_ij_generators_(n_det_input,n_det_input),E_ref
integer(bit_kind), intent(in) :: psi_det_input(N_int,2,n_det_input)
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_bis(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref,psi_det_input,n_det_input)
key_idx = 0
endif
enddo ! ii
! !$OMP ENDDO NOWAIT
enddo ! ispin
call standard_dress_bis(delta_ij_generators_,size_max,Ndet_generators,i_generator,key_idx,keys_out,N_int,iproc,psi_det_generators_input,E_ref,psi_det_input,n_det_input)
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_pertk_single(delta_ij_, Ndet_generators,psi_det_generators_input,E_ref,psi_det_input,n_det_input)
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,n_det_input
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
integer(bit_kind), intent(in) :: psi_det_input(N_int,2,n_det_input)
double precision, intent(inout) :: delta_ij_(n_det_input,n_det_input),E_ref
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
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_, Ndet_generators,psi_det_generators_input,E_ref,n_det_input,psi_det_input)
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(.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_, Ndet_generators,psi_det_generators_input,E_ref,n_det_input,psi_det_input)
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
subroutine standard_dress_bis(delta_ij_generators_,size_buffer,Ndet_generators,i_generator,n_selected,det_buffer,Nint,iproc,psi_det_generators_input,E_ref,psi_det_input,n_det_input)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint, iproc,n_det_input
integer, intent(in) :: Ndet_generators,size_buffer
double precision, intent(inout) :: delta_ij_generators_(n_det_input,n_det_input),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(bit_kind), intent(in) :: psi_det_input(N_int,2,n_det_input)
integer :: i,j,k,m
integer :: new_size
integer :: degree(n_det_input)
integer :: idx(0:n_det_input)
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(n_det_input)
double precision :: contrib,lambda_i,accu
integer :: number_of_holes,n_h, number_of_particles,n_p
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
print*
n_h = number_of_holes(det_buffer(1,1,i))
n_p = number_of_particles(det_buffer(1,1,i))
print*,'n_h,n_p = ',n_h,n_p
call get_excitation_degree_vector(psi_det_input,det_buffer(1,1,i),degree,N_int,n_det_input,idx)
H_array = 0.d0
do k=1,idx(0)
call i_h_j(det_buffer(1,1,i),psi_det_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)
lambda_i = f
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
enddo
end

2
plugins/FOBOCI/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Generators_restart Selectors_no_sorted
Perturbation Selectors_no_sorted Hartree_Fock

333
plugins/FOBOCI/all_singles.irp.f
View File

@ -8,7 +8,7 @@ subroutine all_single
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
threshold_davidson = 1.d-9
soft_touch threshold_davidson davidson_criterion
i = 0
print*,'Doing all the mono excitations !'
@ -52,10 +52,12 @@ subroutine all_single
enddo
endif
E_before = CI_energy
!!!!!!!!!!!!!!!!!!!!!!!!!!! DOING ONLY ONE ITERATION OF SELECTION AS THE SELECTION CRITERION IS SET TO ZERO
exit
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
! threshold_davidson = 1.d-8
! soft_touch threshold_davidson davidson_criterion
! call diagonalize_CI
print*,'Final Step '
print*,'N_det = ',N_det
do i = 1, N_states_diag
@ -67,10 +69,250 @@ subroutine all_single
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
! call save_wavefunction
deallocate(pt2,norm_pert,E_before)
end
subroutine all_1h2p
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 1h2P 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_only_1h2p(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
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
deallocate(pt2,norm_pert,E_before)
end
subroutine all_2h2p
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 2h2P 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_only_2h2p(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
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
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_only_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
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
deallocate(pt2,norm_pert,E_before)
do i = 1, 2
print*,'psi_coef = ',psi_coef(i,1)
enddo
end
subroutine all_single_no_1h_or_1p
implicit none
integer :: i,k
@ -126,7 +368,7 @@ subroutine all_single_no_1h_or_1p
endif
E_before = CI_energy
enddo
threshold_davidson = 1.d-10
threshold_davidson = 1.d-16
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
print*,'Final Step '
@ -217,85 +459,6 @@ subroutine all_single_no_1h_or_1p_or_2p
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

452
plugins/FOBOCI/all_singles_split.irp.f
View File

@ -5,7 +5,7 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N
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
integer :: i,i_hole,j
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
do i = 1, n_inact_orb
@ -22,16 +22,170 @@ subroutine all_single_split(psi_det_generators_input,psi_coef_generators_input,N
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 diagonalize_CI_SC2
! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2)
call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
enddo
do i = 1, n_act_orb
i_hole = list_act(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
! call diagonalize_CI_SC2
! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2)
call provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det_generators_input)
enddo
do i = 1, n_virt_orb
i_hole = list_virt(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
! call diagonalize_CI_SC2
! call update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,Diag_H_elements_SC2)
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_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p)
implicit none
use bitmasks
integer, intent(in) :: i_particl
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)
double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
integer :: i,j
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
call all_single
threshold_davidson = 1.d-12
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
double precision, allocatable :: matrix_ref_1h_1p(:,:)
double precision, allocatable :: matrix_ref_1h_1p_dressing_1h1p(:,:)
double precision, allocatable :: matrix_ref_1h_1p_dressing_1h2p(:,:)
double precision, allocatable :: psi_coef_ref_1h_1p(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_1h2p(:,:)
integer(bit_kind), allocatable :: psi_det_1h2p(:,:,:)
integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:)
integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:)
integer :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
double precision :: hka
double precision,allocatable :: eigenvectors(:,:), eigenvalues(:)
call give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p)
allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
allocate(matrix_ref_1h_1p_dressing_1h2p(n_det_ref_1h_1p,n_det_ref_1h_1p))
allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states))
allocate(psi_det_1h2p(N_int,2,n_det_1h2p), psi_coef_1h2p(n_det_1h2p,N_states))
allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states))
call give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p)
do i = 1, n_det_ref_1h_1p
do j = 1, n_det_ref_1h_1p
call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka)
matrix_ref_1h_1p(i,j) = hka
enddo
enddo
matrix_ref_1h_1p_dressing_1h1p = 0.d0
matrix_ref_1h_1p_dressing_1h2p = 0.d0
call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h2p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
psi_det_1h2p,psi_coef_1h2p,n_det_1h2p)
call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
psi_det_1h1p,psi_coef_1h1p,n_det_1h1p)
do i = 1, n_det_ref_1h_1p
do j = 1, n_det_ref_1h_1p
matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_1h2p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j)
enddo
enddo
allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p))
call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p)
!do j = 1, n_det_ref_1h_1p
! print*,'coef = ',eigenvectors(j,1)
!enddo
print*,''
print*,'-----------------------'
print*,'-----------------------'
print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion
print*,'-----------------------'
! Extract the
integer, allocatable :: index_generator(:)
integer :: n_det_generators_tmp,degree
n_det_generators_tmp = 0
allocate(index_generator(n_det_ref_1h_1p))
do i = 1, n_det_ref_1h_1p
do j = 1, N_det_generators
call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(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_generators_tmp .ne. n_det_generators)then
print*,'PB !!!'
print*,'if(n_det_generators_tmp .ne. n_det_genrators)then'
stop
endif
do i = 1, N_det_generators
print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1)
do j = 1, N_det_generators
dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j))
dressing_matrix_1h2p(i,j) += matrix_ref_1h_1p_dressing_1h2p(index_generator(i),index_generator(j))
enddo
enddo
print*,'-----------------------'
print*,'-----------------------'
deallocate(matrix_ref_1h_1p)
deallocate(matrix_ref_1h_1p_dressing_1h1p)
deallocate(matrix_ref_1h_1p_dressing_1h2p)
deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p)
deallocate(psi_det_1h2p, psi_coef_1h2p)
deallocate(psi_det_1h1p, psi_coef_1h1p)
deallocate(eigenvectors,eigenvalues)
deallocate(index_generator)
end
subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dressing_matrix_extra_1h_or_1p)
implicit none
use bitmasks
@ -39,50 +193,168 @@ subroutine all_single_for_1h(i_hole,dressing_matrix_1h1p,dressing_matrix_2h1p,dr
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)
double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
integer :: i
integer :: i,j
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
integer :: n_det_1h1p,n_det_2h1p,n_det_extra_1h_or_1p
integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
integer(bit_kind), allocatable :: psi_2h1p(:,:,:)
integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:)
double precision, allocatable :: psi_ref_coef_out(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_2h1p(:,:)
double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:)
!call all_single_no_1h_or_1p
call all_single
threshold_davidson = 1.d-12
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p))
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))
allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states))
call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
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)
call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p)
deallocate(psi_ref_out)
deallocate(psi_1h1p)
deallocate(psi_2h1p)
deallocate(psi_extra_1h_or_1p)
deallocate(psi_ref_coef_out)
deallocate(psi_coef_1h1p)
deallocate(psi_coef_2h1p)
deallocate(psi_coef_extra_1h_or_1p)
double precision, allocatable :: matrix_ref_1h_1p(:,:)
double precision, allocatable :: matrix_ref_1h_1p_dressing_1h1p(:,:)
double precision, allocatable :: matrix_ref_1h_1p_dressing_2h1p(:,:)
double precision, allocatable :: psi_coef_ref_1h_1p(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_2h1p(:,:)
integer(bit_kind), allocatable :: psi_det_2h1p(:,:,:)
integer(bit_kind), allocatable :: psi_det_ref_1h_1p(:,:,:)
integer(bit_kind), allocatable :: psi_det_1h1p(:,:,:)
integer :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
double precision :: hka
double precision,allocatable :: eigenvectors(:,:), eigenvalues(:)
call give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p)
allocate(matrix_ref_1h_1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
allocate(matrix_ref_1h_1p_dressing_1h1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
allocate(matrix_ref_1h_1p_dressing_2h1p(n_det_ref_1h_1p,n_det_ref_1h_1p))
allocate(psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p), psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states))
allocate(psi_det_2h1p(N_int,2,n_det_2h1p), psi_coef_2h1p(n_det_2h1p,N_states))
allocate(psi_det_1h1p(N_int,2,n_det_1h1p), psi_coef_1h1p(n_det_1h1p,N_states))
call give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p)
do i = 1, n_det_ref_1h_1p
do j = 1, n_det_ref_1h_1p
call i_h_j(psi_det_ref_1h_1p(1,1,i),psi_det_ref_1h_1p(1,1,j),N_int,hka)
matrix_ref_1h_1p(i,j) = hka
enddo
enddo
matrix_ref_1h_1p_dressing_1h1p = 0.d0
matrix_ref_1h_1p_dressing_2h1p = 0.d0
call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_2h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
psi_det_2h1p,psi_coef_2h1p,n_det_2h1p)
call provide_matrix_dressing_general(matrix_ref_1h_1p_dressing_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,n_det_ref_1h_1p, &
psi_det_1h1p,psi_coef_1h1p,n_det_1h1p)
do i = 1, n_det_ref_1h_1p
do j = 1, n_det_ref_1h_1p
matrix_ref_1h_1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(i,j) + matrix_ref_1h_1p_dressing_1h1p(i,j)
enddo
enddo
allocate(eigenvectors(n_det_ref_1h_1p,n_det_ref_1h_1p), eigenvalues(n_det_ref_1h_1p))
call lapack_diag(eigenvalues,eigenvectors,matrix_ref_1h_1p,n_det_ref_1h_1p,n_det_ref_1h_1p)
!do j = 1, n_det_ref_1h_1p
! print*,'coef = ',eigenvectors(j,1)
!enddo
print*,''
print*,'-----------------------'
print*,'-----------------------'
print*,'e_dressed = ',eigenvalues(1)+nuclear_repulsion
print*,'-----------------------'
! Extract the
integer, allocatable :: index_generator(:)
integer :: n_det_generators_tmp,degree
n_det_generators_tmp = 0
allocate(index_generator(n_det_ref_1h_1p))
do i = 1, n_det_ref_1h_1p
do j = 1, N_det_generators
call get_excitation_degree(psi_det_generators(1,1,j),psi_det_ref_1h_1p(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_generators_tmp .ne. n_det_generators)then
print*,'PB !!!'
print*,'if(n_det_generators_tmp .ne. n_det_genrators)then'
stop
endif
do i = 1, N_det_generators
print*,'psi_coef_dressed = ',eigenvectors(index_generator(i),1)
do j = 1, N_det_generators
dressing_matrix_1h1p(i,j) += matrix_ref_1h_1p_dressing_1h1p(index_generator(i),index_generator(j))
dressing_matrix_2h1p(i,j) += matrix_ref_1h_1p_dressing_2h1p(index_generator(i),index_generator(j))
enddo
enddo
print*,'-----------------------'
print*,'-----------------------'
deallocate(matrix_ref_1h_1p)
deallocate(matrix_ref_1h_1p_dressing_1h1p)
deallocate(matrix_ref_1h_1p_dressing_2h1p)
deallocate(psi_det_ref_1h_1p, psi_coef_ref_1h_1p)
deallocate(psi_det_2h1p, psi_coef_2h1p)
deallocate(psi_det_1h1p, psi_coef_1h1p)
deallocate(eigenvectors,eigenvalues)
deallocate(index_generator)
!return
!
!integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
!integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
!integer(bit_kind), allocatable :: psi_2h1p(:,:,:)
!integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:)
!double precision, allocatable :: psi_ref_coef_out(:,:)
!double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:)
!call all_single_no_1h_or_1p
!call give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p))
!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))
!allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states))
!call split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_det_1h1p,n_det_2h1p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_2h1p,psi_coef_2h1p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
!do i = 1, n_det_extra_1h_or_1p
! print*,'----'
! print*,'c = ',psi_coef_extra_1h_or_1p(i,1)
! call debug_det(psi_extra_1h_or_1p(1,1,i),N_int)
! print*,'----'
!enddo
!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)
!print*,'Dressing 1h1p '
!do j =1, N_det_generators
! print*,' dressing ',dressing_matrix_1h1p(j,:)
!enddo
!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)
!print*,'Dressing 2h1p '
!do j =1, N_det_generators
! print*,' dressing ',dressing_matrix_2h1p(j,:)
!enddo
!call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
! psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p)
!print*,',dressing_matrix_extra_1h_or_1p'
!do j =1, N_det_generators
! print*,' dressing ',dressing_matrix_extra_1h_or_1p(j,:)
!enddo
!deallocate(psi_ref_out)
!deallocate(psi_1h1p)
!deallocate(psi_2h1p)
!deallocate(psi_extra_1h_or_1p)
!deallocate(psi_ref_coef_out)
!deallocate(psi_coef_1h1p)
!deallocate(psi_coef_2h1p)
!deallocate(psi_coef_extra_1h_or_1p)
end
@ -208,56 +480,56 @@ subroutine all_single_split_for_1p(dressing_matrix_1h1p,dressing_matrix_1h2p)
soft_touch n_det_max_jacobi
end
subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p)
implicit none
use bitmasks
integer, intent(in ) :: i_particl
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)
double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
integer :: i
n_det_max_jacobi = 50
soft_touch n_det_max_jacobi
integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p
integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
integer(bit_kind), allocatable :: psi_1h2p(:,:,:)
integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:)
double precision, allocatable :: psi_ref_coef_out(:,:)
double precision, allocatable :: psi_coef_1h1p(:,:)
double precision, allocatable :: psi_coef_1h2p(:,:)
double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:)
!call all_single_no_1h_or_1p_or_2p
call all_single
threshold_davidson = 1.d-12
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p))
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))
allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states))
call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
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)
call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p)
deallocate(psi_ref_out)
deallocate(psi_1h1p)
deallocate(psi_1h2p)
deallocate(psi_ref_coef_out)
deallocate(psi_coef_1h1p)
deallocate(psi_coef_1h2p)
end
! subroutine all_single_for_1p(i_particl,dressing_matrix_1h1p,dressing_matrix_1h2p,dressing_matrix_extra_1h_or_1p)
! implicit none
! use bitmasks
! integer, intent(in ) :: i_particl
! 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)
! double precision, intent(inout) :: dressing_matrix_extra_1h_or_1p(N_det_generators,N_det_generators)
! integer :: i
! n_det_max_jacobi = 50
! soft_touch n_det_max_jacobi
!
! integer :: n_det_1h1p,n_det_1h2p,n_det_extra_1h_or_1p
! integer(bit_kind), allocatable :: psi_ref_out(:,:,:)
! integer(bit_kind), allocatable :: psi_1h1p(:,:,:)
! integer(bit_kind), allocatable :: psi_1h2p(:,:,:)
! integer(bit_kind), allocatable :: psi_extra_1h_or_1p(:,:,:)
! double precision, allocatable :: psi_ref_coef_out(:,:)
! double precision, allocatable :: psi_coef_1h1p(:,:)
! double precision, allocatable :: psi_coef_1h2p(:,:)
! double precision, allocatable :: psi_coef_extra_1h_or_1p(:,:)
!!!!call all_single_no_1h_or_1p_or_2p
! call all_single
!
! threshold_davidson = 1.d-12
! soft_touch threshold_davidson davidson_criterion
! call diagonalize_CI
! call give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,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_extra_1h_or_1p(N_int,2,n_det_extra_1h_or_1p))
! 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))
! allocate(psi_coef_extra_1h_or_1p(n_det_extra_1h_or_1p,N_states))
! call split_wf_generators_and_1h1p_and_1h2p(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p,psi_ref_out,psi_ref_coef_out,psi_1h1p,psi_coef_1h1p,psi_1h2p,psi_coef_1h2p,psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p)
! 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)
! call provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix_extra_1h_or_1p,psi_ref_out,psi_ref_coef_out,N_det_generators, &
! psi_extra_1h_or_1p,psi_coef_extra_1h_or_1p,n_det_extra_1h_or_1p)
!
! deallocate(psi_ref_out)
! deallocate(psi_1h1p)
! deallocate(psi_1h2p)
! deallocate(psi_ref_coef_out)
! deallocate(psi_coef_1h1p)
! deallocate(psi_coef_1h2p)
!
! end

436
plugins/FOBOCI/collect_all_lmct.irp.f Normal file
View File

@ -0,0 +1,436 @@
use bitmasks
subroutine collect_lmct(hole_particle,n_couples)
implicit none
integer, intent(out) :: hole_particle(1000,2), n_couples
BEGIN_DOC
! Collect all the couple holes/particles of the important LMCT
! hole_particle(i,1) = ith hole
! hole_particle(i,2) = ith particle
! n_couples is the number of important excitations
END_DOC
print*,'COLLECTING THE PERTINENT LMCT (1h)'
double precision, allocatable :: tmp(:,:)
allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
integer :: i,j,iorb,jorb
n_couples = 0
do i = 1,n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.1.d-2)then
n_couples +=1
hole_particle(n_couples,1) = jorb
hole_particle(n_couples,2) = iorb
print*,'DM'
print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
endif
enddo
enddo
deallocate(tmp)
print*,'number of meaning full couples of holes/particles '
print*,'n_couples = ',n_couples
end
subroutine collect_mlct(hole_particle,n_couples)
implicit none
integer, intent(out) :: hole_particle(1000,2), n_couples
BEGIN_DOC
! Collect all the couple holes/particles of the important LMCT
! hole_particle(i,1) = ith hole
! hole_particle(i,2) = ith particle
! n_couples is the number of important excitations
END_DOC
print*,'COLLECTING THE PERTINENT MLCT (1p)'
double precision, allocatable :: tmp(:,:)
allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
integer :: i,j,iorb,jorb
n_couples = 0
do i = 1,n_act_orb
iorb = list_act(i)
do j = 1, n_virt_orb
jorb = list_virt(j)
if(dabs(tmp(iorb,jorb)).gt.1.d-3)then
n_couples +=1
hole_particle(n_couples,1) = iorb
hole_particle(n_couples,2) = jorb
print*,'DM'
print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
endif
enddo
enddo
deallocate(tmp)
print*,'number of meaning full couples of holes/particles '
print*,'n_couples = ',n_couples
end
subroutine collect_lmct_mlct(hole_particle,n_couples)
implicit none
integer, intent(out) :: hole_particle(1000,2), n_couples
BEGIN_DOC
! Collect all the couple holes/particles of the important LMCT
! hole_particle(i,1) = ith hole
! hole_particle(i,2) = ith particle
! n_couples is the number of important excitations
END_DOC
double precision, allocatable :: tmp(:,:)
print*,'COLLECTING THE PERTINENT LMCT (1h)'
print*,'AND THE PERTINENT MLCT (1p)'
allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
integer :: i,j,iorb,jorb
n_couples = 0
do i = 1,n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
n_couples +=1
hole_particle(n_couples,1) = jorb
hole_particle(n_couples,2) = iorb
print*,'DM'
print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
endif
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
n_couples +=1
hole_particle(n_couples,1) = iorb
hole_particle(n_couples,2) = jorb
print*,'DM'
print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
endif
enddo
enddo
deallocate(tmp)
print*,'number of meaning full couples of holes/particles '
print*,'n_couples = ',n_couples
end
subroutine collect_1h1p(hole_particle,n_couples)
implicit none
integer, intent(out) :: hole_particle(1000,2), n_couples
BEGIN_DOC
! Collect all the couple holes/particles of the important LMCT
! hole_particle(i,1) = ith hole
! hole_particle(i,2) = ith particle
! n_couples is the number of important excitations
END_DOC
double precision, allocatable :: tmp(:,:)
print*,'COLLECTING THE PERTINENT 1h1p'
allocate(tmp(size(one_body_dm_mo_alpha_osoci,1),size(one_body_dm_mo_alpha_osoci,2)))
tmp = one_body_dm_mo_alpha_osoci + one_body_dm_mo_beta_osoci
integer :: i,j,iorb,jorb
n_couples = 0
do i = 1,n_virt_orb
iorb = list_virt(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.1.d-2)then
n_couples +=1
hole_particle(n_couples,1) = jorb
hole_particle(n_couples,2) = iorb
print*,'DM'
print*,hole_particle(n_couples,1),hole_particle(n_couples,2),tmp(iorb,jorb)
endif
enddo
enddo
deallocate(tmp)
print*,'number of meaning full couples of holes/particles '
print*,'n_couples = ',n_couples
end
subroutine set_lmct_to_generators_restart
implicit none
integer :: i,j,m,n,i_hole,i_particle
integer :: hole_particle(1000,2), n_couples
integer(bit_kind) :: key_tmp(N_int,2)
integer :: N_det_total,i_ok
call collect_lmct(hole_particle,n_couples)
call set_generators_to_generators_restart
N_det_total = N_det_generators_restart
do i = 1, n_couples
i_hole = hole_particle(i,1)
i_particle = hole_particle(i,2)
do m = 1, N_det_cas
do n = 1, N_int
key_tmp(n,1) = psi_cas(n,1,m)
key_tmp(n,2) = psi_cas(n,2,m)
enddo
! You excite the beta electron from i_hole to i_particle
print*,'i_hole,i_particle 2 = ',i_hole,i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
print*,'i_ok = ',i_ok
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
do n = 1, N_int
key_tmp(n,1) = psi_cas(n,1,m)
key_tmp(n,2) = psi_cas(n,2,m)
enddo
! You excite the alpha electron from i_hole to i_particle
print*,'i_hole,i_particle 1 = ',i_hole,i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
print*,'i_ok = ',i_ok
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
enddo
enddo
N_det_generators = N_det_total
do i = 1, N_det_generators
psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total))
enddo
print*,'number of generators in total = ',N_det_generators
touch N_det_generators psi_coef_generators psi_det_generators
end
subroutine set_mlct_to_generators_restart
implicit none
integer :: i,j,m,n,i_hole,i_particle
integer :: hole_particle(1000,2), n_couples
integer(bit_kind) :: key_tmp(N_int,2)
integer :: N_det_total,i_ok
call collect_mlct(hole_particle,n_couples)
call set_generators_to_generators_restart
N_det_total = N_det_generators_restart
do i = 1, n_couples
i_hole = hole_particle(i,1)
i_particle = hole_particle(i,2)
do m = 1, N_det_cas
do n = 1, N_int
key_tmp(n,1) = psi_cas(n,1,m)
key_tmp(n,2) = psi_cas(n,2,m)
enddo
! You excite the beta electron from i_hole to i_particle
print*,'i_hole,i_particle 2 = ',i_hole,i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
print*,'i_ok = ',i_ok
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
do n = 1, N_int
key_tmp(n,1) = psi_cas(n,1,m)
key_tmp(n,2) = psi_cas(n,2,m)
enddo
! You excite the alpha electron from i_hole to i_particle
print*,'i_hole,i_particle 1 = ',i_hole,i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
print*,'i_ok = ',i_ok
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
enddo
enddo
N_det_generators = N_det_total
do i = 1, N_det_generators
psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total))
enddo
print*,'number of generators in total = ',N_det_generators
touch N_det_generators psi_coef_generators psi_det_generators
end
subroutine set_lmct_mlct_to_generators_restart
implicit none
integer :: i,j,m,n,i_hole,i_particle
integer :: hole_particle(1000,2), n_couples
integer(bit_kind) :: key_tmp(N_int,2)
integer :: N_det_total,i_ok
call collect_lmct_mlct(hole_particle,n_couples)
call set_generators_to_generators_restart
N_det_total = N_det_generators_restart
do i = 1, n_couples
i_hole = hole_particle(i,1)
i_particle = hole_particle(i,2)
do m = 1, N_det_cas
do n = 1, N_int
key_tmp(n,1) = psi_cas(n,1,m)
key_tmp(n,2) = psi_cas(n,2,m)
enddo
! You excite the beta electron from i_hole to i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
do n = 1, N_int
key_tmp(n,1) = psi_cas(n,1,m)
key_tmp(n,2) = psi_cas(n,2,m)
enddo
! You excite the alpha electron from i_hole to i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det_generators(n,1,N_det_total) = key_tmp(n,1)
psi_det_generators(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
enddo
enddo
N_det_generators = N_det_total
do i = 1, N_det_generators
psi_coef_generators(i,1) = 1.d0/dsqrt(dble(N_det_total))
enddo
print*,'number of generators in total = ',N_det_generators
touch N_det_generators psi_coef_generators psi_det_generators
end
subroutine set_lmct_mlct_to_psi_det
implicit none
integer :: i,j,m,n,i_hole,i_particle
integer :: hole_particle(1000,2), n_couples
integer(bit_kind) :: key_tmp(N_int,2)
integer :: N_det_total,i_ok
call collect_lmct_mlct(hole_particle,n_couples)
call set_psi_det_to_generators_restart
N_det_total = N_det_generators_restart
do i = 1, n_couples
i_hole = hole_particle(i,1)
i_particle = hole_particle(i,2)
do m = 1, N_det_generators_restart
do n = 1, N_int
key_tmp(n,1) = psi_det_generators_restart(n,1,m)
key_tmp(n,2) = psi_det_generators_restart(n,2,m)
enddo
! You excite the beta electron from i_hole to i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det(n,1,N_det_total) = key_tmp(n,1)
psi_det(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
do n = 1, N_int
key_tmp(n,1) = psi_det_generators_restart(n,1,m)
key_tmp(n,2) = psi_det_generators_restart(n,2,m)
enddo
! You excite the alpha electron from i_hole to i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det(n,1,N_det_total) = key_tmp(n,1)
psi_det(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
enddo
enddo
N_det = N_det_total
integer :: k
do k = 1, N_states
do i = 1, N_det
psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total))
enddo
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
call remove_duplicates_in_psi_det(found_duplicates)
end
subroutine set_1h1p_to_psi_det
implicit none
integer :: i,j,m,n,i_hole,i_particle
integer :: hole_particle(1000,2), n_couples
integer(bit_kind) :: key_tmp(N_int,2)
integer :: N_det_total,i_ok
call collect_1h1p(hole_particle,n_couples)
call set_psi_det_to_generators_restart
N_det_total = N_det_generators_restart
do i = 1, n_couples
i_hole = hole_particle(i,1)
i_particle = hole_particle(i,2)
do m = 1, N_det_generators_restart
do n = 1, N_int
key_tmp(n,1) = psi_det_generators_restart(n,1,m)
key_tmp(n,2) = psi_det_generators_restart(n,2,m)
enddo
! You excite the beta electron from i_hole to i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,2,i_ok)
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det(n,1,N_det_total) = key_tmp(n,1)
psi_det(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
do n = 1, N_int
key_tmp(n,1) = psi_det_generators_restart(n,1,m)
key_tmp(n,2) = psi_det_generators_restart(n,2,m)
enddo
! You excite the alpha electron from i_hole to i_particle
call do_mono_excitation(key_tmp,i_hole,i_particle,1,i_ok)
if(i_ok==1)then
N_det_total +=1
do n = 1, N_int
psi_det(n,1,N_det_total) = key_tmp(n,1)
psi_det(n,2,N_det_total) = key_tmp(n,2)
enddo
endif
enddo
enddo
N_det = N_det_total
integer :: k
do k = 1, N_states
do i = 1, N_det
psi_coef(i,k) = 1.d0/dsqrt(dble(N_det_total))
enddo
enddo
SOFT_TOUCH N_det psi_det psi_coef
logical :: found_duplicates
call remove_duplicates_in_psi_det(found_duplicates)
end

363
plugins/FOBOCI/corr_energy_2h2p.irp.f Normal file
View File

@ -0,0 +1,363 @@
BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_ab, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_aa, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_bb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_2h2p]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_2h2p = 0.d0
corr_energy_2h2p_per_orb_ab = 0.d0
corr_energy_2h2p_per_orb_aa = 0.d0
corr_energy_2h2p_per_orb_bb = 0.d0
do i = 1, n_inact_orb ! beta
i_hole = list_inact(i)
do k = 1, n_virt_orb ! beta
k_part = list_virt(k)
do j = 1, n_inact_orb ! alpha
j_hole = list_inact(j)
do l = 1, n_virt_orb ! alpha
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = (ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = hij*hij/delta_e
total_corr_e_2h2p += contrib
corr_energy_2h2p_per_orb_ab(i_hole) += contrib
corr_energy_2h2p_per_orb_ab(k_part) += contrib
enddo
enddo
enddo
enddo
! alpha alpha correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_virt_orb
k_part = list_virt(k)
do l = k+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 1
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h2p += contrib
corr_energy_2h2p_per_orb_aa(i_hole) += contrib
corr_energy_2h2p_per_orb_aa(k_part) += contrib
enddo
enddo
enddo
enddo
! beta beta correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_virt_orb
k_part = list_virt(k)
do l = k+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 2
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h2p += contrib
corr_energy_2h2p_per_orb_bb(i_hole) += contrib
corr_energy_2h2p_per_orb_bb(k_part) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_aa, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_bb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_2h1p]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_2h1p = 0.d0
corr_energy_2h1p_per_orb_ab = 0.d0
corr_energy_2h1p_per_orb_aa = 0.d0
corr_energy_2h1p_per_orb_bb = 0.d0
do i = 1, n_inact_orb
i_hole = list_inact(i)
do k = 1, n_act_orb
k_part = list_act(k)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do l = 1, n_virt_orb
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h1p += contrib
corr_energy_2h1p_per_orb_ab(i_hole) += contrib
corr_energy_2h1p_per_orb_ab(l_part) += contrib
enddo
enddo
enddo
enddo
! Alpha Alpha spin correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = 1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 1
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h1p += contrib
corr_energy_2h1p_per_orb_aa(i_hole) += contrib
corr_energy_2h1p_per_orb_aa(l_part) += contrib
enddo
enddo
enddo
enddo
! Beta Beta correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = 1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 2
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h1p += contrib
corr_energy_2h1p_per_orb_bb(i_hole) += contrib
corr_energy_2h1p_per_orb_aa(l_part) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_ab, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_aa, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_bb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_1h2p]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_1h2p = 0.d0
corr_energy_1h2p_per_orb_ab = 0.d0
corr_energy_1h2p_per_orb_aa = 0.d0
corr_energy_1h2p_per_orb_bb = 0.d0
do i = 1, n_virt_orb
i_hole = list_virt(i)
do k = 1, n_act_orb
k_part = list_act(k)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do l = 1, n_virt_orb
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h2p += contrib
corr_energy_1h2p_per_orb_ab(i_hole) += contrib
corr_energy_1h2p_per_orb_ab(j_hole) += contrib
enddo
enddo
enddo
enddo
! Alpha Alpha correlation energy
do i = 1, n_virt_orb
i_hole = list_virt(i)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = i+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 1
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h2p += contrib
corr_energy_1h2p_per_orb_aa(i_hole) += contrib
corr_energy_1h2p_per_orb_ab(j_hole) += contrib
enddo
enddo
enddo
enddo
! Beta Beta correlation energy
do i = 1, n_virt_orb
i_hole = list_virt(i)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = i+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 2
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h2p += contrib
corr_energy_1h2p_per_orb_bb(i_hole) += contrib
corr_energy_1h2p_per_orb_ab(j_hole) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, corr_energy_1h1p_spin_flip_per_orb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_1h1p_spin_flip]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_1h1p_spin_flip = 0.d0
corr_energy_1h1p_spin_flip_per_orb = 0.d0
do i = 1, n_inact_orb
i_hole = list_inact(i)
do k = 1, n_act_orb
k_part = list_act(k)
do j = 1, n_act_orb
j_hole = list_act(j)
do l = 1, n_virt_orb
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h1p_spin_flip += contrib
corr_energy_1h1p_spin_flip_per_orb(i_hole) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER

48
plugins/FOBOCI/diag_fock_inactiv_virt.irp.f
View File

@ -3,6 +3,7 @@ subroutine diag_inactive_virt_and_update_mos
integer :: i,j,i_inact,j_inact,i_virt,j_virt
double precision :: tmp(mo_tot_num_align,mo_tot_num)
character*(64) :: label
print*,'Diagonalizing the occ and virt Fock operator'
tmp = 0.d0
do i = 1, mo_tot_num
tmp(i,i) = Fock_matrix_mo(i,i)
@ -33,3 +34,50 @@ subroutine diag_inactive_virt_and_update_mos
end
subroutine diag_inactive_virt_new_and_update_mos
implicit none
integer :: i,j,i_inact,j_inact,i_virt,j_virt,k,k_act
double precision :: tmp(mo_tot_num_align,mo_tot_num),accu,get_mo_bielec_integral_schwartz
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)
accu =0.d0
do k = 1, n_act_orb
k_act = list_act(k)
accu += get_mo_bielec_integral_schwartz(i_inact,k_act,j_inact,k_act,mo_integrals_map)
accu -= get_mo_bielec_integral_schwartz(i_inact,k_act,k_act,j_inact,mo_integrals_map)
enddo
tmp(i_inact,j_inact) = Fock_matrix_mo(i_inact,j_inact) + accu
tmp(j_inact,i_inact) = Fock_matrix_mo(j_inact,i_inact) + accu
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)
accu =0.d0
do k = 1, n_act_orb
k_act = list_act(k)
accu += get_mo_bielec_integral_schwartz(i_virt,k_act,j_virt,k_act,mo_integrals_map)
enddo
tmp(i_virt,j_virt) = Fock_matrix_mo(i_virt,j_virt) - accu
tmp(j_virt,i_virt) = Fock_matrix_mo(j_virt,i_virt) - accu
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

27
plugins/FOBOCI/dress_simple.irp.f
View File

@ -85,33 +85,6 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen
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

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

@ -0,0 +1,35 @@
program foboscf
implicit none
no_oa_or_av_opt = .True.
touch no_oa_or_av_opt
call run_prepare
call routine_fobo_scf
call save_mos
end
subroutine run_prepare
implicit none
call damping_SCF
call diag_inactive_virt_and_update_mos
end
subroutine routine_fobo_scf
implicit none
integer :: i,j
print*,''
print*,''
character*(64) :: label
label = "Natural"
do i = 1, 5
call FOBOCI_lmct_mlct_old_thr
call save_osoci_natural_mos
call damping_SCF
call diag_inactive_virt_and_update_mos
call clear_mo_map
call provide_properties
enddo
end

97
plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
View File

@ -36,6 +36,10 @@ subroutine FOBOCI_lmct_mlct_old_thr
print*,''
print*,''
print*,'DOING FIRST LMCT !!'
integer(bit_kind) , allocatable :: zero_bitmask(:,:)
integer(bit_kind) , allocatable :: psi_singles(:,:,:)
double precision, allocatable :: psi_singles_coef(:,:)
allocate( zero_bitmask(N_int,2) )
do i = 1, n_inact_orb
integer :: i_hole_osoci
i_hole_osoci = list_inact(i)
@ -54,24 +58,85 @@ subroutine FOBOCI_lmct_mlct_old_thr
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))
dressing_matrix = 0.d0
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)
hkl = dressing_matrix(1,1)
do k = 1, N_det_generators
dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
enddo
print*,'Naked matrix'
do k = 1, N_det_generators
write(*,'(100(F12.5,X))')dressing_matrix(k,:)
enddo
! Do all the single excitations on top of the CAS and 1h determinants
call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask)
call all_single
! ! Change the mask of the holes and particles to perform all the
! ! double excitations that starts from the active space in order
! ! to introduce the Coulomb hole in the active space
! ! These are the 1h2p excitations that have the i_hole_osoci hole in common
! ! and the 2p if there is more than one electron in the active space
! do k = 1, N_int
! zero_bitmask(k,1) = 0_bit_kind
! zero_bitmask(k,2) = 0_bit_kind
! enddo
! ! hole is possible only in the orbital i_hole_osoci
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int)
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int)
! ! and in the active space
! do k = 1, n_act_orb
! call set_bit_to_integer(list_act(k),zero_bitmask(1,1),N_int)
! call set_bit_to_integer(list_act(k),zero_bitmask(1,2),N_int)
! enddo
! call set_bitmask_hole_as_input(zero_bitmask)
! call set_bitmask_particl_as_input(reunion_of_bitmask)
! call all_1h2p
! call diagonalize_CI_SC2
! call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
! ! Change the mask of the holes and particles to perform all the
! ! double excitations that from the orbital i_hole_osoci
! do k = 1, N_int
! zero_bitmask(k,1) = 0_bit_kind
! zero_bitmask(k,2) = 0_bit_kind
! enddo
! ! hole is possible only in the orbital i_hole_osoci
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,1),N_int)
! call set_bit_to_integer(i_hole_osoci,zero_bitmask(1,2),N_int)
! call set_bitmask_hole_as_input(zero_bitmask)
! call set_bitmask_particl_as_input(reunion_of_bitmask)
! call set_psi_det_to_generators
! call all_2h2p
! call diagonalize_CI_SC2
double precision :: hkl
call provide_matrix_dressing(dressing_matrix,n_det_generators,psi_det_generators)
hkl = dressing_matrix(1,1)
do k = 1, N_det_generators
dressing_matrix(k,k) = dressing_matrix(k,k) - hkl
enddo
print*,'Dressed matrix'
do k = 1, N_det_generators
write(*,'(100(F12.5,X))')dressing_matrix(k,:)
enddo
! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
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)
@ -132,24 +197,6 @@ subroutine FOBOCI_lmct_mlct_old_thr
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

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

@ -0,0 +1,18 @@
program osoci_program
implicit none
do_it_perturbative = .True.
touch do_it_perturbative
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 set_lmct_mlct_to_psi_det
call diagonalize_CI
call save_wavefunction
end

158
plugins/FOBOCI/generators_restart_save.irp.f
View File

@ -1,126 +1,74 @@
use bitmasks
use bitmasks
BEGIN_PROVIDER [ integer, N_det_generators_restart ]
implicit none
BEGIN_DOC
! Number of determinants in the wave function
! Read the wave function
END_DOC
logical :: exists
character*64 :: label
integer :: i
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)
double precision :: norm
if(ifirst == 0)then
call ezfio_get_determinants_n_det(N_det_generators_restart)
ifirst = 1
!endif
else
print*,'PB in generators_restart restart !!!'
endif
call write_int(output_determinants,N_det_generators_restart,'Number of generators_restart')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,psi_det_size) ]
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators_restart, (N_int,2,N_det_generators_restart) ]
&BEGIN_PROVIDER [ integer(bit_kind), ref_generators_restart, (N_int,2) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators_restart, (N_det_generators_restart,N_states) ]
implicit none
BEGIN_DOC
! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
! is empty
! read wf
!
END_DOC
integer :: i
logical :: exists
character*64 :: label
integer :: i, k
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
double precision, allocatable :: psi_coef_read(:,:)
if(ifirst == 0)then
call read_dets(psi_det_generators_restart,N_int,N_det_generators_restart)
do k = 1, N_int
ref_generators_restart(k,1) = psi_det_generators_restart(k,1,1)
ref_generators_restart(k,2) = psi_det_generators_restart(k,2,1)
enddo
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
ifirst = 1
!endif
deallocate(psi_coef_read)
else
print*,'PB in generators_restart restart !!!'
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
BEGIN_PROVIDER [ integer, size_select_max]
implicit none
BEGIN_DOC
! Size of the select_max array
END_DOC
size_select_max = 10000
END_PROVIDER
BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
implicit none
BEGIN_DOC
! Memo to skip useless selectors
END_DOC
select_max = huge(1.d0)
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_generators ]
&BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,10000) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators, (10000,N_states) ]
END_PROVIDER

83
plugins/FOBOCI/hcc_1h1p.irp.f Normal file
View File

@ -0,0 +1,83 @@
program test_sc2
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
double precision, allocatable :: energies(:),diag_H_elements(:)
double precision, allocatable :: H_matrix(:,:)
allocate(energies(N_states),diag_H_elements(N_det))
call diagonalize_CI
call test_hcc
call test_mulliken
! call SC2_1h1p(psi_det,psi_coef,energies, &
! diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
allocate(H_matrix(N_det,N_det))
call SC2_1h1p_full(psi_det,psi_coef,energies, &
H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
deallocate(H_matrix)
integer :: i,j
double precision :: accu,coef_hf
! coef_hf = 1.d0/psi_coef(1,1)
! do i = 1, N_det
! psi_coef(i,1) *= coef_hf
! enddo
touch psi_coef
call pouet
end
subroutine pouet
implicit none
double precision :: accu,coef_hf
! provide one_body_dm_mo_alpha one_body_dm_mo_beta
! call density_matrix_1h1p(psi_det,psi_coef,one_body_dm_mo_alpha,one_body_dm_mo_beta,accu,size(psi_coef,1),N_det,N_states_diag,N_int)
! touch one_body_dm_mo_alpha one_body_dm_mo_beta
call test_hcc
call test_mulliken
! call save_wavefunction
end
subroutine test_hcc
implicit none
double precision :: accu
integer :: i,j
print*,'Z AU GAUSS MHZ cm^-1'
do i = 1, nucl_num
write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
enddo
end
subroutine test_mulliken
double precision :: accu
integer :: i
integer :: j
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),mulliken_spin_densities(i)
accu += mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
!print*,'AO SPIN POPULATIONS'
accu = 0.d0
!do i = 1, ao_num
! accu += spin_gross_orbital_product(i)
! write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
!enddo
!print*,'sum = ',accu
!accu = 0.d0
!print*,'Angular momentum analysis'
!do i = 0, ao_l_max
! accu += spin_population_angular_momentum(i)
! print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
!print*,'sum = ',accu
!enddo
end

1
plugins/FOBOCI/modify_generators.irp.f
View File

@ -6,6 +6,7 @@ subroutine set_generators_to_psi_det
END_DOC
N_det_generators = N_det
integer :: i,k
print*,'N_det = ',N_det
do i=1,N_det_generators
do k=1,N_int
psi_det_generators(k,1,i) = psi_det(k,1,i)

316
plugins/FOBOCI/new_approach.irp.f
View File

@ -213,6 +213,8 @@ subroutine new_approach
deallocate(dressing_matrix_1h2p)
deallocate(dressing_matrix_extra_1h_or_1p)
enddo
double precision, allocatable :: H_matrix_total(:,:)
integer :: n_det_total
n_det_total = N_det_generators_restart + n_good_det
@ -251,25 +253,79 @@ subroutine new_approach
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
! Adding the correlation energy
logical :: orb_taken_good_det(mo_tot_num)
double precision :: phase
integer :: n_h,n_p,number_of_holes,number_of_particles
integer :: exc(0:2,2,2)
integer :: degree
integer :: h1,h2,p1,p2,s1,s2
logical, allocatable :: one_hole_or_one_p(:)
integer, allocatable :: holes_or_particle(:)
allocate(one_hole_or_one_p(n_good_det), holes_or_particle(n_good_det))
orb_taken_good_det = .False.
do i = 1, n_good_det
n_h = number_of_holes(psi_good_det(1,1,i))
n_p = number_of_particles(psi_good_det(1,1,i))
call get_excitation(ref_bitmask,psi_good_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
if(n_h == 0 .and. n_p == 1)then
orb_taken_good_det(h1) = .True.
one_hole_or_one_p(i) = .True.
holes_or_particle(i) = h1
endif
if(n_h == 1 .and. n_p == 0)then
orb_taken_good_det(p1) = .True.
one_hole_or_one_p(i) = .False.
holes_or_particle(i) = p1
endif
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
do i = 1, N_det_generators_restart
! Add the 2h2p, 2h1p and 1h2p correlation energy
H_matrix_total(i,i) += total_corr_e_2h2p + total_corr_e_2h1p + total_corr_e_1h2p + total_corr_e_1h1p_spin_flip
! Substract the 2h1p part that have already been taken into account
do j = 1, n_inact_orb
iorb = list_inact(j)
if(.not.orb_taken_good_det(iorb))cycle
H_matrix_total(i,i) -= corr_energy_2h1p_per_orb_ab(iorb) - corr_energy_2h1p_per_orb_bb(iorb) - corr_energy_1h1p_spin_flip_per_orb(iorb)
enddo
! Substract the 1h2p part that have already been taken into account
do j = 1, n_virt_orb
iorb = list_virt(j)
if(.not.orb_taken_good_det(iorb))cycle
H_matrix_total(i,i) -= corr_energy_1h2p_per_orb_ab(iorb) - corr_energy_1h2p_per_orb_aa(iorb)
enddo
enddo
do i = 1, N_good_det
! Repeat the 2h2p correlation energy
H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_2h2p
! Substract the part that can not be repeated
! If it is a 1h
if(one_hole_or_one_p(i))then
! 2h2p
H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) &
-corr_energy_2h2p_per_orb_bb(holes_or_particle(i))
! You can repeat a certain part of the 1h2p correlation energy
! that is everything except the part that involves the hole of the 1h
H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += total_corr_e_1h2p
H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_1h2p_per_orb_ab(holes_or_particle(i)) &
-corr_energy_1h2p_per_orb_bb(holes_or_particle(i))
else
! 2h2p
H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h2p_per_orb_ab(holes_or_particle(i)) &
-corr_energy_2h2p_per_orb_aa(holes_or_particle(i))
! You can repeat a certain part of the 2h1p correlation energy
! that is everything except the part that involves the hole of the 1p
! 2h1p
H_matrix_total(N_det_generators_restart+i,N_det_generators_restart+i) += -corr_energy_2h1p_per_orb_ab(holes_or_particle(i)) &
-corr_energy_2h1p_per_orb_aa(holes_or_particle(i))
endif
enddo
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
@ -284,22 +340,222 @@ subroutine new_approach
psi_det_final(j,2,n_det_generators_restart+i) = psi_good_det(j,2,i)
enddo
enddo
norm = 0.d0
double precision :: href
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
integer(bit_kind), allocatable :: psi_det_final(:,:,:)
double precision, allocatable :: psi_coef_final(:,:)
double precision :: norm
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*,''
print*,''
print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1)
print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion
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)
print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1)
enddo
print*,'norm = ',norm
integer(bit_kind), allocatable :: psi_det_remaining_1h_or_1p(:,:,:)
integer(bit_kind), allocatable :: key_tmp(:,:)
integer :: n_det_remaining_1h_or_1p
integer :: ispin,i_ok
allocate(key_tmp(N_int,2),psi_det_remaining_1h_or_1p(N_int,2,n_inact_orb*n_act_orb+n_virt_orb*n_act_orb))
logical :: is_already_present
logical, allocatable :: one_hole_or_one_p_bis(:)
integer, allocatable :: holes_or_particle_bis(:)
double precision,allocatable :: H_array(:)
allocate(one_hole_or_one_p_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb), holes_or_particle_bis(n_inact_orb*n_act_orb+n_virt_orb*n_act_orb))
allocate(H_array(n_det_total))
! Dressing with the remaining 1h determinants
print*,''
print*,''
print*,'Dressing with the remaining 1h determinants'
n_det_remaining_1h_or_1p = 0
do i = 1, n_inact_orb
iorb = list_inact(i)
if(orb_taken_good_det(iorb))cycle
do j = 1, n_act_orb
jorb = list_act(j)
ispin = 2
key_tmp = ref_bitmask
call do_mono_excitation(key_tmp,iorb,jorb,ispin,i_ok)
if(i_ok .ne.1)cycle
is_already_present = .False.
H_array = 0.d0
call i_h_j(key_tmp,key_tmp,N_int,hij)
href = ref_bitmask_energy - hij
href = 1.d0/href
do k = 1, n_det_total
call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int)
if(degree == 0)then
is_already_present = .True.
exit
endif
enddo
if(is_already_present)cycle
n_det_remaining_1h_or_1p +=1
one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .True.
holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb
do k = 1, N_int
psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1)
psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2)
enddo
! do k = 1, n_det_total
! call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij)
! H_array(k) = hij
! enddo
! do k = 1, n_det_total
! do l = 1, n_det_total
! H_matrix_total(k,l) += H_array(k) * H_array(l) * href
! enddo
! enddo
enddo
enddo
! Dressing with the remaining 1p determinants
print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p
print*,'Dressing with the remaining 1p determinants'
do i = 1, n_virt_orb
iorb = list_virt(i)
if(orb_taken_good_det(iorb))cycle
do j = 1, n_act_orb
jorb = list_act(j)
ispin = 1
key_tmp = ref_bitmask
call do_mono_excitation(key_tmp,jorb,iorb,ispin,i_ok)
if(i_ok .ne.1)cycle
is_already_present = .False.
H_array = 0.d0
call i_h_j(key_tmp,key_tmp,N_int,hij)
href = ref_bitmask_energy - hij
href = 1.d0/href
do k = 1, n_det_total
call get_excitation_degree(psi_det_final(1,1,k),key_tmp,degree,N_int)
if(degree == 0)then
is_already_present = .True.
exit
endif
enddo
if(is_already_present)cycle
n_det_remaining_1h_or_1p +=1
one_hole_or_one_p_bis(n_det_remaining_1h_or_1p) = .False.
holes_or_particle_bis(n_det_remaining_1h_or_1p) = iorb
do k = 1, N_int
psi_det_remaining_1h_or_1p(k,1,n_det_remaining_1h_or_1p) = key_tmp(k,1)
psi_det_remaining_1h_or_1p(k,2,n_det_remaining_1h_or_1p) = key_tmp(k,2)
enddo
! do k = 1, n_det_total
! call i_h_j(psi_det_final(1,1,k),key_tmp,N_int,hij)
! H_array(k) = hij
! enddo
! do k = 1, n_det_total
! do l = 1, n_det_total
! H_matrix_total(k,l) += H_array(k) * H_array(l) * href
! enddo
! enddo
enddo
enddo
print*,'n_det_remaining_1h_or_1p = ',n_det_remaining_1h_or_1p
deallocate(key_tmp,H_array)
double precision, allocatable :: eigvalues_bis(:),eigvectors_bis(:,:),H_matrix_total_bis(:,:)
integer :: n_det_final
n_det_final = n_det_total + n_det_remaining_1h_or_1p
allocate(eigvalues_bis(n_det_final),eigvectors_bis(n_det_final,n_det_final),H_matrix_total_bis(n_det_final,n_det_final))
print*,'passed the allocate, building the big matrix'
do i = 1, n_det_total
do j = 1, n_det_total
H_matrix_total_bis(i,j) = H_matrix_total(i,j)
enddo
enddo
do i = 1, n_det_remaining_1h_or_1p
do j = 1, n_det_remaining_1h_or_1p
call i_h_j(psi_det_remaining_1h_or_1p(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij)
H_matrix_total_bis(n_det_total+i,n_det_total+j) = hij
enddo
enddo
do i = 1, n_det_total
do j = 1, n_det_remaining_1h_or_1p
call i_h_j(psi_det_final(1,1,i),psi_det_remaining_1h_or_1p(1,1,j),N_int,hij)
H_matrix_total_bis(i,n_det_total+j) = hij
H_matrix_total_bis(n_det_total+j,i) = hij
enddo
enddo
print*,'passed the matrix'
do i = 1, n_det_remaining_1h_or_1p
if(one_hole_or_one_p_bis(i))then
H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) &
-corr_energy_2h2p_per_orb_bb(holes_or_particle_bis(i))
H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_1h2p_per_orb_ab(holes_or_particle_bis(i)) &
-corr_energy_1h2p_per_orb_bb(holes_or_particle_bis(i))
else
H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_2h2p -corr_energy_2h2p_per_orb_ab(holes_or_particle_bis(i)) &
-corr_energy_2h2p_per_orb_aa(holes_or_particle_bis(i))
H_matrix_total_bis(n_det_total+i,n_det_total+i) += total_corr_e_1h2p -corr_energy_2h1p_per_orb_ab(holes_or_particle_bis(i)) &
-corr_energy_2h1p_per_orb_aa(holes_or_particle_bis(i))
endif
enddo
do i = 2, n_det_final
do j = i+1, n_det_final
H_matrix_total_bis(i,j) = 0.d0
H_matrix_total_bis(j,i) = 0.d0
enddo
enddo
do i = 1, n_det_final
write(*,'(500(F10.5,X))')H_matrix_total_bis(i,:)
enddo
call lapack_diag(eigvalues_bis,eigvectors_bis,H_matrix_total_bis,n_det_final,n_det_final)
print*,'e_dressed = ',eigvalues_bis(1) + nuclear_repulsion
do i = 1, n_det_final
print*,'coef = ',eigvectors_bis(i,1),H_matrix_total_bis(i,i) - H_matrix_total_bis(1,1)
enddo
do j = 1, N_states
do i = 1, n_det_total
psi_coef_final(i,j) = eigvectors_bis(i,j)
norm += psi_coef_final(i,j)**2
enddo
norm = 1.d0/dsqrt(norm)
do i = 1, n_det_total
psi_coef_final(i,j) = psi_coef_final(i,j) * norm
enddo
enddo
deallocate(eigvalues_bis,eigvectors_bis,H_matrix_total_bis)
!print*,'H matrix to diagonalize'
!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
!call lapack_diag(eigvalues,eigvectors,H_matrix_total,n_det_total,n_det_total)
!print*,'H_matrix_total(1,1) = ',H_matrix_total(1,1)
!print*,'e_dressed = ',eigvalues(1) + nuclear_repulsion
!do i = 1, n_det_total
! print*,'coef = ',eigvectors(i,1),H_matrix_total(i,i) - H_matrix_total(1,1)
!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
!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
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

132
plugins/FOBOCI/new_new_approach.irp.f Normal file
View File

@ -0,0 +1,132 @@
program test_new_new
implicit none
read_wf = .True.
touch read_wf
call test
end
subroutine test
implicit none
integer :: i,j,k,l
call diagonalize_CI
call set_generators_to_psi_det
print*,'Initial coefficients'
do i = 1, N_det
print*,''
call debug_det(psi_det(1,1,i),N_int)
print*,'psi_coef = ',psi_coef(i,1)
print*,''
enddo
double precision, allocatable :: dressing_matrix(:,:)
double precision :: hij
double precision :: phase
integer :: n_h,n_p,number_of_holes,number_of_particles
integer :: exc(0:2,2,2)
integer :: degree
integer :: h1,h2,p1,p2,s1,s2
allocate(dressing_matrix(N_det_generators,N_det_generators))
do i = 1, N_det_generators
do j = 1, N_det_generators
call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij)
dressing_matrix(i,j) = hij
enddo
enddo
href = dressing_matrix(1,1)
print*,'Diagonal part of the dressing'
do i = 1, N_det_generators
print*,'delta e = ',dressing_matrix(i,i) - href
enddo
call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
double precision :: href
print*,''
! One considers that the following excitation classes are not repeatable on the 1h and 1p determinants :
! + 1h1p spin flip
! + 2h1p
! + 1h2p
! But the 2h2p are correctly taken into account
!dressing_matrix(1,1) += total_corr_e_1h2p + total_corr_e_2h1p + total_corr_e_1h1p_spin_flip
!do i = 1, N_det_generators
! dressing_matrix(i,i) += total_corr_e_2h2p
! 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
!
! call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int)
! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
! print*,''
! print*,' 1h det '
! print*,''
! call debug_det(psi_det_generators(1,1,i),N_int)
! print*,'h1,p1 = ',h1,p1
! print*,'total_corr_e_2h2p ',total_corr_e_2h2p
! print*,'corr_energy_2h2p_per_orb_ab(h1)',corr_energy_2h2p_per_orb_ab(h1)
! print*,'corr_energy_2h2p_per_orb_bb(h1)',corr_energy_2h2p_per_orb_bb(h1)
! dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(h1) - corr_energy_2h2p_per_orb_bb(h1)
! dressing_matrix(1,1) += -corr_energy_2h1p_per_orb_aa(h1) - corr_energy_2h1p_per_orb_ab(h1) -corr_energy_2h1p_per_orb_bb(h1) &
! -corr_energy_1h1p_spin_flip_per_orb(h1)
! endif
! if(n_h == 0 .and. n_p ==1)then
! call get_excitation(ref_bitmask,psi_det_generators(1,1,i),exc,degree,phase,N_int)
! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
! print*,''
! print*,' 1p det '
! print*,''
! call debug_det(psi_det_generators(1,1,i),N_int)
! print*,'h1,p1 = ',h1,p1
! print*,'total_corr_e_2h2p ',total_corr_e_2h2p
! print*,'corr_energy_2h2p_per_orb_ab(p1)',corr_energy_2h2p_per_orb_ab(p1)
! print*,'corr_energy_2h2p_per_orb_aa(p1)',corr_energy_2h2p_per_orb_aa(p1)
! dressing_matrix(i,i) += -corr_energy_2h2p_per_orb_ab(p1) - corr_energy_2h2p_per_orb_aa(p1)
! dressing_matrix(1,1) += -corr_energy_1h2p_per_orb_aa(p1) - corr_energy_1h2p_per_orb_ab(p1) -corr_energy_1h2p_per_orb_bb(p1)
! endif
!enddo
!href = dressing_matrix(1,1)
!print*,'Diagonal part of the dressing'
!do i = 1, N_det_generators
! print*,'delta e = ',dressing_matrix(i,i) - href
!enddo
call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
print*,'After dressing matrix'
print*,''
print*,''
do i = 1, N_det
print*,'psi_coef = ',psi_coef(i,1)
enddo
!print*,''
!print*,''
!print*,'Canceling the dressing part of the interaction between 1h and 1p'
!do i = 2, N_det_generators
! do j = i+1, N_det_generators
! call i_h_j(psi_det_generators(1,1,i),psi_det_generators(1,1,j),N_int,hij)
! dressing_matrix(i,j) = hij
! dressing_matrix(j,i) = hij
! enddo
!enddo
!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
!print*,''
!print*,''
!do i = 1, N_det
! print*,'psi_coef = ',psi_coef(i,1)
!enddo
!print*,''
!print*,''
!print*,'Canceling the interaction between 1h and 1p'
!print*,''
!print*,''
!do i = 2, N_det_generators
! do j = i+1, N_det_generators
! dressing_matrix(i,j) = 0.d0
! dressing_matrix(j,i) = 0.d0
! enddo
!enddo
!call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
!do i = 1, N_det
! print*,'psi_coef = ',psi_coef(i,1)
!enddo
call save_natural_mos
deallocate(dressing_matrix)
end

314
plugins/FOBOCI/routines_dressing.irp.f
View File

@ -55,15 +55,11 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det
i_pert = 0
endif
do j = 1, ndet_generators_input
if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
if(dabs(H_array(j)*lambda_i).gt.0.1d0)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
@ -79,8 +75,52 @@ subroutine provide_matrix_dressing(dressing_matrix,ndet_generators_input,psi_det
enddo
enddo
enddo
href = dressing_matrix(1,1)
print*,'Diagonal part of the dressing'
do i = 1, ndet_generators_input
print*,'delta e = ',dressing_matrix(i,i) - href
enddo
!print*,'i_pert_count = ',i_pert_count
end
subroutine update_matrix_dressing_sc2(dressing_matrix,ndet_generators_input,psi_det_generators_input,H_jj_in)
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(in) :: H_jj_in(N_det)
double precision, intent(inout) :: dressing_matrix(ndet_generators_input,ndet_generators_input)
integer :: i,j,n_det_ref_tmp,degree
double precision :: href
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
dressing_matrix(j,j) += H_jj_in(i)
n_det_ref_tmp +=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
href = dressing_matrix(1,1)
print*,''
print*,'Update with the SC2 dressing'
print*,''
print*,'Diagonal part of the dressing'
do i = 1, ndet_generators_input
print*,'delta e = ',dressing_matrix(i,i) - href
enddo
end
subroutine provide_matrix_dressing_for_extra_1h_or_1p(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
@ -125,11 +165,14 @@ subroutine provide_matrix_dressing_for_extra_1h_or_1p(dressing_matrix,psi_det_re
i_pert = 0
endif
do j = 1, n_det_ref_input
if(dabs(H_array(j)*lambda_i).gt.0.3d0)then
if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
i_pert = 1
exit
endif
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! i_pert = 0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(i_pert==1)then
lambda_i = f
i_pert_count +=1
@ -178,16 +221,17 @@ subroutine provide_matrix_dressing_general(dressing_matrix,psi_det_ref_input,psi
accu += psi_coef_ref_input(j,1) * hka
enddo
lambda_i = psi_coef_outer_input(i,1)/accu
i_pert = 1
i_pert = 0
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
if(dabs(H_array(j)*lambda_i).gt.0.5d0)then
i_pert = 1
exit
endif
enddo
! i_pert = 0
if(i_pert==1)then
lambda_i = f
i_pert_count +=1
@ -275,19 +319,257 @@ subroutine give_n_1h1p_and_n_2h1p_in_psi_det(i_hole,n_det_extra_1h_or_1p,n_det_1
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
if(n_det_2h1p + n_det_1h1p + n_det_extra_1h_or_1p + n_det_generators .ne. N_det)then
print*,'PB !!!!'
print*,'You have forgotten something in your generators ... '
stop
endif
end
subroutine give_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p)
use bitmasks
implicit none
integer, intent(out) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
integer :: i
integer :: n_det_ref_restart_tmp,n_det_1h
integer :: number_of_holes,n_h, number_of_particles,n_p
logical :: is_the_hole_in_det
n_det_ref_1h_1p = 0
n_det_2h1p = 0
n_det_1h1p = 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_1h_1p +=1
else if (n_h ==1 .and. n_p==0)then
n_det_ref_1h_1p +=1
else if (n_h ==0 .and. n_p==1)then
n_det_ref_1h_1p +=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, 1p, 1h1p or 2h1p'
print*,'n_h,n_p = ',n_h,n_p
call debug_det(psi_det(1,1,i),N_int)
stop
endif
enddo
end
subroutine give_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p)
use bitmasks
implicit none
integer, intent(out) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
integer :: i
integer :: n_det_ref_restart_tmp,n_det_1h
integer :: number_of_holes,n_h, number_of_particles,n_p
logical :: is_the_hole_in_det
n_det_ref_1h_1p = 0
n_det_1h2p = 0
n_det_1h1p = 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_1h_1p +=1
else if (n_h ==1 .and. n_p==0)then
n_det_ref_1h_1p +=1
else if (n_h ==0 .and. n_p==1)then
n_det_ref_1h_1p +=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 1h, 1p, 1h1p or 1h2p'
print*,'n_h,n_p = ',n_h,n_p
call debug_det(psi_det(1,1,i),N_int)
stop
endif
enddo
end
subroutine give_wf_n_ref_1h_1p_and_n_2h1p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
psi_det_2h1p,psi_coef_2h1p,psi_det_1h1p,psi_coef_1h1p)
use bitmasks
implicit none
integer, intent(in) :: n_det_ref_1h_1p,n_det_2h1p,n_det_1h1p
integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p)
integer(bit_kind), intent(out) :: psi_det_2h1p(N_int,2,n_det_2h1p)
integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p)
double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)
double precision, intent(out) :: psi_coef_2h1p(n_det_2h1p,N_states)
double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states)
integer :: n_det_ref_1h_1p_tmp,n_det_2h1p_tmp,n_det_1h1p_tmp
integer :: i,j
integer :: n_det_ref_restart_tmp,n_det_1h
integer :: number_of_holes,n_h, number_of_particles,n_p
logical :: is_the_hole_in_det
integer, allocatable :: index_ref_1h_1p(:)
integer, allocatable :: index_2h1p(:)
integer, allocatable :: index_1h1p(:)
allocate(index_ref_1h_1p(n_det))
allocate(index_2h1p(n_det))
allocate(index_1h1p(n_det))
n_det_ref_1h_1p_tmp = 0
n_det_2h1p_tmp = 0
n_det_1h1p_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 == 0 .and. n_p == 0)then
n_det_ref_1h_1p_tmp +=1
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
else if (n_h ==1 .and. n_p==0)then
n_det_ref_1h_1p_tmp +=1
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
else if (n_h ==0 .and. n_p==1)then
n_det_ref_1h_1p_tmp +=1
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
else 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
else
print*,'PB !!!!'
print*,'You have something else than a 1h, 1p, 1h1p or 2h1p'
print*,'n_h,n_p = ',n_h,n_p
call debug_det(psi_det(1,1,i),N_int)
stop
endif
enddo
do i = 1, n_det_2h1p
do j = 1, N_int
psi_det_2h1p(j,1,i) = psi_det(j,1,index_2h1p(i))
psi_det_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
do i = 1, n_det_1h1p
do j = 1, N_int
psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
psi_det_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_ref_1h_1p
do j = 1, N_int
psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i))
psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i))
enddo
do j = 1, N_states
psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j)
enddo
enddo
end
subroutine give_wf_n_ref_1h_1p_and_n_1h2p_1h1p_in_psi_det(n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p,psi_det_ref_1h_1p,psi_coef_ref_1h_1p,&
psi_det_1h2p,psi_coef_1h2p,psi_det_1h1p,psi_coef_1h1p)
use bitmasks
implicit none
integer, intent(in) :: n_det_ref_1h_1p,n_det_1h2p,n_det_1h1p
integer(bit_kind), intent(out) :: psi_det_ref_1h_1p(N_int,2,n_det_ref_1h_1p)
integer(bit_kind), intent(out) :: psi_det_1h2p(N_int,2,n_det_1h2p)
integer(bit_kind), intent(out) :: psi_det_1h1p(N_int,2,n_det_1h1p)
double precision, intent(out) :: psi_coef_ref_1h_1p(n_det_ref_1h_1p,N_states)
double precision, intent(out) :: psi_coef_1h2p(n_det_1h2p,N_states)
double precision, intent(out) :: psi_coef_1h1p(n_det_1h1p,N_states)
integer :: n_det_ref_1h_1p_tmp,n_det_1h2p_tmp,n_det_1h1p_tmp
integer :: i,j
integer :: n_det_ref_restart_tmp,n_det_1h
integer :: number_of_holes,n_h, number_of_particles,n_p
logical :: is_the_hole_in_det
integer, allocatable :: index_ref_1h_1p(:)
integer, allocatable :: index_1h2p(:)
integer, allocatable :: index_1h1p(:)
allocate(index_ref_1h_1p(n_det))
allocate(index_1h2p(n_det))
allocate(index_1h1p(n_det))
n_det_ref_1h_1p_tmp = 0
n_det_1h2p_tmp = 0
n_det_1h1p_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 == 0 .and. n_p == 0)then
n_det_ref_1h_1p_tmp +=1
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
else if (n_h ==1 .and. n_p==0)then
n_det_ref_1h_1p_tmp +=1
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
else if (n_h ==0 .and. n_p==1)then
n_det_ref_1h_1p_tmp +=1
index_ref_1h_1p(n_det_ref_1h_1p_tmp) = i
else 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
else
print*,'PB !!!!'
print*,'You have something else than a 1h, 1p, 1h1p or 1h2p'
print*,'n_h,n_p = ',n_h,n_p
call debug_det(psi_det(1,1,i),N_int)
stop
endif
enddo
do i = 1, n_det_1h2p
do j = 1, N_int
psi_det_1h2p(j,1,i) = psi_det(j,1,index_1h2p(i))
psi_det_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
do i = 1, n_det_1h1p
do j = 1, N_int
psi_det_1h1p(j,1,i) = psi_det(j,1,index_1h1p(i))
psi_det_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_ref_1h_1p
do j = 1, N_int
psi_det_ref_1h_1p(j,1,i) = psi_det(j,1,index_ref_1h_1p(i))
psi_det_ref_1h_1p(j,2,i) = psi_det(j,2,index_ref_1h_1p(i))
enddo
do j = 1, N_states
psi_coef_ref_1h_1p(i,j) = psi_coef(index_ref_1h_1p(i),j)
enddo
enddo
end
subroutine give_n_1h1p_and_n_1h2p_in_psi_det(i_particl,n_det_extra_1h_or_1p,n_det_1h1p,n_det_1h2p)
use bitmasks
implicit none
@ -353,7 +635,7 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_d
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,n_det_extra_1h_or_1p_tmp
integer :: n_det_extra_1h_tmp
integer :: n_det_1h_tmp
integer, allocatable :: index_generator(:)
integer, allocatable :: index_1h1p(:)
integer, allocatable :: index_2h1p(:)
@ -370,7 +652,7 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_d
n_det_1h1p_tmp = 0
n_det_2h1p_tmp = 0
n_det_extra_1h_or_1p_tmp = 0
n_det_extra_1h_tmp = 0
n_det_1h_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))
@ -385,7 +667,7 @@ subroutine split_wf_generators_and_1h1p_and_2h1p(i_hole,n_det_extra_1h_or_1p,n_d
index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i
else 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_det_extra_1h_tmp +=1
n_det_1h_tmp +=1
else
n_det_extra_1h_or_1p_tmp +=1
index_extra_1h_or_1p(n_det_extra_1h_or_1p_tmp) = i

51
plugins/FOBOCI/routines_foboci.irp.f
View File

@ -332,20 +332,20 @@ subroutine save_osoci_natural_mos
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-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
@ -389,14 +389,14 @@ subroutine save_osoci_natural_mos
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
print*,'INACTIVE '
print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb))
print*,'DM ',iorb,jorb,(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))
print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
endif
enddo
enddo
@ -410,8 +410,9 @@ subroutine save_osoci_natural_mos
enddo
label = "Natural"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
soft_touch mo_coef
!soft_touch mo_coef
deallocate(tmp,occ)
@ -520,14 +521,14 @@ subroutine set_osoci_natural_mos
jorb = list_inact(j)
if(dabs(tmp(iorb,jorb)).gt.threshold_singles)then
print*,'INACTIVE '
print*,'DM ',iorb,jorb,dabs(tmp(iorb,jorb))
print*,'DM ',iorb,jorb,(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))
print*,'DM ',iorb,jorb,(tmp(iorb,jorb))
endif
enddo
enddo
@ -602,15 +603,7 @@ 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
call print_mulliken_sd
call print_hcc
end

28
plugins/Generators_restart/generators.irp.f
View File

@ -1,5 +1,5 @@
use bitmasks
BEGIN_PROVIDER [ integer, N_det_generators ]
implicit none
BEGIN_DOC
@ -8,17 +8,18 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
integer :: i
integer, save :: ifirst = 0
double precision :: norm
read_wf = .True.
if(ifirst == 0)then
N_det_generators = N_det
call ezfio_get_determinants_n_det(N_det_generators)
ifirst = 1
else
print*,'PB in generators restart !!!'
endif
call write_int(output_determinants,N_det_generators,'Number of generators')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ]
implicit none
BEGIN_DOC
! read wf
@ -26,17 +27,20 @@ END_PROVIDER
END_DOC
integer :: i, k
integer, save :: ifirst = 0
double precision, allocatable :: psi_coef_read(:,:)
if(ifirst == 0)then
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
call read_dets(psi_det_generators,N_int,N_det_generators)
allocate (psi_coef_read(N_det_generators,N_states))
call ezfio_get_determinants_psi_coef(psi_coef_read)
do k = 1, N_states
psi_coef_generators(i,k) = psi_coef(i,k)
do i = 1, N_det_generators
psi_coef_generators(i,k) = psi_coef_read(i,k)
enddo
enddo
enddo
ifirst = 1
deallocate(psi_coef_read)
else
print*,'PB in generators restart !!!'
endif
END_PROVIDER

4
plugins/Hartree_Fock/damping_SCF.irp.f
View File

@ -119,7 +119,9 @@ subroutine damping_SCF
write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '===='
write(output_hartree_fock,*)
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
if(.not.no_oa_or_av_opt)then
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
endif
call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')
call ezfio_set_hartree_fock_energy(E_min)

2
plugins/Perturbation/pt2_equations.irp.f
View File

@ -126,6 +126,8 @@ subroutine pt2_moller_plesset ($arguments)
delta_e = (Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)) + &
(Fock_matrix_diag_mo(h2) - Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
! print*,'h1,p1',h1,p1
! print*,'h2,p2',h2,p2
else if (degree == 1) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)

13
plugins/Properties/hyperfine_constants.irp.f
View File

@ -133,3 +133,16 @@ END_PROVIDER
enddo
END_PROVIDER
subroutine print_hcc
implicit none
double precision :: accu
integer :: i,j
print*,'Z AU GAUSS MHZ cm^-1'
do i = 1, nucl_num
write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
enddo
end

31
plugins/Properties/mulliken.irp.f
View File

@ -105,3 +105,34 @@ END_PROVIDER
enddo
END_PROVIDER
subroutine print_mulliken_sd
implicit none
double precision :: accu
integer :: i
integer :: j
print*,'Mulliken spin densities'
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),mulliken_spin_densities(i)
accu += mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
print*,'AO SPIN POPULATIONS'
accu = 0.d0
do i = 1, ao_num
accu += spin_gross_orbital_product(i)
write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
enddo
print*,'sum = ',accu
accu = 0.d0
print*,'Angular momentum analysis'
do i = 0, ao_l_max
accu += spin_population_angular_momentum(i)
print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
print*,'sum = ',accu
enddo
end

15
plugins/Properties/print_hcc.irp.f
View File

@ -1,17 +1,6 @@
program print_hcc
program print_hcc_main
implicit none
read_wf = .True.
touch read_wf
call test
call print_hcc
end
subroutine test
implicit none
double precision :: accu
integer :: i,j
print*,'Z AU GAUSS MHZ cm^-1'
do i = 1, nucl_num
write(*,'(I2,X,F3.1,X,4(F16.6,X))')i,nucl_charge(i),spin_density_at_nucleous(i),iso_hcc_gauss(i),iso_hcc_mhz(i),iso_hcc_cm_1(i)
enddo
end

31
plugins/Properties/print_mulliken.irp.f
View File

@ -2,34 +2,5 @@ program print_mulliken
implicit none
read_wf = .True.
touch read_wf
print*,'Mulliken spin densities'
call test
call print_mulliken_sd
end
subroutine test
double precision :: accu
integer :: i
integer :: j
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),mulliken_spin_densities(i)
accu += mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
print*,'AO SPIN POPULATIONS'
accu = 0.d0
do i = 1, ao_num
accu += spin_gross_orbital_product(i)
write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
enddo
print*,'sum = ',accu
accu = 0.d0
print*,'Angular momentum analysis'
do i = 0, ao_l_max
accu += spin_population_angular_momentum(i)
print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
print*,'sum = ',accu
enddo
end

33
scripts/generate_h_apply.py
View File

@ -24,13 +24,18 @@ skip
init_main
filter_integrals
filter2p
filter2h2p
filter2h2p_double
filter2h2p_single
filter1h
filter1p
only_2p_single
only_2p_double
filter_only_1h1p_single
filter_only_1h1p_double
filter_only_1h2p_single
filter_only_1h2p_double
filter_only_2h2p_single
filter_only_2h2p_double
filterhole
filterparticle
do_double_excitations
@ -194,6 +199,27 @@ class H_apply(object):
if (is_a_1h1p(key).eqv..False.) cycle
"""
def filter_only_2h2p(self):
self["filter_only_2h2p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_two_holes_two_particles(hole).eqv..False.) cycle
"""
self["filter_only_1h1p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_two_holes_two_particles(key).eqv..False.) cycle
"""
def filter_only_1h2p(self):
self["filter_only_1h2p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h2p(hole).eqv..False.) cycle
"""
self["filter_only_1h2p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h2p(key).eqv..False.) cycle
"""
def unset_skip(self):
self["skip"] = """
@ -201,9 +227,12 @@ class H_apply(object):
def set_filter_2h_2p(self):
self["filter2h2p"] = """
self["filter2h2p_double"] = """
if (is_a_two_holes_two_particles(key)) cycle
"""
self["filter2h2p_single"] = """
if (is_a_two_holes_two_particles(hole)) cycle
"""
def set_perturbation(self,pert):

17
src/Bitmask/bitmask_cas_routines.irp.f
View File

@ -212,6 +212,12 @@ logical function is_a_two_holes_two_particles(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: i,i_diff
integer :: number_of_holes, number_of_particles
is_a_two_holes_two_particles = .False.
if(number_of_holes(key_in) == 2 .and. number_of_particles(key_in) == 2)then
is_a_two_holes_two_particles = .True.
return
endif
i_diff = 0
if(N_int == 1)then
i_diff = i_diff &
@ -456,6 +462,17 @@ logical function is_a_1h1p(key_in)
end
logical function is_a_1h2p(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: number_of_particles, number_of_holes
is_a_1h2p = .False.
if(number_of_holes(key_in).eq.1 .and. number_of_particles(key_in).eq.2)then
is_a_1h2p = .True.
endif
end
logical function is_a_1h(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)

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

@ -98,9 +98,40 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ]
END_PROVIDER
BEGIN_PROVIDER [ integer, N_generators_bitmask_restart ]
implicit none
BEGIN_DOC
! Number of bitmasks for generators
END_DOC
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_bitmasks_N_mask_gen(exists)
if (exists) then
call ezfio_get_bitmasks_N_mask_gen(N_generators_bitmask_restart)
integer :: N_int_check
integer :: bit_kind_check
call ezfio_get_bitmasks_bit_kind(bit_kind_check)
if (bit_kind_check /= bit_kind) then
print *, bit_kind_check, bit_kind
print *, 'Error: bit_kind is not correct in EZFIO file'
endif
call ezfio_get_bitmasks_N_int(N_int_check)
if (N_int_check /= N_int) then
print *, N_int_check, N_int
print *, 'Error: N_int is not correct in EZFIO file'
endif
else
N_generators_bitmask_restart = 1
endif
ASSERT (N_generators_bitmask_restart > 0)
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask) ]
BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask_restart, (N_int,2,6,N_generators_bitmask_restart) ]
implicit none
BEGIN_DOC
! Bitmasks for generator determinants.
@ -258,7 +289,7 @@ BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ]
call ezfio_get_bitmasks_cas(cas_bitmask)
print*,'---------------------'
else
if(N_generators_bitmask == 1)then
if(N_generators_bitmask_restart == 1)then
do i=1,N_cas_bitmask
cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
enddo
@ -302,7 +333,7 @@ END_PROVIDER
n_inact_orb = 0
n_virt_orb = 0
if(N_generators_bitmask == 1)then
if(N_generators_bitmask_restart == 1)then
do j = 1, N_int
inact_bitmask(j,1) = xor(generators_bitmask_restart(j,1,1,1),cas_bitmask(j,1,1))
inact_bitmask(j,2) = xor(generators_bitmask_restart(j,2,1,1),cas_bitmask(j,2,1))
@ -315,15 +346,15 @@ END_PROVIDER
i_hole = 1
i_gen = 1
do i = 1, N_int
inact_bitmask(i,1) = generators_bitmask(i,1,i_hole,i_gen)
inact_bitmask(i,2) = generators_bitmask(i,2,i_hole,i_gen)
inact_bitmask(i,1) = generators_bitmask_restart(i,1,i_hole,i_gen)
inact_bitmask(i,2) = generators_bitmask_restart(i,2,i_hole,i_gen)
n_inact_orb += popcnt(inact_bitmask(i,1))
enddo
i_part = 2
i_gen = 3
do i = 1, N_int
virt_bitmask(i,1) = generators_bitmask(i,1,i_part,i_gen)
virt_bitmask(i,2) = generators_bitmask(i,2,i_part,i_gen)
virt_bitmask(i,1) = generators_bitmask_restart(i,1,i_part,i_gen)
virt_bitmask(i,2) = generators_bitmask_restart(i,2,i_part,i_gen)
n_virt_orb += popcnt(virt_bitmask(i,1))
enddo
endif

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

@ -166,6 +166,7 @@ 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_1h2p
logical :: is_a_1h
logical :: is_a_1p
logical :: is_a_2p
@ -301,8 +302,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
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)
$filter2h2p
$filter2h2p_double
$filter_only_1h1p_double
$filter_only_1h2p_double
$filter_only_2h2p_double
$only_2p_double
key_idx += 1
do k=1,N_int
@ -353,8 +356,10 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
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)
$filter2h2p
$filter2h2p_double
$filter_only_1h1p_double
$filter_only_1h2p_double
$filter_only_2h2p_double
$only_2p_double
key_idx += 1
do k=1,N_int
@ -427,6 +432,7 @@ 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_1h2p
logical :: is_a_1h
logical :: is_a_1p
logical :: is_a_2p
@ -505,8 +511,10 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
$filter1h
$filter1p
$filter2p
$filter2h2p
$filter2h2p_single
$filter_only_1h1p_single
$filter_only_1h2p_single
$filter_only_2h2p_single
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
@ -566,7 +574,6 @@ subroutine $subroutine($params_main)
iproc = 0
allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) )
do i_generator=1,nmax
progress_bar(1) = i_generator
if (abort_here) then
@ -600,6 +607,16 @@ subroutine $subroutine($params_main)
not(psi_det_generators(k,ispin,i_generator)) )
enddo
enddo
! print*,'generator in '
! call debug_det(psi_det_generators(1,1,i_generator),N_int)
! print*,'hole 1'
! call debug_det(mask(1,1,d_hole1),N_int)
! print*,'hole 2'
! call debug_det(mask(1,1,d_hole2),N_int)
! print*,'part 1'
! call debug_det(mask(1,1,d_part1),N_int)
! print*,'part 2'
! call debug_det(mask(1,1,d_part2),N_int)
if($do_double_excitations)then
call $subroutine_diexc(psi_det_generators(1,1,i_generator), &
psi_det_generators(1,1,1), &

6
src/Determinants/SC2.irp.f
View File

@ -1,4 +1,4 @@
subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
subroutine CISD_SC2(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
use bitmasks
implicit none
BEGIN_DOC
@ -21,6 +21,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st)
double precision, intent(out) :: diag_H_elements(dim_in)
double precision, intent(in) :: convergence
ASSERT (N_st > 0)
ASSERT (sze > 0)
@ -200,6 +201,9 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
converged = dabs(e_corr_double - e_corr_double_before) < convergence
converged = converged .or. abort_here
if (converged) then
do i = 1, dim_in
diag_H_elements(i) = H_jj_dressed(i) - H_jj_ref(i)
enddo
exit
endif
e_corr_double_before = e_corr_double

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

@ -58,7 +58,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
else
psi_det_size = 1
endif
psi_det_size = max(psi_det_size,10000)
psi_det_size = max(psi_det_size,100000)
call write_int(output_determinants,psi_det_size,'Dimension of the psi arrays')
END_PROVIDER

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

@ -23,8 +23,10 @@ END_PROVIDER
threshold_convergence_SC2 = 1.d-10
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
@ -39,7 +41,8 @@ END_PROVIDER
enddo
call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
! size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
END_PROVIDER
subroutine diagonalize_CI_SC2
@ -54,5 +57,6 @@ subroutine diagonalize_CI_SC2
psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
enddo
enddo
SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors
SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2
! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors
end

1
src/Determinants/save_natorb.irp.f
View File

@ -2,5 +2,6 @@ program save_natorb
read_wf = .True.
touch read_wf
call save_natural_mos
call save_ref_determinant
end

1
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -230,7 +230,6 @@ subroutine clear_ao_map
end
!! MO Map
!! ======

13
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -72,7 +72,7 @@ subroutine add_integrals_to_map(mask_ijkl)
integer :: i2,i3,i4
double precision,parameter :: thr_coef = 1.d-10
PROVIDE ao_bielec_integrals_in_map
PROVIDE ao_bielec_integrals_in_map mo_coef
!Get list of MOs for i,j,k and l
!-------------------------------
@ -341,7 +341,7 @@ end
double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
if (.not.do_direct_integrals) then
PROVIDE ao_bielec_integrals_in_map
PROVIDE ao_bielec_integrals_in_map mo_coef
endif
mo_bielec_integral_jj_from_ao = 0.d0
@ -513,4 +513,13 @@ subroutine clear_mo_map
call map_deinit(mo_integrals_map)
FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
end
subroutine provide_all_mo_integrals
implicit none
provide mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
provide mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
end

1
src/Integrals_Monoelec/mo_mono_ints.irp.f
View File

@ -5,6 +5,7 @@ BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_to
! array of the mono electronic hamiltonian on the MOs basis
! : sum of the kinetic and nuclear electronic potential
END_DOC
print*,'Providing the mono electronic integrals'
do j = 1, mo_tot_num
do i = 1, mo_tot_num
mo_mono_elec_integral(i,j) = mo_nucl_elec_integral(i,j) + mo_kinetic_integral(i,j) + mo_pseudo_integral(i,j)