mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-25 05:51:46 +02:00
no vvv integrals is ok
This commit is contained in:
parent
50d1f364e0
commit
82a29d5603
@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
|
|||||||
# 0 : Deactivate
|
# 0 : Deactivate
|
||||||
#
|
#
|
||||||
[OPTION]
|
[OPTION]
|
||||||
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
|
||||||
CACHE : 1 ; Enable cache_compile.py
|
CACHE : 1 ; Enable cache_compile.py
|
||||||
OPENMP : 1 ; Append OpenMP flags
|
OPENMP : 1 ; Append OpenMP flags
|
||||||
|
|
||||||
|
@ -89,11 +89,11 @@ subroutine standard_dress(delta_ij_generators_,size_buffer,Ndet_generators,i_gen
|
|||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
subroutine is_a_good_candidate(threshold,is_ok,verbose)
|
subroutine is_a_good_candidate(threshold,is_ok,verbose,exit_loop)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
double precision, intent(in) :: threshold
|
double precision, intent(in) :: threshold
|
||||||
logical, intent(out) :: is_ok
|
logical, intent(out) :: is_ok,exit_loop
|
||||||
logical, intent(in) :: verbose
|
logical, intent(in) :: verbose
|
||||||
|
|
||||||
integer :: l,k,m
|
integer :: l,k,m
|
||||||
@ -111,7 +111,7 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose)
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input)
|
!call H_apply_dressed_pert(dressed_H_matrix,N_det_generators,psi_det_generators_input)
|
||||||
call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose)
|
call dress_H_matrix_from_psi_det_input(psi_det_generators_input,N_det_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop)
|
||||||
if(do_it_perturbative)then
|
if(do_it_perturbative)then
|
||||||
if(is_ok)then
|
if(is_ok)then
|
||||||
N_det = N_det_generators
|
N_det = N_det_generators
|
||||||
@ -135,14 +135,14 @@ subroutine is_a_good_candidate(threshold,is_ok,verbose)
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose)
|
subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_generators,is_ok,psi_coef_diagonalized_tmp, dressed_H_matrix,threshold,verbose,exit_loop)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
|
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
|
||||||
integer, intent(in) :: Ndet_generators
|
integer, intent(in) :: Ndet_generators
|
||||||
double precision, intent(in) :: threshold
|
double precision, intent(in) :: threshold
|
||||||
logical, intent(in) :: verbose
|
logical, intent(in) :: verbose
|
||||||
logical, intent(out) :: is_ok
|
logical, intent(out) :: is_ok,exit_loop
|
||||||
double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states)
|
double precision, intent(out) :: psi_coef_diagonalized_tmp(Ndet_generators,N_states)
|
||||||
double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators)
|
double precision, intent(inout) :: dressed_H_matrix(Ndet_generators, Ndet_generators)
|
||||||
|
|
||||||
@ -151,6 +151,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
|
|||||||
double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij
|
double precision :: eigvalues(Ndet_generators), eigvectors(Ndet_generators,Ndet_generators),hij
|
||||||
double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average
|
double precision :: psi_coef_ref(Ndet_generators,N_states),diag_h_mat_average,diag_h_mat_no_ref_average
|
||||||
logical :: is_a_ref_det(Ndet_generators)
|
logical :: is_a_ref_det(Ndet_generators)
|
||||||
|
exit_loop = .False.
|
||||||
|
|
||||||
is_a_ref_det = .False.
|
is_a_ref_det = .False.
|
||||||
do i = 1, N_det_generators
|
do i = 1, N_det_generators
|
||||||
@ -191,6 +192,7 @@ subroutine dress_H_matrix_from_psi_det_input(psi_det_generators_input,Ndet_gener
|
|||||||
if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then
|
if(number_of_holes(psi_det_generators_input(1,1,i)).eq.0 .and. number_of_particles(psi_det_generators_input(1,1,i)).eq.1)then
|
||||||
if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then
|
if(diag_h_mat_average - dressed_H_matrix(index_ref_generators_restart,index_ref_generators_restart) .gt.2.d0)then
|
||||||
is_ok = .False.
|
is_ok = .False.
|
||||||
|
exit_loop = .True.
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
endif
|
endif
|
||||||
|
@ -1,6 +1,6 @@
|
|||||||
program foboscf
|
program foboscf
|
||||||
implicit none
|
implicit none
|
||||||
call run_prepare
|
!call run_prepare
|
||||||
no_oa_or_av_opt = .True.
|
no_oa_or_av_opt = .True.
|
||||||
touch no_oa_or_av_opt
|
touch no_oa_or_av_opt
|
||||||
call routine_fobo_scf
|
call routine_fobo_scf
|
||||||
|
@ -38,6 +38,7 @@ subroutine FOBOCI_lmct_mlct_old_thr
|
|||||||
integer(bit_kind) , allocatable :: psi_singles(:,:,:)
|
integer(bit_kind) , allocatable :: psi_singles(:,:,:)
|
||||||
logical :: lmct
|
logical :: lmct
|
||||||
double precision, allocatable :: psi_singles_coef(:,:)
|
double precision, allocatable :: psi_singles_coef(:,:)
|
||||||
|
logical :: exit_loop
|
||||||
allocate( zero_bitmask(N_int,2) )
|
allocate( zero_bitmask(N_int,2) )
|
||||||
do i = 1, n_inact_orb
|
do i = 1, n_inact_orb
|
||||||
lmct = .True.
|
lmct = .True.
|
||||||
@ -55,7 +56,7 @@ subroutine FOBOCI_lmct_mlct_old_thr
|
|||||||
print*,'Passed set generators'
|
print*,'Passed set generators'
|
||||||
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
||||||
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
||||||
call is_a_good_candidate(threshold_lmct,is_ok,verbose)
|
call is_a_good_candidate(threshold_lmct,is_ok,verbose,exit_loop)
|
||||||
print*,'is_ok = ',is_ok
|
print*,'is_ok = ',is_ok
|
||||||
if(.not.is_ok)cycle
|
if(.not.is_ok)cycle
|
||||||
allocate(dressing_matrix(N_det_generators,N_det_generators))
|
allocate(dressing_matrix(N_det_generators,N_det_generators))
|
||||||
@ -161,9 +162,9 @@ subroutine FOBOCI_lmct_mlct_old_thr
|
|||||||
|
|
||||||
print*,'--------------------------'
|
print*,'--------------------------'
|
||||||
! First set the current generators to the one of restart
|
! First set the current generators to the one of restart
|
||||||
|
call check_symetry(i_particl_osoci,thr,test_sym)
|
||||||
call set_generators_to_generators_restart
|
call set_generators_to_generators_restart
|
||||||
call set_psi_det_to_generators
|
call set_psi_det_to_generators
|
||||||
call check_symetry(i_particl_osoci,thr,test_sym)
|
|
||||||
if(.not.test_sym)cycle
|
if(.not.test_sym)cycle
|
||||||
print*,'i_particl_osoci= ',i_particl_osoci
|
print*,'i_particl_osoci= ',i_particl_osoci
|
||||||
! Initialize the bitmask to the restart ones
|
! Initialize the bitmask to the restart ones
|
||||||
@ -180,9 +181,15 @@ subroutine FOBOCI_lmct_mlct_old_thr
|
|||||||
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
||||||
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
||||||
!! ! so all the mono excitation on the new generators
|
!! ! so all the mono excitation on the new generators
|
||||||
call is_a_good_candidate(threshold_mlct,is_ok,verbose)
|
call is_a_good_candidate(threshold_mlct,is_ok,verbose,exit_loop)
|
||||||
print*,'is_ok = ',is_ok
|
print*,'is_ok = ',is_ok
|
||||||
if(.not.is_ok)cycle
|
if(.not. is_ok)then
|
||||||
|
if(exit_loop)then
|
||||||
|
exit
|
||||||
|
else
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
endif
|
||||||
allocate(dressing_matrix(N_det_generators,N_det_generators))
|
allocate(dressing_matrix(N_det_generators,N_det_generators))
|
||||||
if(.not.do_it_perturbative)then
|
if(.not.do_it_perturbative)then
|
||||||
dressing_matrix = 0.d0
|
dressing_matrix = 0.d0
|
||||||
@ -234,7 +241,7 @@ subroutine FOBOCI_mlct_old
|
|||||||
double precision :: norm_tmp,norm_total
|
double precision :: norm_tmp,norm_total
|
||||||
logical :: test_sym
|
logical :: test_sym
|
||||||
double precision :: thr
|
double precision :: thr
|
||||||
logical :: verbose,is_ok
|
logical :: verbose,is_ok,exit_loop
|
||||||
verbose = .False.
|
verbose = .False.
|
||||||
thr = 1.d-12
|
thr = 1.d-12
|
||||||
allocate(unpaired_bitmask(N_int,2))
|
allocate(unpaired_bitmask(N_int,2))
|
||||||
@ -274,7 +281,7 @@ subroutine FOBOCI_mlct_old
|
|||||||
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
||||||
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
||||||
! ! so all the mono excitation on the new generators
|
! ! so all the mono excitation on the new generators
|
||||||
call is_a_good_candidate(threshold_mlct,is_ok,verbose)
|
call is_a_good_candidate(threshold_mlct,is_ok,verbose,exit_loop)
|
||||||
print*,'is_ok = ',is_ok
|
print*,'is_ok = ',is_ok
|
||||||
is_ok =.True.
|
is_ok =.True.
|
||||||
if(.not.is_ok)cycle
|
if(.not.is_ok)cycle
|
||||||
@ -308,7 +315,7 @@ subroutine FOBOCI_lmct_old
|
|||||||
double precision :: norm_tmp,norm_total
|
double precision :: norm_tmp,norm_total
|
||||||
logical :: test_sym
|
logical :: test_sym
|
||||||
double precision :: thr
|
double precision :: thr
|
||||||
logical :: verbose,is_ok
|
logical :: verbose,is_ok,exit_loop
|
||||||
verbose = .False.
|
verbose = .False.
|
||||||
thr = 1.d-12
|
thr = 1.d-12
|
||||||
allocate(unpaired_bitmask(N_int,2))
|
allocate(unpaired_bitmask(N_int,2))
|
||||||
@ -346,7 +353,7 @@ subroutine FOBOCI_lmct_old
|
|||||||
call set_generators_to_psi_det
|
call set_generators_to_psi_det
|
||||||
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
call set_bitmask_particl_as_input(reunion_of_bitmask)
|
||||||
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
call set_bitmask_hole_as_input(reunion_of_bitmask)
|
||||||
call is_a_good_candidate(threshold_lmct,is_ok,verbose)
|
call is_a_good_candidate(threshold_lmct,is_ok,verbose,exit_loop)
|
||||||
print*,'is_ok = ',is_ok
|
print*,'is_ok = ',is_ok
|
||||||
if(.not.is_ok)cycle
|
if(.not.is_ok)cycle
|
||||||
! ! so all the mono excitation on the new generators
|
! ! so all the mono excitation on the new generators
|
||||||
|
@ -212,8 +212,8 @@ subroutine update_density_matrix_osoci
|
|||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
do i = 1, mo_tot_num
|
do i = 1, mo_tot_num
|
||||||
do j = 1, mo_tot_num
|
do j = 1, mo_tot_num
|
||||||
one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha(i,j) - one_body_dm_mo_alpha_generators_restart(i,j))
|
one_body_dm_mo_alpha_osoci(i,j) = one_body_dm_mo_alpha_osoci(i,j) + (one_body_dm_mo_alpha_average(i,j) - one_body_dm_mo_alpha_generators_restart(i,j))
|
||||||
one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta(i,j) - one_body_dm_mo_beta_generators_restart(i,j))
|
one_body_dm_mo_beta_osoci(i,j) = one_body_dm_mo_beta_osoci(i,j) + (one_body_dm_mo_beta_average(i,j) - one_body_dm_mo_beta_generators_restart(i,j))
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
@ -588,14 +588,14 @@ end
|
|||||||
integer :: i
|
integer :: i
|
||||||
double precision :: accu_tot,accu_sd
|
double precision :: accu_tot,accu_sd
|
||||||
print*,'touched the one_body_dm_mo_beta'
|
print*,'touched the one_body_dm_mo_beta'
|
||||||
one_body_dm_mo_alpha = one_body_dm_mo_alpha_osoci
|
one_body_dm_mo_alpha_average = one_body_dm_mo_alpha_osoci
|
||||||
one_body_dm_mo_beta = one_body_dm_mo_beta_osoci
|
one_body_dm_mo_beta_average = one_body_dm_mo_beta_osoci
|
||||||
touch one_body_dm_mo_alpha one_body_dm_mo_beta
|
touch one_body_dm_mo_alpha one_body_dm_mo_beta
|
||||||
accu_tot = 0.d0
|
accu_tot = 0.d0
|
||||||
accu_sd = 0.d0
|
accu_sd = 0.d0
|
||||||
do i = 1, mo_tot_num
|
do i = 1, mo_tot_num
|
||||||
accu_tot += one_body_dm_mo_alpha(i,i) + one_body_dm_mo_beta(i,i)
|
accu_tot += one_body_dm_mo_alpha_average(i,i) + one_body_dm_mo_beta_average(i,i)
|
||||||
accu_sd += one_body_dm_mo_alpha(i,i) - one_body_dm_mo_beta(i,i)
|
accu_sd += one_body_dm_mo_alpha_average(i,i) - one_body_dm_mo_beta_average(i,i)
|
||||||
enddo
|
enddo
|
||||||
print*,'accu_tot = ',accu_tot
|
print*,'accu_tot = ',accu_tot
|
||||||
print*,'accu_sdt = ',accu_sd
|
print*,'accu_sdt = ',accu_sd
|
||||||
|
@ -10,14 +10,28 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)]
|
|||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)]
|
||||||
|
implicit none
|
||||||
|
integer :: i
|
||||||
|
double precision :: energies(N_states_diag)
|
||||||
|
do i = 1, N_states
|
||||||
|
call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states_diag,i)
|
||||||
|
energy_cas_dyall_no_exchange(i) = energies(i)
|
||||||
|
print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i)
|
||||||
|
enddo
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
|
BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
integer :: ispin
|
integer :: ispin
|
||||||
integer :: orb, hole_particle,spin_exc
|
integer :: orb, hole_particle,spin_exc
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,N_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(N_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
|
||||||
integer :: iorb
|
integer :: iorb
|
||||||
@ -45,6 +59,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)]
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -54,9 +69,10 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
|
|||||||
integer :: ispin
|
integer :: ispin
|
||||||
integer :: orb, hole_particle,spin_exc
|
integer :: orb, hole_particle,spin_exc
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb
|
integer :: iorb
|
||||||
integer :: state_target
|
integer :: state_target
|
||||||
@ -83,6 +99,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)]
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -93,9 +110,10 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
|
|||||||
integer :: orb_i, hole_particle_i,spin_exc_i
|
integer :: orb_i, hole_particle_i,spin_exc_i
|
||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: state_target
|
integer :: state_target
|
||||||
@ -131,6 +149,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states)
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -141,9 +160,10 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
|
|||||||
integer :: orb_i, hole_particle_i,spin_exc_i
|
integer :: orb_i, hole_particle_i,spin_exc_i
|
||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: state_target
|
integer :: state_target
|
||||||
@ -178,6 +198,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states)
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -188,10 +209,11 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
|
|||||||
integer :: orb_i, hole_particle_i,spin_exc_i
|
integer :: orb_i, hole_particle_i,spin_exc_i
|
||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: state_target
|
integer :: state_target
|
||||||
double precision :: energies(n_states_diag)
|
double precision :: energies(n_states_diag)
|
||||||
@ -219,8 +241,13 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
|
|||||||
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
||||||
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, &
|
||||||
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states_diag)
|
||||||
|
if(orb_i == orb_j .and. ispin .ne. jspin)then
|
||||||
|
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
||||||
|
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target)
|
||||||
|
else
|
||||||
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
||||||
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target)
|
||||||
|
endif
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -229,6 +256,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2
|
|||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
|
BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_act_orb,2,2,2,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
@ -237,9 +265,10 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
|
|||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
integer :: orb_k, hole_particle_k,spin_exc_k
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: korb
|
integer :: korb
|
||||||
@ -286,6 +315,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -297,9 +327,10 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
|
|||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
integer :: orb_k, hole_particle_k,spin_exc_k
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: korb
|
integer :: korb
|
||||||
@ -345,6 +376,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -356,9 +388,10 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
|
|||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
integer :: orb_k, hole_particle_k,spin_exc_k
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: korb
|
integer :: korb
|
||||||
@ -404,6 +437,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -415,9 +449,10 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
|
|||||||
integer :: orb_j, hole_particle_j,spin_exc_j
|
integer :: orb_j, hole_particle_j,spin_exc_j
|
||||||
integer :: orb_k, hole_particle_k,spin_exc_k
|
integer :: orb_k, hole_particle_k,spin_exc_k
|
||||||
double precision :: norm_out(N_states_diag)
|
double precision :: norm_out(N_states_diag)
|
||||||
integer(bit_kind) :: psi_in_out(N_int,2,n_det)
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
double precision :: psi_in_out_coef(n_det,N_states_diag)
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
integer :: iorb,jorb
|
integer :: iorb,jorb
|
||||||
integer :: korb
|
integer :: korb
|
||||||
@ -463,5 +498,617 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt, (n_inact_orb,n_virt_orb,N_States)]
|
||||||
|
&BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt_norm, (n_inact_orb,n_virt_orb,N_States,2)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,vorb,j
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: orb_i, hole_particle_i
|
||||||
|
integer :: orb_v
|
||||||
|
double precision :: norm_out(N_states_diag)
|
||||||
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
|
integer :: iorb,jorb,i_ok
|
||||||
|
integer :: state_target
|
||||||
|
double precision :: energies(n_states_diag)
|
||||||
|
double precision :: hij
|
||||||
|
double precision :: norm(N_states,2),norm_no_inv(N_states,2),norm_bis(N_states,2)
|
||||||
|
double precision :: energies_alpha_beta(N_states,2)
|
||||||
|
|
||||||
|
|
||||||
|
double precision :: thresh_norm
|
||||||
|
|
||||||
|
thresh_norm = 1.d-10
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
do vorb = 1,n_virt_orb
|
||||||
|
orb_v = list_virt(vorb)
|
||||||
|
do iorb = 1, n_inact_orb
|
||||||
|
orb_i = list_inact(iorb)
|
||||||
|
norm = 0.d0
|
||||||
|
norm_bis = 0.d0
|
||||||
|
do ispin = 1,2
|
||||||
|
do state_target =1 , N_states
|
||||||
|
one_anhil_one_creat_inact_virt_norm(iorb,vorb,state_target,ispin) = 0.d0
|
||||||
|
enddo
|
||||||
|
do i = 1, n_det
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||||
|
enddo
|
||||||
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
|
||||||
|
if(i_ok.ne.1)then
|
||||||
|
print*, orb_i,orb_v
|
||||||
|
call debug_det(psi_in_out,N_int)
|
||||||
|
print*, 'pb, i_ok ne 0 !!!'
|
||||||
|
endif
|
||||||
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
||||||
|
do j = 1, n_states
|
||||||
|
double precision :: coef,contrib
|
||||||
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
||||||
|
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
||||||
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do j = 1, N_states
|
||||||
|
if (dabs(norm(j,ispin)) .le. thresh_norm)then
|
||||||
|
norm(j,ispin) = 0.d0
|
||||||
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
||||||
|
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 0.d0
|
||||||
|
else
|
||||||
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
||||||
|
one_anhil_one_creat_inact_virt_norm(iorb,vorb,j,ispin) = 1.d0 / norm(j,ispin)
|
||||||
|
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do i = 1, N_det
|
||||||
|
do j = 1, N_states
|
||||||
|
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
|
||||||
|
norm_bis(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
||||||
|
enddo
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_active(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_active(j,2,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
energies_alpha_beta(state_target, ispin) = - mo_bielec_integral_jj_exchange(orb_i,orb_v)
|
||||||
|
! energies_alpha_beta(state_target, ispin) = 0.d0
|
||||||
|
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
|
||||||
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
||||||
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo ! ispin
|
||||||
|
do state_target = 1, N_states
|
||||||
|
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
||||||
|
! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.5d0 * &
|
||||||
|
! ( energy_cas_dyall(state_target) - energies_alpha_beta(state_target,1) + &
|
||||||
|
! energy_cas_dyall(state_target) - energies_alpha_beta(state_target,2) )
|
||||||
|
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
|
||||||
|
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
|
||||||
|
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall(state_target) - &
|
||||||
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
||||||
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
||||||
|
else
|
||||||
|
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = 0.d0
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_States)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,iorb,j
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: orb_i, hole_particle_i
|
||||||
|
double precision :: norm_out(N_states_diag)
|
||||||
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
|
integer :: jorb,i_ok,aorb,orb_a
|
||||||
|
integer :: state_target
|
||||||
|
double precision :: energies(n_states_diag)
|
||||||
|
double precision :: hij
|
||||||
|
double precision :: norm(N_states,2),norm_no_inv(N_states,2)
|
||||||
|
double precision :: energies_alpha_beta(N_states,2)
|
||||||
|
double precision :: norm_alpha_beta(N_states,2)
|
||||||
|
|
||||||
|
double precision :: thresh_norm
|
||||||
|
|
||||||
|
thresh_norm = 1.d-10
|
||||||
|
|
||||||
|
do aorb = 1,n_act_orb
|
||||||
|
orb_a = list_act(aorb)
|
||||||
|
do iorb = 1, n_inact_orb
|
||||||
|
orb_i = list_inact(iorb)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
one_anhil_inact(iorb,aorb,state_target) = 0.d0
|
||||||
|
enddo
|
||||||
|
norm_alpha_beta = 0.d0
|
||||||
|
norm = 0.d0
|
||||||
|
norm_bis = 0.d0
|
||||||
|
do ispin = 1,2
|
||||||
|
do i = 1, n_det
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||||
|
enddo
|
||||||
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_a,ispin,i_ok)
|
||||||
|
if(i_ok.ne.1)then
|
||||||
|
do j = 1, n_states
|
||||||
|
psi_in_out_coef(i,j) = 0.d0
|
||||||
|
enddo
|
||||||
|
else
|
||||||
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
||||||
|
do j = 1, n_states
|
||||||
|
double precision :: coef,contrib
|
||||||
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
||||||
|
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
||||||
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do j = 1, N_states
|
||||||
|
if (dabs(norm(j,ispin)) .le. thresh_norm)then
|
||||||
|
norm(j,ispin) = 0.d0
|
||||||
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
||||||
|
else
|
||||||
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
||||||
|
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
double precision :: norm_bis(N_states,2)
|
||||||
|
do i = 1, N_det
|
||||||
|
do j = 1, N_states
|
||||||
|
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
|
||||||
|
norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j)
|
||||||
|
enddo
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
|
||||||
|
psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
energies_alpha_beta(state_target, ispin) = 0.d0
|
||||||
|
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
|
||||||
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
||||||
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo ! ispin
|
||||||
|
do state_target = 1, N_states
|
||||||
|
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
||||||
|
one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall(state_target) - &
|
||||||
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
||||||
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
||||||
|
else
|
||||||
|
one_anhil_inact(iorb,aorb,state_target) = 0.d0
|
||||||
|
endif
|
||||||
|
! print*, '********'
|
||||||
|
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
|
||||||
|
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
|
||||||
|
! print*, one_anhil_inact(iorb,aorb,state_target)
|
||||||
|
! print*, one_creat(aorb,1,state_target)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_States)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,vorb,j
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: orb_i, hole_particle_i
|
||||||
|
integer :: orb_v
|
||||||
|
double precision :: norm_out(N_states_diag)
|
||||||
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag))
|
||||||
|
|
||||||
|
integer :: iorb,jorb,i_ok,aorb,orb_a
|
||||||
|
integer :: state_target
|
||||||
|
double precision :: energies(n_states_diag)
|
||||||
|
double precision :: hij
|
||||||
|
double precision :: norm(N_states,2),norm_no_inv(N_states,2)
|
||||||
|
double precision :: energies_alpha_beta(N_states,2)
|
||||||
|
double precision :: norm_alpha_beta(N_states,2)
|
||||||
|
|
||||||
|
double precision :: thresh_norm
|
||||||
|
|
||||||
|
thresh_norm = 1.d-10
|
||||||
|
|
||||||
|
do aorb = 1,n_act_orb
|
||||||
|
orb_a = list_act(aorb)
|
||||||
|
do vorb = 1, n_virt_orb
|
||||||
|
orb_v = list_virt(vorb)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
one_creat_virt(aorb,vorb,state_target) = 0.d0
|
||||||
|
enddo
|
||||||
|
norm_alpha_beta = 0.d0
|
||||||
|
norm = 0.d0
|
||||||
|
norm_bis = 0.d0
|
||||||
|
do ispin = 1,2
|
||||||
|
do i = 1, n_det
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||||
|
enddo
|
||||||
|
call do_mono_excitation(psi_in_out(1,1,i),orb_a,orb_v,ispin,i_ok)
|
||||||
|
if(i_ok.ne.1)then
|
||||||
|
do j = 1, n_states
|
||||||
|
psi_in_out_coef(i,j) = 0.d0
|
||||||
|
enddo
|
||||||
|
else
|
||||||
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
||||||
|
do j = 1, n_states
|
||||||
|
double precision :: coef,contrib
|
||||||
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
||||||
|
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
||||||
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do j = 1, N_states
|
||||||
|
if (dabs(norm(j,ispin)) .le. thresh_norm)then
|
||||||
|
norm(j,ispin) = 0.d0
|
||||||
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
||||||
|
else
|
||||||
|
norm_no_inv(j,ispin) = norm(j,ispin)
|
||||||
|
norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin))
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
double precision :: norm_bis(N_states,2)
|
||||||
|
do i = 1, N_det
|
||||||
|
do j = 1, N_states
|
||||||
|
psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin)
|
||||||
|
norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j)
|
||||||
|
enddo
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = iand(psi_in_out(j,1,i),cas_bitmask(j,1,1))
|
||||||
|
psi_in_out(j,2,i) = iand(psi_in_out(j,2,i),cas_bitmask(j,1,1))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
energies_alpha_beta(state_target, ispin) = 0.d0
|
||||||
|
if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then
|
||||||
|
call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states_diag,state_target)
|
||||||
|
! print*, energies(state_target)
|
||||||
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo ! ispin
|
||||||
|
do state_target = 1, N_states
|
||||||
|
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
||||||
|
one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall(state_target) - &
|
||||||
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
||||||
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
||||||
|
else
|
||||||
|
one_creat_virt(aorb,vorb,state_target) = 0.d0
|
||||||
|
endif
|
||||||
|
! print*, '********'
|
||||||
|
! print*, energies_alpha_beta(state_target,1) , energies_alpha_beta(state_target,2)
|
||||||
|
! print*, norm_bis(state_target,1) , norm_bis(state_target,2)
|
||||||
|
! print*, one_creat_virt(aorb,vorb,state_target)
|
||||||
|
! print*, one_anhil(aorb,1,state_target)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef)
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, one_anhil_one_creat_inact_virt_bis, (n_inact_orb,n_virt_orb,N_det,N_States)]
|
||||||
|
&BEGIN_PROVIDER [ double precision, corr_e_from_1h1p, (N_States)]
|
||||||
|
implicit none
|
||||||
|
integer :: i,vorb,j
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: orb_i, hole_particle_i
|
||||||
|
integer :: orb_v
|
||||||
|
double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det)
|
||||||
|
double precision :: delta_e_inact_virt(N_states)
|
||||||
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
|
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:)
|
||||||
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1))
|
||||||
|
allocate (eigenvectors(size(H_matrix,1),N_det+1))
|
||||||
|
allocate (eigenvalues(N_det+1))
|
||||||
|
|
||||||
|
integer :: iorb,jorb,i_ok
|
||||||
|
integer :: state_target
|
||||||
|
double precision :: energies(n_states_diag)
|
||||||
|
double precision :: hij
|
||||||
|
double precision :: energies_alpha_beta(N_states,2)
|
||||||
|
|
||||||
|
|
||||||
|
double precision :: accu(N_states),norm
|
||||||
|
double precision :: amplitudes_alpha_beta(N_det,2)
|
||||||
|
double precision :: delta_e_alpha_beta(N_det,2)
|
||||||
|
|
||||||
|
corr_e_from_1h1p = 0.d0
|
||||||
|
do vorb = 1,n_virt_orb
|
||||||
|
orb_v = list_virt(vorb)
|
||||||
|
do iorb = 1, n_inact_orb
|
||||||
|
orb_i = list_inact(iorb)
|
||||||
|
! print*, '---------------------------------------------------------------------------'
|
||||||
|
do j = 1, N_states
|
||||||
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(orb_i,j) &
|
||||||
|
- fock_virt_total_spin_trace(orb_v,j)
|
||||||
|
enddo
|
||||||
|
do ispin = 1,2
|
||||||
|
do i = 1, n_det
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||||
|
enddo
|
||||||
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
|
||||||
|
if(i_ok.ne.1)then
|
||||||
|
print*, orb_i,orb_v
|
||||||
|
call debug_det(psi_in_out,N_int)
|
||||||
|
print*, 'pb, i_ok ne 0 !!!'
|
||||||
|
endif
|
||||||
|
interact_psi0(i) = 0.d0
|
||||||
|
do j = 1 , N_det
|
||||||
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
|
||||||
|
interact_psi0(i) += hij * psi_coef(j,1)
|
||||||
|
enddo
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_active(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_active(j,2,i)
|
||||||
|
enddo
|
||||||
|
call i_H_j_dyall(psi_active(1,1,i),psi_active(1,1,i),N_int,hij)
|
||||||
|
diag_elem(i) = hij
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
! Building the Hamiltonian matrix
|
||||||
|
H_matrix(1,1) = energy_cas_dyall(state_target)
|
||||||
|
do i = 1, N_det
|
||||||
|
! interaction with psi0
|
||||||
|
H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target)
|
||||||
|
H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target)
|
||||||
|
! diagonal elements
|
||||||
|
H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target)
|
||||||
|
! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1)
|
||||||
|
do j = i+1, N_det
|
||||||
|
call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij)
|
||||||
|
H_matrix(i+1,j+1) = hij !0.d0 !
|
||||||
|
H_matrix(j+1,i+1) = hij !0.d0 !
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
print*, '***'
|
||||||
|
do i = 1, N_det+1
|
||||||
|
write(*,'(100(F16.10,X))')H_matrix(i,:)
|
||||||
|
enddo
|
||||||
|
call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1)
|
||||||
|
corr_e_from_1h1p(state_target) += eigenvalues(1) - energy_cas_dyall(state_target)
|
||||||
|
norm = 0.d0
|
||||||
|
do i = 1, N_det
|
||||||
|
psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1)
|
||||||
|
!! if(dabs(psi_coef(i,state_target)*) .gt. 1.d-8)then
|
||||||
|
if(dabs(psi_in_out_coef(i,state_target)) .gt. 1.d-8)then
|
||||||
|
! if(dabs(interact_psi0(i)) .gt. 1.d-8)then
|
||||||
|
delta_e_alpha_beta(i,ispin) = H_matrix(1,i+1) / psi_in_out_coef(i,state_target)
|
||||||
|
! delta_e_alpha_beta(i,ispin) = interact_psi0(i) / psi_in_out_coef(i,state_target)
|
||||||
|
amplitudes_alpha_beta(i,ispin) = psi_in_out_coef(i,state_target) / psi_coef(i,state_target)
|
||||||
|
else
|
||||||
|
amplitudes_alpha_beta(i,ispin) = 0.d0
|
||||||
|
delta_e_alpha_beta(i,ispin) = delta_e_inact_virt(state_target)
|
||||||
|
endif
|
||||||
|
!! one_anhil_one_creat_inact_virt_bis(iorb,vorb,i,ispin,state_target) = amplitudes_alpha_beta(i,ispin)
|
||||||
|
norm += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target)
|
||||||
|
enddo
|
||||||
|
print*, 'Coef '
|
||||||
|
write(*,'(100(X,F16.10))')psi_coef(1:N_det,state_target)
|
||||||
|
write(*,'(100(X,F16.10))')psi_in_out_coef(:,state_target)
|
||||||
|
double precision :: coef_tmp(N_det)
|
||||||
|
do i = 1, N_det
|
||||||
|
coef_tmp(i) = psi_coef(i,1) * interact_psi0(i) / delta_e_alpha_beta(i,ispin)
|
||||||
|
enddo
|
||||||
|
write(*,'(100(X,F16.10))')coef_tmp(:)
|
||||||
|
print*, 'naked interactions'
|
||||||
|
write(*,'(100(X,F16.10))')interact_psi0(:)
|
||||||
|
print*, ''
|
||||||
|
|
||||||
|
print*, 'norm ',norm
|
||||||
|
norm = 1.d0/(norm)
|
||||||
|
accu(state_target) = 0.d0
|
||||||
|
do i = 1, N_det
|
||||||
|
accu(state_target) += psi_in_out_coef(i,state_target) * psi_in_out_coef(i,state_target) * H_matrix(i+1,i+1)
|
||||||
|
do j = i+1, N_det
|
||||||
|
accu(state_target) += 2.d0 * psi_in_out_coef(i,state_target) * psi_in_out_coef(j,state_target) * H_matrix(i+1,j+1)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
accu(state_target) = accu(state_target) * norm
|
||||||
|
print*, delta_e_inact_virt(state_target)
|
||||||
|
print*, eigenvalues(1),accu(state_target),eigenvectors(1,1)
|
||||||
|
print*, energy_cas_dyall(state_target) - accu(state_target), one_anhil_one_creat_inact_virt(iorb,vorb,state_target) + delta_e_inact_virt(state_target)
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo ! ispin
|
||||||
|
do state_target = 1, N_states
|
||||||
|
do i = 1, N_det
|
||||||
|
one_anhil_one_creat_inact_virt_bis(iorb,vorb,i,state_target) = 0.5d0 * &
|
||||||
|
( delta_e_alpha_beta(i,1) + delta_e_alpha_beta(i,1))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
print*, '***'
|
||||||
|
write(*,'(100(X,F16.10))')
|
||||||
|
write(*,'(100(X,F16.10))')delta_e_alpha_beta(:,2)
|
||||||
|
! write(*,'(100(X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,1,:)
|
||||||
|
! write(*,'(100(X,F16.10))')one_anhil_one_creat_inact_virt_bis(iorb,vorb,:,2,:)
|
||||||
|
print*, '---------------------------------------------------------------------------'
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef,H_matrix,eigenvectors,eigenvalues)
|
||||||
|
print*, 'corr_e_from_1h1p,',corr_e_from_1h1p(:)
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from_1h1p_singles)
|
||||||
|
implicit none
|
||||||
|
double precision , intent(inout) :: matrix_1h1p(N_det,N_det,N_states)
|
||||||
|
double precision , intent(out) :: e_corr_from_1h1p_singles(N_states)
|
||||||
|
integer :: i,vorb,j
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: orb_i, hole_particle_i
|
||||||
|
integer :: orb_v
|
||||||
|
double precision :: norm_out(N_states_diag),diag_elem(N_det),interact_psi0(N_det)
|
||||||
|
double precision :: delta_e_inact_virt(N_states)
|
||||||
|
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
|
||||||
|
double precision, allocatable :: psi_in_out_coef(:,:)
|
||||||
|
double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:)
|
||||||
|
double precision, allocatable :: delta_e_det(:,:)
|
||||||
|
use bitmasks
|
||||||
|
allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states_diag),H_matrix(N_det+1,N_det+1))
|
||||||
|
allocate (eigenvectors(size(H_matrix,1),N_det+1))
|
||||||
|
allocate (eigenvalues(N_det+1),interact_cas(N_det,N_det))
|
||||||
|
allocate (delta_e_det(N_det,N_det))
|
||||||
|
|
||||||
|
integer :: iorb,jorb,i_ok
|
||||||
|
integer :: state_target
|
||||||
|
double precision :: energies(n_states_diag)
|
||||||
|
double precision :: hij
|
||||||
|
double precision :: energies_alpha_beta(N_states,2)
|
||||||
|
double precision :: lamda_pt2(N_det)
|
||||||
|
|
||||||
|
|
||||||
|
double precision :: accu(N_states),norm
|
||||||
|
double precision :: amplitudes_alpha_beta(N_det,2)
|
||||||
|
double precision :: delta_e_alpha_beta(N_det,2)
|
||||||
|
double precision :: coef_array(N_states)
|
||||||
|
double precision :: coef_perturb(N_det)
|
||||||
|
double precision :: coef_perturb_bis(N_det)
|
||||||
|
|
||||||
|
do vorb = 1,n_virt_orb
|
||||||
|
orb_v = list_virt(vorb)
|
||||||
|
do iorb = 1, n_inact_orb
|
||||||
|
orb_i = list_inact(iorb)
|
||||||
|
do j = 1, N_states
|
||||||
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(orb_i,j) &
|
||||||
|
- fock_virt_total_spin_trace(orb_v,j)
|
||||||
|
enddo
|
||||||
|
do ispin = 1,2
|
||||||
|
do i = 1, n_det
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_det(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_det(j,2,i)
|
||||||
|
enddo
|
||||||
|
call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok)
|
||||||
|
if(i_ok.ne.1)then
|
||||||
|
print*, orb_i,orb_v
|
||||||
|
call debug_det(psi_in_out,N_int)
|
||||||
|
print*, 'pb, i_ok ne 0 !!!'
|
||||||
|
endif
|
||||||
|
interact_psi0(i) = 0.d0
|
||||||
|
do j = 1 , N_det
|
||||||
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij)
|
||||||
|
call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),coef_array,hij,delta_e_det(i,j))
|
||||||
|
interact_cas(i,j) = hij
|
||||||
|
interact_psi0(i) += hij * psi_coef(j,1)
|
||||||
|
enddo
|
||||||
|
do j = 1, N_int
|
||||||
|
psi_in_out(j,1,i) = psi_active(j,1,i)
|
||||||
|
psi_in_out(j,2,i) = psi_active(j,2,i)
|
||||||
|
enddo
|
||||||
|
call i_H_j_dyall(psi_active(1,1,i),psi_active(1,1,i),N_int,hij)
|
||||||
|
diag_elem(i) = hij
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
! Building the Hamiltonian matrix
|
||||||
|
H_matrix(1,1) = energy_cas_dyall(state_target)
|
||||||
|
do i = 1, N_det
|
||||||
|
! interaction with psi0
|
||||||
|
H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target)
|
||||||
|
H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target)
|
||||||
|
! diagonal elements
|
||||||
|
H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target)
|
||||||
|
! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1)
|
||||||
|
do j = i+1, N_det
|
||||||
|
call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij)
|
||||||
|
H_matrix(i+1,j+1) = hij !0.d0 !
|
||||||
|
H_matrix(j+1,i+1) = hij !0.d0 !
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1)
|
||||||
|
e_corr_from_1h1p_singles(state_target) += eigenvalues(1) - energy_cas_dyall(state_target)
|
||||||
|
|
||||||
|
do i = 1, N_det
|
||||||
|
psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1)
|
||||||
|
coef_perturb(i) = 0.d0
|
||||||
|
do j = 1, N_det
|
||||||
|
coef_perturb(i) += psi_coef(j,state_target) * interact_cas(i,j) *1.d0/delta_e_det(i,j)
|
||||||
|
enddo
|
||||||
|
coef_perturb_bis(i) = interact_psi0(i) / (eigenvalues(1) - H_matrix(i+1,i+1))
|
||||||
|
if(dabs(interact_psi0(i)) .gt. 1.d-12)then
|
||||||
|
lamda_pt2(i) = psi_in_out_coef(i,state_target) / interact_psi0(i)
|
||||||
|
else
|
||||||
|
lamda_pt2(i) =energy_cas_dyall(state_target) - H_matrix(i+1,i+1)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
if(dabs(eigenvalues(1) - energy_cas_dyall(state_target)).gt.1.d-10)then
|
||||||
|
print*, ''
|
||||||
|
do i = 1, N_det+1
|
||||||
|
write(*,'(100(F16.10))') H_matrix(i,:)
|
||||||
|
enddo
|
||||||
|
accu = 0.d0
|
||||||
|
do i = 1, N_det
|
||||||
|
accu(state_target) += psi_in_out_coef(i,state_target) * interact_psi0(i)
|
||||||
|
enddo
|
||||||
|
print*, ''
|
||||||
|
print*, 'e corr diagonal ',accu(state_target)
|
||||||
|
accu = 0.d0
|
||||||
|
do i = 1, N_det
|
||||||
|
accu(state_target) += coef_perturb(i) * interact_psi0(i)
|
||||||
|
enddo
|
||||||
|
print*, 'e corr perturb ',accu(state_target)
|
||||||
|
accu = 0.d0
|
||||||
|
do i = 1, N_det
|
||||||
|
accu(state_target) += coef_perturb_bis(i) * interact_psi0(i)
|
||||||
|
enddo
|
||||||
|
print*, 'e corr perturb EN',accu(state_target)
|
||||||
|
print*, ''
|
||||||
|
print*, 'coef diagonalized'
|
||||||
|
write(*,'(100(F16.10,X))')psi_in_out_coef(:,state_target)
|
||||||
|
print*, 'coef_perturb'
|
||||||
|
write(*,'(100(F16.10,X))')coef_perturb(:)
|
||||||
|
print*, 'coef_perturb EN'
|
||||||
|
write(*,'(100(F16.10,X))')coef_perturb_bis(:)
|
||||||
|
endif
|
||||||
|
integer :: k
|
||||||
|
do k = 1, N_det
|
||||||
|
do i = 1, N_det
|
||||||
|
matrix_1h1p(i,i,state_target) += interact_cas(k,i) * interact_cas(k,i) * lamda_pt2(k)
|
||||||
|
do j = i+1, N_det
|
||||||
|
matrix_1h1p(i,j,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k)
|
||||||
|
matrix_1h1p(j,i,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo ! ispin
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
deallocate(psi_in_out,psi_in_out_coef,H_matrix,eigenvectors,eigenvalues,interact_cas)
|
||||||
|
deallocate(delta_e_det)
|
||||||
|
end
|
||||||
|
@ -382,8 +382,6 @@ subroutine i_H_j_dyall(key_i,key_j,Nint,hij)
|
|||||||
|
|
||||||
endif
|
endif
|
||||||
hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) )
|
hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) )
|
||||||
! hij = phase*(mo_mono_elec_integral(m,p) ) ! + fock_operator_active_from_core_inact(m,p) )
|
|
||||||
! hij = 0.d0
|
|
||||||
|
|
||||||
case (0)
|
case (0)
|
||||||
hij = diag_H_mat_elem_no_elec_check(key_i,Nint)
|
hij = diag_H_mat_elem_no_elec_check(key_i,Nint)
|
||||||
@ -422,3 +420,289 @@ subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coe
|
|||||||
energies(state_target) = accu
|
energies(state_target) = accu
|
||||||
deallocate(psi_coef_tmp)
|
deallocate(psi_coef_tmp)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
double precision function coulomb_value_no_check(det_in,Nint)
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Computes <i|H|i>
|
||||||
|
END_DOC
|
||||||
|
integer,intent(in) :: Nint
|
||||||
|
integer(bit_kind),intent(in) :: det_in(Nint,2)
|
||||||
|
|
||||||
|
integer :: i, j, iorb, jorb
|
||||||
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
|
integer :: elec_num_tab_local(2)
|
||||||
|
|
||||||
|
double precision :: core_act
|
||||||
|
double precision :: alpha_alpha
|
||||||
|
double precision :: alpha_beta
|
||||||
|
double precision :: beta_beta
|
||||||
|
double precision :: mono_elec
|
||||||
|
core_act = 0.d0
|
||||||
|
alpha_alpha = 0.d0
|
||||||
|
alpha_beta = 0.d0
|
||||||
|
beta_beta = 0.d0
|
||||||
|
mono_elec = 0.d0
|
||||||
|
|
||||||
|
coulomb_value_no_check = 0.d0
|
||||||
|
call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), N_int)
|
||||||
|
call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), N_int)
|
||||||
|
! alpha - alpha
|
||||||
|
do i = 1, elec_num_tab_local(1)
|
||||||
|
iorb = occ(i,1)
|
||||||
|
do j = i+1, elec_num_tab_local(1)
|
||||||
|
jorb = occ(j,1)
|
||||||
|
coulomb_value_no_check += mo_bielec_integral_jj_anti(jorb,iorb)
|
||||||
|
alpha_alpha += mo_bielec_integral_jj_anti(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! beta - beta
|
||||||
|
do i = 1, elec_num_tab_local(2)
|
||||||
|
iorb = occ(i,2)
|
||||||
|
do j = i+1, elec_num_tab_local(2)
|
||||||
|
jorb = occ(j,2)
|
||||||
|
coulomb_value_no_check += mo_bielec_integral_jj_anti(jorb,iorb)
|
||||||
|
beta_beta += mo_bielec_integral_jj_anti(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
! alpha - beta
|
||||||
|
do i = 1, elec_num_tab_local(2)
|
||||||
|
iorb = occ(i,2)
|
||||||
|
do j = 1, elec_num_tab_local(1)
|
||||||
|
jorb = occ(j,1)
|
||||||
|
coulomb_value_no_check += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
alpha_beta += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine i_H_j_dyall_no_exchange(key_i,key_j,Nint,hij)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Returns <i|H|j> where i and j are determinants
|
||||||
|
END_DOC
|
||||||
|
integer, intent(in) :: Nint
|
||||||
|
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
|
||||||
|
double precision, intent(out) :: hij
|
||||||
|
|
||||||
|
integer :: exc(0:2,2,2)
|
||||||
|
integer :: degree
|
||||||
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
|
integer :: m,n,p,q
|
||||||
|
integer :: i,j,k
|
||||||
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
|
double precision :: diag_H_mat_elem_no_elec_check_no_exchange, phase,phase_2
|
||||||
|
integer :: n_occ_ab(2)
|
||||||
|
logical :: has_mipi(Nint*bit_kind_size)
|
||||||
|
double precision :: mipi(Nint*bit_kind_size)
|
||||||
|
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
|
||||||
|
|
||||||
|
ASSERT (Nint > 0)
|
||||||
|
ASSERT (Nint == N_int)
|
||||||
|
|
||||||
|
hij = 0.d0
|
||||||
|
!DIR$ FORCEINLINE
|
||||||
|
call get_excitation_degree(key_i,key_j,degree,Nint)
|
||||||
|
select case (degree)
|
||||||
|
case (2)
|
||||||
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
|
if (exc(0,1,1) == 1) then
|
||||||
|
! Mono alpha, mono beta
|
||||||
|
if(exc(1,1,1) == exc(1,2,2) .and. exc(1,2,1) == exc(1,1,2))then
|
||||||
|
hij = 0.d0
|
||||||
|
else
|
||||||
|
hij = phase*get_mo_bielec_integral_schwartz( &
|
||||||
|
exc(1,1,1), &
|
||||||
|
exc(1,1,2), &
|
||||||
|
exc(1,2,1), &
|
||||||
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
|
endif
|
||||||
|
else if (exc(0,1,1) == 2) then
|
||||||
|
! Double alpha
|
||||||
|
hij = phase*get_mo_bielec_integral_schwartz( &
|
||||||
|
exc(1,1,1), &
|
||||||
|
exc(2,1,1), &
|
||||||
|
exc(1,2,1), &
|
||||||
|
exc(2,2,1) ,mo_integrals_map)
|
||||||
|
else if (exc(0,1,2) == 2) then
|
||||||
|
! Double beta
|
||||||
|
hij = phase*get_mo_bielec_integral_schwartz( &
|
||||||
|
exc(1,1,2), &
|
||||||
|
exc(2,1,2), &
|
||||||
|
exc(1,2,2), &
|
||||||
|
exc(2,2,2) ,mo_integrals_map)
|
||||||
|
endif
|
||||||
|
case (1)
|
||||||
|
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
|
!DIR$ FORCEINLINE
|
||||||
|
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
|
||||||
|
has_mipi = .False.
|
||||||
|
if (exc(0,1,1) == 1) then
|
||||||
|
! Mono alpha
|
||||||
|
m = exc(1,1,1)
|
||||||
|
p = exc(1,2,1)
|
||||||
|
do k = 1, n_occ_ab(1)
|
||||||
|
i = occ(k,1)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do k = 1, n_occ_ab(2)
|
||||||
|
i = occ(k,2)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do k = 1, n_occ_ab(1)
|
||||||
|
hij = hij + mipi(occ(k,1))
|
||||||
|
enddo
|
||||||
|
do k = 1, n_occ_ab(2)
|
||||||
|
hij = hij + mipi(occ(k,2))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
else
|
||||||
|
! Mono beta
|
||||||
|
m = exc(1,1,2)
|
||||||
|
p = exc(1,2,2)
|
||||||
|
do k = 1, n_occ_ab(2)
|
||||||
|
i = occ(k,2)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do k = 1, n_occ_ab(1)
|
||||||
|
i = occ(k,1)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do k = 1, n_occ_ab(1)
|
||||||
|
hij = hij + mipi(occ(k,1))
|
||||||
|
enddo
|
||||||
|
do k = 1, n_occ_ab(2)
|
||||||
|
hij = hij + mipi(occ(k,2))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
endif
|
||||||
|
hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) )
|
||||||
|
|
||||||
|
case (0)
|
||||||
|
hij = diag_H_mat_elem_no_elec_check_no_exchange(key_i,Nint)
|
||||||
|
end select
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
double precision function diag_H_mat_elem_no_elec_check_no_exchange(det_in,Nint)
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Computes <i|H|i>
|
||||||
|
END_DOC
|
||||||
|
integer,intent(in) :: Nint
|
||||||
|
integer(bit_kind),intent(in) :: det_in(Nint,2)
|
||||||
|
|
||||||
|
integer :: i, j, iorb, jorb
|
||||||
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
|
integer :: elec_num_tab_local(2)
|
||||||
|
|
||||||
|
double precision :: core_act_exchange(2)
|
||||||
|
core_act_exchange = 0.d0
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange = 0.d0
|
||||||
|
call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), N_int)
|
||||||
|
call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), N_int)
|
||||||
|
! alpha - alpha
|
||||||
|
do i = 1, elec_num_tab_local(1)
|
||||||
|
iorb = occ(i,1)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb)
|
||||||
|
do j = i+1, elec_num_tab_local(1)
|
||||||
|
jorb = occ(j,1)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! beta - beta
|
||||||
|
do i = 1, elec_num_tab_local(2)
|
||||||
|
iorb = occ(i,2)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += mo_mono_elec_integral(iorb,iorb)
|
||||||
|
do j = i+1, elec_num_tab_local(2)
|
||||||
|
jorb = occ(j,2)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
! alpha - beta
|
||||||
|
do i = 1, elec_num_tab_local(2)
|
||||||
|
iorb = occ(i,2)
|
||||||
|
do j = 1, elec_num_tab_local(1)
|
||||||
|
jorb = occ(j,1)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
! alpha - core-act
|
||||||
|
do i = 1, elec_num_tab_local(1)
|
||||||
|
iorb = occ(i,1)
|
||||||
|
do j = 1, n_core_inact_orb
|
||||||
|
jorb = list_core_inact(j)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
core_act_exchange(1) += - mo_bielec_integral_jj_exchange(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! beta - core-act
|
||||||
|
do i = 1, elec_num_tab_local(2)
|
||||||
|
iorb = occ(i,2)
|
||||||
|
do j = 1, n_core_inact_orb
|
||||||
|
jorb = list_core_inact(j)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_exchange += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
core_act_exchange(2) += - mo_bielec_integral_jj_exchange(jorb,iorb)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine u0_H_dyall_u0_no_exchange(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: N_states_in,ndet,dim_psi_in,dim_psi_coef,state_target
|
||||||
|
integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in)
|
||||||
|
double precision, intent(in) :: psi_in_coef(dim_psi_coef,N_states_in)
|
||||||
|
double precision, intent(out) :: energies(N_states_in)
|
||||||
|
|
||||||
|
integer :: i,j
|
||||||
|
double precision :: hij,accu
|
||||||
|
energies = 0.d0
|
||||||
|
accu = 0.d0
|
||||||
|
double precision, allocatable :: psi_coef_tmp(:)
|
||||||
|
allocate(psi_coef_tmp(ndet))
|
||||||
|
|
||||||
|
do i = 1, ndet
|
||||||
|
psi_coef_tmp(i) = psi_in_coef(i,state_target)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
double precision :: hij_bis
|
||||||
|
do i = 1, ndet
|
||||||
|
if(psi_coef_tmp(i)==0.d0)cycle
|
||||||
|
do j = 1, ndet
|
||||||
|
if(psi_coef_tmp(j)==0.d0)cycle
|
||||||
|
call i_H_j_dyall_no_exchange(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
|
||||||
|
accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
energies(state_target) = accu
|
||||||
|
deallocate(psi_coef_tmp)
|
||||||
|
end
|
||||||
|
@ -79,8 +79,12 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip
|
|||||||
phase_array =0.d0
|
phase_array =0.d0
|
||||||
do i = 1,idx_alpha(0)
|
do i = 1,idx_alpha(0)
|
||||||
index_i = idx_alpha(i)
|
index_i = idx_alpha(i)
|
||||||
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e)
|
|
||||||
call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha)
|
call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha)
|
||||||
|
double precision :: coef_array(N_states)
|
||||||
|
do i_state = 1, N_states
|
||||||
|
coef_array(i_state) = psi_coef(index_i,i_state)
|
||||||
|
enddo
|
||||||
|
call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e)
|
||||||
hij_array(index_i) = hialpha
|
hij_array(index_i) = hialpha
|
||||||
call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
|
call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int)
|
||||||
! phase_array(index_i) = phase
|
! phase_array(index_i) = phase
|
||||||
|
@ -57,6 +57,9 @@
|
|||||||
! 1h1p
|
! 1h1p
|
||||||
delta_ij_tmp = 0.d0
|
delta_ij_tmp = 0.d0
|
||||||
call H_apply_mrpt_1h1p(delta_ij_tmp,N_det)
|
call H_apply_mrpt_1h1p(delta_ij_tmp,N_det)
|
||||||
|
double precision :: e_corr_from_1h1p_singles(N_states)
|
||||||
|
!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles)
|
||||||
|
!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp)
|
||||||
accu = 0.d0
|
accu = 0.d0
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
do i = 1, N_det
|
do i = 1, N_det
|
||||||
@ -69,6 +72,23 @@
|
|||||||
enddo
|
enddo
|
||||||
print*, '1h1p = ',accu
|
print*, '1h1p = ',accu
|
||||||
|
|
||||||
|
! 1h1p third order
|
||||||
|
delta_ij_tmp = 0.d0
|
||||||
|
call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
|
||||||
|
!call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles)
|
||||||
|
!call give_1h1p_only_doubles_spin_cross(delta_ij_tmp)
|
||||||
|
accu = 0.d0
|
||||||
|
do i_state = 1, N_states
|
||||||
|
do i = 1, N_det
|
||||||
|
do j = 1, N_det
|
||||||
|
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
|
||||||
|
delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
second_order_pt_new_1h1p(i_state) = accu(i_state)
|
||||||
|
enddo
|
||||||
|
print*, '1h1p(3)',accu
|
||||||
|
|
||||||
! 2h
|
! 2h
|
||||||
delta_ij_tmp = 0.d0
|
delta_ij_tmp = 0.d0
|
||||||
call H_apply_mrpt_2h(delta_ij_tmp,N_det)
|
call H_apply_mrpt_2h(delta_ij_tmp,N_det)
|
||||||
@ -101,6 +121,7 @@
|
|||||||
|
|
||||||
! 1h2p
|
! 1h2p
|
||||||
delta_ij_tmp = 0.d0
|
delta_ij_tmp = 0.d0
|
||||||
|
!call give_1h2p_contrib(delta_ij_tmp)
|
||||||
call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
|
call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
|
||||||
accu = 0.d0
|
accu = 0.d0
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
@ -116,6 +137,7 @@
|
|||||||
|
|
||||||
! 2h1p
|
! 2h1p
|
||||||
delta_ij_tmp = 0.d0
|
delta_ij_tmp = 0.d0
|
||||||
|
!call give_2h1p_contrib(delta_ij_tmp)
|
||||||
call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
|
call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
|
||||||
accu = 0.d0
|
accu = 0.d0
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
@ -159,7 +181,7 @@
|
|||||||
accu = 0.d0
|
accu = 0.d0
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
do i = 1, N_det
|
do i = 1, N_det
|
||||||
write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
|
! write(*,'(1000(F16.10,x))')delta_ij(i,:,:)
|
||||||
do j = i_state, N_det
|
do j = i_state, N_det
|
||||||
accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
|
accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
|
||||||
enddo
|
enddo
|
||||||
|
@ -391,3 +391,568 @@ subroutine give_1h2p_contrib(matrix_1h2p)
|
|||||||
|
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine give_1h1p_contrib(matrix_1h1p)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*)
|
||||||
|
integer :: i,j,r,a,b
|
||||||
|
integer :: iorb, jorb, rorb, aorb, borb
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: idet,jdet
|
||||||
|
integer :: inint
|
||||||
|
integer :: elec_num_tab_local(2),acu_elec
|
||||||
|
integer(bit_kind) :: det_tmp(N_int,2)
|
||||||
|
integer :: exc(0:2,2,2)
|
||||||
|
integer :: accu_elec
|
||||||
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
|
double precision :: active_int(n_act_orb,2)
|
||||||
|
double precision :: hij,phase
|
||||||
|
integer :: degree(N_det)
|
||||||
|
integer :: idx(0:N_det)
|
||||||
|
integer :: istate
|
||||||
|
double precision :: hja,delta_e_inact_virt(N_states)
|
||||||
|
integer :: kspin,degree_scalar
|
||||||
|
!matrix_1h1p = 0.d0
|
||||||
|
|
||||||
|
elec_num_tab_local = 0
|
||||||
|
do inint = 1, N_int
|
||||||
|
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
|
||||||
|
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
|
||||||
|
enddo
|
||||||
|
do i = 1, n_inact_orb ! First inactive
|
||||||
|
iorb = list_inact(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
do j = 1, N_states
|
||||||
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(iorb,j) &
|
||||||
|
- fock_virt_total_spin_trace(rorb,j)
|
||||||
|
enddo
|
||||||
|
do idet = 1, N_det
|
||||||
|
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations
|
||||||
|
do jdet = 1, idx(0)
|
||||||
|
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = psi_det(inint,1,idet)
|
||||||
|
det_tmp(inint,2) = psi_det(inint,2,idet)
|
||||||
|
enddo
|
||||||
|
! Do the excitation inactive -- > virtual
|
||||||
|
double precision :: himono,delta_e(N_states),coef_mono(N_states)
|
||||||
|
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
|
||||||
|
call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin
|
||||||
|
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono)
|
||||||
|
|
||||||
|
do state_target = 1, N_states
|
||||||
|
! delta_e(state_target) = one_anhil_one_creat_inact_virt(i,r,state_target) + delta_e_inact_virt(state_target)
|
||||||
|
delta_e(state_target) = one_anhil_one_creat_inact_virt_bis(i,r,idet,state_target)
|
||||||
|
coef_mono(state_target) = himono / delta_e(state_target)
|
||||||
|
enddo
|
||||||
|
if(idx(jdet).ne.idet)then
|
||||||
|
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
|
||||||
|
if (exc(0,1,1) == 1) then
|
||||||
|
! Mono alpha
|
||||||
|
aorb = (exc(1,2,1)) !!! a^{\dagger}_a
|
||||||
|
borb = (exc(1,1,1)) !!! a_{b}
|
||||||
|
jspin = 1
|
||||||
|
else
|
||||||
|
! Mono beta
|
||||||
|
aorb = (exc(1,2,2)) !!! a^{\dagger}_a
|
||||||
|
borb = (exc(1,1,2)) !!! a_{b}
|
||||||
|
jspin = 2
|
||||||
|
endif
|
||||||
|
|
||||||
|
call get_excitation_degree(psi_det(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
|
||||||
|
if(degree_scalar .ne. 2)then
|
||||||
|
print*, 'pb !!!'
|
||||||
|
print*, degree_scalar
|
||||||
|
call debug_det(psi_det(1,1,idx(jdet)),N_int)
|
||||||
|
call debug_det(det_tmp,N_int)
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
|
||||||
|
if(ispin == jspin )then
|
||||||
|
hij = -get_mo_bielec_integral_schwartz(iorb,aorb,rorb,borb,mo_integrals_map) &
|
||||||
|
+ get_mo_bielec_integral_schwartz(iorb,aorb,borb,rorb,mo_integrals_map)
|
||||||
|
else
|
||||||
|
hij = get_mo_bielec_integral_schwartz(iorb,borb,rorb,aorb,mo_integrals_map)
|
||||||
|
endif
|
||||||
|
hij = hij * phase
|
||||||
|
double precision :: hij_test
|
||||||
|
integer :: state_target
|
||||||
|
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test)
|
||||||
|
if(dabs(hij - hij_test).gt.1.d-10)then
|
||||||
|
print*, 'ahah pb !!'
|
||||||
|
print*, 'hij .ne. hij_test'
|
||||||
|
print*, hij,hij_test
|
||||||
|
call debug_det(psi_det(1,1,idx(jdet)),N_int)
|
||||||
|
call debug_det(det_tmp,N_int)
|
||||||
|
print*, ispin, jspin
|
||||||
|
print*,iorb,borb,rorb,aorb
|
||||||
|
print*, phase
|
||||||
|
call i_H_j_verbose(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test)
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
do state_target = 1, N_states
|
||||||
|
matrix_1h1p(idx(jdet),idet,state_target) += hij* coef_mono(state_target)
|
||||||
|
enddo
|
||||||
|
else
|
||||||
|
do state_target = 1, N_states
|
||||||
|
matrix_1h1p(idet,idet,state_target) += himono * coef_mono(state_target)
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*)
|
||||||
|
integer :: i,j,r,a,b
|
||||||
|
integer :: iorb, jorb, rorb, aorb, borb,s,sorb
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: idet,jdet
|
||||||
|
integer :: inint
|
||||||
|
integer :: elec_num_tab_local(2),acu_elec
|
||||||
|
integer(bit_kind) :: det_tmp(N_int,2),det_tmp_bis(N_int,2)
|
||||||
|
integer(bit_kind) :: det_pert(N_int,2,n_inact_orb,n_virt_orb,2)
|
||||||
|
double precision :: coef_det_pert(n_inact_orb,n_virt_orb,2,N_states,2)
|
||||||
|
double precision :: delta_e_det_pert(n_inact_orb,n_virt_orb,2,N_states)
|
||||||
|
double precision :: hij_det_pert(n_inact_orb,n_virt_orb,2,N_states)
|
||||||
|
integer :: exc(0:2,2,2)
|
||||||
|
integer :: accu_elec
|
||||||
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
|
double precision :: active_int(n_act_orb,2)
|
||||||
|
double precision :: hij,phase
|
||||||
|
integer :: degree(N_det)
|
||||||
|
integer :: idx(0:N_det)
|
||||||
|
integer :: istate
|
||||||
|
double precision :: hja,delta_e_inact_virt(N_states)
|
||||||
|
integer :: kspin,degree_scalar
|
||||||
|
!matrix_1h1p = 0.d0
|
||||||
|
|
||||||
|
elec_num_tab_local = 0
|
||||||
|
do inint = 1, N_int
|
||||||
|
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
|
||||||
|
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
|
||||||
|
enddo
|
||||||
|
double precision :: himono,delta_e(N_states),coef_mono(N_states)
|
||||||
|
integer :: state_target
|
||||||
|
do idet = 1, N_det
|
||||||
|
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
|
||||||
|
do i = 1, n_inact_orb ! First inactive
|
||||||
|
iorb = list_inact(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
coef_det_pert(i,r,ispin,state_target,1) = 0.d0
|
||||||
|
coef_det_pert(i,r,ispin,state_target,2) = 0.d0
|
||||||
|
enddo
|
||||||
|
do j = 1, N_states
|
||||||
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(iorb,j) &
|
||||||
|
- fock_virt_total_spin_trace(rorb,j)
|
||||||
|
enddo
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = psi_det(inint,1,idet)
|
||||||
|
det_tmp(inint,2) = psi_det(inint,2,idet)
|
||||||
|
enddo
|
||||||
|
! Do the excitation inactive -- > virtual
|
||||||
|
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
|
||||||
|
call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin
|
||||||
|
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono)
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_pert(inint,1,i,r,ispin) = det_tmp(inint,1)
|
||||||
|
det_pert(inint,2,i,r,ispin) = det_tmp(inint,2)
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
delta_e_det_pert(i,r,ispin,state_target) = one_anhil_one_creat_inact_virt(i,r,state_target) + delta_e_inact_virt(state_target)
|
||||||
|
coef_det_pert(i,r,ispin,state_target,1) = himono / delta_e_det_pert(i,r,ispin,state_target)
|
||||||
|
enddo
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations
|
||||||
|
enddo ! ispin
|
||||||
|
enddo ! rorb
|
||||||
|
enddo ! iorb
|
||||||
|
|
||||||
|
do i = 1, n_inact_orb ! First inactive
|
||||||
|
iorb = list_inact(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = det_pert(inint,1,i,r,ispin)
|
||||||
|
det_tmp(inint,2) = det_pert(inint,2,i,r,ispin)
|
||||||
|
enddo
|
||||||
|
do j = 1, n_inact_orb ! First inactive
|
||||||
|
jorb = list_inact(j)
|
||||||
|
do s = 1, n_virt_orb ! First virtual
|
||||||
|
sorb = list_virt(s)
|
||||||
|
do jspin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
if(i==j.and.r==s.and.ispin==jspin)cycle
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp_bis(inint,1) = det_pert(inint,1,j,s,jspin)
|
||||||
|
det_tmp_bis(inint,2) = det_pert(inint,2,j,s,jspin)
|
||||||
|
enddo
|
||||||
|
call i_H_j(det_tmp_bis,det_tmp,N_int,himono)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
coef_det_pert(i,r,ispin,state_target,2) += &
|
||||||
|
coef_det_pert(j,s,jspin,state_target,1) * himono / delta_e_det_pert(i,r,ispin,state_target)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo ! ispin
|
||||||
|
enddo ! rorb
|
||||||
|
enddo ! iorb
|
||||||
|
do i = 1, n_inact_orb ! First inactive
|
||||||
|
iorb = list_inact(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
coef_det_pert(i,r,ispin,state_target,1) += coef_det_pert(i,r,ispin,state_target,2)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = det_pert(inint,1,i,r,ispin)
|
||||||
|
det_tmp(inint,2) = det_pert(inint,2,i,r,ispin)
|
||||||
|
enddo
|
||||||
|
do jdet = 1, idx(0)
|
||||||
|
!
|
||||||
|
if(idx(jdet).ne.idet)then
|
||||||
|
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
|
||||||
|
if (exc(0,1,1) == 1) then
|
||||||
|
! Mono alpha
|
||||||
|
aorb = (exc(1,2,1)) !!! a^{\dagger}_a
|
||||||
|
borb = (exc(1,1,1)) !!! a_{b}
|
||||||
|
jspin = 1
|
||||||
|
else
|
||||||
|
aorb = (exc(1,2,2)) !!! a^{\dagger}_a
|
||||||
|
borb = (exc(1,1,2)) !!! a_{b}
|
||||||
|
jspin = 2
|
||||||
|
endif
|
||||||
|
|
||||||
|
call get_excitation_degree(psi_det(1,1,idx(jdet)),det_tmp,degree_scalar,N_int)
|
||||||
|
if(degree_scalar .ne. 2)then
|
||||||
|
print*, 'pb !!!'
|
||||||
|
print*, degree_scalar
|
||||||
|
call debug_det(psi_det(1,1,idx(jdet)),N_int)
|
||||||
|
call debug_det(det_tmp,N_int)
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
|
||||||
|
double precision :: hij_test
|
||||||
|
hij_test = 0.d0
|
||||||
|
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij_test)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
matrix_1h1p(idx(jdet),idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2)
|
||||||
|
enddo
|
||||||
|
else
|
||||||
|
hij_test = 0.d0
|
||||||
|
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hij_test)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
matrix_1h1p(idet,idet,state_target) += hij_test* coef_det_pert(i,r,ispin,state_target,2)
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
enddo ! idet
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine give_1p_sec_order_singles_contrib(matrix_1p)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
double precision , intent(inout) :: matrix_1p(N_det,N_det,*)
|
||||||
|
integer :: i,j,r,a,b
|
||||||
|
integer :: iorb, jorb, rorb, aorb, borb,s,sorb
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: idet,jdet
|
||||||
|
integer :: inint
|
||||||
|
integer :: elec_num_tab_local(2),acu_elec
|
||||||
|
integer(bit_kind) :: det_tmp(N_int,2),det_tmp_bis(N_int,2)
|
||||||
|
integer(bit_kind) :: det_pert(N_int,2,n_act_orb,n_virt_orb,2)
|
||||||
|
double precision :: coef_det_pert(n_act_orb,n_virt_orb,2,N_states,2)
|
||||||
|
double precision :: delta_e_det_pert(n_act_orb,n_virt_orb,2,N_states)
|
||||||
|
double precision :: hij_det_pert(n_act_orb,n_virt_orb,2)
|
||||||
|
integer :: exc(0:2,2,2)
|
||||||
|
integer :: accu_elec
|
||||||
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
|
double precision :: hij,phase
|
||||||
|
integer :: degree(N_det)
|
||||||
|
integer :: idx(0:N_det)
|
||||||
|
integer :: istate
|
||||||
|
double precision :: hja,delta_e_act_virt(N_states)
|
||||||
|
integer :: kspin,degree_scalar
|
||||||
|
!matrix_1p = 0.d0
|
||||||
|
|
||||||
|
elec_num_tab_local = 0
|
||||||
|
do inint = 1, N_int
|
||||||
|
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
|
||||||
|
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
|
||||||
|
enddo
|
||||||
|
double precision :: himono,delta_e(N_states),coef_mono(N_states)
|
||||||
|
integer :: state_target
|
||||||
|
do idet = 1, N_det
|
||||||
|
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
|
||||||
|
do i = 1, n_act_orb ! First active
|
||||||
|
iorb = list_act(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
coef_det_pert(i,r,ispin,state_target,1) = 0.d0
|
||||||
|
coef_det_pert(i,r,ispin,state_target,2) = 0.d0
|
||||||
|
enddo
|
||||||
|
do j = 1, N_states
|
||||||
|
delta_e_act_virt(j) = - fock_virt_total_spin_trace(rorb,j)
|
||||||
|
enddo
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = psi_det(inint,1,idet)
|
||||||
|
det_tmp(inint,2) = psi_det(inint,2,idet)
|
||||||
|
enddo
|
||||||
|
! Do the excitation active -- > virtual
|
||||||
|
call do_mono_excitation(det_tmp,iorb,rorb,ispin,i_ok)
|
||||||
|
integer :: i_ok
|
||||||
|
if(i_ok .ne.1)then
|
||||||
|
do state_target = 1, N_states
|
||||||
|
coef_det_pert(i,r,ispin,state_target,1) = -1.d+10
|
||||||
|
coef_det_pert(i,r,ispin,state_target,2) = -1.d+10
|
||||||
|
hij_det_pert(i,r,ispin) = 0.d0
|
||||||
|
enddo
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_pert(inint,1,i,r,ispin) = 0_bit_kind
|
||||||
|
det_pert(inint,2,i,r,ispin) = 0_bit_kind
|
||||||
|
enddo
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,himono)
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_pert(inint,1,i,r,ispin) = det_tmp(inint,1)
|
||||||
|
det_pert(inint,2,i,r,ispin) = det_tmp(inint,2)
|
||||||
|
enddo
|
||||||
|
do state_target = 1, N_states
|
||||||
|
delta_e_det_pert(i,r,ispin,state_target) = one_creat_virt(i,r,state_target) + delta_e_act_virt(state_target)
|
||||||
|
coef_det_pert(i,r,ispin,state_target,1) = himono / delta_e_det_pert(i,r,ispin,state_target)
|
||||||
|
hij_det_pert(i,r,ispin) = himono
|
||||||
|
enddo
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations
|
||||||
|
enddo ! ispin
|
||||||
|
enddo ! rorb
|
||||||
|
enddo ! iorb
|
||||||
|
|
||||||
|
! do i = 1, n_act_orb ! First active
|
||||||
|
! do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
! if(coef_det_pert(i,1,ispin,1,1) == -1.d+10)cycle
|
||||||
|
! iorb = list_act(i)
|
||||||
|
! do r = 1, n_virt_orb ! First virtual
|
||||||
|
! rorb = list_virt(r)
|
||||||
|
! do inint = 1, N_int
|
||||||
|
! det_tmp(inint,1) = det_pert(inint,1,i,r,ispin)
|
||||||
|
! det_tmp(inint,2) = det_pert(inint,2,i,r,ispin)
|
||||||
|
! enddo
|
||||||
|
! do j = 1, n_act_orb ! First active
|
||||||
|
! do jspin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
! if(coef_det_pert(j,1,jspin,1,1) == -1.d+10)cycle
|
||||||
|
! jorb = list_act(j)
|
||||||
|
! do s = 1, n_virt_orb ! First virtual
|
||||||
|
! sorb = list_virt(s)
|
||||||
|
! if(i==j.and.r==s.and.ispin==jspin)cycle
|
||||||
|
! do inint = 1, N_int
|
||||||
|
! det_tmp_bis(inint,1) = det_pert(inint,1,j,s,jspin)
|
||||||
|
! det_tmp_bis(inint,2) = det_pert(inint,2,j,s,jspin)
|
||||||
|
! enddo
|
||||||
|
! call i_H_j(det_tmp_bis,det_tmp,N_int,himono)
|
||||||
|
! do state_target = 1, N_states
|
||||||
|
! coef_det_pert(i,r,ispin,state_target,2) += &
|
||||||
|
! coef_det_pert(j,s,jspin,state_target,1) * himono / delta_e_det_pert(i,r,ispin,state_target)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo ! ispin
|
||||||
|
! enddo ! rorb
|
||||||
|
! enddo ! iorb
|
||||||
|
|
||||||
|
do i = 1, n_act_orb ! First active
|
||||||
|
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
|
||||||
|
if(coef_det_pert(i,1,ispin,1,1) == -1.d+10)cycle
|
||||||
|
iorb = list_act(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
! do state_target = 1, N_states
|
||||||
|
! coef_det_pert(i,r,ispin,state_target,1) += coef_det_pert(i,r,ispin,state_target,2)
|
||||||
|
! enddo
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = det_pert(inint,1,i,r,ispin)
|
||||||
|
det_tmp(inint,2) = det_pert(inint,2,i,r,ispin)
|
||||||
|
enddo
|
||||||
|
do jdet = 1,N_det
|
||||||
|
double precision :: coef_array(N_states),hij_test
|
||||||
|
call i_H_j(det_tmp,psi_det(1,1,jdet),N_int,himono)
|
||||||
|
call get_delta_e_dyall(psi_det(1,1,jdet),det_tmp,coef_array,hij_test,delta_e)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
! matrix_1p(idet,jdet,state_target) += himono * coef_det_pert(i,r,ispin,state_target,1)
|
||||||
|
matrix_1p(idet,jdet,state_target) += himono * hij_det_pert(i,r,ispin) / delta_e(state_target)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
enddo ! idet
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
double precision , intent(inout) :: matrix_1h1p(N_det,N_det,*)
|
||||||
|
integer :: i,j,r,a,b
|
||||||
|
integer :: iorb, jorb, rorb, aorb, borb
|
||||||
|
integer :: ispin,jspin
|
||||||
|
integer :: idet,jdet
|
||||||
|
integer :: inint
|
||||||
|
integer :: elec_num_tab_local(2),acu_elec
|
||||||
|
integer(bit_kind) :: det_tmp(N_int,2)
|
||||||
|
integer :: exc(0:2,2,2)
|
||||||
|
integer :: accu_elec
|
||||||
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
|
double precision :: active_int(n_act_orb,2)
|
||||||
|
double precision :: hij,phase
|
||||||
|
integer :: degree(N_det)
|
||||||
|
integer :: idx(0:N_det)
|
||||||
|
integer :: istate
|
||||||
|
double precision :: hja,delta_e_inact_virt(N_states)
|
||||||
|
integer(bit_kind) :: pert_det(N_int,2,n_act_orb,n_act_orb,2)
|
||||||
|
double precision :: pert_det_coef(n_act_orb,n_act_orb,2,N_states)
|
||||||
|
integer :: kspin,degree_scalar
|
||||||
|
integer :: other_spin(2)
|
||||||
|
other_spin(1) = 2
|
||||||
|
other_spin(2) = 1
|
||||||
|
double precision :: hidouble,delta_e(N_states)
|
||||||
|
!matrix_1h1p = 0.d0
|
||||||
|
|
||||||
|
elec_num_tab_local = 0
|
||||||
|
do inint = 1, N_int
|
||||||
|
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
|
||||||
|
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
|
||||||
|
enddo
|
||||||
|
do i = 1, n_inact_orb ! First inactive
|
||||||
|
iorb = list_inact(i)
|
||||||
|
do r = 1, n_virt_orb ! First virtual
|
||||||
|
rorb = list_virt(r)
|
||||||
|
do j = 1, N_states
|
||||||
|
delta_e_inact_virt(j) = fock_core_inactive_total_spin_trace(iorb,j) &
|
||||||
|
- fock_virt_total_spin_trace(rorb,j)
|
||||||
|
enddo
|
||||||
|
do idet = 1, N_det
|
||||||
|
call get_excitation_degree_vector_double_alpha_beta(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of the mono excitations
|
||||||
|
do ispin = 1, 2
|
||||||
|
jspin = other_spin(ispin)
|
||||||
|
do a = 1, n_act_orb
|
||||||
|
aorb = list_act(a)
|
||||||
|
do b = 1, n_act_orb
|
||||||
|
borb = list_act(b)
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = psi_det(inint,1,idet)
|
||||||
|
det_tmp(inint,2) = psi_det(inint,2,idet)
|
||||||
|
enddo
|
||||||
|
! Do the excitation (i-->a)(ispin) + (b-->r)(other_spin(ispin))
|
||||||
|
integer :: i_ok,corb,dorb
|
||||||
|
call do_mono_excitation(det_tmp,iorb,aorb,ispin,i_ok)
|
||||||
|
if(i_ok .ne. 1)then
|
||||||
|
do state_target = 1, N_states
|
||||||
|
pert_det_coef(a,b,ispin,state_target) = -100000.d0
|
||||||
|
enddo
|
||||||
|
do inint = 1, N_int
|
||||||
|
pert_det(inint,1,a,b,ispin) = 0_bit_kind
|
||||||
|
pert_det(inint,2,a,b,ispin) = 0_bit_kind
|
||||||
|
enddo
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
call do_mono_excitation(det_tmp,borb,rorb,jspin,i_ok)
|
||||||
|
if(i_ok .ne. 1)then
|
||||||
|
do state_target = 1, N_states
|
||||||
|
pert_det_coef(a,b,ispin,state_target) = -100000.d0
|
||||||
|
enddo
|
||||||
|
do inint = 1, N_int
|
||||||
|
pert_det(inint,1,a,b,ispin) = 0_bit_kind
|
||||||
|
pert_det(inint,2,a,b,ispin) = 0_bit_kind
|
||||||
|
enddo
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
|
do inint = 1, N_int
|
||||||
|
pert_det(inint,1,a,b,ispin) = det_tmp(inint,1)
|
||||||
|
pert_det(inint,2,a,b,ispin) = det_tmp(inint,2)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hidouble)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
delta_e(state_target) = one_anhil_one_creat(a,b,ispin,jspin,state_target) + delta_e_inact_virt(state_target)
|
||||||
|
pert_det_coef(a,b,ispin,state_target) = hidouble / delta_e(state_target)
|
||||||
|
matrix_1h1p(idet,idet,state_target) += hidouble * pert_det_coef(a,b,ispin,state_target)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do jdet = 1, idx(0)
|
||||||
|
if(idx(jdet).ne.idet)then
|
||||||
|
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
|
||||||
|
integer :: c,d,state_target
|
||||||
|
integer(bit_kind) :: det_tmp_bis(N_int,2)
|
||||||
|
! excitation from I --> J
|
||||||
|
! (a->c) (alpha) + (b->d) (beta)
|
||||||
|
aorb = exc(1,1,1)
|
||||||
|
corb = exc(1,2,1)
|
||||||
|
c = list_act_reverse(corb)
|
||||||
|
borb = exc(1,1,2)
|
||||||
|
dorb = exc(1,2,2)
|
||||||
|
d = list_act_reverse(dorb)
|
||||||
|
ispin = 1
|
||||||
|
jspin = 2
|
||||||
|
do inint = 1, N_int
|
||||||
|
det_tmp(inint,1) = pert_det(inint,1,c,d,1)
|
||||||
|
det_tmp(inint,2) = pert_det(inint,2,c,d,1)
|
||||||
|
det_tmp_bis(inint,1) = pert_det(inint,1,c,d,2)
|
||||||
|
det_tmp_bis(inint,2) = pert_det(inint,2,c,d,2)
|
||||||
|
enddo
|
||||||
|
double precision :: hjdouble_1,hjdouble_2
|
||||||
|
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hjdouble_1)
|
||||||
|
call i_H_j(psi_det(1,1,idx(jdet)),det_tmp_bis,N_int,hjdouble_2)
|
||||||
|
do state_target = 1, N_states
|
||||||
|
matrix_1h1p(idx(jdet),idet,state_target) += (pert_det_coef(c,d,1,state_target) * hjdouble_1 + pert_det_coef(c,d,2,state_target) * hjdouble_2 )
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
@ -14,8 +14,6 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)]
|
|||||||
psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1))
|
psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1))
|
||||||
psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1))
|
psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
! call debug_det(psi_active(1,1,i),N_int)
|
|
||||||
enddo
|
enddo
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -154,7 +152,7 @@ subroutine give_particles_in_virt_space(det_1,n_particles_spin,n_particles,parti
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! routine that returns the delta_e with the Moller Plesset and Dyall operators
|
! routine that returns the delta_e with the Moller Plesset and Dyall operators
|
||||||
!
|
!
|
||||||
@ -172,6 +170,7 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
|||||||
use bitmasks
|
use bitmasks
|
||||||
double precision, intent(out) :: delta_e_final(N_states)
|
double precision, intent(out) :: delta_e_final(N_states)
|
||||||
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
|
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
|
||||||
|
double precision, intent(in) :: coef_array(N_states),hij
|
||||||
integer :: i,j,k,l
|
integer :: i,j,k,l
|
||||||
integer :: i_state
|
integer :: i_state
|
||||||
|
|
||||||
@ -292,23 +291,52 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
|||||||
! delta_e_act += one_creat_spin_trace(i_particle_act )
|
! delta_e_act += one_creat_spin_trace(i_particle_act )
|
||||||
ispin = particle_list_practical(1,1)
|
ispin = particle_list_practical(1,1)
|
||||||
i_particle_act = particle_list_practical(2,1)
|
i_particle_act = particle_list_practical(2,1)
|
||||||
|
call get_excitation_degree(det_1,det_2,degree,N_int)
|
||||||
|
if(degree == 1)then
|
||||||
|
call get_excitation(det_1,det_2,exc,degree,phase,N_int)
|
||||||
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||||
|
i_hole = list_inact_reverse(h1)
|
||||||
|
i_part = list_act_reverse(p1)
|
||||||
|
do i_state = 1, N_states
|
||||||
|
delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state)
|
||||||
|
enddo
|
||||||
|
else if (degree == 2)then
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state)
|
delta_e_act(i_state) += one_creat(i_particle_act,ispin,i_state)
|
||||||
enddo
|
enddo
|
||||||
|
endif
|
||||||
|
|
||||||
else if (n_holes_act == 1 .and. n_particles_act == 0) then
|
else if (n_holes_act == 1 .and. n_particles_act == 0) then
|
||||||
! i_hole_act = holes_active_list_spin_traced(1)
|
! i_hole_act = holes_active_list_spin_traced(1)
|
||||||
! delta_e_act += one_anhil_spin_trace(i_hole_act )
|
! delta_e_act += one_anhil_spin_trace(i_hole_act )
|
||||||
ispin = hole_list_practical(1,1)
|
ispin = hole_list_practical(1,1)
|
||||||
i_hole_act = hole_list_practical(2,1)
|
i_hole_act = hole_list_practical(2,1)
|
||||||
|
call get_excitation_degree(det_1,det_2,degree,N_int)
|
||||||
|
if(degree == 1)then
|
||||||
|
call get_excitation(det_1,det_2,exc,degree,phase,N_int)
|
||||||
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||||
|
i_hole = list_act_reverse(h1)
|
||||||
|
i_part = list_virt_reverse(p1)
|
||||||
|
do i_state = 1, N_states
|
||||||
|
if(isnan(one_creat_virt(i_hole,i_part,i_state)))then
|
||||||
|
print*, i_hole,i_part,i_state
|
||||||
|
call debug_det(det_1,N_int)
|
||||||
|
call debug_det(det_2,N_int)
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
delta_e_act(i_state) += one_creat_virt(i_hole,i_part,i_state)
|
||||||
|
enddo
|
||||||
|
else if (degree == 2)then
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state)
|
delta_e_act(i_state) += one_anhil(i_hole_act , ispin,i_state)
|
||||||
enddo
|
enddo
|
||||||
|
endif
|
||||||
|
|
||||||
else if (n_holes_act == 1 .and. n_particles_act == 1) then
|
else if (n_holes_act == 1 .and. n_particles_act == 1) then
|
||||||
! i_hole_act = holes_active_list_spin_traced(1)
|
! i_hole_act = holes_active_list_spin_traced(1)
|
||||||
! i_particle_act = particles_active_list_spin_traced(1)
|
! i_particle_act = particles_active_list_spin_traced(1)
|
||||||
! delta_e_act += one_anhil_one_creat_spin_trace(i_hole_act,i_particle_act)
|
! delta_e_act += one_anhil_one_creat_spin_trace(i_hole_act,i_particle_act)
|
||||||
|
|
||||||
! first hole
|
! first hole
|
||||||
ispin = hole_list_practical(1,1)
|
ispin = hole_list_practical(1,1)
|
||||||
i_hole_act = hole_list_practical(2,1)
|
i_hole_act = hole_list_practical(2,1)
|
||||||
@ -423,6 +451,34 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
|||||||
delta_e_act(i_state) += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin,i_state)
|
delta_e_act(i_state) += three_creat(i_particle_act,j_particle_act,k_particle_act,ispin,jspin,kspin,i_state)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
else if (n_holes_act .eq. 0 .and. n_particles_act .eq.0)then
|
||||||
|
integer :: degree
|
||||||
|
integer(bit_kind) :: det_1_active(N_int,2)
|
||||||
|
integer :: h1,h2,p1,p2,s1,s2
|
||||||
|
integer :: exc(0:2,2,2)
|
||||||
|
integer :: i_hole, i_part
|
||||||
|
double precision :: phase
|
||||||
|
call get_excitation_degree(det_1,det_2,degree,N_int)
|
||||||
|
if(degree == 1)then
|
||||||
|
! call debug_det(det_1,N_int)
|
||||||
|
call get_excitation(det_1,det_2,exc,degree,phase,N_int)
|
||||||
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||||
|
i_hole = list_inact_reverse(h1)
|
||||||
|
i_part = list_virt_reverse(p1)
|
||||||
|
do i_state = 1, N_states
|
||||||
|
! if(one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1).gt.1.d-10)then
|
||||||
|
! print*, hij, one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1)
|
||||||
|
! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state) &
|
||||||
|
! * coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1)
|
||||||
|
! print*, coef_array(i_state)* hij*coef_array(i_state)* hij,one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1), &
|
||||||
|
! coef_array(i_state)* hij*coef_array(i_state)* hij *one_anhil_one_creat_inact_virt_norm(i_hole,i_part,i_state,s1)
|
||||||
|
! else
|
||||||
|
delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state)
|
||||||
|
! endif
|
||||||
|
enddo
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then
|
else if (n_holes_act .ge. 2 .and. n_particles_act .ge.2) then
|
||||||
|
|
||||||
delta_e_act = -10000000.d0
|
delta_e_act = -10000000.d0
|
||||||
@ -438,3 +494,4 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final)
|
|||||||
enddo
|
enddo
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -192,6 +192,11 @@ subroutine pt2_moller_plesset ($arguments)
|
|||||||
endif
|
endif
|
||||||
do i =1,N_st
|
do i =1,N_st
|
||||||
H_pert_diag(i) = h
|
H_pert_diag(i) = h
|
||||||
|
! if(dabs(i_H_psi_array(i)).gt.1.d-8)then
|
||||||
|
! print*, i_H_psi_array(i)
|
||||||
|
! call debug_det(det_pert,N_int)
|
||||||
|
! print*, h1,p1,h2,p2,s1,s2
|
||||||
|
! endif
|
||||||
c_pert(i) = i_H_psi_array(i) *delta_e
|
c_pert(i) = i_H_psi_array(i) *delta_e
|
||||||
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
|
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
|
||||||
enddo
|
enddo
|
||||||
|
@ -102,6 +102,11 @@ END_PROVIDER
|
|||||||
conversion_factor_gauss_hcc(3) = 619.9027742370165d0
|
conversion_factor_gauss_hcc(3) = 619.9027742370165d0
|
||||||
conversion_factor_cm_1_hcc(3) = 579.4924475562677d0
|
conversion_factor_cm_1_hcc(3) = 579.4924475562677d0
|
||||||
|
|
||||||
|
! boron
|
||||||
|
conversion_factor_mhz_hcc(5) = 1434.3655101868d0
|
||||||
|
conversion_factor_gauss_hcc(5) = 511.817264334d0
|
||||||
|
conversion_factor_cm_1_hcc(5) = 478.4528336953d0
|
||||||
|
|
||||||
! carbon
|
! carbon
|
||||||
conversion_factor_mhz_hcc(6) = 1124.18303629792945d0
|
conversion_factor_mhz_hcc(6) = 1124.18303629792945d0
|
||||||
conversion_factor_gauss_hcc(6) = 401.136570647523058d0
|
conversion_factor_gauss_hcc(6) = 401.136570647523058d0
|
||||||
|
@ -98,7 +98,14 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id)
|
|||||||
if (N_det_selectors_read > 0) then
|
if (N_det_selectors_read > 0) then
|
||||||
N_det_selectors = N_det_selectors_read
|
N_det_selectors = N_det_selectors_read
|
||||||
endif
|
endif
|
||||||
SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators
|
SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators psi_coef_generators psi_det_generators
|
||||||
|
! n_det_generators
|
||||||
|
! n_det_selectors
|
||||||
|
! psi_coef
|
||||||
|
! psi_coef_generators
|
||||||
|
! psi_det
|
||||||
|
! psi_det_generators
|
||||||
|
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -541,3 +541,24 @@ use bitmasks
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
logical function is_i_in_virtual(i)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
integer,intent(in) :: i
|
||||||
|
integer(bit_kind) :: key(N_int)
|
||||||
|
integer :: k,j
|
||||||
|
integer :: accu
|
||||||
|
is_i_in_virtual = .False.
|
||||||
|
key= 0_bit_kind
|
||||||
|
k = ishft(i-1,-bit_kind_shift)+1
|
||||||
|
j = i-ishft(k-1,bit_kind_shift)-1
|
||||||
|
key(k) = ibset(key(k),j)
|
||||||
|
accu = 0
|
||||||
|
do k = 1, N_int
|
||||||
|
accu += popcnt(iand(key(k),virt_bitmask(k,1)))
|
||||||
|
enddo
|
||||||
|
if(accu .ne. 0)then
|
||||||
|
is_i_in_virtual = .True.
|
||||||
|
endif
|
||||||
|
|
||||||
|
end
|
||||||
|
@ -473,15 +473,31 @@ END_PROVIDER
|
|||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)]
|
BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)]
|
||||||
|
&BEGIN_PROVIDER [ integer, n_core_inact_act_orb ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Reunion of the core, inactive and active bitmasks
|
! Reunion of the core, inactive and active bitmasks
|
||||||
END_DOC
|
END_DOC
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
|
|
||||||
|
n_core_inact_act_orb = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1))
|
reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1))
|
||||||
reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,1,1))
|
reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,1,1))
|
||||||
|
n_core_inact_act_orb +=popcnt(reunion_of_core_inact_act_bitmask(i,1))
|
||||||
|
enddo
|
||||||
|
END_PROVIDER
|
||||||
|
BEGIN_PROVIDER [ integer, list_core_inact_act, (n_core_inact_act_orb)]
|
||||||
|
&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (mo_tot_num)]
|
||||||
|
implicit none
|
||||||
|
integer :: occ_inact(N_int*bit_kind_size)
|
||||||
|
integer :: itest,i
|
||||||
|
occ_inact = 0
|
||||||
|
call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), occ_inact(1), itest, N_int)
|
||||||
|
list_inact_reverse = 0
|
||||||
|
do i = 1, n_core_inact_act_orb
|
||||||
|
list_core_inact_act(i) = occ_inact(i)
|
||||||
|
list_core_inact_act_reverse(occ_inact(i)) = i
|
||||||
enddo
|
enddo
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -543,8 +559,8 @@ END_PROVIDER
|
|||||||
integer :: i,j
|
integer :: i,j
|
||||||
n_core_orb = 0
|
n_core_orb = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
core_bitmask(i,1) = xor(closed_shell_ref_bitmask(i,1),reunion_of_cas_inact_bitmask(i,1))
|
core_bitmask(i,1) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,1),virt_bitmask(i,1)))
|
||||||
core_bitmask(i,2) = xor(closed_shell_ref_bitmask(i,2),reunion_of_cas_inact_bitmask(i,2))
|
core_bitmask(i,2) = xor(full_ijkl_bitmask(i),ior(reunion_of_cas_inact_bitmask(i,2),virt_bitmask(i,1)))
|
||||||
n_core_orb += popcnt(core_bitmask(i,1))
|
n_core_orb += popcnt(core_bitmask(i,1))
|
||||||
enddo
|
enddo
|
||||||
print*,'n_core_orb = ',n_core_orb
|
print*,'n_core_orb = ',n_core_orb
|
||||||
|
@ -37,6 +37,72 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok)
|
|||||||
endif
|
endif
|
||||||
end
|
end
|
||||||
|
|
||||||
|
subroutine do_spin_flip(key_in,i_flip,ispin,i_ok)
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! flip the spin ispin in the orbital i_flip
|
||||||
|
! on key_in
|
||||||
|
! ispin = 1 == alpha
|
||||||
|
! ispin = 2 == beta
|
||||||
|
! i_ok = 1 == the flip is possible
|
||||||
|
! i_ok = -1 == the flip is not possible
|
||||||
|
END_DOC
|
||||||
|
integer, intent(in) :: i_flip,ispin
|
||||||
|
integer(bit_kind), intent(inout) :: key_in(N_int,2)
|
||||||
|
integer, intent(out) :: i_ok
|
||||||
|
integer :: k,j,i
|
||||||
|
integer(bit_kind) :: key_tmp(N_int,2)
|
||||||
|
i_ok = -1
|
||||||
|
key_tmp = 0_bit_kind
|
||||||
|
k = ishft(i_flip-1,-bit_kind_shift)+1
|
||||||
|
j = i_flip-ishft(k-1,bit_kind_shift)-1
|
||||||
|
key_tmp(k,1) = ibset(key_tmp(k,1),j)
|
||||||
|
integer :: other_spin(2)
|
||||||
|
other_spin(1) = 2
|
||||||
|
other_spin(2) = 1
|
||||||
|
if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then
|
||||||
|
! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip"
|
||||||
|
key_in(k,ispin) = ibclr(key_in(k,ispin),j) ! destroy the electron ispin in the orbital i_flip
|
||||||
|
key_in(k,other_spin(ispin)) = ibset(key_in(k,other_spin(ispin)),j) ! create an electron of spin other_spin in the same orbital
|
||||||
|
i_ok = 1
|
||||||
|
else
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
logical function is_spin_flip_possible(key_in,i_flip,ispin)
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! returns .True. if the spin-flip of spin ispin in the orbital i_flip is possible
|
||||||
|
! on key_in
|
||||||
|
END_DOC
|
||||||
|
integer, intent(in) :: i_flip,ispin
|
||||||
|
integer(bit_kind), intent(in) :: key_in(N_int,2)
|
||||||
|
integer :: k,j,i
|
||||||
|
integer(bit_kind) :: key_tmp(N_int,2)
|
||||||
|
is_spin_flip_possible = .False.
|
||||||
|
key_tmp = 0_bit_kind
|
||||||
|
k = ishft(i_flip-1,-bit_kind_shift)+1
|
||||||
|
j = i_flip-ishft(k-1,bit_kind_shift)-1
|
||||||
|
key_tmp(k,1) = ibset(key_tmp(k,1),j)
|
||||||
|
integer :: other_spin(2)
|
||||||
|
other_spin(1) = 2
|
||||||
|
other_spin(2) = 1
|
||||||
|
if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then
|
||||||
|
! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip"
|
||||||
|
is_spin_flip_possible = .True.
|
||||||
|
return
|
||||||
|
else
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
subroutine set_bit_to_integer(i_physical,key,Nint)
|
subroutine set_bit_to_integer(i_physical,key,Nint)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
|
@ -334,7 +334,7 @@ subroutine diagonalize_s2_betweenstates(keys_tmp,psi_coefs_inout,n,nmax_keys,nma
|
|||||||
accu_precision_of_diag = 0.d0
|
accu_precision_of_diag = 0.d0
|
||||||
do i = 1, nstates
|
do i = 1, nstates
|
||||||
do j = i+1, nstates
|
do j = i+1, nstates
|
||||||
if( ( dabs(s2(i,i) - s2(j,j)) .le.1.d-10 ) .and. (dabs(s2(i,j) + dabs(s2(i,j)))) .le.1.d-10) then
|
if( ( dabs(s2(i,i) - s2(j,j)) .le.0.5d0 ) ) then
|
||||||
s2(i,j) = 0.d0
|
s2(i,j) = 0.d0
|
||||||
s2(j,i) = 0.d0
|
s2(j,i) = 0.d0
|
||||||
endif
|
endif
|
||||||
|
@ -1389,6 +1389,55 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s
|
|||||||
l = l+1
|
l = l+1
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
else
|
||||||
|
|
||||||
|
print*, 'get_excitation_degree_vector_mono_or_exchange not yet implemented for N_int > 1 ...'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
idx(0) = l-1
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
subroutine get_excitation_degree_vector_double_alpha_beta(key1,key2,degree,Nint,sze,idx)
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Applies get_excitation_degree to an array of determinants and return only the mono excitations
|
||||||
|
! and the connections through exchange integrals
|
||||||
|
END_DOC
|
||||||
|
integer, intent(in) :: Nint, sze
|
||||||
|
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
|
||||||
|
integer(bit_kind), intent(in) :: key2(Nint,2)
|
||||||
|
integer, intent(out) :: degree(sze)
|
||||||
|
integer, intent(out) :: idx(0:sze)
|
||||||
|
integer(bit_kind) :: key_tmp(Nint,2)
|
||||||
|
|
||||||
|
integer :: i,l,d,m
|
||||||
|
integer :: degree_alpha, degree_beta
|
||||||
|
|
||||||
|
ASSERT (Nint > 0)
|
||||||
|
ASSERT (sze > 0)
|
||||||
|
|
||||||
|
l=1
|
||||||
|
if (Nint==1) then
|
||||||
|
|
||||||
|
!DIR$ LOOP COUNT (1000)
|
||||||
|
do i=1,sze
|
||||||
|
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
|
||||||
|
popcnt(xor( key1(1,2,i), key2(1,2)))
|
||||||
|
if (d .ne.4)cycle
|
||||||
|
key_tmp(1,1) = xor(key1(1,1,i),key2(1,1))
|
||||||
|
key_tmp(1,2) = xor(key1(1,2,i),key2(1,2))
|
||||||
|
degree_alpha = popcnt(key_tmp(1,1))
|
||||||
|
degree_beta = popcnt(key_tmp(1,2))
|
||||||
|
if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin
|
||||||
|
degree(l) = ishft(d,-1)
|
||||||
|
idx(l) = i
|
||||||
|
l = l+1
|
||||||
|
enddo
|
||||||
else if (Nint==2) then
|
else if (Nint==2) then
|
||||||
|
|
||||||
!DIR$ LOOP COUNT (1000)
|
!DIR$ LOOP COUNT (1000)
|
||||||
@ -1397,29 +1446,17 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s
|
|||||||
popcnt(xor( key1(1,2,i), key2(1,2))) + &
|
popcnt(xor( key1(1,2,i), key2(1,2))) + &
|
||||||
popcnt(xor( key1(2,1,i), key2(2,1))) + &
|
popcnt(xor( key1(2,1,i), key2(2,1))) + &
|
||||||
popcnt(xor( key1(2,2,i), key2(2,2)))
|
popcnt(xor( key1(2,2,i), key2(2,2)))
|
||||||
exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,2),key2(1,2)))) + &
|
if (d .ne.4)cycle
|
||||||
popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,2),key2(2,2))))
|
key_tmp(1,1) = xor(key1(1,1,i),key2(1,1))
|
||||||
exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + &
|
key_tmp(1,2) = xor(key1(1,2,i),key2(1,2))
|
||||||
popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2))))
|
key_tmp(2,1) = xor(key1(2,1,i),key2(2,1))
|
||||||
if (d > 4)cycle
|
key_tmp(2,2) = xor(key1(2,2,i),key2(2,2))
|
||||||
if (d ==4)then
|
degree_alpha = popcnt(key_tmp(1,1)) + popcnt(key_tmp(2,1))
|
||||||
if(exchange_1 .eq. 0 ) then
|
degree_beta = popcnt(key_tmp(1,2)) + popcnt(key_tmp(2,2))
|
||||||
|
if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin
|
||||||
degree(l) = ishft(d,-1)
|
degree(l) = ishft(d,-1)
|
||||||
idx(l) = i
|
idx(l) = i
|
||||||
l = l+1
|
l = l+1
|
||||||
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
|
|
||||||
degree(l) = ishft(d,-1)
|
|
||||||
idx(l) = i
|
|
||||||
l = l+1
|
|
||||||
else
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
! pause
|
|
||||||
else
|
|
||||||
degree(l) = ishft(d,-1)
|
|
||||||
idx(l) = i
|
|
||||||
l = l+1
|
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
else if (Nint==3) then
|
else if (Nint==3) then
|
||||||
@ -1432,31 +1469,19 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s
|
|||||||
popcnt(xor( key1(2,2,i), key2(2,2))) + &
|
popcnt(xor( key1(2,2,i), key2(2,2))) + &
|
||||||
popcnt(xor( key1(3,1,i), key2(3,1))) + &
|
popcnt(xor( key1(3,1,i), key2(3,1))) + &
|
||||||
popcnt(xor( key1(3,2,i), key2(3,2)))
|
popcnt(xor( key1(3,2,i), key2(3,2)))
|
||||||
exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))) + &
|
if (d .ne.4)cycle
|
||||||
popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,1),key2(2,2)))) + &
|
key_tmp(1,1) = xor(key1(1,1,i),key2(1,1))
|
||||||
popcnt(xor(ior(key1(3,1,i),key1(3,2,i)),ior(key2(3,1),key2(3,2))))
|
key_tmp(1,2) = xor(key1(1,2,i),key2(1,2))
|
||||||
exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + &
|
key_tmp(2,1) = xor(key1(2,1,i),key2(2,1))
|
||||||
popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2)))) + &
|
key_tmp(2,2) = xor(key1(2,2,i),key2(2,2))
|
||||||
popcnt(ior(xor(key1(3,1,i),key2(3,1)),xor(key1(3,2,i),key2(3,2))))
|
key_tmp(3,1) = xor(key1(3,1,i),key2(3,1))
|
||||||
if (d > 4)cycle
|
key_tmp(3,2) = xor(key1(3,2,i),key2(3,2))
|
||||||
if (d ==4)then
|
degree_alpha = popcnt(key_tmp(1,1)) + popcnt(key_tmp(2,1)) + popcnt(key_tmp(3,1))
|
||||||
if(exchange_1 .eq. 0 ) then
|
degree_beta = popcnt(key_tmp(1,2)) + popcnt(key_tmp(2,2)) + popcnt(key_tmp(3,2))
|
||||||
|
if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin
|
||||||
degree(l) = ishft(d,-1)
|
degree(l) = ishft(d,-1)
|
||||||
idx(l) = i
|
idx(l) = i
|
||||||
l = l+1
|
l = l+1
|
||||||
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
|
|
||||||
degree(l) = ishft(d,-1)
|
|
||||||
idx(l) = i
|
|
||||||
l = l+1
|
|
||||||
else
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
! pause
|
|
||||||
else
|
|
||||||
degree(l) = ishft(d,-1)
|
|
||||||
idx(l) = i
|
|
||||||
l = l+1
|
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
else
|
else
|
||||||
@ -1464,39 +1489,26 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s
|
|||||||
!DIR$ LOOP COUNT (1000)
|
!DIR$ LOOP COUNT (1000)
|
||||||
do i=1,sze
|
do i=1,sze
|
||||||
d = 0
|
d = 0
|
||||||
exchange_1 = 0
|
|
||||||
!DIR$ LOOP COUNT MIN(4)
|
!DIR$ LOOP COUNT MIN(4)
|
||||||
do m=1,Nint
|
do m=1,Nint
|
||||||
d = d + popcnt(xor( key1(m,1,i), key2(m,1))) &
|
d = d + popcnt(xor( key1(m,1,i), key2(m,1))) &
|
||||||
+ popcnt(xor( key1(m,2,i), key2(m,2)))
|
+ popcnt(xor( key1(m,2,i), key2(m,2)))
|
||||||
exchange_1 = popcnt(xor(ior(key1(m,1,i),key1(m,2,i)),ior(key2(m,1),key2(m,2))))
|
key_tmp(m,1) = xor(key1(m,1,i),key2(m,1))
|
||||||
exchange_2 = popcnt(ior(xor(key1(m,1,i),key2(m,1)),xor(key1(m,2,i),key2(m,2))))
|
key_tmp(m,2) = xor(key1(m,2,i),key2(m,2))
|
||||||
|
degree_alpha = popcnt(key_tmp(m,1))
|
||||||
|
degree_beta = popcnt(key_tmp(m,2))
|
||||||
enddo
|
enddo
|
||||||
if (d > 4)cycle
|
if(degree_alpha .gt.3 .or. degree_beta .gt.3 )cycle !! no double excitations of same spin
|
||||||
if (d ==4)then
|
|
||||||
if(exchange_1 .eq. 0 ) then
|
|
||||||
degree(l) = ishft(d,-1)
|
degree(l) = ishft(d,-1)
|
||||||
idx(l) = i
|
idx(l) = i
|
||||||
l = l+1
|
l = l+1
|
||||||
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
|
|
||||||
degree(l) = ishft(d,-1)
|
|
||||||
idx(l) = i
|
|
||||||
l = l+1
|
|
||||||
else
|
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
! pause
|
|
||||||
else
|
|
||||||
degree(l) = ishft(d,-1)
|
|
||||||
idx(l) = i
|
|
||||||
l = l+1
|
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
endif
|
endif
|
||||||
idx(0) = l-1
|
idx(0) = l-1
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
subroutine get_excitation_degree_vector_mono_or_exchange_verbose(key1,key2,degree,Nint,sze,idx)
|
subroutine get_excitation_degree_vector_mono_or_exchange_verbose(key1,key2,degree,Nint,sze,idx)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
|
@ -10,7 +10,7 @@ type: logical
|
|||||||
doc: If True, do not compute the bielectronic integrals when 4 indices are virtual
|
doc: If True, do not compute the bielectronic integrals when 4 indices are virtual
|
||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: False
|
default: False
|
||||||
ezfio_name: None
|
ezfio_name: no_vvvv_integrals
|
||||||
|
|
||||||
[disk_access_mo_integrals]
|
[disk_access_mo_integrals]
|
||||||
type: Disk_access
|
type: Disk_access
|
||||||
|
File diff suppressed because it is too large
Load Diff
@ -484,6 +484,7 @@ subroutine map_get(map, key, value)
|
|||||||
integer(map_size_kind) :: idx_cache
|
integer(map_size_kind) :: idx_cache
|
||||||
integer(cache_map_size_kind) :: idx
|
integer(cache_map_size_kind) :: idx
|
||||||
|
|
||||||
|
! index in tha pointers array
|
||||||
idx_cache = ishft(key,map_shift)
|
idx_cache = ishft(key,map_shift)
|
||||||
!DIR$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx)
|
call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx)
|
||||||
@ -853,3 +854,25 @@ subroutine get_cache_map(map,map_idx,keys,values,n_elements)
|
|||||||
enddo
|
enddo
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
subroutine get_cache_map_verbose(map,map_idx)
|
||||||
|
use map_module
|
||||||
|
implicit none
|
||||||
|
type (map_type), intent(in) :: map
|
||||||
|
integer(map_size_kind), intent(in) :: map_idx
|
||||||
|
integer(cache_map_size_kind) :: n_elements
|
||||||
|
integer(key_kind) :: keys(2**16)
|
||||||
|
double precision :: values(2**16)
|
||||||
|
integer(cache_map_size_kind) :: i
|
||||||
|
integer(key_kind) :: shift
|
||||||
|
|
||||||
|
shift = ishft(map_idx,-map_shift)
|
||||||
|
|
||||||
|
n_elements = map%map(map_idx)%n_elements
|
||||||
|
do i=1,n_elements
|
||||||
|
keys(i) = map%map(map_idx)%key(i) + shift
|
||||||
|
values(i) = map%map(map_idx)%value(i)
|
||||||
|
print*, ',key,values',keys(i),values(i)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
end
|
||||||
|
Loading…
Reference in New Issue
Block a user