diff --git a/configure b/configure index 7370d540..29c99d24 100755 --- a/configure +++ b/configure @@ -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() diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index fb10ff2f..22464605 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -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 diff --git a/plugins/MP2/H_apply.irp.f b/plugins/MP2/H_apply.irp.f index 2f15391f..a79e3879 100644 --- a/plugins/MP2/H_apply.irp.f +++ b/plugins/MP2/H_apply.irp.f @@ -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 diff --git a/plugins/MP2/mp2_wf.irp.f b/plugins/MP2/mp2_wf.irp.f new file mode 100644 index 00000000..ad068b8a --- /dev/null +++ b/plugins/MP2/mp2_wf.irp.f @@ -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 diff --git a/plugins/Perturbation/Moller_plesset.irp.f b/plugins/Perturbation/Moller_plesset.irp.f index 7435f70c..38fe6610 100644 --- a/plugins/Perturbation/Moller_plesset.irp.f +++ b/plugins/Perturbation/Moller_plesset.irp.f @@ -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) diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 4e95e0f1..26f19d0d 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -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 diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index 4d1fc61d..c708511b 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -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 diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index f90eb0c8..7ee88e28 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -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 diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index f72b337c..5f087642 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -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 diff --git a/src/Determinants/filter_connected.irp.f b/src/Determinants/filter_connected.irp.f index cc304b59..005aa020 100644 --- a/src/Determinants/filter_connected.irp.f +++ b/src/Determinants/filter_connected.irp.f @@ -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 diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index b60bf0f7..46be7c18 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -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