LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
Code Issues Releases Wiki Activity

MRCC acceleration

Browse Source
This commit is contained in:
Anthony Scemama 2016-03-29 23:18:26 +02:00
parent a6911ff85b
commit c504518542
6 changed files with 238 additions and 211 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
plugins/MRCC_Utils/H_apply.irp.f
View File

@ -3,19 +3,19 @@ BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import * from generate_h_apply import *
s = H_apply("mrcc") s = H_apply("mrcc")
s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref" s.data["parameters"] = ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["declarations"] += """ s.data["declarations"] += """
integer, intent(in) :: Ndet_ref,Ndet_non_ref integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) double precision, intent(in) :: delta_ij_(Nstates, Ndet_non_ref, Ndet_ref)
double precision, intent(in) :: delta_ii_(Ndet_ref,*) double precision, intent(in) :: delta_ii_(Nstates, Ndet_ref)
""" """
s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Nstates,Ndet_non_ref,Ndet_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" s.data["params_post"] += ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" s.data["params_main"] += "delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["decls_main"] += """ s.data["decls_main"] += """
integer, intent(in) :: Ndet_ref,Ndet_non_ref integer, intent(in) :: Ndet_ref, Ndet_non_ref, Nstates
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) double precision, intent(in) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
double precision, intent(in) :: delta_ii_(Ndet_ref,*) double precision, intent(in) :: delta_ii_(Nstates,Ndet_ref)
""" """
s.data["finalization"] = "" s.data["finalization"] = ""
s.data["copy_buffer"] = "" s.data["copy_buffer"] = ""
@ -24,27 +24,5 @@ s.data["size_max"] = "3072"
print s print s
s = H_apply("mrcepa")
s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref"
s.data["declarations"] += """
integer, intent(in) :: Ndet_ref,Ndet_non_ref
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
"""
s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
s.data["decls_main"] += """
integer, intent(in) :: Ndet_ref,Ndet_non_ref
double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(in) :: delta_ii_(Ndet_ref,*)
"""
s.data["finalization"] = ""
s.data["copy_buffer"] = ""
s.data["generate_psi_guess"] = ""
s.data["size_max"] = "3072"
# print s
END_SHELL END_SHELL

86
plugins/MRCC_Utils/mrcc_dress.irp.f
View File

@ -14,14 +14,14 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
END_PROVIDER END_PROVIDER
subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
use bitmasks use bitmasks
implicit none implicit none
integer, intent(in) :: i_generator,n_selected, Nint, iproc integer, intent(in) :: i_generator,n_selected, Nint, iproc
integer, intent(in) :: Ndet_ref, Ndet_non_ref integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) double precision, intent(inout) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
double precision, intent(inout) :: delta_ii_(Ndet_ref,*) double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l integer :: i,j,k,l
@ -32,10 +32,10 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
integer(bit_kind) :: tq(Nint,2,n_selected) integer(bit_kind) :: tq(Nint,2,n_selected)
integer :: N_tq, c_ref ,degree integer :: N_tq, c_ref ,degree
double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) double precision :: hIk, hla, hIl, dIk(Nstates), dka(Nstates), dIa(Nstates)
double precision, allocatable :: dIa_hla(:,:) double precision, allocatable :: dIa_hla(:,:)
double precision :: haj, phase, phase2 double precision :: haj, phase, phase2
double precision :: f(N_states), ci_inv(N_states) double precision :: f(Nstates), ci_inv(Nstates)
integer :: exc(0:2,2,2) integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2 integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind) :: tmp_det(Nint,2) integer(bit_kind) :: tmp_det(Nint,2)
@ -46,10 +46,11 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
integer(bit_kind),intent(in) :: key_mask(Nint, 2) integer(bit_kind),intent(in) :: key_mask(Nint, 2)
integer,allocatable :: idx_miniList(:) integer,allocatable :: idx_miniList(:)
integer :: N_miniList, ni, leng integer :: N_miniList, ni, leng
double precision, allocatable :: hij_cache(:)
leng = max(N_det_generators, N_det_non_ref) leng = max(N_det_generators, N_det_non_ref)
allocate(miniList(Nint, 2, leng), idx_miniList(leng)) allocate(miniList(Nint, 2, leng), idx_miniList(leng), hij_cache(N_det_non_ref))
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
@ -61,7 +62,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
allocate (dIa_hla(N_states,Ndet_non_ref)) allocate (dIa_hla(Nstates,Ndet_non_ref))
! |I> ! |I>
@ -80,6 +81,11 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
idx_alpha(j) = idx_miniList(idx_alpha(j)) idx_alpha(j) = idx_miniList(idx_alpha(j))
end do end do
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd))
enddo
! |I> ! |I>
do i_I=1,N_det_ref do i_I=1,N_det_ref
! Find triples and quadruple grand parents ! Find triples and quadruple grand parents
@ -88,20 +94,36 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
cycle cycle
endif endif
do i_state=1,N_states do i_state=1,Nstates
dIa(i_state) = 0.d0 dIa(i_state) = 0.d0
enddo enddo
! <I| <> |alpha> ! <I| <> |alpha>
do k_sd=1,idx_alpha(0) do k_sd=1,idx_alpha(0)
! Loop if lambda == 0
logical :: loop
loop = .True.
do i_state=1,Nstates
if (lambda_mrcc(i_state,idx_alpha(k_sd)) /= 0.d0) then
loop = .False.
exit
endif
enddo
if (loop) then
cycle
endif
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint) call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint)
if (degree > 2) then if (degree > 2) then
cycle cycle
endif endif
! <I| /k\ |alpha> ! <I| /k\ |alpha>
! <I|H|k> ! <I|H|k>
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk) hIk = hij_mrcc(idx_alpha(k_sd),i_I)
do i_state=1,N_states ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
do i_state=1,Nstates
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
enddo enddo
! |l> = Exc(k -> alpha) |I> ! |l> = Exc(k -> alpha) |I>
@ -133,51 +155,63 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
endif endif
! <I| \l/ |alpha> ! <I| \l/ |alpha>
do i_state=1,N_states do i_state=1,Nstates
dka(i_state) = 0.d0 dka(i_state) = 0.d0
enddo enddo
do l_sd=k_sd+1,idx_alpha(0) do l_sd=k_sd+1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
if (degree == 0) then if (degree == 0) then
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) loop = .True.
do i_state=1,N_states do i_state=1,Nstates
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then
enddo loop = .False.
exit exit
endif endif
enddo enddo
do i_state=1,N_states if (.not.loop) then
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
hIl = hij_mrcc(idx_alpha(l_sd),i_I)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
do i_state=1,Nstates
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
enddo
endif
exit
endif
enddo
do i_state=1,Nstates
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo enddo
enddo enddo
do i_state=1,N_states do i_state=1,Nstates
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state) ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
enddo enddo
do l_sd=1,idx_alpha(0) do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd) k_sd = idx_alpha(l_sd)
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla) hla = hij_cache(k_sd)
do i_state=1,N_states ! call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
do i_state=1,Nstates
dIa_hla(i_state,k_sd) = dIa(i_state) * hla dIa_hla(i_state,k_sd) = dIa(i_state) * hla
enddo enddo
enddo enddo
call omp_set_lock( psi_ref_lock(i_I) ) call omp_set_lock( psi_ref_lock(i_I) )
do l_sd=1,idx_alpha(0) do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd) k_sd = idx_alpha(l_sd)
do i_state=1,N_states do i_state=1,Nstates
delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd) delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state) delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
else else
delta_ii_(i_I,i_state) = 0.d0 delta_ii_(i_state,i_I) = 0.d0
endif endif
enddo enddo
enddo enddo
call omp_unset_lock( psi_ref_lock(i_I) ) call omp_unset_lock( psi_ref_lock(i_I) )
enddo enddo
enddo enddo
deallocate (dIa_hla) deallocate (dIa_hla,hij_cache)
deallocate(miniList, idx_miniList) deallocate(miniList, idx_miniList)
end end

52
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -41,8 +41,22 @@ END_PROVIDER
!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref) !call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref)
!END_PROVIDER !END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ] BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ] implicit none
BEGIN_DOC
! < ref | H | Non-ref > matrix
END_DOC
integer :: i_I, k_sd
do i_I=1,N_det_ref
do k_sd=1,N_det_non_ref
call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,k_sd),N_int,hij_mrcc(k_sd,i_I))
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states,N_det_ref) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Dressing matrix in N_det basis ! Dressing matrix in N_det basis
@ -50,32 +64,7 @@ END_PROVIDER
integer :: i,j,m integer :: i,j,m
delta_ij = 0.d0 delta_ij = 0.d0
delta_ii = 0.d0 delta_ii = 0.d0
call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref) call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref)
double precision :: max_delta
double precision :: accu
integer :: imax,jmax
max_delta = 0.d0
accu = 0.d0
do i = 1, N_det_ref
do j = 1, N_det_non_ref
accu += psi_non_ref_coef(j,1) * psi_ref_coef(i,1) * delta_ij(i,j,1)
if(dabs(delta_ij(i,j,1)).gt.max_delta)then
max_delta = dabs(delta_ij(i,j,1))
imax = i
jmax = j
endif
enddo
enddo
print*,''
print*,''
print*,'<psi| Delta H |psi> = ',accu
print*,'MAX VAL OF DRESING = ',delta_ij(imax,jmax,1)
print*,'imax,jmax = ',imax,jmax
print*,'psi_ref_coef(imax,1) = ',psi_ref_coef(imax,1)
print*,'psi_non_ref_coef(jmax,1) = ',psi_non_ref_coef(jmax,1)
do i = 1, N_det_ref
print*,'delta_ii(i,1) = ',delta_ii(i,1)
enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
@ -92,11 +81,11 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
enddo enddo
do ii = 1, N_det_ref do ii = 1, N_det_ref
i =idx_ref(ii) i =idx_ref(ii)
h_matrix_dressed(i,i,istate) += delta_ii(ii,istate) h_matrix_dressed(i,i,istate) += delta_ii(istate,ii)
do jj = 1, N_det_non_ref do jj = 1, N_det_non_ref
j =idx_non_ref(jj) j =idx_non_ref(jj)
h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate) h_matrix_dressed(i,j,istate) += delta_ij(istate,jj,ii)
h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate) h_matrix_dressed(j,i,istate) += delta_ij(istate,jj,ii)
enddo enddo
enddo enddo
enddo enddo
@ -200,3 +189,4 @@ subroutine diagonalize_CI_dressed
SOFT_TOUCH psi_coef SOFT_TOUCH psi_coef
end end

25
plugins/Psiref_Utils/psi_ref_utils.irp.f
View File

@ -14,6 +14,31 @@ use bitmasks
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_ref_coef_transp, (n_states,psi_det_size) ]
implicit none
BEGIN_DOC
! Transposed psi_ref_coef
END_DOC
integer :: i,j
do j=1,N_det_ref
do i=1, n_states
psi_ref_coef_transp(i,j) = psi_ref_coef(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_non_ref_coef_transp, (n_states,psi_det_size) ]
implicit none
BEGIN_DOC
! Transposed psi_non_ref_coef
END_DOC
integer :: i,j
do j=1,N_det_non_ref
do i=1, n_states
psi_non_ref_coef_transp(i,j) = psi_non_ref_coef(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ] BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ] &BEGIN_PROVIDER [ double precision, psi_non_ref_coef, (psi_det_size,n_states) ]

72
src/Determinants/slater_rules.irp.f
View File

@ -443,7 +443,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
integer :: exc(0:2,2,2) integer :: exc(0:2,2,2)
integer :: degree integer :: degree
double precision :: get_mo_bielec_integral_schwartz double precision :: get_mo_bielec_integral
integer :: m,n,p,q integer :: m,n,p,q
integer :: i,j,k integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2) integer :: occ(Nint*bit_kind_size,2)
@ -468,31 +468,31 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
call get_double_excitation(key_i,key_j,exc,phase,Nint) call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then if (exc(0,1,1) == 1) then
! Mono alpha, mono beta ! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral_schwartz( & hij = phase*get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(1,1,2), & exc(1,1,2), &
exc(1,2,1), & exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map) exc(1,2,2) ,mo_integrals_map)
else if (exc(0,1,1) == 2) then else if (exc(0,1,1) == 2) then
! Double alpha ! Double alpha
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
exc(1,2,1), & exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map) - & exc(2,2,1) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( & get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
exc(2,2,1), & exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) ) exc(1,2,1) ,mo_integrals_map) )
else if (exc(0,1,2) == 2) then else if (exc(0,1,2) == 2) then
! Double beta ! Double beta
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
exc(1,2,2), & exc(1,2,2), &
exc(2,2,2) ,mo_integrals_map) - & exc(2,2,2) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( & get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
exc(2,2,2), & exc(2,2,2), &
@ -510,15 +510,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
@ -537,15 +537,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
@ -579,7 +579,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
integer,intent(out) :: exc(0:2,2,2) integer,intent(out) :: exc(0:2,2,2)
integer,intent(out) :: degree integer,intent(out) :: degree
double precision :: get_mo_bielec_integral_schwartz double precision :: get_mo_bielec_integral
integer :: m,n,p,q integer :: m,n,p,q
integer :: i,j,k integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2) integer :: occ(Nint*bit_kind_size,2)
@ -604,31 +604,31 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
call get_double_excitation(key_i,key_j,exc,phase,Nint) call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then if (exc(0,1,1) == 1) then
! Mono alpha, mono beta ! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral_schwartz( & hij = phase*get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(1,1,2), & exc(1,1,2), &
exc(1,2,1), & exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map) exc(1,2,2) ,mo_integrals_map)
else if (exc(0,1,1) == 2) then else if (exc(0,1,1) == 2) then
! Double alpha ! Double alpha
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
exc(1,2,1), & exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map) - & exc(2,2,1) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( & get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
exc(2,2,1), & exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) ) exc(1,2,1) ,mo_integrals_map) )
else if (exc(0,1,2) == 2) then else if (exc(0,1,2) == 2) then
! Double beta ! Double beta
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
exc(1,2,2), & exc(1,2,2), &
exc(2,2,2) ,mo_integrals_map) - & exc(2,2,2) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( & get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
exc(2,2,2), & exc(2,2,2), &
@ -646,15 +646,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
@ -673,15 +673,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
@ -715,7 +715,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
integer :: exc(0:2,2,2) integer :: exc(0:2,2,2)
integer :: degree integer :: degree
double precision :: get_mo_bielec_integral_schwartz double precision :: get_mo_bielec_integral
integer :: m,n,p,q integer :: m,n,p,q
integer :: i,j,k integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2) integer :: occ(Nint*bit_kind_size,2)
@ -742,31 +742,31 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
call get_double_excitation(key_i,key_j,exc,phase,Nint) call get_double_excitation(key_i,key_j,exc,phase,Nint)
if (exc(0,1,1) == 1) then if (exc(0,1,1) == 1) then
! Mono alpha, mono beta ! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral_schwartz( & hij = phase*get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(1,1,2), & exc(1,1,2), &
exc(1,2,1), & exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map) exc(1,2,2) ,mo_integrals_map)
else if (exc(0,1,1) == 2) then else if (exc(0,1,1) == 2) then
! Double alpha ! Double alpha
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
exc(1,2,1), & exc(1,2,1), &
exc(2,2,1) ,mo_integrals_map) - & exc(2,2,1) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( & get_mo_bielec_integral( &
exc(1,1,1), & exc(1,1,1), &
exc(2,1,1), & exc(2,1,1), &
exc(2,2,1), & exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) ) exc(1,2,1) ,mo_integrals_map) )
else if (exc(0,1,2) == 2) then else if (exc(0,1,2) == 2) then
! Double beta ! Double beta
hij = phase*(get_mo_bielec_integral_schwartz( & hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
exc(1,2,2), & exc(1,2,2), &
exc(2,2,2) ,mo_integrals_map) - & exc(2,2,2) ,mo_integrals_map) - &
get_mo_bielec_integral_schwartz( & get_mo_bielec_integral( &
exc(1,1,2), & exc(1,1,2), &
exc(2,1,2), & exc(2,1,2), &
exc(2,2,2), & exc(2,2,2), &
@ -784,15 +784,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
@ -811,15 +811,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
do k = 1, elec_beta_num do k = 1, elec_beta_num
i = occ(k,2) i = occ(k,2)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo
do k = 1, elec_alpha_num do k = 1, elec_alpha_num
i = occ(k,1) i = occ(k,1)
if (.not.has_mipi(i)) then if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map)
has_mipi(i) = .True. has_mipi(i) = .True.
endif endif
enddo enddo

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

@ -324,9 +324,9 @@ double precision function mo_bielec_integral(i,j,k,l)
! Returns one integral <ij|kl> in the MO basis ! Returns one integral <ij|kl> in the MO basis
END_DOC END_DOC
integer, intent(in) :: i,j,k,l integer, intent(in) :: i,j,k,l
double precision :: get_mo_bielec_integral_schwartz double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map PROVIDE mo_bielec_integrals_in_map
mo_bielec_integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map) mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
return return
end end