mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-23 12:56:14 +01:00
now eigenfunction of S^2
This commit is contained in:
Emmanuel Giner
parent
90042a19f4
commit
8a91b293bf
@ -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
|
||||||
|
|
||||||
|
|||||||
@ -520,7 +520,7 @@ END_PROVIDER
|
|||||||
integer :: iorb,jorb,i_ok
|
integer :: iorb,jorb,i_ok
|
||||||
integer :: state_target
|
integer :: state_target
|
||||||
double precision :: energies(n_states)
|
double precision :: energies(n_states)
|
||||||
double precision :: hij
|
double precision :: hij,hij_test
|
||||||
double precision :: norm(N_states,2),norm_no_inv(N_states,2),norm_bis(N_states,2)
|
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 :: energies_alpha_beta(N_states,2)
|
||||||
|
|
||||||
@ -552,11 +552,18 @@ END_PROVIDER
|
|||||||
call debug_det(psi_in_out,N_int)
|
call debug_det(psi_in_out,N_int)
|
||||||
print*, 'pb, i_ok ne 0 !!!'
|
print*, 'pb, i_ok ne 0 !!!'
|
||||||
endif
|
endif
|
||||||
|
! call i_H_j_no_k_operators_from_act(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij_test)
|
||||||
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij)
|
||||||
|
! if(i==1.and.dabs(hij)>1.d-8)then
|
||||||
|
! if(dabs(hij)>1.d-8)then
|
||||||
|
! print*, ispin,vorb,iorb
|
||||||
|
! print*, i,hij,hij_test
|
||||||
|
! pause
|
||||||
|
! endif
|
||||||
do j = 1, n_states
|
do j = 1, n_states
|
||||||
double precision :: coef,contrib
|
double precision :: coef,contrib
|
||||||
coef = psi_coef(i,j) !* psi_coef(i,j)
|
coef = psi_coef(i,j) !* psi_coef(i,j)
|
||||||
psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij
|
psi_in_out_coef(i,j) = coef * hij
|
||||||
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -582,22 +589,18 @@ END_PROVIDER
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
do state_target = 1, N_states
|
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
|
||||||
! 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
|
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,state_target)
|
! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
|
||||||
|
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
|
||||||
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
enddo ! ispin
|
enddo ! ispin
|
||||||
do state_target = 1, N_states
|
do state_target = 1, N_states
|
||||||
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
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 * &
|
! one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall(state_target) - &
|
||||||
! ( energy_cas_dyall(state_target) - energies_alpha_beta(state_target,1) + &
|
one_anhil_one_creat_inact_virt(iorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - &
|
||||||
! 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) ) &
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
||||||
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
||||||
else
|
else
|
||||||
@ -688,24 +691,20 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta
|
|||||||
do state_target = 1, N_states
|
do state_target = 1, N_states
|
||||||
energies_alpha_beta(state_target, ispin) = 0.d0
|
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
|
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,state_target)
|
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
|
||||||
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
enddo ! ispin
|
enddo ! ispin
|
||||||
do state_target = 1, N_states
|
do state_target = 1, N_states
|
||||||
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
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) - &
|
! one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall(state_target) - &
|
||||||
|
one_anhil_inact(iorb,aorb,state_target) = energy_cas_dyall_no_exchange(state_target) - &
|
||||||
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
||||||
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
||||||
else
|
else
|
||||||
one_anhil_inact(iorb,aorb,state_target) = 0.d0
|
one_anhil_inact(iorb,aorb,state_target) = 0.d0
|
||||||
endif
|
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
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -735,7 +734,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
|
|||||||
|
|
||||||
double precision :: thresh_norm
|
double precision :: thresh_norm
|
||||||
|
|
||||||
thresh_norm = 1.d-10
|
thresh_norm = 1.d-20
|
||||||
|
|
||||||
do aorb = 1,n_act_orb
|
do aorb = 1,n_act_orb
|
||||||
orb_a = list_act(aorb)
|
orb_a = list_act(aorb)
|
||||||
@ -791,15 +790,16 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
|
|||||||
do state_target = 1, N_states
|
do state_target = 1, N_states
|
||||||
energies_alpha_beta(state_target, ispin) = 0.d0
|
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
|
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,state_target)
|
! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
|
||||||
! print*, energies(state_target)
|
call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target)
|
||||||
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
energies_alpha_beta(state_target, ispin) += energies(state_target)
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
enddo ! ispin
|
enddo ! ispin
|
||||||
do state_target = 1, N_states
|
do state_target = 1, N_states
|
||||||
if((norm_no_inv(state_target,1) + norm_no_inv(state_target,2)) .ne. 0.d0)then
|
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) - &
|
! one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall(state_target) - &
|
||||||
|
one_creat_virt(aorb,vorb,state_target) = energy_cas_dyall_no_exchange(state_target) - &
|
||||||
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
( energies_alpha_beta(state_target,1) + energies_alpha_beta(state_target,2) ) &
|
||||||
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
/( norm_bis(state_target,1) + norm_bis(state_target,2) )
|
||||||
else
|
else
|
||||||
@ -815,154 +815,6 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State
|
|||||||
enddo
|
enddo
|
||||||
deallocate(psi_in_out,psi_in_out_coef)
|
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_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),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)
|
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from_1h1p_singles)
|
subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from_1h1p_singles)
|
||||||
|
|||||||
@ -212,52 +212,97 @@ double precision function diag_H_mat_elem_no_elec_check(det_in,Nint)
|
|||||||
core_act += 2.d0 * mo_bielec_integral_jj(jorb,iorb) - mo_bielec_integral_jj_exchange(jorb,iorb)
|
core_act += 2.d0 * mo_bielec_integral_jj(jorb,iorb) - mo_bielec_integral_jj_exchange(jorb,iorb)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
! print*,'core_act = ',core_act
|
|
||||||
! print*,'alpha_alpha = ',alpha_alpha
|
|
||||||
! print*,'alpha_beta = ',alpha_beta
|
|
||||||
! print*,'beta_beta = ',beta_beta
|
|
||||||
! print*,'mono_elec = ',mono_elec
|
|
||||||
|
|
||||||
! do i = 1, n_core_inact_orb
|
end
|
||||||
! iorb = list_core_inact(i)
|
|
||||||
! diag_H_mat_elem_no_elec_check += 2.d0 * fock_core_inactive_total_spin_trace(iorb,1)
|
|
||||||
! enddo
|
|
||||||
|
|
||||||
|
|
||||||
!!!!!!!!!!!!
|
|
||||||
return
|
|
||||||
!!!!!!!!!!!!
|
|
||||||
|
|
||||||
|
|
||||||
|
double precision function diag_H_mat_elem_no_elec_check_no_spin(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
|
||||||
|
|
||||||
|
diag_H_mat_elem_no_elec_check_no_spin = 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
|
! alpha - alpha
|
||||||
do i = 1, n_core_inact_orb
|
do i = 1, elec_num_tab_local(1)
|
||||||
iorb = list_core_inact(i)
|
iorb = occ(i,1)
|
||||||
diag_H_mat_elem_no_elec_check += 1.d0 * mo_mono_elec_integral(iorb,iorb)
|
diag_H_mat_elem_no_elec_check_no_spin += mo_mono_elec_integral(iorb,iorb)
|
||||||
do j = i+1, n_core_inact_orb
|
mono_elec += mo_mono_elec_integral(iorb,iorb)
|
||||||
jorb = list_core_inact(j)
|
do j = i+1, elec_num_tab_local(1)
|
||||||
diag_H_mat_elem_no_elec_check += 1.d0 * mo_bielec_integral_jj(jorb,iorb) - 1.d0 * mo_bielec_integral_jj_exchange(jorb,iorb)
|
jorb = occ(j,1)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_spin += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
alpha_alpha += mo_bielec_integral_jj(jorb,iorb)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do i = 1, n_core_inact_orb
|
! beta - beta
|
||||||
iorb = list_core_inact(i)
|
do i = 1, elec_num_tab_local(2)
|
||||||
diag_H_mat_elem_no_elec_check += 1.d0 * mo_mono_elec_integral(iorb,iorb)
|
iorb = occ(i,2)
|
||||||
do j = i+1, n_core_inact_orb
|
diag_H_mat_elem_no_elec_check_no_spin += mo_mono_elec_integral(iorb,iorb)
|
||||||
jorb = list_core_inact(j)
|
mono_elec += mo_mono_elec_integral(iorb,iorb)
|
||||||
diag_H_mat_elem_no_elec_check += 1.d0 * mo_bielec_integral_jj(jorb,iorb) - 1.d0 * mo_bielec_integral_jj_exchange(jorb,iorb)
|
do j = i+1, elec_num_tab_local(2)
|
||||||
|
jorb = occ(j,2)
|
||||||
|
diag_H_mat_elem_no_elec_check_no_spin += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
beta_beta += mo_bielec_integral_jj(jorb,iorb)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do i = 1, n_core_inact_orb
|
|
||||||
iorb = list_core_inact(i)
|
! 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_spin += mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
alpha_beta += 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
|
do j = 1, n_core_inact_orb
|
||||||
jorb = list_core_inact(j)
|
jorb = list_core_inact(j)
|
||||||
diag_H_mat_elem_no_elec_check += 1.d0 * mo_bielec_integral_jj(jorb,iorb)
|
diag_H_mat_elem_no_elec_check_no_spin += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
core_act += 2.d0 * mo_bielec_integral_jj(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_spin += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
|
||||||
|
core_act += 2.d0 * mo_bielec_integral_jj(jorb,iorb)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
subroutine i_H_j_dyall(key_i,key_j,Nint,hij)
|
subroutine i_H_j_dyall(key_i,key_j,Nint,hij)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
@ -389,6 +434,133 @@ subroutine i_H_j_dyall(key_i,key_j,Nint,hij)
|
|||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine i_H_j_dyall_no_spin(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
|
||||||
|
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, phase,phase_2
|
||||||
|
integer :: n_occ_ab(2)
|
||||||
|
logical :: has_mipi(Nint*bit_kind_size)
|
||||||
|
double precision :: mipi(Nint*bit_kind_size), miip(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,1,2) .and. exc(1,1,2) == exc(1,2,1) )then
|
||||||
|
hij = 0.d0
|
||||||
|
else
|
||||||
|
hij = phase*get_mo_bielec_integral( &
|
||||||
|
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( &
|
||||||
|
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( &
|
||||||
|
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(m,i,p,i,mo_integrals_map)
|
||||||
|
miip(i) = get_mo_bielec_integral(m,i,i,p,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(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)) !- miip(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(m,i,p,i,mo_integrals_map)
|
||||||
|
miip(i) = get_mo_bielec_integral(m,i,i,p,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(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)) !- miip(occ(k,2))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
endif
|
||||||
|
hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) )
|
||||||
|
|
||||||
|
case (0)
|
||||||
|
double precision :: diag_H_mat_elem_no_elec_check_no_spin
|
||||||
|
hij = diag_H_mat_elem_no_elec_check_no_spin(key_i,Nint)
|
||||||
|
end select
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
|
subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
@ -414,6 +586,7 @@ subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coe
|
|||||||
do j = 1, ndet
|
do j = 1, ndet
|
||||||
if(psi_coef_tmp(j)==0.d0)cycle
|
if(psi_coef_tmp(j)==0.d0)cycle
|
||||||
call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
|
call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
|
||||||
|
! call i_H_j_dyall_no_spin(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
|
||||||
accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij
|
accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
@ -393,126 +393,6 @@ 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
|
|
||||||
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(iorb,aorb,rorb,borb,mo_integrals_map) &
|
|
||||||
+ get_mo_bielec_integral(iorb,aorb,borb,rorb,mo_integrals_map)
|
|
||||||
else
|
|
||||||
hij = get_mo_bielec_integral(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)
|
subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
|
|||||||
@ -2,7 +2,17 @@ program print_1h2p
|
|||||||
implicit none
|
implicit none
|
||||||
read_wf = .True.
|
read_wf = .True.
|
||||||
touch read_wf
|
touch read_wf
|
||||||
call routine
|
call routine_2
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine routine_2
|
||||||
|
implicit none
|
||||||
|
integer :: i,j
|
||||||
|
do i =1, n_inact_orb
|
||||||
|
write(*,'(100(F16.10,X))')one_anhil_one_creat_inact_virt(i,:,1)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
subroutine routine
|
subroutine routine
|
||||||
|
|||||||
@ -293,27 +293,38 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
|
|||||||
if (n_holes_act == 0 .and. n_particles_act == 1) then
|
if (n_holes_act == 0 .and. n_particles_act == 1) then
|
||||||
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)
|
call get_excitation_degree(det_1,det_2,degree,N_int)
|
||||||
! if(degree == 1)then
|
if(degree == 1)then
|
||||||
! call get_excitation(det_1,det_2,exc,degree,phase,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)
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||||
! i_hole = list_inact_reverse(h1)
|
i_hole = list_inact_reverse(h1)
|
||||||
! i_part = list_act_reverse(p1)
|
i_part = list_act_reverse(p1)
|
||||||
! do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
! delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state)
|
delta_e_act(i_state) += one_anhil_inact(i_hole,i_part,i_state)
|
||||||
! enddo
|
enddo
|
||||||
! else if (degree == 2)then
|
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
|
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
|
||||||
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
|
||||||
|
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
|
||||||
! first hole
|
! first hole
|
||||||
@ -415,7 +426,7 @@ subroutine get_delta_e_dyall(det_1,det_2,coef_array,hij,delta_e_final)
|
|||||||
i_hole = list_inact_reverse(h1)
|
i_hole = list_inact_reverse(h1)
|
||||||
i_part = list_virt_reverse(p1)
|
i_part = list_virt_reverse(p1)
|
||||||
do i_state = 1, N_states
|
do i_state = 1, N_states
|
||||||
! delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state)
|
delta_e_act(i_state) += one_anhil_one_creat_inact_virt(i_hole,i_part,i_state)
|
||||||
enddo
|
enddo
|
||||||
endif
|
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
|
||||||
|
|||||||
183
plugins/MRPT_Utils/special_hij.irp.f
Normal file
183
plugins/MRPT_Utils/special_hij.irp.f
Normal file
@ -0,0 +1,183 @@
|
|||||||
|
|
||||||
|
|
||||||
|
subroutine i_H_j_no_k_operators_from_act(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, phase
|
||||||
|
integer :: m,n,p,q
|
||||||
|
integer :: i,j,k
|
||||||
|
integer :: occ(Nint*bit_kind_size,2)
|
||||||
|
double precision :: diag_H_mat_elem
|
||||||
|
integer :: n_occ_ab(2)
|
||||||
|
logical :: has_mipi(Nint*bit_kind_size)
|
||||||
|
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size), miip_other(Nint*bit_kind_size)
|
||||||
|
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
|
||||||
|
|
||||||
|
ASSERT (Nint > 0)
|
||||||
|
ASSERT (Nint == N_int)
|
||||||
|
ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
|
||||||
|
ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
|
||||||
|
ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
|
||||||
|
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
|
||||||
|
|
||||||
|
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
|
||||||
|
hij = phase*get_mo_bielec_integral( &
|
||||||
|
exc(1,1,1), &
|
||||||
|
exc(1,1,2), &
|
||||||
|
exc(1,2,1), &
|
||||||
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
|
else if (exc(0,1,1) == 2) then
|
||||||
|
! Double alpha
|
||||||
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
|
exc(1,1,1), &
|
||||||
|
exc(2,1,1), &
|
||||||
|
exc(1,2,1), &
|
||||||
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
|
get_mo_bielec_integral( &
|
||||||
|
exc(1,1,1), &
|
||||||
|
exc(2,1,1), &
|
||||||
|
exc(2,2,1), &
|
||||||
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
|
else if (exc(0,1,2) == 2) then
|
||||||
|
! Double beta
|
||||||
|
hij = phase*(get_mo_bielec_integral( &
|
||||||
|
exc(1,1,2), &
|
||||||
|
exc(2,1,2), &
|
||||||
|
exc(1,2,2), &
|
||||||
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
|
get_mo_bielec_integral( &
|
||||||
|
exc(1,1,2), &
|
||||||
|
exc(2,1,2), &
|
||||||
|
exc(2,2,2), &
|
||||||
|
exc(1,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.
|
||||||
|
logical :: is_i_in_active
|
||||||
|
double precision :: accu_a, accu_b, accu_core
|
||||||
|
accu_a = 0.d0
|
||||||
|
accu_b = 0.d0
|
||||||
|
if (exc(0,1,1) == 1) then
|
||||||
|
! Mono alpha
|
||||||
|
m = exc(1,1,1)
|
||||||
|
p = exc(1,2,1)
|
||||||
|
do k = 1, elec_alpha_num
|
||||||
|
i = occ(k,1)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
if(.not.is_i_in_active(i))then
|
||||||
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
else
|
||||||
|
! print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
accu_a += miip(i)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do k = 1, elec_beta_num
|
||||||
|
i = occ(k,2)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
if(.not.is_i_in_active(i))then
|
||||||
|
miip_other(i) = 0.d0
|
||||||
|
else
|
||||||
|
! print*, i,get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
accu_b += miip(i)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
! print*, accu_a,accu_b,accu_a + accu_b
|
||||||
|
accu_a = 0.d0
|
||||||
|
accu_b = 0.d0
|
||||||
|
accu_core = mo_mono_elec_integral(m,p)
|
||||||
|
|
||||||
|
do k = 1, elec_alpha_num
|
||||||
|
hij = hij + mipi(occ(k,1)) - miip(occ(k,1))
|
||||||
|
accu_Core += mipi(occ(k,1))
|
||||||
|
if(is_i_in_active(occ(k,1)))then
|
||||||
|
accu_a += miip(occ(k,1))
|
||||||
|
else
|
||||||
|
accu_Core -= miip(occ(k,1))
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
! print*, hij,accu_core
|
||||||
|
do k = 1, elec_beta_num
|
||||||
|
hij = hij + mipi(occ(k,2)) - miip_other(occ(k,2))
|
||||||
|
accu_Core += mipi(occ(k,2))
|
||||||
|
if(is_i_in_active(occ(k,2)))then
|
||||||
|
accu_b += miip_other(occ(k,2))
|
||||||
|
else
|
||||||
|
accu_Core -= miip_other(occ(k,2))
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
! print*, hij,accu_core,accu_core - accu_a - accu_b
|
||||||
|
! print*, accu_a,accu_b,accu_a + accu_b
|
||||||
|
! pause
|
||||||
|
|
||||||
|
else
|
||||||
|
! Mono beta
|
||||||
|
m = exc(1,1,2)
|
||||||
|
p = exc(1,2,2)
|
||||||
|
do k = 1, elec_beta_num
|
||||||
|
i = occ(k,2)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
if(.not.is_i_in_active(i))then
|
||||||
|
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
else
|
||||||
|
miip(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do k = 1, elec_alpha_num
|
||||||
|
i = occ(k,1)
|
||||||
|
if (.not.has_mipi(i)) then
|
||||||
|
mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
|
||||||
|
has_mipi(i) = .True.
|
||||||
|
endif
|
||||||
|
if(.not.is_i_in_active(i))then
|
||||||
|
miip_other(i) = 0.d0
|
||||||
|
else
|
||||||
|
miip_other(i) = 1.0d0 * get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
||||||
|
endif
|
||||||
|
enddo
|
||||||
|
do k = 1, elec_alpha_num
|
||||||
|
hij = hij + mipi(occ(k,1)) - miip_other(occ(k,1))
|
||||||
|
enddo
|
||||||
|
do k = 1, elec_beta_num
|
||||||
|
hij = hij + mipi(occ(k,2)) - miip(occ(k,2))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
endif
|
||||||
|
hij = phase*(hij + mo_mono_elec_integral(m,p))
|
||||||
|
|
||||||
|
case (0)
|
||||||
|
hij = diag_H_mat_elem(key_i,Nint)
|
||||||
|
end select
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -560,3 +560,24 @@ logical function is_i_in_virtual(i)
|
|||||||
endif
|
endif
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
logical function is_i_in_active(i)
|
||||||
|
implicit none
|
||||||
|
integer,intent(in) :: i
|
||||||
|
integer(bit_kind) :: key(N_int)
|
||||||
|
integer :: k,j
|
||||||
|
integer :: accu
|
||||||
|
is_i_in_active = .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),cas_bitmask(k,1,1)))
|
||||||
|
enddo
|
||||||
|
if(accu .ne. 0)then
|
||||||
|
is_i_in_active= .True.
|
||||||
|
endif
|
||||||
|
|
||||||
|
end
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user