LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-21 20:52:18 +02:00
Code Issues Releases Wiki Activity

Merge pull request #4 from LCPQ/master

Browse Source 
merge with LCPQ
This commit is contained in:
garniron 2015-11-23 09:16:53 +01:00
commit 4f501c0a98
11 changed files with 211 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
configure vendored
View File

@ -33,8 +33,8 @@ import sys
from os.path import join
if not any(i in ["--production", "--development"] for i in sys.argv):
print __doc__
sys.exit()
sys.argv += ["--development"]
if len(sys.argv) != 3:
print __doc__
sys.exit()

52
ocaml/qp_create_ezfio_from_xyz.ml
View File

@ -100,26 +100,21 @@ let run ?o b c d m p xyz_file =
match Hashtbl.find basis_table key with
| Some in_channel ->
in_channel
| None ->
begin
Printf.printf "%s is not defined in basis %s.\nEnter alternate basis : %!"
(Element.to_long_string element) b ;
let bas =
match In_channel.input_line stdin with
| Some line -> String.strip line |> String.lowercase
| None -> failwith "Aborted"
in
let new_channel = In_channel.create
(Qpackage.root ^ "/data/basis/" ^ bas)
in
Hashtbl.add_exn basis_table ~key:key ~data:new_channel;
new_channel
end
| None ->
let msg =
Printf.sprintf "%s is not defined in basis %s.%!"
(Element.to_long_string element) b ;
in
failwith msg
in
let temp_filename =
Filename.temp_file "qp_create_" ".basis"
in
let () =
Sys.remove temp_filename
in
let rec build_basis = function
| [] -> ()
| elem_and_basis_name :: rest ->
@ -130,10 +125,10 @@ let run ?o b c d m p xyz_file =
let command =
if (p) then
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "\" \"" ^ basis ^"\" pseudo"
^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo"
else
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "\" \"" ^ basis ^"\""
^ "." ^ basis ^ "\" \"" ^ basis ^"\""
in
begin
let filename =
@ -163,8 +158,8 @@ let run ?o b c d m p xyz_file =
Element.to_string elem
in
let command =
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^
"\" \"" ^ basis ^ "\" " ^ key
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename
^ "." ^ basis ^ "\" \"" ^ basis ^ "\" "
in
begin
let filename =
@ -232,11 +227,15 @@ let run ?o b c d m p xyz_file =
(* Write Basis set *)
let basis =
let nmax = Nucl_number.get_max () in
let nmax =
Nucl_number.get_max ()
in
let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
| [] -> accu
| e::tail ->
let new_accu = (e,(Nucl_number.of_int ~max:nmax n))::accu in
let new_accu =
(e,(Nucl_number.of_int ~max:nmax n))::accu
in
do_work new_accu (n+1) tail
in
let result = do_work [] 1 nuclei
@ -250,15 +249,8 @@ let run ?o b c d m p xyz_file =
in
Basis.read_element (basis_channel x.Atom.element) i e
with
| End_of_file ->
begin
let alt_channel = basis_channel x.Atom.element in
try
Basis.read_element alt_channel i x.Atom.element
with
End_of_file -> failwith
("Element "^(Element.to_string x.Atom.element)^" not found")
end
| End_of_file -> failwith
("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
)
|> List.concat
in

4
plugins/MP2/H_apply.irp.f
View File

@ -6,5 +6,9 @@ from perturbation import perturbations
s = H_apply("mp2")
s.set_perturbation("Moller_plesset")
print s
s = H_apply("mp2_selection")
s.set_selection_pt2("Moller_plesset")
print s
END_SHELL

31
plugins/MP2/mp2_wf.irp.f Normal file
View File

@ -0,0 +1,31 @@
program mp2_wf
implicit none
BEGIN_DOC
! Save the MP2 wave function
END_DOC
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, iter
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
pt2 = 1.d0
selection_criterion = 1.e-12
selection_criterion_min = 1.e-12
TOUCH selection_criterion_min selection_criterion selection_criterion_factor
call H_apply_mp2_selection(pt2, norm_pert, H_pert_diag, N_st)
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
touch N_det psi_det psi_coef
print*,'N_det = ',N_det
print*,'-----'
print *, 'PT2 = ', pt2(1)
print *, 'E = ', HF_energy
print *, 'E_before +PT2 = ', HF_energy+pt2(1)
N_det = min(N_det,N_det_max)
call save_wavefunction
call ezfio_set_mp2_energy(HF_energy+pt2(1))
deallocate(pt2,norm_pert,H_pert_diag)
end

16
plugins/Perturbation/Moller_plesset.irp.f
View File

@ -26,10 +26,18 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
if (degree == 2) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
else if (degree == 1) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)
delta_e = 1.d0/delta_e
else
delta_e = 0.d0
endif
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)

6
scripts/generate_h_apply.py
View File

@ -60,7 +60,7 @@ class H_apply(object):
s["omp_master"] = "!$OMP MASTER"
s["omp_end_master"] = "!$OMP END MASTER"
s["omp_barrier"] = "!$OMP BARRIER"
s["omp_do"] = "!$OMP DO SCHEDULE (static)"
s["omp_do"] = "!$OMP DO SCHEDULE (static,1)"
s["omp_enddo"] = "!$OMP ENDDO NOWAIT"
d = { True : '.True.', False : '.False.'}
@ -201,7 +201,7 @@ class H_apply(object):
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
"""
self.data["keys_work"] = """
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
@ -219,7 +219,7 @@ class H_apply(object):
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag(N_st)
double precision :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
PROVIDE CI_electronic_energy N_det_generators
PROVIDE N_det_generators
do k=1,N_st
pt2(k) = 0.d0
norm_pert(k) = 0.d0

2
scripts/get_basis.sh
View File

@ -8,6 +8,8 @@
# Prints in stdout the name of a temporary file containing the basis set.
#
#DEBUG:
#echo $0 $@ 1>&2
if [[ -z ${QP_ROOT} ]]
then

21
src/Determinants/H_apply.template.f
View File

@ -10,6 +10,9 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
integer :: highest, p1,p2,sp,ni,i,mi,nt,ns
integer(bit_kind), intent(in) :: key_prev(N_int, 2, *)
PROVIDE N_int
PROVIDE N_det
$declarations
@ -183,10 +186,8 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1))
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole(1,1),occ_hole(1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole(1,2),occ_hole(1,2),N_elec_in_key_hole_1(2),N_int)
call bitstring_to_list_ab(particle,occ_particle,N_elec_in_key_part_1,N_int)
call bitstring_to_list_ab(hole,occ_hole,N_elec_in_key_hole_1,N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2), &
ib_jb_pairs(2,0:(elec_alpha_num)*mo_tot_num))
@ -249,10 +250,8 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
particle_tmp(j,2) = iand(xor(particl_2(j,2),hole(j,2)),particl_2(j,2))
enddo
call bitstring_to_list(particle_tmp(1,1),occ_particle_tmp(1,1),N_elec_in_key_part_2(1),N_int)
call bitstring_to_list(particle_tmp(1,2),occ_particle_tmp(1,2),N_elec_in_key_part_2(2),N_int)
call bitstring_to_list(hole_tmp (1,1),occ_hole_tmp (1,1),N_elec_in_key_hole_2(1),N_int)
call bitstring_to_list(hole_tmp (1,2),occ_hole_tmp (1,2),N_elec_in_key_hole_2(2),N_int)
call bitstring_to_list_ab(particle_tmp,occ_particle_tmp,N_elec_in_key_part_2,N_int)
call bitstring_to_list_ab(hole_tmp,occ_hole_tmp,N_elec_in_key_hole_2,N_int)
! hole = a^(+)_j_a(ispin) a_i_a(ispin)|key_in> : mono exc :: orb(i_a,ispin) --> orb(j_a,ispin)
hole_save = hole
@ -444,10 +443,8 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int)
call bitstring_to_list_ab(particle,occ_particle,N_elec_in_key_part_1,N_int)
call bitstring_to_list_ab(hole,occ_hole,N_elec_in_key_hole_1,N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
do ispin=1,2

2
src/Determinants/density_matrix.irp.f
View File

@ -71,7 +71,6 @@
one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b
!$OMP END CRITICAL
deallocate(tmp_a,tmp_b)
!$OMP BARRIER
!$OMP END PARALLEL
endif
@ -157,7 +156,6 @@ END_PROVIDER
one_body_single_double_dm_mo_beta = one_body_single_double_dm_mo_beta + tmp_b
!$OMP END CRITICAL
deallocate(tmp_a,tmp_b)
!$OMP BARRIER
!$OMP END PARALLEL
END_PROVIDER

29
src/Determinants/filter_connected.irp.f
View File

@ -132,7 +132,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
else if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
@ -148,7 +148,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
else if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
@ -166,7 +166,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
else if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
@ -175,23 +175,28 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
outer: do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do m=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) +&
popcnt(xor( key1(m,2,i), key2(m,2)))
if (degree_x2 > 4) then
exit
if ( key1(m,1,i) /= key2(m,1)) then
degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1)))
if (degree_x2 > 4) then
cycle outer
endif
endif
if ( key1(m,2,i) /= key2(m,2)) then
degree_x2 = degree_x2+ popcnt(xor( key1(m,2,i), key2(m,2)))
if (degree_x2 > 4) then
cycle outer
endif
endif
enddo
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
enddo
enddo outer
endif
idx(0) = l-1

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

@ -349,6 +349,80 @@ subroutine get_mono_excitation(det1,det2,exc,phase,Nint)
enddo
end
subroutine bitstring_to_list_ab( string, list, n_elements, Nint)
use bitmasks
implicit none
BEGIN_DOC
! Gives the inidices(+1) of the bits set to 1 in the bit string
! For alpha/beta determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: string(Nint,2)
integer, intent(out) :: list(Nint*bit_kind_size,2)
integer, intent(out) :: n_elements(2)
integer :: i, j, ishift
integer(bit_kind) :: l
n_elements(1) = 0
n_elements(2) = 0
ishift = 1
do i=1,Nint
l = string(i,1)
do while (l /= 0_bit_kind)
j = trailz(l)
n_elements(1) = n_elements(1)+1
l = ibclr(l,j)
list(n_elements(1),1) = ishift+j
enddo
l = string(i,2)
do while (l /= 0_bit_kind)
j = trailz(l)
n_elements(2) = n_elements(2)+1
l = ibclr(l,j)
list(n_elements(2),2) = ishift+j
enddo
ishift = ishift + bit_kind_size
enddo
end
subroutine bitstring_to_list_ab_old( string, list, n_elements, Nint)
use bitmasks
implicit none
BEGIN_DOC
! Gives the inidices(+1) of the bits set to 1 in the bit string
! For alpha/beta determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: string(Nint,2)
integer, intent(out) :: list(Nint*bit_kind_size,2)
integer, intent(out) :: n_elements(2)
integer :: i, ishift
integer(bit_kind) :: l
n_elements(1) = 0
n_elements(2) = 0
ishift = 2
do i=1,Nint
l = string(i,1)
do while (l /= 0_bit_kind)
n_elements(1) = n_elements(1)+1
list(n_elements(1),1) = ishift+popcnt(l-1_bit_kind) - popcnt(l)
l = iand(l,l-1_bit_kind)
enddo
l = string(i,2)
do while (l /= 0_bit_kind)
n_elements(2) = n_elements(2)+1
list(n_elements(2),2) = ishift+popcnt(l-1_bit_kind) - popcnt(l)
l = iand(l,l-1_bit_kind)
enddo
ishift = ishift + bit_kind_size
enddo
end
@ -370,7 +444,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_alpha, n_occ_beta
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
@ -422,8 +496,8 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
@ -506,7 +580,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem
integer :: n_occ_alpha, n_occ_beta
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
@ -558,8 +632,8 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
@ -642,7 +716,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_alpha, n_occ_beta
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
@ -696,8 +770,8 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
@ -1229,15 +1303,15 @@ double precision function diag_H_mat_elem(det_in,Nint)
endif
!call debug_det(det_in,Nint)
integer :: tmp
call bitstring_to_list(particle(1,1), occ_particle(1,1), tmp, Nint)
ASSERT (tmp == nexc(1))
call bitstring_to_list(particle(1,2), occ_particle(1,2), tmp, Nint)
ASSERT (tmp == nexc(2))
call bitstring_to_list(hole(1,1), occ_hole(1,1), tmp, Nint)
ASSERT (tmp == nexc(1))
call bitstring_to_list(hole(1,2), occ_hole(1,2), tmp, Nint)
ASSERT (tmp == nexc(2))
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
ASSERT (tmp(1) == nexc(1))
ASSERT (tmp(2) == nexc(2))
!DIR$ FORCEINLINE
call bitstring_to_list_ab(hole, occ_hole, tmp, Nint)
ASSERT (tmp(1) == nexc(1))
ASSERT (tmp(2) == nexc(2))
det_tmp = ref_bitmask
do ispin=1,2
@ -1266,6 +1340,7 @@ subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb)
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i
integer :: tmp(2)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
@ -1279,19 +1354,19 @@ subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb)
other_spin = iand(ispin,1)+1
!DIR$ FORCEINLINE
call get_occ_from_key(key,occ,Nint)
na -= 1
call bitstring_to_list_ab(key, occ, tmp, Nint)
na = na-1
hjj -= mo_mono_elec_integral(iorb,iorb)
hjj = hjj - mo_mono_elec_integral(iorb,iorb)
! Same spin
do i=1,na
hjj -= mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
hjj = hjj - mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
hjj -= mo_bielec_integral_jj(occ(i,other_spin),iorb)
hjj = hjj - mo_bielec_integral_jj(occ(i,other_spin),iorb)
enddo
end
@ -1317,13 +1392,11 @@ subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
integer :: tmp
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint)
ASSERT (tmp == elec_alpha_num)
!DIR$ FORCEINLINE
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
ASSERT (tmp == elec_beta_num)
call bitstring_to_list_ab(key, occ, tmp, Nint)
ASSERT (tmp(1) == elec_alpha_num)
ASSERT (tmp(2) == elec_beta_num)
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k > 0)
@ -1331,18 +1404,18 @@ subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)
key(k,ispin) = ibset(key(k,ispin),l)
other_spin = iand(ispin,1)+1
hjj += mo_mono_elec_integral(iorb,iorb)
hjj = hjj + mo_mono_elec_integral(iorb,iorb)
! Same spin
do i=1,na
hjj += mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
hjj = hjj + mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
hjj += mo_bielec_integral_jj(occ(i,other_spin),iorb)
hjj = hjj + mo_bielec_integral_jj(occ(i,other_spin),iorb)
enddo
na += 1
na = na+1
end
subroutine get_occ_from_key(key,occ,Nint)
@ -1354,10 +1427,10 @@ subroutine get_occ_from_key(key,occ,Nint)
integer(bit_kind), intent(in) :: key(Nint,2)
integer , intent(in) :: Nint
integer , intent(out) :: occ(Nint*bit_kind_size,2)
integer :: tmp
integer :: tmp(2)
call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint)
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
end