10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-14 01:53:55 +01:00
quantum_package/plugins/Perturbation/pert_single.irp.f

50 lines
1.4 KiB
Fortran
Raw Normal View History

2015-03-19 21:14:52 +01:00
subroutine pt2_h_core(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem, h
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
integer :: exc(0:2,2,2)
integer :: degree
double precision :: phase
call get_excitation(ref_bitmask,det_pert,exc,degree,phase,N_int)
h = diag_H_mat_elem(det_pert,N_int)
print*,'delta E = ',h-ref_bitmask_energy
if(h<ref_bitmask_energy)then
c_pert = 0.d0
e_2_pert = 0.d0
H_pert_diag = 0.d0
return
endif
if(degree>1)then
c_pert = 0.d0
e_2_pert = 0.d0
H_pert_diag = 0.d0
return
endif
integer :: h1,p1,h2,p2,s1,s2
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
c_pert = phase * mo_mono_elec_integral(h1,p1)
e_2_pert = -dabs(mo_mono_elec_integral(h1,p1)+1.d0)
end