890 lines
31 KiB
Fortran
890 lines
31 KiB
Fortran
subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
|
double precision, intent(out) :: diag_H_elements(dim_in)
|
|
double precision, intent(in) :: convergence
|
|
|
|
integer :: i,j,k,l
|
|
integer :: n_singles
|
|
integer :: index_singles(sze),hole_particles_singles(sze,3)
|
|
integer :: n_doubles
|
|
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
|
|
integer :: index_hf
|
|
double precision :: e_corr_singles(mo_tot_num,2)
|
|
double precision :: e_corr_doubles(mo_tot_num)
|
|
double precision :: e_corr_singles_total(2)
|
|
double precision :: e_corr_doubles_1h1p
|
|
|
|
integer :: exc(0:2,2,2),degree
|
|
integer :: h1,h2,p1,p2,s1,s2
|
|
integer :: other_spin(2)
|
|
double precision :: phase
|
|
integer(bit_kind) :: key_tmp(N_int,2)
|
|
integer :: i_ok
|
|
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral
|
|
double precision :: hij,c_ref,contrib
|
|
integer :: iorb
|
|
|
|
other_spin(1) = 2
|
|
other_spin(2) = 1
|
|
|
|
n_singles = 0
|
|
n_doubles = 0
|
|
do i = 1,sze
|
|
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
call i_H_j(dets_in(1,1,i),dets_in(1,1,i),N_int,hij)
|
|
diag_H_elements(i) = hij
|
|
if(degree == 0)then
|
|
index_hf = i
|
|
else if (degree == 1)then
|
|
n_singles +=1
|
|
index_singles(n_singles) = i
|
|
! h1 = inactive orbital of the hole
|
|
hole_particles_singles(n_singles,1) = h1
|
|
! p1 = virtual orbital of the particle
|
|
hole_particles_singles(n_singles,2) = p1
|
|
! s1 = spin of the electron excited
|
|
hole_particles_singles(n_singles,3) = s1
|
|
else if (degree == 2)then
|
|
n_doubles +=1
|
|
index_doubles(n_doubles) = i
|
|
! h1 = inactive orbital of the hole (beta of course)
|
|
hole_particles_doubles(n_doubles,1) = h1
|
|
! p1 = virtual orbital of the particle (alpha of course)
|
|
hole_particles_doubles(n_doubles,2) = p2
|
|
else
|
|
print*,'PB !! found out other thing than a single or double'
|
|
print*,'stopping ..'
|
|
stop
|
|
endif
|
|
enddo
|
|
|
|
e_corr_singles = 0.d0
|
|
e_corr_doubles = 0.d0
|
|
e_corr_singles_total = 0.d0
|
|
e_corr_doubles_1h1p = 0.d0
|
|
c_ref = 1.d0/u_in(index_hf,1)
|
|
print*,'c_ref = ',c_ref
|
|
do i = 1,sze
|
|
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
call i_H_j(ref_bitmask,dets_in(1,1,i),N_int,hij)
|
|
contrib = hij * u_in(i,1) * c_ref
|
|
if (degree == 1)then
|
|
e_corr_singles(h1,s1) += contrib
|
|
e_corr_singles(p1,s1) += contrib
|
|
e_corr_singles_total(s1)+= contrib
|
|
else if (degree == 2)then
|
|
e_corr_doubles_1h1p += contrib
|
|
e_corr_doubles(h1) += contrib
|
|
e_corr_doubles(p2) += contrib
|
|
endif
|
|
enddo
|
|
print*,'e_corr_singles alpha = ',e_corr_singles_total(1)
|
|
print*,'e_corr_singles beta = ',e_corr_singles_total(2)
|
|
print*,'e_corr_doubles_1h1p = ',e_corr_doubles_1h1p
|
|
|
|
! repeat all the correlation energy on the singles
|
|
do i = 1,n_singles
|
|
! you can repeat all the correlation energy of the single excitation of the other spin
|
|
diag_H_elements(index_singles(i)) += e_corr_singles_total(other_spin(hole_particles_singles(i,3)))
|
|
|
|
! you can repeat all the correlation energy of the single excitation of the same spin
|
|
do j = 1, n_inact_orb
|
|
iorb = list_inact(j)
|
|
! except the one of the hole
|
|
if(iorb == hole_particles_singles(i,1))cycle
|
|
! ispin = hole_particles_singles(i,3)
|
|
diag_H_elements(index_singles(i)) += e_corr_singles(iorb,hole_particles_singles(i,3))
|
|
enddo
|
|
! also exclude all the energy coming from the virtual orbital
|
|
diag_H_elements(index_singles(i)) -= e_corr_singles(hole_particles_singles(i,2),hole_particles_singles(i,3))
|
|
|
|
! If it is a single excitation alpha, you can repeat :
|
|
! +) all the double excitation 1h1p, appart the part involving the virtual orbital "r"
|
|
! If it is a single excitation alpha, you can repeat :
|
|
! +) all the double excitation 1h1p, appart the part involving the inactive orbital "i"
|
|
diag_H_elements(index_singles(i)) += e_corr_doubles_1h1p
|
|
if(hole_particles_singles(i,3) == 1)then ! alpha single excitation
|
|
diag_H_elements(index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,2))
|
|
else ! beta single exctitation
|
|
diag_H_elements(index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,1))
|
|
endif
|
|
enddo
|
|
|
|
! repeat all the correlation energy on the doubles
|
|
! as all the doubles involve the active space, you cannot repeat any of them one on another
|
|
do i = 1, n_doubles
|
|
! on a given double, you can repeat all the correlation energy of the singles alpha
|
|
do j = 1, n_inact_orb
|
|
iorb = list_inact(j)
|
|
! ispin = hole_particles_singles(i,3)
|
|
diag_H_elements(index_doubles(i)) += e_corr_singles(iorb,1)
|
|
enddo
|
|
! except the part involving the virtual orbital "hole_particles_doubles(i,2)"
|
|
diag_H_elements(index_doubles(i)) -= e_corr_singles(hole_particles_doubles(i,2),1)
|
|
! on a given double, you can repeat all the correlation energy of the singles beta
|
|
do j = 1, n_inact_orb
|
|
iorb = list_inact(j)
|
|
! except the one of the hole
|
|
if(iorb == hole_particles_doubles(i,1))cycle
|
|
! ispin = hole_particles_singles(i,3)
|
|
diag_H_elements(index_doubles(i)) += e_corr_singles(iorb,2)
|
|
enddo
|
|
enddo
|
|
|
|
|
|
! Taking into account the connected part of the 2h2p on the HF determinant
|
|
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
|
|
|
|
! diag_H_elements(index_hf) += total_corr_e_2h2p
|
|
return
|
|
c_ref = c_ref * c_ref
|
|
print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf)
|
|
do i = 1, n_singles
|
|
! start on the single excitation "|i>"
|
|
h1 = hole_particles_singles(i,1)
|
|
p1 = hole_particles_singles(i,2)
|
|
do j = 1, n_singles
|
|
do k = 1, N_int
|
|
key_tmp(k,1) = dets_in(k,1,index_singles(i))
|
|
key_tmp(k,2) = dets_in(k,2,index_singles(i))
|
|
enddo
|
|
h2 = hole_particles_singles(j,1)
|
|
p2 = hole_particles_singles(j,2)
|
|
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
|
|
! apply the excitation operator from the single excitation "|j>"
|
|
if(i_ok .ne. 1)cycle
|
|
double precision :: phase_ref_other_single,diag_H_mat_elem,hijj,contrib_e2,coef_1
|
|
call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_single_double,N_int)
|
|
call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_ref_other_single,N_int)
|
|
call i_H_j(ref_bitmask,key_tmp,N_int,hij)
|
|
diag_H_elements(index_hf) += u_in(index_singles(i),1) * u_in(index_singles(j),1) * c_ref * hij &
|
|
* phase_single_double * phase_ref_other_single
|
|
enddo
|
|
enddo
|
|
print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf)
|
|
|
|
end
|
|
|
|
|
|
subroutine dressing_1h1p_by_2h2p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
|
double precision, intent(out) :: diag_H_elements(0:dim_in)
|
|
double precision, intent(in) :: convergence
|
|
|
|
integer :: i,j,k,l
|
|
integer :: r,s,i0,j0,r0,s0
|
|
integer :: n_singles
|
|
integer :: index_singles(sze),hole_particles_singles(sze,3)
|
|
integer :: n_doubles
|
|
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
|
|
integer :: index_hf
|
|
double precision :: e_corr_singles(mo_tot_num,2)
|
|
double precision :: e_corr_doubles(mo_tot_num)
|
|
double precision :: e_corr_singles_total(2)
|
|
double precision :: e_corr_doubles_1h1p
|
|
|
|
integer :: exc(0:2,2,2),degree
|
|
integer :: h1,h2,p1,p2,s1,s2
|
|
integer :: other_spin(2)
|
|
double precision :: phase
|
|
integer(bit_kind) :: key_tmp(N_int,2)
|
|
integer :: i_ok
|
|
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral
|
|
double precision :: hij,c_ref,contrib
|
|
integer :: iorb
|
|
|
|
other_spin(1) = 2
|
|
other_spin(2) = 1
|
|
|
|
n_singles = 0
|
|
n_doubles = 0
|
|
do i = 1,sze
|
|
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
call i_H_j(dets_in(1,1,i),dets_in(1,1,i),N_int,hij)
|
|
diag_H_elements(i) = hij
|
|
if(degree == 0)then
|
|
index_hf = i
|
|
else if (degree == 1)then
|
|
n_singles +=1
|
|
index_singles(n_singles) = i
|
|
! h1 = inactive orbital of the hole
|
|
hole_particles_singles(n_singles,1) = h1
|
|
! p1 = virtual orbital of the particle
|
|
hole_particles_singles(n_singles,2) = p1
|
|
! s1 = spin of the electron excited
|
|
hole_particles_singles(n_singles,3) = s1
|
|
else if (degree == 2)then
|
|
n_doubles +=1
|
|
index_doubles(n_doubles) = i
|
|
! h1 = inactive orbital of the hole (beta of course)
|
|
hole_particles_doubles(n_doubles,1) = h1
|
|
! p1 = virtual orbital of the particle (alpha of course)
|
|
hole_particles_doubles(n_doubles,2) = p2
|
|
else
|
|
print*,'PB !! found out other thing than a single or double'
|
|
print*,'stopping ..'
|
|
stop
|
|
endif
|
|
enddo
|
|
double precision :: delta_e
|
|
double precision :: coef_ijrs
|
|
diag_H_elements = 0.d0
|
|
do i0 = 1, n_core_inact_orb
|
|
i= list_core_inact(i0)
|
|
do j0 = i0+1, n_core_inact_orb
|
|
j = list_core_inact(j0)
|
|
print*, i,j
|
|
do r0 = 1, n_virt_orb
|
|
r = list_virt(r0)
|
|
do s0 = r0+1, n_virt_orb
|
|
s = list_virt(s0)
|
|
!!! alpha (i-->r) / beta (j-->s)
|
|
s1 = 1
|
|
s2 = 2
|
|
key_tmp = ref_bitmask
|
|
call do_mono_excitation(key_tmp,i,r,s1,i_ok)
|
|
if(i_ok .ne.1)then
|
|
print*, 'pb !!'
|
|
stop
|
|
endif
|
|
call do_mono_excitation(key_tmp,j,s,s2,i_ok)
|
|
if(i_ok .ne.1)then
|
|
print*, 'pb !!'
|
|
stop
|
|
endif
|
|
call i_H_j(ref_bitmask, key_tmp, N_int,hij)
|
|
delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s)
|
|
coef_ijrs = hij/delta_e
|
|
do k = 1, n_singles
|
|
l = index_singles(k)
|
|
call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij)
|
|
diag_H_elements(l) += coef_ijrs * hij
|
|
enddo
|
|
!if(i>j.and.r>s)then
|
|
!! alpha (i-->r) / alpha (j-->s)
|
|
s1 = 1
|
|
s2 = 1
|
|
key_tmp = ref_bitmask
|
|
call do_mono_excitation(key_tmp,i,r,s1,i_ok)
|
|
if(i_ok .ne.1)then
|
|
print*, 'pb !!'
|
|
stop
|
|
endif
|
|
call do_mono_excitation(key_tmp,j,s,s2,i_ok)
|
|
if(i_ok .ne.1)then
|
|
print*, 'pb !!'
|
|
stop
|
|
endif
|
|
call i_H_j(ref_bitmask, key_tmp, N_int,hij)
|
|
delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s)
|
|
coef_ijrs = hij/delta_e
|
|
do k = 1, n_singles
|
|
l = index_singles(k)
|
|
call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij)
|
|
diag_H_elements(l) += coef_ijrs * hij
|
|
enddo
|
|
!! beta (i-->r) / beta (j-->s)
|
|
s1 = 2
|
|
s2 = 2
|
|
key_tmp = ref_bitmask
|
|
call do_mono_excitation(key_tmp,i,r,s1,i_ok)
|
|
if(i_ok .ne.1)then
|
|
print*, 'pb !!'
|
|
stop
|
|
endif
|
|
call do_mono_excitation(key_tmp,j,s,s2,i_ok)
|
|
if(i_ok .ne.1)then
|
|
print*, 'pb !!'
|
|
stop
|
|
endif
|
|
call i_H_j(ref_bitmask, key_tmp, N_int,hij)
|
|
delta_e = Fock_matrix_diag_mo(i) + Fock_matrix_diag_mo(j) - Fock_matrix_diag_mo(r) - Fock_matrix_diag_mo(s)
|
|
coef_ijrs = hij/delta_e
|
|
do k = 1, n_singles
|
|
l = index_singles(k)
|
|
call i_H_j(dets_in(1,1,l), key_tmp, N_int,hij)
|
|
diag_H_elements(l) += coef_ijrs * hij
|
|
enddo
|
|
!endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
c_ref = 1.d0/u_in(index_hf,1)
|
|
do k = 1, n_singles
|
|
l = index_singles(k)
|
|
diag_H_elements(0) -= diag_H_elements(l)
|
|
enddo
|
|
! do k = 1, n_doubles
|
|
! l = index_doubles(k)
|
|
! diag_H_elements(0) += diag_H_elements(l)
|
|
! enddo
|
|
|
|
|
|
end
|
|
|
|
|
|
subroutine dressing_1h1p_full(dets_in,u_in,H_matrix,dim_in,sze,N_st,Nint,convergence)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(in) :: u_in(dim_in,N_st)
|
|
double precision, intent(inout) :: H_matrix(sze,sze)
|
|
double precision, intent(in) :: convergence
|
|
|
|
integer :: i,j,k,l
|
|
integer :: n_singles
|
|
integer :: index_singles(sze),hole_particles_singles(sze,3)
|
|
integer :: n_doubles
|
|
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
|
|
integer :: index_hf
|
|
double precision :: e_corr_singles(mo_tot_num,2)
|
|
double precision :: e_corr_doubles(mo_tot_num)
|
|
double precision :: e_corr_singles_total(2)
|
|
double precision :: e_corr_doubles_1h1p
|
|
|
|
integer :: exc(0:2,2,2),degree
|
|
integer :: h1,h2,p1,p2,s1,s2
|
|
integer :: other_spin(2)
|
|
double precision :: phase
|
|
integer(bit_kind) :: key_tmp(N_int,2)
|
|
integer :: i_ok
|
|
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral
|
|
double precision :: hij,c_ref,contrib
|
|
integer :: iorb
|
|
|
|
other_spin(1) = 2
|
|
other_spin(2) = 1
|
|
|
|
n_singles = 0
|
|
n_doubles = 0
|
|
do i = 1,sze
|
|
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
if(degree == 0)then
|
|
index_hf = i
|
|
else if (degree == 1)then
|
|
n_singles +=1
|
|
index_singles(n_singles) = i
|
|
! h1 = inactive orbital of the hole
|
|
hole_particles_singles(n_singles,1) = h1
|
|
! p1 = virtual orbital of the particle
|
|
hole_particles_singles(n_singles,2) = p1
|
|
! s1 = spin of the electron excited
|
|
hole_particles_singles(n_singles,3) = s1
|
|
else if (degree == 2)then
|
|
n_doubles +=1
|
|
index_doubles(n_doubles) = i
|
|
! h1 = inactive orbital of the hole (beta of course)
|
|
hole_particles_doubles(n_doubles,1) = h1
|
|
! p1 = virtual orbital of the particle (alpha of course)
|
|
hole_particles_doubles(n_doubles,2) = p2
|
|
else
|
|
print*,'PB !! found out other thing than a single or double'
|
|
print*,'stopping ..'
|
|
stop
|
|
endif
|
|
enddo
|
|
double precision, allocatable :: dressing_H_mat_elem(:)
|
|
allocate(dressing_H_mat_elem(N_det))
|
|
logical :: lmct
|
|
dressing_H_mat_elem = 0.d0
|
|
call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det)
|
|
lmct = .False.
|
|
call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,1000)
|
|
lmct = .true.
|
|
call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,1000)
|
|
do i = 1, N_det
|
|
H_matrix(i,i) += dressing_H_mat_elem(i)
|
|
enddo
|
|
|
|
e_corr_singles = 0.d0
|
|
e_corr_doubles = 0.d0
|
|
e_corr_singles_total = 0.d0
|
|
e_corr_doubles_1h1p = 0.d0
|
|
c_ref = 1.d0/u_in(index_hf,1)
|
|
print*,'c_ref = ',c_ref
|
|
do i = 1,sze
|
|
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
call i_H_j(ref_bitmask,dets_in(1,1,i),N_int,hij)
|
|
contrib = hij * u_in(i,1) * c_ref
|
|
if (degree == 1)then
|
|
e_corr_singles(h1,s1) += contrib
|
|
e_corr_singles(p1,s1) += contrib
|
|
e_corr_singles_total(s1)+= contrib
|
|
else if (degree == 2)then
|
|
e_corr_doubles_1h1p += contrib
|
|
e_corr_doubles(h1) += contrib
|
|
e_corr_doubles(p2) += contrib
|
|
endif
|
|
enddo
|
|
print*,'e_corr_singles alpha = ',e_corr_singles_total(1)
|
|
print*,'e_corr_singles beta = ',e_corr_singles_total(2)
|
|
print*,'e_corr_doubles_1h1p = ',e_corr_doubles_1h1p
|
|
|
|
|
|
! repeat all the correlation energy on the singles
|
|
! do i = 1,n_singles
|
|
! ! you can repeat all the correlation energy of the single excitation of the other spin
|
|
! H_matrix(index_singles(i),index_singles(i)) += e_corr_singles_total(other_spin(hole_particles_singles(i,3)))
|
|
|
|
! ! you can repeat all the correlation energy of the single excitation of the same spin
|
|
! do j = 1, n_inact_orb
|
|
! iorb = list_inact(j)
|
|
! ! except the one of the hole
|
|
! if(iorb == hole_particles_singles(i,1))cycle
|
|
! ! ispin = hole_particles_singles(i,3)
|
|
! H_matrix(index_singles(i),index_singles(i)) += e_corr_singles(iorb,hole_particles_singles(i,3))
|
|
! enddo
|
|
! ! also exclude all the energy coming from the virtual orbital
|
|
! H_matrix(index_singles(i),index_singles(i)) -= e_corr_singles(hole_particles_singles(i,2),hole_particles_singles(i,3))
|
|
!
|
|
! ! If it is a single excitation alpha, you can repeat :
|
|
! ! +) all the double excitation 1h1p, appart the part involving the virtual orbital "r"
|
|
! ! If it is a single excitation alpha, you can repeat :
|
|
! ! +) all the double excitation 1h1p, appart the part involving the inactive orbital "i"
|
|
! H_matrix(index_singles(i),index_singles(i)) += e_corr_doubles_1h1p
|
|
! if(hole_particles_singles(i,3) == 1)then ! alpha single excitation
|
|
! H_matrix(index_singles(i),index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,2))
|
|
! else ! beta single exctitation
|
|
! H_matrix(index_singles(i),index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,1))
|
|
! endif
|
|
! enddo
|
|
|
|
! ! repeat all the correlation energy on the doubles
|
|
! ! as all the doubles involve the active space, you cannot repeat any of them one on another
|
|
! do i = 1, n_doubles
|
|
! ! on a given double, you can repeat all the correlation energy of the singles alpha
|
|
! do j = 1, n_inact_orb
|
|
! iorb = list_inact(j)
|
|
! ! ispin = hole_particles_singles(i,3)
|
|
! H_matrix(index_doubles(i),index_doubles(i)) += e_corr_singles(iorb,1)
|
|
! enddo
|
|
! ! except the part involving the virtual orbital "hole_particles_doubles(i,2)"
|
|
! H_matrix(index_doubles(i),index_doubles(i)) -= e_corr_singles(hole_particles_doubles(i,2),1)
|
|
! ! on a given double, you can repeat all the correlation energy of the singles beta
|
|
! do j = 1, n_inact_orb
|
|
! iorb = list_inact(j)
|
|
! ! except the one of the hole
|
|
! if(iorb == hole_particles_doubles(i,1))cycle
|
|
! ! ispin = hole_particles_singles(i,3)
|
|
! H_matrix(index_doubles(i),index_doubles(i)) += e_corr_singles(iorb,2)
|
|
! enddo
|
|
! enddo
|
|
|
|
|
|
! Taking into account the connected part of the 2h2p on the HF determinant
|
|
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
|
|
|
|
! H_matrix(index_hf) += total_corr_e_2h2p
|
|
print*,'H_matrix(index_hf,index_hf) = ',H_matrix(index_hf,index_hf)
|
|
do i = 1, n_singles
|
|
! start on the single excitation "|i>"
|
|
h1 = hole_particles_singles(i,1)
|
|
p1 = hole_particles_singles(i,2)
|
|
print*,'i = ',i
|
|
do j = i+1, n_singles
|
|
do k = 1, N_int
|
|
key_tmp(k,1) = dets_in(k,1,index_singles(i))
|
|
key_tmp(k,2) = dets_in(k,2,index_singles(i))
|
|
enddo
|
|
h2 = hole_particles_singles(j,1)
|
|
p2 = hole_particles_singles(j,2)
|
|
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
|
|
! apply the excitation operator from the single excitation "|j>"
|
|
if(i_ok .ne. 1)cycle
|
|
double precision :: H_array(sze),diag_H_mat_elem,hjj
|
|
do k = 1, sze
|
|
call get_excitation_degree(dets_in(1,1,k),key_tmp,degree,N_int)
|
|
H_array(k) = 0.d0
|
|
if(degree > 2)cycle
|
|
call i_H_j(dets_in(1,1,k),key_tmp,N_int,hij)
|
|
H_array(k) = hij
|
|
enddo
|
|
hjj = 1.d0/(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
|
|
! contrib_e2 = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * hij * hij))
|
|
do l = 2, sze
|
|
! pause
|
|
H_matrix(l,l) += H_array(l) * H_array(l) * hjj
|
|
! H_matrix(1,l) += H_array(1) * H_array(l) * hjj
|
|
! H_matrix(l,1) += H_array(1) * H_array(l) * hjj
|
|
enddo
|
|
enddo
|
|
enddo
|
|
print*,'H_matrix(index_hf,index_hf) = ',H_matrix(index_hf,index_hf)
|
|
|
|
end
|
|
|
|
subroutine SC2_1h1p_full(dets_in,u_in,energies,H_matrix,dim_in,sze,N_st,Nint,convergence)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
|
double precision, intent(out) :: energies(N_st)
|
|
double precision, intent(out) :: H_matrix(sze,sze)
|
|
double precision, intent(in) :: convergence
|
|
integer :: i,j,iter
|
|
print*,'sze = ',sze
|
|
H_matrix = 0.d0
|
|
do iter = 1, 1
|
|
! if(sze<=N_det_max_jacobi)then
|
|
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
|
|
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze))
|
|
H_matrix_tmp = 0.d0
|
|
call dressing_1h1p_full(dets_in,u_in,H_matrix_tmp,dim_in,sze,N_st,Nint,convergence)
|
|
do j=1,sze
|
|
do i=1,sze
|
|
H_matrix_tmp(i,j) += H_matrix_all_dets(i,j)
|
|
enddo
|
|
enddo
|
|
print*,'passed the dressing'
|
|
call lapack_diag(eigenvalues,eigenvectors, &
|
|
H_matrix_tmp,size(H_matrix_all_dets,1),sze)
|
|
do j=1,min(N_states_diag,sze)
|
|
do i=1,sze
|
|
u_in(i,j) = eigenvectors(i,j)
|
|
enddo
|
|
energies(j) = eigenvalues(j)
|
|
enddo
|
|
deallocate (H_matrix_tmp, eigenvalues, eigenvectors)
|
|
! else
|
|
! call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants)
|
|
! endif
|
|
print*,'E = ',energies(1) + nuclear_repulsion
|
|
|
|
enddo
|
|
|
|
|
|
end
|
|
|
|
|
|
subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
|
double precision, intent(out) :: energies(N_st)
|
|
double precision, intent(out) :: diag_H_elements(dim_in)
|
|
double precision :: extra_diag_H_elements(dim_in)
|
|
double precision, intent(in) :: convergence
|
|
integer :: i,j,iter
|
|
DIAG_H_ELEMENTS = 0.d0
|
|
do iter = 1, 1
|
|
! call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
|
|
call dressing_1h1p_by_2h2p(dets_in,u_in,extra_diag_H_elements,dim_in,sze,N_st,Nint,convergence)
|
|
! if(sze<=N_det_max_jacobi)then
|
|
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
|
|
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze))
|
|
do j=1,sze
|
|
do i=1,sze
|
|
H_matrix_tmp(i,j) = H_matrix_all_dets(i,j)
|
|
enddo
|
|
enddo
|
|
H_matrix_tmp(1,1) += extra_diag_H_elements(1)
|
|
do i = 2,sze
|
|
H_matrix_tmp(1,i) += extra_diag_H_elements(i)
|
|
H_matrix_tmp(i,1) += extra_diag_H_elements(i)
|
|
enddo
|
|
!do i = 1,sze
|
|
! H_matrix_tmp(i,i) = diag_H_elements(i)
|
|
!enddo
|
|
call lapack_diag(eigenvalues,eigenvectors, &
|
|
H_matrix_tmp,size(H_matrix_all_dets,1),sze)
|
|
do j=1,min(N_states_diag,sze)
|
|
do i=1,sze
|
|
u_in(i,j) = eigenvectors(i,j)
|
|
enddo
|
|
energies(j) = eigenvalues(j)
|
|
enddo
|
|
deallocate (H_matrix_tmp, eigenvalues, eigenvectors)
|
|
! else
|
|
! call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants)
|
|
! endif
|
|
print*,'E = ',energies(1) + nuclear_repulsion
|
|
|
|
enddo
|
|
|
|
|
|
end
|
|
|
|
|
|
subroutine density_matrix_1h1p(dets_in,u_in,density_matrix_alpha,density_matrix_beta,norm,dim_in,sze,N_st,Nint)
|
|
use bitmasks
|
|
implicit none
|
|
BEGIN_DOC
|
|
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
|
|
!
|
|
! dets_in : bitmasks corresponding to determinants
|
|
!
|
|
! u_in : guess coefficients on the various states. Overwritten
|
|
! on exit
|
|
!
|
|
! dim_in : leftmost dimension of u_in
|
|
!
|
|
! sze : Number of determinants
|
|
!
|
|
! N_st : Number of eigenstates
|
|
!
|
|
! Initial guess vectors are not necessarily orthonormal
|
|
END_DOC
|
|
integer, intent(in) :: dim_in, sze, N_st, Nint
|
|
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
|
|
double precision, intent(inout) :: u_in(dim_in,N_st)
|
|
double precision, intent(inout) :: density_matrix_alpha(mo_tot_num_align,mo_tot_num)
|
|
double precision, intent(inout) :: density_matrix_beta(mo_tot_num_align,mo_tot_num)
|
|
double precision, intent(inout) :: norm
|
|
|
|
integer :: i,j,k,l
|
|
integer :: n_singles
|
|
integer :: index_singles(sze),hole_particles_singles(sze,3)
|
|
integer :: n_doubles
|
|
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
|
|
integer :: index_hf
|
|
|
|
integer :: exc(0:2,2,2),degree
|
|
integer :: h1,h2,p1,p2,s1,s2
|
|
integer :: other_spin(2)
|
|
double precision :: phase
|
|
integer(bit_kind) :: key_tmp(N_int,2)
|
|
integer :: i_ok
|
|
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral
|
|
double precision :: hij,c_ref,contrib
|
|
integer :: iorb
|
|
|
|
other_spin(1) = 2
|
|
other_spin(2) = 1
|
|
|
|
n_singles = 0
|
|
n_doubles = 0
|
|
norm = 0.d0
|
|
do i = 1,sze
|
|
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
|
|
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
|
norm += u_in(i,1)* u_in(i,1)
|
|
if(degree == 0)then
|
|
index_hf = i
|
|
c_ref = 1.d0/psi_coef(i,1)
|
|
else if (degree == 1)then
|
|
n_singles +=1
|
|
index_singles(n_singles) = i
|
|
! h1 = inactive orbital of the hole
|
|
hole_particles_singles(n_singles,1) = h1
|
|
! p1 = virtual orbital of the particle
|
|
hole_particles_singles(n_singles,2) = p1
|
|
! s1 = spin of the electron excited
|
|
hole_particles_singles(n_singles,3) = s1
|
|
else if (degree == 2)then
|
|
n_doubles +=1
|
|
index_doubles(n_doubles) = i
|
|
! h1 = inactive orbital of the hole (beta of course)
|
|
hole_particles_doubles(n_doubles,1) = h1
|
|
! p1 = virtual orbital of the particle (alpha of course)
|
|
hole_particles_doubles(n_doubles,2) = p2
|
|
else
|
|
print*,'PB !! found out other thing than a single or double'
|
|
print*,'stopping ..'
|
|
stop
|
|
endif
|
|
enddo
|
|
print*,'norm = ',norm
|
|
|
|
! Taking into account the connected part of the 2h2p on the HF determinant
|
|
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
|
|
|
|
do i = 1, n_singles
|
|
! start on the single excitation "|i>"
|
|
h1 = hole_particles_singles(i,1)
|
|
p1 = hole_particles_singles(i,2)
|
|
do j = 1, n_singles
|
|
do k = 1, N_int
|
|
key_tmp(k,1) = dets_in(k,1,index_singles(i))
|
|
key_tmp(k,2) = dets_in(k,2,index_singles(i))
|
|
enddo
|
|
h2 = hole_particles_singles(j,1)
|
|
p2 = hole_particles_singles(j,2)
|
|
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
|
|
! apply the excitation operator from the single excitation "|j>"
|
|
if(i_ok .ne. 1)cycle
|
|
double precision :: coef_ijrs,phase_other_single_ref
|
|
integer :: occ(N_int*bit_kind_size,2),n_occ(2)
|
|
call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_single_double,N_int)
|
|
call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
|
|
call get_excitation(key_tmp,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
|
|
coef_ijrs = u_in(index_singles(i),1) * u_in(index_singles(j),1) * c_ref * c_ref &
|
|
* phase_single_double * phase_other_single_ref
|
|
call bitstring_to_list_ab(key_tmp, occ, n_occ, N_int)
|
|
do k=1,elec_alpha_num
|
|
l = occ(k,1)
|
|
density_matrix_alpha(l,l) += coef_ijrs*coef_ijrs
|
|
enddo
|
|
do k=1,elec_beta_num
|
|
l = occ(k,1)
|
|
density_matrix_beta(l,l) += coef_ijrs*coef_ijrs
|
|
enddo
|
|
norm += coef_ijrs* coef_ijrs
|
|
if(hole_particles_singles(j,3) == 1)then ! single alpha
|
|
density_matrix_alpha(h2,p2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
|
|
density_matrix_alpha(p2,h2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
|
|
else
|
|
density_matrix_beta(h2,p2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
|
|
density_matrix_beta(p2,h2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
|
|
do i = 1, n_doubles
|
|
! start on the double excitation "|i>"
|
|
h1 = hole_particles_doubles(i,1)
|
|
p1 = hole_particles_doubles(i,2)
|
|
do j = 1, n_singles
|
|
do k = 1, N_int
|
|
key_tmp(k,1) = dets_in(k,1,index_doubles(i))
|
|
key_tmp(k,2) = dets_in(k,2,index_doubles(i))
|
|
enddo
|
|
h2 = hole_particles_singles(j,1)
|
|
p2 = hole_particles_singles(j,2)
|
|
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
|
|
! apply the excitation operator from the single excitation "|j>"
|
|
if(i_ok .ne. 1)cycle
|
|
double precision :: coef_ijrs_kv,phase_double_triple
|
|
call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_double_triple,N_int)
|
|
call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
|
|
call get_excitation(key_tmp,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
|
|
coef_ijrs_kv = u_in(index_doubles(i),1) * u_in(index_singles(j),1) * c_ref * c_ref &
|
|
* phase_double_triple * phase_other_single_ref
|
|
call bitstring_to_list_ab(key_tmp, occ, n_occ, N_int)
|
|
do k=1,elec_alpha_num
|
|
l = occ(k,1)
|
|
density_matrix_alpha(l,l) += coef_ijrs_kv*coef_ijrs_kv
|
|
enddo
|
|
do k=1,elec_beta_num
|
|
l = occ(k,1)
|
|
density_matrix_beta(l,l) += coef_ijrs_kv*coef_ijrs_kv
|
|
enddo
|
|
norm += coef_ijrs_kv* coef_ijrs_kv
|
|
if(hole_particles_singles(j,3) == 1)then ! single alpha
|
|
density_matrix_alpha(h2,p2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
|
|
density_matrix_alpha(p2,h2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
|
|
else
|
|
density_matrix_beta(h2,p2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
|
|
density_matrix_beta(p2,h2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
|
|
|
|
|
|
print*,'norm = ',norm
|
|
norm = 1.d0/norm
|
|
do i = 1, mo_tot_num
|
|
do j = 1, mo_tot_num
|
|
density_matrix_alpha(i,j) *= norm
|
|
density_matrix_beta(i,j) *= norm
|
|
enddo
|
|
enddo
|
|
coef_ijrs = 0.d0
|
|
do i = 1, mo_tot_num
|
|
coef_ijrs += density_matrix_beta(i,i) + density_matrix_beta(i,i)
|
|
enddo
|
|
print*,'accu = ',coef_ijrs
|
|
|
|
end
|
|
|