diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 22464605..b2eaaae0 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -91,7 +91,10 @@ let run ?o b c d m p xyz_file = in - let basis_table = Hashtbl.Poly.create () in + let basis_table = + Hashtbl.Poly.create () + in + (* Open basis set channels *) let basis_channel element = let key = @@ -115,31 +118,46 @@ let run ?o b c d m p xyz_file = Sys.remove temp_filename in + let fetch_channel basis = + let command = + if (p) then + Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename + ^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo" + else + Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename + ^ "." ^ basis ^ "\" \"" ^ basis ^"\"" + in + match Sys.is_file basis with + | `Yes -> + In_channel.create basis + | _ -> + begin + let filename = + Unix.open_process_in command + |> In_channel.input_all + |> String.strip + in + let new_channel = + In_channel.create filename + in + Unix.unlink filename; + new_channel + end + in + let rec build_basis = function | [] -> () | elem_and_basis_name :: rest -> begin match (String.lsplit2 ~on:':' elem_and_basis_name) with | None -> (* Principal basis *) - let basis = elem_and_basis_name in - let command = - if (p) then - Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename - ^ "." ^ basis ^ "\" \"" ^ basis ^"\" pseudo" - else - Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename - ^ "." ^ basis ^ "\" \"" ^ basis ^"\"" - in begin - let filename = - Unix.open_process_in command - |> In_channel.input_all - |> String.strip + let basis = + elem_and_basis_name in let new_channel = - In_channel.create filename + fetch_channel basis in - Unix.unlink filename; List.iter nuclei ~f:(fun elem-> let key = Element.to_string elem.Atom.element @@ -151,26 +169,18 @@ let run ?o b c d m p xyz_file = end | Some (key, basis) -> (*Aux basis *) begin - let elem = Element.of_string key - and basis = String.lowercase basis + let elem = + Element.of_string key + and basis = + String.lowercase basis in let key = Element.to_string elem in - let command = - Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename - ^ "." ^ basis ^ "\" \"" ^ basis ^ "\" " + let new_channel = + fetch_channel basis in begin - let filename = - Unix.open_process_in command - |> In_channel.input_all - |> String.strip - in - let new_channel = - In_channel.create filename - in - Unix.unlink filename; match Hashtbl.add basis_table ~key:key ~data:new_channel with | `Ok -> () | `Duplicate -> failwith ("Duplicate definition of basis for "^(Element.to_long_string elem)) @@ -335,13 +345,15 @@ let command = Command.basic ~summary: "Quantum Package command" ~readme:(fun () -> " -Creates an EZFIO directory from a standard xyz file. -The basis set is defined as a single string if all the -atoms are taken from the same basis set, otherwise specific -elements can be defined as follows: +Creates an EZFIO directory from a standard xyz file. The basis set is defined +as a single string if all the atoms are taken from the same basis set, +otherwise specific elements can be defined as follows: -b \"cc-pcvdz | H:cc-pvdz | C:6-31g\" +If a file with the same name as the basis set exists, this file will be read. +Otherwise, the basis set is obtained from the database. + ") spec (fun o b c d m p xyz_file () -> diff --git a/plugins/Hartree_Fock/README.rst b/plugins/Hartree_Fock/README.rst index ffe80c75..26976226 100644 --- a/plugins/Hartree_Fock/README.rst +++ b/plugins/Hartree_Fock/README.rst @@ -32,11 +32,11 @@ Documentation .. by the `update_README.py` script. -`ao_bi_elec_integral_alpha `_ +`ao_bi_elec_integral_alpha `_ Alpha Fock matrix in AO basis set -`ao_bi_elec_integral_beta `_ +`ao_bi_elec_integral_beta `_ Alpha Fock matrix in AO basis set @@ -62,23 +62,23 @@ Documentation Diagonal Fock matrix in the MO basis -`fock_matrix_alpha_ao `_ +`fock_matrix_alpha_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_ao `_ +`fock_matrix_ao `_ Fock matrix in AO basis set -`fock_matrix_beta_ao `_ +`fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis @@ -114,7 +114,7 @@ Documentation .br -`fock_mo_to_ao `_ +`fock_mo_to_ao `_ Undocumented @@ -134,7 +134,7 @@ Documentation S^-1 Beta density matrix in the AO basis x S^-1 -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy @@ -142,7 +142,11 @@ Documentation Build the MOs using the extended Huckel model -`mo_guess_type `_ +`level_shift `_ + Energy shift on the virtual MOs to improve SCF convergence + + +`mo_guess_type `_ Initial MO guess. Can be [ Huckel | HCore ] @@ -161,6 +165,6 @@ Documentation optional: mo_basis.mo_coef -`thresh_scf `_ +`thresh_scf `_ Threshold on the convergence of the Hartree Fock energy diff --git a/plugins/MRCC_Utils/README.rst b/plugins/MRCC_Utils/README.rst index 7db472a6..8b97bfbe 100644 --- a/plugins/MRCC_Utils/README.rst +++ b/plugins/MRCC_Utils/README.rst @@ -90,10 +90,6 @@ Documentation N_states lowest eigenvalues of the dressed CI matrix -`create_minilist `_ - Undocumented - - `davidson_diag_hjj_mrcc `_ Davidson diagonalization with specific diagonal elements of the H matrix .br @@ -206,7 +202,7 @@ Documentation Find A.C = B -`find_triples_and_quadruples `_ +`find_triples_and_quadruples `_ Undocumented @@ -232,19 +228,19 @@ Documentation Undocumented -`gen_det_idx `_ +`gen_det_idx `_ Undocumented -`gen_det_shortcut `_ +`gen_det_shortcut `_ Undocumented -`gen_det_sorted `_ +`gen_det_sorted `_ Undocumented -`gen_det_version `_ +`gen_det_version `_ Undocumented @@ -312,7 +308,7 @@ h_apply_mrcc_monoexc Dressed H with Delta_ij -`h_u_0_mrcc `_ +`h_u_0_mrcc `_ Computes v_0 = H|u_0> .br n : number of determinants diff --git a/plugins/Perturbation/Moller_plesset.irp.f b/plugins/Perturbation/Moller_plesset.irp.f deleted file mode 100644 index 38fe6610..00000000 --- a/plugins/Perturbation/Moller_plesset.irp.f +++ /dev/null @@ -1,50 +0,0 @@ -subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,n_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - - BEGIN_DOC - ! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution - ! - ! for the various n_st states. - ! - ! c_pert(i) = /(difference of orbital energies) - ! - ! e_2_pert(i) = ^2/(difference of orbital energies) - ! - END_DOC - - integer :: i,j - double precision :: diag_H_mat_elem - integer :: exc(0:2,2,2) - integer :: degree - double precision :: phase,delta_e,h - integer :: h1,h2,p1,p2,s1,s2 - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint) - 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) - do i =1,n_st - H_pert_diag(i) = h - c_pert(i) = i_H_psi_array(i) *delta_e - e_2_pert(i) = c_pert(i) * i_H_psi_array(i) - enddo - -end diff --git a/plugins/Perturbation/README.rst b/plugins/Perturbation/README.rst index 55f2077e..4bf62a2a 100644 --- a/plugins/Perturbation/README.rst +++ b/plugins/Perturbation/README.rst @@ -239,7 +239,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet `_ +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states. @@ -250,7 +250,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various N_st states. @@ -261,7 +261,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet_sc2 `_ +`pt2_epstein_nesbet_sc2 `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, but with the CISD_SC2 energies and coefficients @@ -272,7 +272,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet_sc2_no_projected `_ +`pt2_epstein_nesbet_sc2_no_projected `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, @@ -296,7 +296,7 @@ perturb_buffer_moller_plesset H_pert_diag = c_pert -`pt2_epstein_nesbet_sc2_projected `_ +`pt2_epstein_nesbet_sc2_projected `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, @@ -336,7 +336,7 @@ perturb_buffer_moller_plesset than pt2_max in absolute value -`pt2_moller_plesset `_ +`pt2_moller_plesset `_ compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -352,7 +352,7 @@ perturb_buffer_moller_plesset provided. -`repeat_all_e_corr `_ +`repeat_all_e_corr `_ Undocumented diff --git a/plugins/Perturbation/epstein_nesbet.irp.f b/plugins/Perturbation/epstein_nesbet.irp.f deleted file mode 100644 index ede75df9..00000000 --- a/plugins/Perturbation/epstein_nesbet.irp.f +++ /dev/null @@ -1,106 +0,0 @@ -subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states. - ! - ! c_pert(i) = /( E(i) - ) - ! - ! e_2_pert(i) = ^2/( E(i) - ) - ! - END_DOC - - integer :: i,j - double precision :: diag_H_mat_elem, h - PROVIDE selection_criterion - - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) - - - h = diag_H_mat_elem(det_pert,Nint) - do i =1,N_st - if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then - c_pert(i) = -1.d0 - e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) - H_pert_diag(i) = h*c_pert(i)*c_pert(i) - e_2_pert(i) = c_pert(i) * i_H_psi_array(i) - else - c_pert(i) = -1.d0 - e_2_pert(i) = -dabs(i_H_psi_array(i)) - H_pert_diag(i) = h - endif - enddo - -end - -subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution - ! - ! for the various N_st states. - ! - ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) - ! - ! c_pert(i) = e_2_pert(i)/ - ! - END_DOC - - integer :: i,j - double precision :: diag_H_mat_elem,delta_e, h - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - PROVIDE CI_electronic_energy - - !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) - - h = diag_H_mat_elem(det_pert,Nint) - do i =1,N_st - if (i_H_psi_array(i) /= 0.d0) then - delta_e = h - CI_electronic_energy(i) - if (delta_e > 0.d0) then - e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) - else - e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) - endif - if (dabs(i_H_psi_array(i)) > 1.d-6) then - c_pert(i) = e_2_pert(i)/i_H_psi_array(i) - else - c_pert(i) = 0.d0 - endif - H_pert_diag(i) = h*c_pert(i)*c_pert(i) - else - e_2_pert(i) = 0.d0 - c_pert(i) = 0.d0 - H_pert_diag(i) = 0.d0 - endif - enddo - -end diff --git a/plugins/Perturbation/pert_sc2.irp.f b/plugins/Perturbation/pert_sc2.irp.f index 1ae6ce73..f225e757 100644 --- a/plugins/Perturbation/pert_sc2.irp.f +++ b/plugins/Perturbation/pert_sc2.irp.f @@ -1,166 +1,3 @@ - -subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(0:ndet) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states, - ! - ! but with the correction in the denominator - ! - ! comming from the interaction of that determinant with all the others determinants - ! - ! that can be repeated by repeating all the double excitations - ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) - ! - ! that could be repeated to this determinant. - ! - ! In addition, for the perturbative energetic contribution you have the standard second order - ! - ! e_2_pert = ^2/(Delta_E) - ! - ! and also the purely projected contribution - ! - ! H_pert_diag = c_pert - END_DOC - - integer :: i,j,degree,l - double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E - double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda - - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) - accu_e_corr = 0.d0 - !$IVDEP - do i = 1, idx_repeat(0) - accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i)) - enddo - h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) - - - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) - e_2_pert(1) = i_H_psi_array(1) * c_pert(1) - - do i =2,N_st - H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) - e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) - else - c_pert(i) = i_H_psi_array(i) - e_2_pert(i) = -dabs(i_H_psi_array(i)) - endif - enddo - - degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + & - popcnt(xor( ref_bitmask(1,2), det_pert(1,2))) - !DEC$ NOUNROLL - do l=2,Nint - degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + & - popcnt(xor( ref_bitmask(l,2), det_pert(l,2))) - enddo - if(degree==4)then - ! - e_2_pert_fonda = e_2_pert(1) - H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1) - do i = 1, N_st - do j = 1, idx_repeat(0) - e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i) - enddo - enddo - endif - -end - - -subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(0:ndet) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states, - ! - ! but with the correction in the denominator - ! - ! comming from the interaction of that determinant with all the others determinants - ! - ! that can be repeated by repeating all the double excitations - ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) - ! - ! that could be repeated to this determinant. - ! - ! In addition, for the perturbative energetic contribution you have the standard second order - ! - ! e_2_pert = ^2/(Delta_E) - ! - ! and also the purely projected contribution - ! - ! H_pert_diag = c_pert - END_DOC - - integer :: i,j,degree,l - double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E - double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda - - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) - accu_e_corr = 0.d0 - !$IVDEP - do i = 1, idx_repeat(0) - accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i)) - enddo - h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) - - - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) - e_2_pert(1) = i_H_psi_array(1) * c_pert(1) - - do i =2,N_st - H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) - e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) - else - c_pert(i) = i_H_psi_array(i) - e_2_pert(i) = -dabs(i_H_psi_array(i)) - endif - enddo -end - - - - - double precision function repeat_all_e_corr(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) @@ -190,54 +27,3 @@ double precision function repeat_all_e_corr(key_in) end - -subroutine pt2_epstein_nesbet_sc2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist) - use bitmasks - implicit none - integer, intent(in) :: Nint,ndet,N_st - integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) - double precision :: i_H_psi_array(N_st) - - integer, intent(in) :: N_minilist - integer, intent(in) :: idx_minilist(0:N_det_selectors) - integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) - - BEGIN_DOC - ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - ! - ! for the various N_st states, but with the CISD_SC2 energies and coefficients - ! - ! c_pert(i) = /( E(i) - ) - ! - ! e_2_pert(i) = ^2/( E(i) - ) - ! - END_DOC - - integer :: i,j - double precision :: diag_H_mat_elem, h - PROVIDE selection_criterion - - ASSERT (Nint == N_int) - ASSERT (Nint > 0) - !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) - - - h = diag_H_mat_elem(det_pert,Nint) - do i =1,N_st - if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then - c_pert(i) = -1.d0 - e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) - H_pert_diag(i) = h*c_pert(i)*c_pert(i) - e_2_pert(i) = c_pert(i) * i_H_psi_array(i) - else - c_pert(i) = -1.d0 - e_2_pert(i) = -dabs(i_H_psi_array(i)) - H_pert_diag(i) = h - endif - enddo - -end diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index 05176fe6..fa01cc99 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -2,7 +2,7 @@ BEGIN_SHELL [ /usr/bin/env python ] import perturbation END_SHELL -subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask) +subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -12,6 +12,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c integer, intent(in) :: Nint, N_st, buffer_size, i_generator integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) integer(bit_kind),intent(in) :: key_mask(Nint,2) + double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) @@ -43,24 +44,8 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c end if - buffer_loop : do i = 1,buffer_size + do i=1,buffer_size -! do k=1,N_minilist_gen -! ex = 0 -! do ni=1,Nint -! ex += popcnt(xor(minilist_gen(ni,1,k), buffer(ni,1,i))) + popcnt(xor(minilist_gen(ni,2,k), buffer(ni,2,i))) -! end do -! if(ex <= 4) then -! cycle buffer_loop -! end if -! end do - -! c_ref = connected_to_ref(buffer(1,1,i),miniList_gen,Nint,N_minilist_gen+1,N_minilist_gen) -! -! if (c_ref /= 0) then -! cycle -! endif - if(is_connected_to(buffer(1,1,i), miniList_gen, Nint, N_minilist_gen)) then cycle end if @@ -71,20 +56,18 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c integer :: degree call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int) -! call pt2_$PERT(buffer(1,1,i), & -! c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st,minilist,idx_minilist) - call pt2_$PERT(buffer(1,1,i), & - c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) !!!!!!!!!!!!!!!!! MAUVAISE SIGNATURE PR LES AUTRES PT2_* !!!!! + call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st - e_2_pert_buffer(k,i) = e_2_pert(k) - coef_pert_buffer(k,i) = c_pert(k) - sum_norm_pert(k) += c_pert(k) * c_pert(k) - sum_e_2_pert(k) += e_2_pert(k) - sum_H_pert_diag(k) += H_pert_diag(k) + e_2_pert_buffer(k,i) = e_2_pert(k) + coef_pert_buffer(k,i) = c_pert(k) + sum_norm_pert(k) = sum_norm_pert(k) + c_pert(k) * c_pert(k) + sum_e_2_pert(k) = sum_e_2_pert(k) + e_2_pert(k) + sum_H_pert_diag(k) = sum_H_pert_diag(k) + H_pert_diag(k) enddo - enddo buffer_loop + enddo end diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f new file mode 100644 index 00000000..f49ee2ff --- /dev/null +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -0,0 +1,367 @@ +BEGIN_TEMPLATE + +subroutine pt2_epstein_nesbet ($arguments) + use bitmasks + implicit none + $declarations + + BEGIN_DOC + ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states. + ! + ! c_pert(i) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem_fock, h + double precision :: i_H_psi_array(N_st) + PROVIDE selection_criterion + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + + + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + do i =1,N_st + if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then + c_pert(i) = -1.d0 + e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 + else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) + H_pert_diag(i) = h*c_pert(i)*c_pert(i) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + else + c_pert(i) = -1.d0 + e_2_pert(i) = -dabs(i_H_psi_array(i)) + H_pert_diag(i) = h + endif + enddo + +end + +subroutine pt2_epstein_nesbet_2x2 ($arguments) + use bitmasks + implicit none + $declarations + + BEGIN_DOC + ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + ! + ! for the various N_st states. + ! + ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + ! + ! c_pert(i) = e_2_pert(i)/ + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem_fock,delta_e, h + double precision :: i_H_psi_array(N_st) + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + PROVIDE CI_electronic_energy + + !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + do i =1,N_st + if (i_H_psi_array(i) /= 0.d0) then + delta_e = h - CI_electronic_energy(i) + if (delta_e > 0.d0) then + e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) + else + e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) + endif + if (dabs(i_H_psi_array(i)) > 1.d-6) then + c_pert(i) = e_2_pert(i)/i_H_psi_array(i) + else + c_pert(i) = 0.d0 + endif + H_pert_diag(i) = h*c_pert(i)*c_pert(i) + else + e_2_pert(i) = 0.d0 + c_pert(i) = 0.d0 + H_pert_diag(i) = 0.d0 + endif + enddo + +end + +subroutine pt2_moller_plesset ($arguments) + use bitmasks + implicit none + $declarations + + BEGIN_DOC + ! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states. + ! + ! c_pert(i) = /(difference of orbital energies) + ! + ! e_2_pert(i) = ^2/(difference of orbital energies) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem_fock + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase,delta_e,h + double precision :: i_H_psi_array(N_st) + integer :: h1,h2,p1,p2,s1,s2 + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint) + 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_fock(det_ref,det_pert,fock_diag_tmp,Nint) + do i =1,n_st + H_pert_diag(i) = h + c_pert(i) = i_H_psi_array(i) *delta_e + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end + + +subroutine pt2_epstein_nesbet_SC2_projected ($arguments) + use bitmasks + implicit none + $declarations + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! In addition, for the perturbative energetic contribution you have the standard second order + ! + ! e_2_pert = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = c_pert + END_DOC + + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(0:ndet) + integer :: i,j,degree,l + double precision :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E + double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + !$IVDEP + do i = 1, idx_repeat(0) + accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i)) + enddo + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + h = h + accu_e_corr + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + + + c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + e_2_pert(1) = i_H_psi_array(1) * c_pert(1) + + do i =2,N_st + H_pert_diag(i) = h + if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) + else + c_pert(i) = i_H_psi_array(i) + e_2_pert(i) = -dabs(i_H_psi_array(i)) + endif + enddo + + degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + & + popcnt(xor( ref_bitmask(1,2), det_pert(1,2))) + !DEC$ NOUNROLL + do l=2,Nint + degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + & + popcnt(xor( ref_bitmask(l,2), det_pert(l,2))) + enddo + if(degree==4)then + ! + e_2_pert_fonda = e_2_pert(1) + H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1) + do i = 1, N_st + do j = 1, idx_repeat(0) + e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i) + enddo + enddo + endif + +end + + +subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments) + use bitmasks + implicit none + $declarations + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! In addition, for the perturbative energetic contribution you have the standard second order + ! + ! e_2_pert = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = c_pert + END_DOC + + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(0:ndet) + integer :: i,j,degree,l + double precision :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E + double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + !$IVDEP + do i = 1, idx_repeat(0) + accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i)) + enddo + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + h = h + accu_e_corr + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + + + c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + e_2_pert(1) = i_H_psi_array(1) * c_pert(1) + + do i =2,N_st + H_pert_diag(i) = h + if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) + else + c_pert(i) = i_H_psi_array(i) + e_2_pert(i) = -dabs(i_H_psi_array(i)) + endif + enddo +end + + + + + +subroutine pt2_epstein_nesbet_sc2 ($arguments) + use bitmasks + implicit none + $declarations + BEGIN_DOC + ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, but with the CISD_SC2 energies and coefficients + ! + ! c_pert(i) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + END_DOC + + integer :: i,j + double precision :: i_H_psi_array(N_st) + double precision :: diag_H_mat_elem_fock, h + PROVIDE selection_criterion + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + + + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + do i =1,N_st + if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then + c_pert(i) = -1.d0 + e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 + else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) + H_pert_diag(i) = h*c_pert(i)*c_pert(i) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + else + c_pert(i) = -1.d0 + e_2_pert(i) = -dabs(i_H_psi_array(i)) + H_pert_diag(i) = h + endif + enddo + +end + + + +SUBST [ arguments, declarations ] + +det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ; + + integer, intent(in) :: Nint + integer, intent(in) :: ndet + integer, intent(in) :: N_st + integer, intent(in) :: N_minilist + integer(bit_kind), intent(in) :: det_ref (Nint,2) + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1) + double precision , intent(out) :: c_pert(N_st) + double precision , intent(out) :: e_2_pert(N_st) + double precision, intent(out) :: H_pert_diag(N_st) + integer, intent(in) :: idx_minilist(0:N_det_selectors) + integer(bit_kind), intent(in) :: minilist(Nint,2,N_det_selectors) +;; + + +END_TEMPLATE + +! Note : If the arguments are changed here, they should also be changed accordingly in +! the perturbation.template.f file. + diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 26f19d0d..e1c915bc 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -131,10 +131,10 @@ class H_apply(object): def filter_vvvv_excitation(self): self["filter_vvvv_excitation"] = """ key_union_hole_part = 0_bit_kind - call set_bite_to_integer(i_a,key_union_hole_part,N_int) - call set_bite_to_integer(j_a,key_union_hole_part,N_int) - call set_bite_to_integer(i_b,key_union_hole_part,N_int) - call set_bite_to_integer(j_b,key_union_hole_part,N_int) + call set_bit_to_integer(i_a,key_union_hole_part,N_int) + call set_bit_to_integer(j_a,key_union_hole_part,N_int) + call set_bit_to_integer(i_b,key_union_hole_part,N_int) + call set_bit_to_integer(j_b,key_union_hole_part,N_int) do jtest_vvvv = 1, N_int if(iand(key_union_hole_part(jtest_vvvv),virt_bitmask(jtest_vvvv,1).ne.key_union_hole_part(jtest_vvvv)))then b_cycle = .False. @@ -157,7 +157,6 @@ class H_apply(object): def set_filter_2h_2p(self): self["filter2h2p"] = """ -! ! DIR$ FORCEINLINE if (is_a_two_holes_two_particles(key)) cycle """ @@ -205,7 +204,7 @@ class H_apply(object): """ 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, & - sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask) + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) """%(pert,) self.data["finalization"] = """ """ diff --git a/src/Determinants/Fock_diag.irp.f b/src/Determinants/Fock_diag.irp.f new file mode 100644 index 00000000..a99bbcad --- /dev/null +++ b/src/Determinants/Fock_diag.irp.f @@ -0,0 +1,84 @@ +subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Build the diagonal of the Fock matrix corresponding to a generator +! determinant. F_00 is = E0. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det_ref(Nint,2) + double precision, intent(out) :: fock_diag_tmp(2,mo_tot_num+1) + + integer :: occ(Nint*bit_kind_size,2) + integer :: ne(2), i, j, ii, jj + double precision :: E0 + + ! Compute Fock matrix diagonal elements + call bitstring_to_list_ab(det_ref,occ,Ne,Nint) + + fock_diag_tmp = 0.d0 + E0 = 0.d0 + + ! Occupied MOs + do ii=1,elec_alpha_num + i = occ(ii,1) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i) + E0 = E0 + mo_mono_elec_integral(i,i) + do jj=1,elec_alpha_num + j = occ(jj,1) + if (i==j) cycle + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j) + enddo + do jj=1,elec_beta_num + j = occ(jj,2) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j) + E0 = E0 + mo_bielec_integral_jj(i,j) + enddo + enddo + do ii=1,elec_beta_num + i = occ(ii,2) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i) + E0 = E0 + mo_mono_elec_integral(i,i) + do jj=1,elec_beta_num + j = occ(jj,2) + if (i==j) cycle + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_bielec_integral_jj_anti(i,j) + enddo + do jj=1,elec_alpha_num + j = occ(jj,1) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j) + enddo + enddo + + ! Virtual MOs + do i=1,mo_tot_num + if (fock_diag_tmp(1,i) /= 0.d0) cycle + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_mono_elec_integral(i,i) + do jj=1,elec_alpha_num + j = occ(jj,1) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj_anti(i,j) + enddo + do jj=1,elec_beta_num + j = occ(jj,2) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bielec_integral_jj(i,j) + enddo + enddo + do i=1,mo_tot_num + if (fock_diag_tmp(2,i) /= 0.d0) cycle + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_mono_elec_integral(i,i) + do jj=1,elec_beta_num + j = occ(jj,2) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj_anti(i,j) + enddo + do jj=1,elec_alpha_num + j = occ(jj,1) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bielec_integral_jj(i,j) + enddo + enddo + + fock_diag_tmp(1,mo_tot_num+1) = E0 + fock_diag_tmp(2,mo_tot_num+1) = E0 + +end diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 7ee88e28..58ae8b08 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -1,6 +1,6 @@ -subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) +subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters ) integer(bit_kind), intent(in) :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2) integer(bit_kind), intent(in) :: particl_1(N_int, 2), particl_2(N_int, 2) @@ -8,8 +8,8 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl integer,intent(in) :: i_generator,iproc_in integer(bit_kind) :: status(N_int*bit_kind_size, 2) integer :: highest, p1,p2,sp,ni,i,mi,nt,ns - - integer(bit_kind), intent(in) :: key_prev(N_int, 2, *) + double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1) + integer(bit_kind), intent(in) :: key_prev(N_int, 2, *) PROVIDE N_int PROVIDE N_det @@ -72,7 +72,7 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl if((status(p1, sp) == 1 .and. status(p2, sp) > 1) .or. & (status(p1, sp) == 2 .and. status(p2, sp) == 3) .or. & (status(p1, sp) == 3 .and. status(p2, sp) == 3 .and. p2 > p1)) then - call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, i_generator, iproc_in $parameters ) + call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters ) end if end do end do @@ -89,16 +89,17 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl (status(p1, 1) == 1 .and. status(p2, 2) >= 2) .or. & (status(p1, 1) == 2 .and. status(p2, 2) /= 2)) then - call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, i_generator, iproc_in $parameters ) + call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters ) end if end do end do end subroutine -subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, i_generator, iproc_in $parameters ) +subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters ) integer(bit_kind), intent(in) :: key_in(N_int, 2), particl_1(N_int, 2), particl_2(N_int, 2) + double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1) integer(bit_kind) :: p1_mask(N_int, 2), p2_mask(N_int, 2), key_mask(N_int, 2) integer,intent(in) :: fh1,fh2,fs1,fs2,i_generator,iproc_in integer(bit_kind) :: miniList(N_int, 2, N_det) @@ -115,11 +116,11 @@ subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, key_mask(ishft(fh1,-bit_kind_shift) + 1, fs1) -= ishft(1,iand(fh1-1,bit_kind_size-1)) key_mask(ishft(fh2,-bit_kind_shift) + 1, fs2) -= ishft(1,iand(fh2-1,bit_kind_size-1)) - call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters ) + call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters ) end subroutine -subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, i_generator, iproc_in $parameters ) +subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, fock_diag_tmp, i_generator, iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -136,6 +137,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2) integer, intent(in) :: iproc_in + double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1) integer(bit_kind), allocatable :: hole_save(:,:) integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:) integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:) @@ -175,6 +177,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),& occ_hole_tmp(N_int*bit_kind_size,2),key_union_hole_part(N_int)) + $init_thread @@ -362,17 +365,17 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl enddo ! ispin $keys_work $deinit_thread - deallocate (ia_ja_pairs, ib_jb_pairs, & - keys_out, hole_save, & - key,hole, particle, hole_tmp,& - particle_tmp, occ_particle, & - occ_hole, occ_particle_tmp,& + deallocate (ia_ja_pairs, ib_jb_pairs, & + keys_out, hole_save, & + key,hole, particle, hole_tmp, & + particle_tmp, occ_particle, & + occ_hole, occ_particle_tmp, & occ_hole_tmp,array_pairs,key_union_hole_part) $omp_end_parallel $finalization end -subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $parameters ) +subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generator,iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -387,6 +390,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa integer(bit_kind),intent(in) :: key_in(N_int,2) integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer, intent(in) :: iproc_in + double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1) integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind),allocatable :: hole_save(:,:) integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:) @@ -526,6 +530,7 @@ subroutine $subroutine($params_main) integer(bit_kind), allocatable :: mask(:,:,:) integer :: ispin, k integer :: iproc + double precision, allocatable :: fock_diag_tmp(:,:) $initialization PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators @@ -539,7 +544,7 @@ subroutine $subroutine($params_main) call wall_time(wall_0) iproc = 0 - allocate( mask(N_int,2,6) ) + allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) ) do i_generator=1,nmax progress_bar(1) = i_generator @@ -549,6 +554,9 @@ subroutine $subroutine($params_main) endif $skip + ! Compute diagonal of the Fock matrix + call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int) + ! Create bit masks for holes and particles do ispin=1,2 do k=1,N_int @@ -577,12 +585,12 @@ subroutine $subroutine($params_main) psi_det_generators(1,1,1), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & - i_generator, iproc $params_post) + fock_diag_tmp, i_generator, iproc $params_post) endif if($do_mono_excitations)then call $subroutine_monoexc(psi_det_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & - i_generator, iproc $params_post) + fock_diag_tmp, i_generator, iproc $params_post) endif call wall_time(wall_1) $printout_always @@ -592,13 +600,13 @@ subroutine $subroutine($params_main) endif enddo - deallocate( mask ) + deallocate( mask, fock_diag_tmp ) !$OMP PARALLEL DEFAULT(SHARED) & - !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc) + !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc,fock_diag_tmp) call wall_time(wall_0) !$ iproc = omp_get_thread_num() - allocate( mask(N_int,2,6) ) + allocate( mask(N_int,2,6), fock_diag_tmp(2,mo_tot_num+1) ) !$OMP DO SCHEDULE(dynamic,1) do i_generator=nmax+1,N_det_generators if (iproc == 0) then @@ -609,6 +617,9 @@ subroutine $subroutine($params_main) endif $skip + ! Compute diagonal of the Fock matrix + call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int) + ! Create bit masks for holes and particles do ispin=1,2 do k=1,N_int @@ -638,12 +649,12 @@ subroutine $subroutine($params_main) psi_det_generators(1,1,1), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & - i_generator, iproc $params_post) + fock_diag_tmp, i_generator, iproc $params_post) endif if($do_mono_excitations)then call $subroutine_monoexc(psi_det_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & - i_generator, iproc $params_post) + fock_diag_tmp, i_generator, iproc $params_post) endif !$ call omp_set_lock(lck) call wall_time(wall_1) @@ -655,7 +666,7 @@ subroutine $subroutine($params_main) !$ call omp_unset_lock(lck) enddo !$OMP END DO - deallocate( mask ) + deallocate( mask, fock_diag_tmp ) !$OMP END PARALLEL !$ call omp_destroy_lock(lck) diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index 735a8d63..19eec306 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -15,19 +15,19 @@ Documentation .. by the `update_README.py` script. -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem -`abs_psi_coef_max `_ +`abs_psi_coef_max `_ Max and min values of the coefficients -`abs_psi_coef_min `_ +`abs_psi_coef_min `_ Max and min values of the coefficients -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem @@ -39,6 +39,21 @@ Documentation Energy of the reference bitmask used in Slater rules +`bitstring_to_list_ab `_ + Gives the inidices(+1) of the bits set to 1 in the bit string + For alpha/beta determinants + + +`bitstring_to_list_ab_old `_ + Gives the inidices(+1) of the bits set to 1 in the bit string + For alpha/beta determinants + + +`build_fock_tmp `_ + Build the diagonal of the Fock matrix corresponding to a generator + determinant. F_00 is = E0. + + `ci_eigenvectors `_ Eigenvectors/values of the CI matrix @@ -100,29 +115,37 @@ Documentation Initial guess vectors are not necessarily orthonormal -`connected_to_ref `_ +`connected_to_ref `_ Undocumented -`connected_to_ref_by_mono `_ +`connected_to_ref_by_mono `_ Undocumented -`copy_h_apply_buffer_to_wf `_ +`copy_h_apply_buffer_to_wf `_ Copies the H_apply buffer to psi_coef. After calling this subroutine, N_det, psi_det and psi_coef need to be touched +`create_minilist `_ + Undocumented + + +`create_minilist_find_previous `_ + Undocumented + + `create_wf_of_psi_bilinear_matrix `_ Generate a wave function containing all possible products of alpha and beta determinants -`davidson_converged `_ +`davidson_converged `_ True if the Davidson algorithm is converged -`davidson_criterion `_ +`davidson_criterion `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] @@ -145,7 +168,7 @@ Documentation Initial guess vectors are not necessarily orthonormal -`davidson_diag_hjj `_ +`davidson_diag_hjj `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -174,10 +197,6 @@ Documentation Max number of Davidson sizes -`davidson_threshold `_ - Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] - - `decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes @@ -206,7 +225,7 @@ Documentation Undocumented -`det_occ `_ +`det_occ `_ det_occ @@ -222,10 +241,15 @@ Documentation Diagonalization algorithm (Davidson or Lapack) -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes +`diag_h_mat_elem_fock `_ + Computes when i is at most a double excitation from + a reference. + + `diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix @@ -262,11 +286,11 @@ Documentation Expected value of S2 : S*(S+1) -`fill_h_apply_buffer_no_selection `_ +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD -`filter_3_highest_electrons `_ +`filter_3_highest_electrons `_ Returns a determinant with only the 3 highest electrons @@ -282,17 +306,7 @@ Documentation idx(0) is the number of determinants that interact with key1 -`filter_connected_davidson_warp `_ - Filters out the determinants that are not connected by H - returns the array idx which contains the index of the - determinants in the array key1 that interact - via the H operator with key2. - .br - idx(0) is the number of determinants that interact with key1 - key1 should come from psi_det_sorted_ab. - - -`filter_connected_i_h_psi0 `_ +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -302,7 +316,7 @@ Documentation idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`filter_connected_i_h_psi0_sc2 `_ standard filter_connected_i_H_psi but returns in addition .br the array of the index of the non connected determinants to key1 @@ -314,7 +328,11 @@ Documentation to repeat the excitations -`generate_all_alpha_beta_det_products `_ +`first_guess `_ + Select all the determinants with the lowest energy as a starting point. + + +`generate_all_alpha_beta_det_products `_ Create a wave function from all possible alpha x beta determinants @@ -330,7 +348,7 @@ Documentation Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants @@ -350,7 +368,7 @@ Documentation Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring @@ -384,7 +402,7 @@ Documentation Undocumented -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -392,23 +410,35 @@ Documentation H_jj : array of -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_j_phase_out `_ +`i_h_j_phase_out `_ Returns where i and j are determinants -`i_h_j_verbose `_ +`i_h_j_verbose `_ Returns where i and j are determinants -`i_h_psi `_ - for the various Nstates +`i_h_psi `_ + Computes = \sum_J c_J . + .br + Uses filter_connected_i_H_psi0 to get all the |J> to which |i> + is connected. + The i_H_psi_minilist is much faster but requires to build the + minilists -`i_h_psi_sc2 `_ +`i_h_psi_minilist `_ + Computes = \sum_J c_J . + .br + Uses filter_connected_i_H_psi0 to get all the |J> to which |i> + is connected. The |J> are searched in short pre-computed lists. + + +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -422,7 +452,7 @@ Documentation to repeat the excitations -`i_h_psi_sc2_verbose `_ +`i_h_psi_sc2_verbose `_ for the various Nstate .br returns in addition @@ -436,7 +466,7 @@ Documentation to repeat the excitations -`i_h_psi_sec_ord `_ +`i_h_psi_sec_ord `_ for the various Nstates @@ -451,7 +481,7 @@ Documentation idx_non_cas gives the indice of the determinant in psi_det. -`int_of_3_highest_electrons `_ +`int_of_3_highest_electrons `_ Returns an integer*8 as : .br |_<--- 21 bits ---><--- 21 bits ---><--- 21 bits --->| @@ -463,6 +493,10 @@ Documentation .br +`is_connected_to `_ + Undocumented + + `is_in_wavefunction `_ True if the determinant ``det`` is in the wave function @@ -475,7 +509,7 @@ Documentation Undocumented -`max_degree_exc `_ +`max_degree_exc `_ Maximum degree of excitation in the wf @@ -508,7 +542,7 @@ Documentation Maximum number of determinants diagonalized by Jacobi -`n_det_max_property `_ +`n_det_max_property `_ Max number of determinants in the wave function when you select for a given property @@ -556,7 +590,7 @@ Documentation Number of possible determinants for a given occ_pattern -`one_body_dm_mo `_ +`one_body_dm_mo `_ One-body density matrix @@ -568,15 +602,15 @@ Documentation Alpha and beta one-body density matrix for each state -`one_body_single_double_dm_mo_alpha `_ +`one_body_single_double_dm_mo_alpha `_ Alpha and beta one-body density matrix for each state -`one_body_single_double_dm_mo_beta `_ +`one_body_single_double_dm_mo_beta `_ Alpha and beta one-body density matrix for each state -`one_body_spin_density_mo `_ +`one_body_spin_density_mo `_ rho(alpha) - rho(beta) @@ -584,15 +618,11 @@ Documentation If true, The One body DM is calculated with ignoring the Double<->Doubles extra diag elements -`pouet `_ - Undocumented - - -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) @@ -644,7 +674,7 @@ Documentation function. -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty @@ -653,33 +683,33 @@ Documentation Undocumented -`psi_coef_max `_ +`psi_coef_max `_ Max and min values of the coefficients -`psi_coef_min `_ +`psi_coef_min `_ Max and min values of the coefficients -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_ab `_ +`psi_coef_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave function. -`psi_det `_ +`psi_det `_ The wave function determinants. Initialized with Hartree-Fock if the EZFIO file is empty @@ -700,29 +730,29 @@ Documentation Unique beta determinants -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_ab `_ +`psi_det_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave function. -`psi_det_sorted_next_ab `_ +`psi_det_sorted_next_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate @@ -761,7 +791,7 @@ Documentation Undocumented -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file @@ -777,20 +807,16 @@ Documentation Energy of the reference bitmask used in Slater rules -`remove_duplicates_in_psi_det `_ +`remove_duplicates_in_psi_det `_ Removes duplicate determinants in the wave function. -`resize_h_apply_buffer `_ +`resize_h_apply_buffer `_ Resizes the H_apply buffer of proc iproc. The buffer lock should be set before calling this function. -`routine `_ - Undocumented - - -`s2_eig `_ +`s2_eig `_ Force the wave function to be an eigenfunction of S^2 @@ -810,31 +836,31 @@ Documentation Undocumented -`save_natural_mos `_ +`save_natural_mos `_ Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`save_wavefunction_general `_ +`save_wavefunction_general `_ Save the wave function into the EZFIO file -`save_wavefunction_specified `_ +`save_wavefunction_specified `_ Save the wave function into the EZFIO file -`save_wavefunction_unsorted `_ +`save_wavefunction_unsorted `_ Save the wave function into the EZFIO file -`set_bite_to_integer `_ +`set_bit_to_integer `_ Undocumented -`set_natural_mos `_ +`set_natural_mos `_ Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis @@ -844,26 +870,26 @@ Documentation for a given couple of hole/particle excitations i. -`sort_dets_ab `_ - Undocumented +`sort_dets_ab `_ + Uncodumented : TODO -`sort_dets_ab_v `_ - Undocumented +`sort_dets_ab_v `_ + Uncodumented : TODO -`sort_dets_ba_v `_ - Undocumented +`sort_dets_ba_v `_ + Uncodumented : TODO -`sort_dets_by_3_highest_electrons `_ +`sort_dets_by_3_highest_electrons `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`sort_dets_by_det_search_key `_ +`sort_dets_by_det_search_key `_ Determinants are sorted are sorted according to their det_search_key. Useful to accelerate the search of a random determinant in the wave function. @@ -873,12 +899,12 @@ Documentation Return an integer*8 corresponding to a determinant index for searching -`state_average_weight `_ +`state_average_weight `_ Weights in the state-average calculation of the density matrix `tamiser `_ - Undocumented + Uncodumented : TODO `target_energy `_ @@ -889,7 +915,11 @@ Documentation convergence of the correlation energy of SC2 iterations -`threshold_generators `_ +`threshold_davidson `_ + Thresholds of Davidson's algorithm + + +`threshold_generators `_ Thresholds on generators (fraction of the norm) diff --git a/src/Determinants/create_excitations.irp.f b/src/Determinants/create_excitations.irp.f index a2acc8df..b7233beb 100644 --- a/src/Determinants/create_excitations.irp.f +++ b/src/Determinants/create_excitations.irp.f @@ -35,7 +35,7 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) endif end -subroutine set_bite_to_integer(i_physical,key,Nint) +subroutine set_bit_to_integer(i_physical,key,Nint) use bitmasks implicit none integer, intent(in) :: i_physical,Nint diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 6e980402..32e84532 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1264,6 +1264,75 @@ end +double precision function diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes when i is at most a double excitation from + ! a reference. + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_ref(Nint,2), det_pert(Nint,2) + double precision, intent(in) :: fock_diag_tmp(2,mo_tot_num+1) + + integer :: degree + double precision :: phase, E0 + integer :: exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + + call get_excitation_degree(det_ref,det_pert,degree,Nint) + E0 = fock_diag_tmp(1,mo_tot_num+1) + if (degree == 2) then + call get_double_excitation(det_ref,det_pert,exc,phase,Nint) + call decode_exc(exc,2,h1,p1,h2,p2,s1,s2) + + if ( (s1 == 1).and.(s2 == 1) ) then ! alpha/alpha + diag_H_mat_elem_fock = E0 & + - fock_diag_tmp(1,h1) & + + ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) & + - ( fock_diag_tmp(1,h2) - mo_bielec_integral_jj_anti(h1,h2) & + + mo_bielec_integral_jj_anti(p1,h2) ) & + + ( fock_diag_tmp(1,p2) - mo_bielec_integral_jj_anti(h1,p2) & + + mo_bielec_integral_jj_anti(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) ) + + else if ( (s1 == 2).and.(s2 == 2) ) then ! beta/beta + diag_H_mat_elem_fock = E0 & + - fock_diag_tmp(2,h1) & + + ( fock_diag_tmp(2,p1) - mo_bielec_integral_jj_anti(h1,p1) ) & + - ( fock_diag_tmp(2,h2) - mo_bielec_integral_jj_anti(h1,h2) & + + mo_bielec_integral_jj_anti(p1,h2) ) & + + ( fock_diag_tmp(2,p2) - mo_bielec_integral_jj_anti(h1,p2) & + + mo_bielec_integral_jj_anti(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) ) + + else ! alpha/beta + diag_H_mat_elem_fock = E0 & + - fock_diag_tmp(1,h1) & + + ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) & + - ( fock_diag_tmp(2,h2) - mo_bielec_integral_jj(h1,h2) & + + mo_bielec_integral_jj(p1,h2) ) & + + ( fock_diag_tmp(2,p2) - mo_bielec_integral_jj(h1,p2) & + + mo_bielec_integral_jj(p1,p2) - mo_bielec_integral_jj_anti(h2,p2) ) + + endif + + else if (degree == 1) then + call get_mono_excitation(det_ref,det_pert,exc,phase,Nint) + call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) + if (s1 == 1) then + diag_H_mat_elem_fock = E0 - fock_diag_tmp(1,h1) & + + ( fock_diag_tmp(1,p1) - mo_bielec_integral_jj_anti(h1,p1) ) + else + diag_H_mat_elem_fock = E0 - fock_diag_tmp(2,h1) & + + ( fock_diag_tmp(2,p1) - mo_bielec_integral_jj_anti(h1,p1) ) + endif + + else if (degree == 0) then + diag_H_mat_elem_fock = E0 + else + STOP 'Bug in diag_H_mat_elem_fock' + endif +end + double precision function diag_H_mat_elem(det_in,Nint) implicit none BEGIN_DOC diff --git a/src/Ezfio_files/README.rst b/src/Ezfio_files/README.rst index e5726d91..c7183a64 100644 --- a/src/Ezfio_files/README.rst +++ b/src/Ezfio_files/README.rst @@ -203,34 +203,10 @@ output_bitmask Output file for Bitmask -output_cas_sd - Output file for CAS_SD - - -output_cid - Output file for CID - - -output_cisd - Output file for CISD - - -output_cisd_sc2_selected - Output file for CISD_SC2_selected - - -output_cisd_selected - Output file for CISD_selected - - `output_cpu_time_0 `_ Initial CPU and wall times when printing in the output files -output_ddci_selected - Output file for DDCI_selected - - output_determinants Output file for Determinants @@ -243,18 +219,10 @@ output_ezfio_files Output file for Ezfio_files -output_fcidump - Output file for FCIdump - - output_full_ci Output file for Full_CI -output_generators_cas - Output file for Generators_CAS - - output_generators_full Output file for Generators_full @@ -271,10 +239,6 @@ output_integrals_monoelec Output file for Integrals_Monoelec -output_loc_cele - Output file for loc_cele - - output_mo_basis Output file for MO_Basis @@ -283,14 +247,6 @@ output_moguess Output file for MOGuess -output_molden - Output file for Molden - - -output_mp2 - Output file for MP2 - - output_mrcc_cassd Output file for MRCC_CASSD @@ -323,18 +279,10 @@ output_psiref_utils Output file for Psiref_Utils -output_qmcchem - Output file for QmcChem - - output_selectors_full Output file for Selectors_full -output_singlerefmethod - Output file for SingleRefMethod - - output_utils Output file for Utils diff --git a/src/Integrals_Bielec/.gitignore b/src/Integrals_Bielec/.gitignore index 54da4aed..f4bdeaca 100644 --- a/src/Integrals_Bielec/.gitignore +++ b/src/Integrals_Bielec/.gitignore @@ -1,20 +1,20 @@ -# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py -IRPF90_temp +# Automatically created by $QP_ROOT/scripts/module/module_handler.py +.ninja_deps +.ninja_log +AO_Basis +Bitmask +Electrons +Ezfio_files IRPF90_man -irpf90_entities -tags -irpf90.make +IRPF90_temp +MO_Basis Makefile Makefile.depend -build.ninja -.ninja_log -.ninja_deps -ezfio_interface.irp.f -Ezfio_files -MO_Basis -Utils +Nuclei Pseudo -Bitmask -AO_Basis -Electrons -Nuclei \ No newline at end of file +Utils +ezfio_interface.irp.f +irpf90.make +irpf90_entities +tags +test_integrals \ No newline at end of file diff --git a/src/Integrals_Bielec/README.rst b/src/Integrals_Bielec/README.rst index b71d9c0d..f09d1a9c 100644 --- a/src/Integrals_Bielec/README.rst +++ b/src/Integrals_Bielec/README.rst @@ -47,7 +47,7 @@ Documentation i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integral_schwartz `_ +`ao_bielec_integral_schwartz `_ Needed to compute Schwartz inequalities @@ -56,7 +56,7 @@ Documentation i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integrals_in_map `_ +`ao_bielec_integrals_in_map `_ Map of Atomic integrals i(r1) j(r2) 1/r12 k(r1) l(r2) @@ -73,6 +73,10 @@ Documentation Computes the product of l values of i,j,k,and l +`bench_maps `_ + Performs timing benchmarks on integral access + + `bielec_integrals_index `_ Undocumented @@ -85,7 +89,7 @@ Documentation Frees the memory of the AO map -`clear_mo_map `_ +`clear_mo_map `_ Frees the memory of the MO map @@ -105,15 +109,15 @@ Documentation Compute integrals on the fly -`dump_ao_integrals `_ +`dump_ao_integrals `_ Save to disk the $ao integrals -`dump_mo_integrals `_ +`dump_mo_integrals `_ Save to disk the $ao integrals -`eri `_ +`eri `_ ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) @@ -135,7 +139,7 @@ Documentation t_w(i,2,k) = t(i) -`general_primitive_integral `_ +`general_primitive_integral `_ Computes the integral where p,q,r,s are Gaussian primitives @@ -161,52 +165,56 @@ Documentation Returns one integral in the MO basis -`get_mo_bielec_integrals `_ +`get_mo_bielec_integral_schwartz `_ + Returns one integral in the MO basis + + +`get_mo_bielec_integrals `_ Returns multiple integrals in the MO basis, all i for j,k,l fixed. -`get_mo_bielec_integrals_existing_ik `_ +`get_mo_bielec_integrals_existing_ik `_ Returns multiple integrals in the MO basis, all i(1)j(1) 1/r12 k(2)l(2) i for j,k,l fixed. -`get_mo_map_size `_ +`get_mo_map_size `_ Return the number of elements in the MO map -`give_polynom_mult_center_x `_ +`give_polynom_mult_center_x `_ subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw : I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) -`i_x1_new `_ +`i_x1_new `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult `_ +`i_x1_pol_mult `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a1 `_ +`i_x1_pol_mult_a1 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a2 `_ +`i_x1_pol_mult_a2 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_recurs `_ +`i_x1_pol_mult_recurs `_ recursive function involved in the bielectronic integral -`i_x2_new `_ +`i_x2_new `_ recursive function involved in the bielectronic integral -`i_x2_pol_mult `_ +`i_x2_pol_mult `_ recursive function involved in the bielectronic integral @@ -218,21 +226,21 @@ Documentation Create new entry into MO map, or accumulate in an existing entry -`integrale_new `_ +`integrale_new `_ calculate the integral of the polynom :: I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) between ( 0 ; 1) -`load_ao_integrals `_ +`load_ao_integrals `_ Read from disk the $ao integrals -`load_mo_integrals `_ +`load_mo_integrals `_ Read from disk the $ao integrals -`mo_bielec_integral `_ +`mo_bielec_integral `_ Returns one integral in the MO basis @@ -272,6 +280,10 @@ Documentation mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij +`mo_bielec_integral_schwartz `_ + Needed to compute Schwartz inequalities + + `mo_bielec_integrals_in_map `_ If True, the map of MO bielectronic integrals is provided @@ -292,7 +304,7 @@ Documentation Aligned n_pt_max_integrals -`n_pt_sup `_ +`n_pt_sup `_ Returns the upper boundary of the degree of the polynomial involved in the bielctronic integral : Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) diff --git a/src/Integrals_Monoelec/README.rst b/src/Integrals_Monoelec/README.rst index 16230cbf..1d2d158b 100644 --- a/src/Integrals_Monoelec/README.rst +++ b/src/Integrals_Monoelec/README.rst @@ -111,11 +111,11 @@ Documentation Pseudo-potential -`ao_pseudo_integral_local `_ +`ao_pseudo_integral_local `_ Local pseudo-potential -`ao_pseudo_integral_non_local `_ +`ao_pseudo_integral_non_local `_ Local pseudo-potential