LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-01 19:05:25 +02:00
Code Issues Releases Wiki Activity

Merge branch 'master' into develop

Browse Source
This commit is contained in:
Anthony Scemama 2016-03-30 01:21:48 +02:00
commit 7587b2bff9
13 changed files with 278 additions and 678 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -31,7 +31,7 @@ OPENMP : 1 ; Append OpenMP flags
# -ftz : Flushes denormal results to zero
#
[OPT]
FCFLAGS : -xHost -O2 -ip -ftz -g
FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g
# Profiling flags
#################

10
plugins/Full_CI/H_apply.irp.f
View File

@ -2,20 +2,20 @@ use bitmasks
BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import *
s = H_apply_zmq("FCI")
s = H_apply("FCI")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_openmp()
#s.unset_openmp()
print s
s = H_apply_zmq("FCI_PT2")
s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp()
#s.unset_openmp()
print s
s = H_apply_zmq("FCI_no_skip")
s = H_apply("FCI_no_skip")
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.unset_openmp()
#s.unset_openmp()
print s
s = H_apply("FCI_mono")

13
plugins/MRCC_CASSD/EZFIO.cfg
View File

@ -2,3 +2,16 @@
type: double precision
doc: Calculated energy
interface: ezfio
[thresh_mrcc]
type: Threshold
doc: Threshold on the convergence of the MRCC energy
interface: ezfio,provider,ocaml
default: 1.e-7
[n_it_mrcc_max]
type: Strictly_positive_int
doc: Maximum number of MRCC iterations
interface: ezfio,provider,ocaml
default: 20

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 *
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"] += """
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,*)
integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
double precision, intent(in) :: delta_ij_(Nstates, Ndet_non_ref, 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["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["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_, Nstates, Ndet_non_ref, Ndet_ref"
s.data["params_main"] += "delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_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,*)
integer, intent(in) :: Ndet_ref, Ndet_non_ref, Nstates
double precision, intent(in) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
double precision, intent(in) :: delta_ii_(Nstates,Ndet_ref)
"""
s.data["finalization"] = ""
s.data["copy_buffer"] = ""
@ -24,27 +24,5 @@ s.data["size_max"] = "3072"
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

254
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
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
implicit none
integer, intent(in) :: i_generator,n_selected, Nint, iproc
integer, intent(in) :: Ndet_ref, Ndet_non_ref
double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref
double precision, intent(inout) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref)
double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
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 :: 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 :: 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 :: h1,h2,p1,p2,s1,s2
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,allocatable :: idx_miniList(:)
integer :: N_miniList, ni, leng
double precision, allocatable :: hij_cache(:)
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)
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
@ -60,124 +61,157 @@ 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)
allocate (dIa_hla(N_states,Ndet_non_ref))
allocate (dIa_hla(Nstates,Ndet_non_ref))
! |I>
! |alpha>
if(N_tq > 0) then
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
end if
if(N_tq > 0) then
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
end if
do i_alpha=1,N_tq
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
do j=1,idx_alpha(0)
idx_alpha(j) = idx_miniList(idx_alpha(j))
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>
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
if (degree > 4) then
cycle
endif
do i_state=1,N_states
dIa(i_state) = 0.d0
enddo
! <I| <> |alpha>
do k_sd=1,idx_alpha(0)
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
cycle
endif
! <I| /k\ |alpha>
! <I|H|k>
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,N_states
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
enddo
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do k=1,N_int
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree_alpha(k_sd) == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
! <I| \l/ |alpha>
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
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)
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)
do i_state=1,N_states
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
enddo
exit
endif
enddo
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo
enddo
do i_state=1,N_states
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
enddo
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,hla)
do i_state=1,N_states
dIa_hla(i_state,k_sd) = dIa(i_state) * hla
enddo
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
do i_state=1,N_states
delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd)
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)
else
delta_ii_(i_I,i_state) = 0.d0
endif
enddo
enddo
call omp_unset_lock( psi_ref_lock(i_I) )
! Find triples and quadruple grand parents
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
if (degree > 4) then
cycle
endif
do i_state=1,Nstates
dIa(i_state) = 0.d0
enddo
! <I| <> |alpha>
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)
if (degree > 2) then
cycle
endif
! <I| /k\ |alpha>
! <I|H|k>
hIk = hij_mrcc(idx_alpha(k_sd),i_I)
! 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))
enddo
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do k=1,N_int
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree_alpha(k_sd) == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
! <I| \l/ |alpha>
do i_state=1,Nstates
dka(i_state) = 0.d0
enddo
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)
if (degree == 0) then
loop = .True.
do i_state=1,Nstates
if (lambda_mrcc(i_state,idx_alpha(l_sd)) /= 0.d0) then
loop = .False.
exit
endif
enddo
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)
enddo
enddo
do i_state=1,Nstates
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
enddo
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
hla = hij_cache(k_sd)
! 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
enddo
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
do i_state=1,Nstates
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
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
delta_ii_(i_state,i_I) = 0.d0
endif
enddo
enddo
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
deallocate (dIa_hla)
deallocate (dIa_hla,hij_cache)
deallocate(miniList, idx_miniList)
end

16
plugins/MRCC_Utils/mrcc_general.irp.f
View File

@ -31,23 +31,7 @@ subroutine mrcc_iterations
E_past(j) = E_new
j +=1
if(j>4)then
j=1
endif
if(iteration > 4) then
if(delta_E > 1.d-10)then
if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then
print*,'OSCILLATIONS !!!'
oscillations = .True.
i_oscillations +=1
lambda_mrcc_tmp = lambda_mrcc
endif
endif
endif
call save_wavefunction
! if (i_oscillations > 5) then
! exit
! endif
if (iteration > 200) then
exit
endif

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

@ -1,102 +1,35 @@
BEGIN_PROVIDER [integer, pert_determinants, (N_states, psi_det_size) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, lambda_pert, (N_states,psi_det_size) ]
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
implicit none
BEGIN_DOC
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
END_DOC
integer :: i,k,j
double precision :: ihpsi(N_states), hii,delta_e_eff,ihpsi_current(N_states),hij
integer :: i_ok,i_pert,i_pert_count
i_ok = 0
double precision :: phase_restart(N_states),tmp
do k = 1, N_states
phase_restart(k) = dsign(1.d0,psi_ref_coef_restart(1,k)/psi_ref_coef(1,k))
enddo
i_pert_count = 0
integer :: i,k
double precision :: ihpsi(N_states),ihpsi_current(N_states)
integer :: i_pert_count
do i=1,N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref_restart, psi_ref_coef_restart, N_int, N_det_ref,&
size(psi_ref_coef_restart,1), n_states, ihpsi)
call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
! TODO --- Test perturbatif ------
do k=1,N_states
lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
! TODO : i_h_psi peut sortir de la boucle?
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
if (ihpsi_current(k) == 0.d0) then
ihpsi_current(k) = 1.d-32
endif
tmp = psi_non_ref_coef(i,k)/ihpsi_current(k)
i_pert = 0
! Perturbation only if 1st order < 0.5 x second order
if((ihpsi(k) * lambda_pert(k,i)) < 0.5d0 * psi_non_ref_coef_restart(i,k) )then
i_pert = 1
else
do j = 1, N_det_ref
call i_H_j(psi_non_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
! Perturbation diverges when hij*tmp > 0.5
if(dabs(hij * tmp).ge.0.5d0)then
i_pert_count +=1
i_pert = 1
exit
endif
enddo
endif
if( i_pert == 1)then
pert_determinants(k,i) = i_pert
endif
if(pert_determinants(k,i) == 1)then
i_ok +=1
lambda_mrcc(k,i) = lambda_pert(k,i)
else
lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
endif
enddo
! TODO --- Fin test perturbatif ------
i_pert_count = 0
lambda_mrcc = 0.d0
do i=1,N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
do k=1,N_states
if (ihpsi_current(k) == 0.d0) then
ihpsi_current(k) = 1.d-32
endif
if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then
i_pert_count +=1
else
lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
endif
enddo
enddo
!if(oscillations)then
! print*,'AVERAGING the lambda_mrcc with those of the previous iterations'
! do i = 1, N_det_non_ref
! do k = 1, N_states
! double precision :: tmp
! tmp = lambda_mrcc(k,i)
! lambda_mrcc(k,i) += lambda_mrcc_tmp(k,i)
! lambda_mrcc(k,i) = lambda_mrcc(k,i) * 0.5d0
! if(dabs(tmp - lambda_mrcc(k,i)).ge.1.d-9)then
! print*,''
! print*,'i = ',i
! print*,'psi_non_ref_coef(i,k) = ',psi_non_ref_coef(i,k)
! print*,'lambda_mrcc(k,i) = ',lambda_mrcc(k,i)
! print*,' tmp = ',tmp
! endif
! enddo
! enddo
!endif
print*,'N_det_non_ref = ',N_det_non_ref
print*,'Number of Perturbatively treated determinants = ',i_ok
print*,'i_pert_count = ',i_pert_count
print*,'Number of ignored determinants = ',i_pert_count
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
END_PROVIDER
BEGIN_PROVIDER [ double precision, lambda_mrcc_tmp, (N_states,psi_det_size) ]
implicit none
lambda_mrcc_tmp = 0.d0
END_PROVIDER
BEGIN_PROVIDER [ logical, oscillations ]
implicit none
oscillations = .False.
END_PROVIDER
!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ]
@ -108,8 +41,22 @@ END_PROVIDER
!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref)
!END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
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
BEGIN_DOC
! Dressing matrix in N_det basis
@ -117,32 +64,7 @@ END_PROVIDER
integer :: i,j,m
delta_ij = 0.d0
delta_ii = 0.d0
call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_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
call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref)
END_PROVIDER
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
@ -159,11 +81,11 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
enddo
do ii = 1, N_det_ref
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
j =idx_non_ref(jj)
h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate)
h_matrix_dressed(j,i,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(istate,jj,ii)
enddo
enddo
enddo
@ -267,3 +189,4 @@ subroutine diagonalize_CI_dressed
SOFT_TOUCH psi_coef
end

260
plugins/MRCC_Utils/mrcepa_dress.irp.f
View File

@ -1,260 +0,0 @@
use omp_lib
use bitmasks
subroutine mrcepa_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint, iproc
integer, intent(in) :: Ndet_ref, Ndet_non_ref
double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l
integer :: degree_alpha(psi_det_size)
integer :: idx_alpha(0:psi_det_size)
logical :: good, fullMatch
integer(bit_kind) :: tq(Nint,2,n_selected)
integer :: N_tq, c_ref ,degree
double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states)
double precision, allocatable :: dIa_hia(:,:)
double precision :: haj, phase, phase2
double precision :: f(N_states), ci_inv(N_states)
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind) :: tmp_det(Nint,2)
integer(bit_kind) :: tmp_det_0(Nint,2)
integer :: iint, ipos
integer :: i_state, i_sd, k_sd, l_sd, i_I, i_alpha
integer(bit_kind),allocatable :: miniList(:,:,:)
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
integer,allocatable :: idx_miniList(:)
integer :: N_miniList, ni, leng
integer(bit_kind) :: isum
double precision :: hia
integer, allocatable :: index_sorted(:)
leng = max(N_det_generators, N_det_non_ref)
allocate(miniList(Nint, 2, leng), idx_miniList(leng), index_sorted(N_det))
!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)
if(fullMatch) then
return
end if
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
allocate (dIa_hia(N_states,Ndet_non_ref))
! |I>
! |alpha>
if(N_tq > 0) then
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
end if
do i_alpha=1,N_tq
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
integer, external :: get_index_in_psi_det_sorted_bit
index_sorted = huge(-1)
do j=1,idx_alpha(0)
idx_alpha(j) = idx_miniList(idx_alpha(j))
index_sorted( get_index_in_psi_det_sorted_bit( psi_non_ref(1,1,idx_alpha(j)), N_int ) ) = idx_alpha(j)
end do
! |I>
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
call get_excitation_degree(tq(1,1,i_alpha),psi_ref(1,1,i_I),degree,Nint)
if (degree > 4) then
cycle
endif
do i_state=1,N_states
dIa(i_state) = 0.d0
enddo
!TODO: MR
do i_sd=1,idx_alpha(0)
call get_excitation_degree(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),degree,Nint)
if (degree > 2) then
cycle
endif
call get_excitation(psi_non_ref(1,1,idx_alpha(i_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
tmp_det_0 = 0_bit_kind
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s1) = ibset(tmp_det_0(iint,s1),ipos)
if (degree == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det_0(iint,s2) = ibset(tmp_det_0(iint,s2),ipos)
endif
call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(i_sd)),Nint,hia)
! <I| <> |alpha>
do k_sd=1,idx_alpha(0)
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
cycle
endif
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
tmp_det = 0_bit_kind
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
isum = 0_bit_kind
do iint = 1,N_int
isum = isum + iand(tmp_det(iint,1), tmp_det_0(iint,1)) &
+ iand(tmp_det(iint,2), tmp_det_0(iint,2))
enddo
if (isum /= 0_bit_kind) then
cycle
endif
! <I| /k\ |alpha>
! <I|H|k>
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,N_states
dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
enddo
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
do k=1,N_int
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
! <I| \l/ |alpha>
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
! l_sd = index_sorted( get_index_in_psi_det_sorted_bit( tmp_det, N_int ) )
! call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,Nint)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),Nint,hIl)
! do i_state=1,N_states
! dka(i_state) = hIl * lambda_mrcc(i_state,l_sd) * phase * phase2
! enddo
do l_sd=1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
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)
do i_state=1,N_states
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
enddo
exit
endif
enddo
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo
enddo
do i_state=1,N_states
ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
enddo
k_sd = idx_alpha(i_sd)
do i_state=1,N_states
dIa_hia(i_state,k_sd) = dIa(i_state) * hia
enddo
call omp_set_lock( psi_ref_lock(i_I) )
do i_state=1,N_states
delta_ij_(i_I,k_sd,i_state) += dIa_hia(i_state,k_sd)
if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
delta_ii_(i_I,i_state) -= dIa_hia(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef(k_sd,i_state)
else
delta_ii_(i_I,i_state) = 0.d0
endif
enddo
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
enddo
deallocate (dIa_hia,index_sorted)
deallocate(miniList, idx_miniList)
end

97
plugins/MRCC_Utils/mrcepa_general.irp.f
View File

@ -1,97 +0,0 @@
subroutine run_mrcepa
implicit none
call set_generators_bitmasks_as_holes_and_particles
call mrcepa_iterations
end
subroutine mrcepa_iterations
implicit none
integer :: i,j
double precision :: E_new, E_old, delta_e
integer :: iteration,i_oscillations
double precision :: E_past(4)
E_new = 0.d0
delta_E = 1.d0
iteration = 0
j = 1
i_oscillations = 0
do while (delta_E > 1.d-7)
iteration += 1
print *, '==========================='
print *, 'MRCEPA Iteration', iteration
print *, '==========================='
print *, ''
E_old = sum(ci_energy_dressed)
call write_double(6,ci_energy_dressed(1),"MRCEPA energy")
call diagonalize_ci_dressed
E_new = sum(ci_energy_dressed)
delta_E = dabs(E_new - E_old)
E_past(j) = E_new
j +=1
if(j>4)then
j=1
endif
if(iteration > 4) then
if(delta_E > 1.d-10)then
if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then
print*,'OSCILLATIONS !!!'
oscillations = .True.
i_oscillations +=1
lambda_mrcc_tmp = lambda_mrcc
endif
endif
endif
call save_wavefunction
! if (i_oscillations > 5) then
! exit
! endif
if (iteration > 200) then
exit
endif
print*,'------------'
print*,'VECTOR'
do i = 1, N_det_ref
print*,''
print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
print*,'delta_ii(i,1) = ',delta_ii(i,1)
enddo
print*,'------------'
enddo
call write_double(6,ci_energy_dressed(1),"Final MRCEPA energy")
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
call save_wavefunction
end
subroutine set_generators_bitmasks_as_holes_and_particles
implicit none
integer :: i,k
do k = 1, N_generators_bitmask
do i = 1, N_int
! Pure single part
generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha
generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha
generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta
generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta
! Double excitation
generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha
generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha
generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta
generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta
generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha
generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha
generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta
generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta
enddo
enddo
touch generators_bitmask
end

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

@ -14,6 +14,31 @@ use bitmasks
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 [ 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 :: degree
double precision :: get_mo_bielec_integral_schwartz
double precision :: get_mo_bielec_integral
integer :: m,n,p,q
integer :: i,j,k
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)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral_schwartz( &
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_schwartz( &
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_schwartz( &
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_schwartz( &
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_schwartz( &
get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,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
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
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, elec_beta_num
i = occ(k,2)
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.
endif
enddo
@ -537,15 +537,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
do k = 1, elec_beta_num
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
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, elec_alpha_num
i = occ(k,1)
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.
endif
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) :: degree
double precision :: get_mo_bielec_integral_schwartz
double precision :: get_mo_bielec_integral
integer :: m,n,p,q
integer :: i,j,k
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)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral_schwartz( &
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_schwartz( &
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_schwartz( &
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_schwartz( &
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_schwartz( &
get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,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
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
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, elec_beta_num
i = occ(k,2)
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.
endif
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
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
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, elec_alpha_num
i = occ(k,1)
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.
endif
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 :: degree
double precision :: get_mo_bielec_integral_schwartz
double precision :: get_mo_bielec_integral
integer :: m,n,p,q
integer :: i,j,k
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)
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
hij = phase*get_mo_bielec_integral_schwartz( &
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_schwartz( &
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_schwartz( &
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_schwartz( &
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_schwartz( &
get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,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
i = occ(k,1)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
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, elec_beta_num
i = occ(k,2)
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.
endif
enddo
@ -811,15 +811,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
do k = 1, elec_beta_num
i = occ(k,2)
if (.not.has_mipi(i)) then
mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map)
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
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, elec_alpha_num
i = occ(k,1)
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.
endif
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
END_DOC
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
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
end

2
tests/bats/qp.bats
View File

@ -155,7 +155,7 @@ function run_all_1h_1p() {
ezfio set determinants read_wf True
qp_run mrcc_cassd $INPUT
energy="$(ezfio get mrcc_cassd energy)"
eq $energy -0.762303253805911E+02 1.E-3
eq $energy -76.2289109271715 1.E-3
}