Merge pull request #2 from LCPQ/master

merge with LCPQ

configure

@@ -62,8 +62,9 @@ d_dependency = {
     "resultsFile": ["python"],
     "emsl": ["python"],
     "gcc": [],
+    "g++": [],
     "python": [],
-    "ninja": ["gcc", "python"],
+    "ninja": ["g++", "python"],
     "make": [],
     "p_graphviz": ["python"]
 }

ocaml/Input_determinants_by_hand.ml

@@ -7,8 +7,6 @@ module Determinants_by_hand : sig
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
-      n_states               : States_number.t;
-      n_states_diag          : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -23,8 +21,6 @@ end = struct
     { n_int                  : N_int_number.t;
       bit_kind               : Bit_kind.t;
       n_det                  : Det_number.t;
-      n_states               : States_number.t;
-      n_states_diag          : States_number.t;
       expected_s2            : Positive_float.t;
       psi_coef               : Det_coef.t array;
       psi_det                : Determinant.t array;
@@ -146,11 +142,12 @@ end = struct
     |> Array.map ~f:Det_coef.of_float
   ;;
 
-  let write_psi_coef ~n_det ~n_states c =
+  let write_psi_coef ~n_det c =
     let n_det = Det_number.to_int n_det
     and c = Array.to_list c
             |> List.map ~f:Det_coef.to_float
-    and n_states = States_number.to_int n_states
+    and n_states = 
+      read_n_states () |> States_number.to_int
     in
     Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
     |> Ezfio.set_determinants_psi_coef 
@@ -214,8 +211,6 @@ end = struct
       { n_int                  = read_n_int ()                ;
         bit_kind               = read_bit_kind ()             ;
         n_det                  = read_n_det ()                ;
-        n_states               = read_n_states ()             ;
-        n_states_diag          = read_n_states_diag ()        ;
         expected_s2            = read_expected_s2 ()          ;
         psi_coef               = read_psi_coef ()             ;
         psi_det                = read_psi_det ()              ;
@@ -227,8 +222,6 @@ end = struct
   let write { n_int                ;
               bit_kind             ;
               n_det                ;
-              n_states             ;
-              n_states_diag        ;
               expected_s2          ;
               psi_coef             ;
               psi_det              ;
@@ -236,10 +229,8 @@ end = struct
      write_n_int n_int ;
      write_bit_kind bit_kind;
      write_n_det n_det;
-     write_n_states n_states;
-     write_n_states_diag ~n_states:n_states n_states_diag;
      write_expected_s2 expected_s2;
-     write_psi_coef ~n_det:n_det psi_coef ~n_states:n_states;
+     write_psi_coef ~n_det:n_det psi_coef ;
      write_psi_det ~n_int:n_int ~n_det:n_det psi_det;
   ;;
 
@@ -249,7 +240,7 @@ end = struct
     let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in
     let det_text = 
       let nstates =
-        States_number.to_int b.n_states
+        read_n_states () |> States_number.to_int
       and ndet =
         Det_number.to_int b.n_det
       in
@@ -284,12 +275,6 @@ If true, input the expected value of S^2 ::
 
   expected_s2 = %s
 
-Number of requested states, and number of states used for the
-Davidson diagonalization ::
-
-  n_states      = %s
-  n_states_diag = %s
-
 Number of determinants ::
 
   n_det = %s
@@ -299,8 +284,6 @@ Determinants ::
 %s
 "
      (b.expected_s2   |> Positive_float.to_string)
-     (b.n_states      |> States_number.to_string)
-     (b.n_states_diag |> States_number.to_string)
      (b.n_det         |> Det_number.to_string)
      det_text
      |> Rst_string.of_string
@@ -313,8 +296,6 @@ Determinants ::
 n_int                  = %s
 bit_kind               = %s
 n_det                  = %s
-n_states               = %s
-n_states_diag          = %s
 expected_s2            = %s
 psi_coef               = %s
 psi_det                = %s
@@ -322,8 +303,6 @@ psi_det                = %s
      (b.n_int         |> N_int_number.to_string)
      (b.bit_kind      |> Bit_kind.to_string)
      (b.n_det         |> Det_number.to_string)
-     (b.n_states      |> States_number.to_string)
-     (b.n_states_diag |> States_number.to_string)
      (b.expected_s2   |> Positive_float.to_string)
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")

ocaml/Qputils.ml

@@ -12,7 +12,6 @@ let rec transpose = function
 ;;
 *)
 
-
 let input_to_sexp s =
     let result = 
       String.split_lines s 

ocaml/qptypes_generator.ml

@@ -83,6 +83,12 @@ let input_data = "
   assert (x >= 0.) ;
   assert (x <= 1.) ;
 
+* Energy : float
+  assert (x <=0.) ;
+
+* S2 : float
+  assert (x >=0.) ;
+
 * PT2_energy : float
   assert (x >=0.) ;
 

plugins/CAS_SD/README.rst

@@ -118,101 +118,101 @@ Documentation
   Undocumented
 
 
-`h_apply_cas_sd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L537>`_
+h_apply_cas_sd
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L3>`_
+h_apply_cas_sd_diexc
   Undocumented
 
 
-`h_apply_cas_sd_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L119>`_
+h_apply_cas_sd_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L96>`_
+h_apply_cas_sd_diexcp
   Undocumented
 
 
-`h_apply_cas_sd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L387>`_
+h_apply_cas_sd_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L3175>`_
+h_apply_cas_sd_pt2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2540>`_
+h_apply_cas_sd_pt2_diexc
   Undocumented
 
 
-`h_apply_cas_sd_pt2_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2678>`_
+h_apply_cas_sd_pt2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_pt2_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2644>`_
+h_apply_cas_sd_pt2_diexcp
   Undocumented
 
 
-`h_apply_cas_sd_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2987>`_
+h_apply_cas_sd_pt2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2291>`_
+h_apply_cas_sd_selected
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_selected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1620>`_
+h_apply_cas_sd_selected_diexc
   Undocumented
 
 
-`h_apply_cas_sd_selected_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1760>`_
+h_apply_cas_sd_selected_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1725>`_
+h_apply_cas_sd_selected_diexcp
   Undocumented
 
 
-`h_apply_cas_sd_selected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L2089>`_
+h_apply_cas_sd_selected_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1399>`_
+h_apply_cas_sd_selected_no_skip
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cas_sd_selected_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L728>`_
+h_apply_cas_sd_selected_no_skip_diexc
   Undocumented
 
 
-`h_apply_cas_sd_selected_no_skip_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L868>`_
+h_apply_cas_sd_selected_no_skip_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cas_sd_selected_no_skip_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L833>`_
+h_apply_cas_sd_selected_no_skip_diexcp
   Undocumented
 
 
-`h_apply_cas_sd_selected_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CAS_SD/H_apply.irp.f_shell_22#L1197>`_
+h_apply_cas_sd_selected_no_skip_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/CID/README.rst

@@ -203,3 +203,48 @@ Documentation
   particles.
   Assume N_int is already provided.
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/plugins/SingleRefMethod>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`cid <http://github.com/LCPQ/quantum_package/tree/master/plugins/CID/cid_lapack.irp.f#L1>`_
+  Undocumented
+
+
+h_apply_cid
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cid_diexc
+  Undocumented
+
+
+h_apply_cid_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cid_diexcp
+  Undocumented
+
+
+h_apply_cid_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+

plugins/CISD/.gitignore

@@ -20,6 +20,7 @@ Pseudo
 Selectors_full
 SingleRefMethod
 Utils
+cisd
 cisd_lapack
 ezfio_interface.irp.f
 irpf90.make

plugins/CISD/EZFIO.cfg

@@ -0,0 +1,10 @@
+[energy]
+type: double precision
+doc: Variational CISD energy
+interface: ezfio
+
+[energy_pt2]
+type: double precision
+doc: Estimated CISD energy (including PT2)
+interface: ezfio
+

plugins/CISD/README.rst

@@ -59,22 +59,26 @@ Documentation
 .. by the `update_README.py` script.
 
 
-`cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/cisd_lapack.irp.f#L1>`_
+h_apply_cisd
-  Undocumented
-
-
-`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/H_apply.irp.f_shell_8#L414>`_
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/H_apply.irp.f_shell_8#L1>`_
+h_apply_cisd_diexc
+  Undocumented
+
+
+h_apply_cisd_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD/H_apply.irp.f_shell_8#L269>`_
+h_apply_cisd_diexcp
+  Undocumented
+
+
+h_apply_cisd_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/CISD_SC2_selected/README.rst

@@ -72,3 +72,122 @@ Needed Modules
 
 * `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`cisd_sc2_selected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_SC2_selected/cisd_sc2_selection.irp.f#L1>`_
+  Undocumented
+
+
+h_apply_cisd
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_cisd_diexc
+  Undocumented
+
+
+h_apply_cisd_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_cisd_diexcp
+  Undocumented
+
+
+h_apply_cisd_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_pt2_diexc
+  Undocumented
+
+
+h_apply_pt2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2_diexcp
+  Undocumented
+
+
+h_apply_pt2_en_sc2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_pt2_en_sc2_diexc
+  Undocumented
+
+
+h_apply_pt2_en_sc2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2_en_sc2_diexcp
+  Undocumented
+
+
+h_apply_pt2_en_sc2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_pt2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_sc2_selected
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_sc2_selected_diexc
+  Undocumented
+
+
+h_apply_sc2_selected_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_sc2_selected_diexcp
+  Undocumented
+
+
+h_apply_sc2_selected_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+

plugins/CISD_selected/README.rst

@@ -196,22 +196,26 @@ Documentation
 .. by the `update_README.py` script.
 
 
-`cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/cisd_selection.irp.f#L1>`_
+h_apply_cisd
-  Undocumented
-
-
-`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_8#L414>`_
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_8#L1>`_
+h_apply_cisd_diexc
+  Undocumented
+
+
+h_apply_cisd_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_8#L269>`_
+h_apply_cisd_diexcp
+  Undocumented
+
+
+h_apply_cisd_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
@@ -221,154 +225,226 @@ Documentation
   Undocumented
 
 
-`h_apply_cisd_selection_delta_rho_one_point <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5931>`_
+h_apply_cisd_selection_delta_rho_one_point
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_delta_rho_one_point_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5405>`_
+h_apply_cisd_selection_delta_rho_one_point_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_delta_rho_one_point_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_delta_rho_one_point_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5734>`_
+h_apply_cisd_selection_delta_rho_one_point_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_delta_rho_one_point_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L5159>`_
+h_apply_cisd_selection_dipole_moment_z
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_dipole_moment_z_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4633>`_
+h_apply_cisd_selection_dipole_moment_z_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_dipole_moment_z_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_dipole_moment_z_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4962>`_
+h_apply_cisd_selection_dipole_moment_z_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_dipole_moment_z_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3615>`_
+h_apply_cisd_selection_epstein_nesbet
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4387>`_
+h_apply_cisd_selection_epstein_nesbet_2x2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_2x2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3861>`_
+h_apply_cisd_selection_epstein_nesbet_2x2_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_2x2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_2x2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L4190>`_
+h_apply_cisd_selection_epstein_nesbet_2x2_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_2x2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3089>`_
+h_apply_cisd_selection_epstein_nesbet_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L3418>`_
+h_apply_cisd_selection_epstein_nesbet_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2843>`_
+h_apply_cisd_selection_epstein_nesbet_sc2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2317>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2646>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L2071>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1545>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1874>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_no_projected_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1299>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_projected
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L773>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_epstein_nesbet_sc2_projected_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1102>`_
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_epstein_nesbet_sc2_projected_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L527>`_
+h_apply_cisd_selection_h_core
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_h_core_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L1>`_
+h_apply_cisd_selection_h_core_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_h_core_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_h_core_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L330>`_
+h_apply_cisd_selection_h_core_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_h_core_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L6703>`_
+h_apply_cisd_selection_moller_plesset
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_cisd_selection_moller_plesset_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L6177>`_
+h_apply_cisd_selection_moller_plesset_diexc
+  Undocumented
+
+
+h_apply_cisd_selection_moller_plesset_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_cisd_selection_moller_plesset_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/CISD_selected/H_apply.irp.f_shell_10#L6506>`_
+h_apply_cisd_selection_moller_plesset_diexcp
+  Undocumented
+
+
+h_apply_cisd_selection_moller_plesset_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/DDCI_selected/README.rst

@@ -57,3 +57,74 @@ Needed Modules
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `Generators_CAS <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_CAS>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
+* `Generators_CAS <http://github.com/LCPQ/quantum_package/tree/master/plugins/Generators_CAS>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`ddci <http://github.com/LCPQ/quantum_package/tree/master/plugins/DDCI_selected/ddci.irp.f#L1>`_
+  Undocumented
+
+
+h_apply_ddci_pt2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_ddci_pt2_diexc
+  Undocumented
+
+
+h_apply_ddci_pt2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_ddci_pt2_diexcp
+  Undocumented
+
+
+h_apply_ddci_pt2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_ddci_selection
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_ddci_selection_diexc
+  Undocumented
+
+
+h_apply_ddci_selection_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_ddci_selection_diexcp
+  Undocumented
+
+
+h_apply_ddci_selection_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+

plugins/Full_CI/.gitignore

@@ -0,0 +1,32 @@
+# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
+.ninja_deps
+.ninja_log
+AO_Basis
+Bitmask
+Determinants
+Electrons
+Ezfio_files
+Generators_full
+Hartree_Fock
+IRPF90_man
+IRPF90_temp
+Integrals_Bielec
+Integrals_Monoelec
+MOGuess
+MO_Basis
+Makefile
+Makefile.depend
+Nuclei
+Perturbation
+Properties
+Pseudo
+Selectors_full
+Utils
+ezfio_interface.irp.f
+full_ci
+full_ci_no_skip
+irpf90.make
+irpf90_entities
+tags
+target_pt2
+var_pt2_ratio

plugins/Full_CI/README.rst

@@ -27,201 +27,201 @@ Documentation
   Undocumented
 
 
-`h_apply_fci <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L674>`_
+h_apply_fci
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3>`_
+h_apply_fci_diexc
   Undocumented
 
 
-`h_apply_fci_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L143>`_
+h_apply_fci_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L108>`_
+h_apply_fci_diexcp
   Undocumented
 
 
-`h_apply_fci_mono <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3333>`_
+h_apply_fci_mono
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_mono_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2660>`_
+h_apply_fci_mono_diexc
   Undocumented
 
 
-`h_apply_fci_mono_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2800>`_
+h_apply_fci_mono_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_mono_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2765>`_
+h_apply_fci_mono_diexcp
   Undocumented
 
 
-`h_apply_fci_mono_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3129>`_
+h_apply_fci_mono_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L472>`_
+h_apply_fci_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_no_skip <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2439>`_
+h_apply_fci_no_skip
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_no_skip_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1768>`_
+h_apply_fci_no_skip_diexc
   Undocumented
 
 
-`h_apply_fci_no_skip_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1908>`_
+h_apply_fci_no_skip_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_no_skip_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1873>`_
+h_apply_fci_no_skip_diexcp
   Undocumented
 
 
-`h_apply_fci_no_skip_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L2237>`_
+h_apply_fci_no_skip_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1558>`_
+h_apply_fci_pt2
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_fci_pt2_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L923>`_
+h_apply_fci_pt2_diexc
   Undocumented
 
 
-`h_apply_fci_pt2_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1061>`_
+h_apply_fci_pt2_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_fci_pt2_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1027>`_
+h_apply_fci_pt2_diexcp
   Undocumented
 
 
-`h_apply_fci_pt2_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L1370>`_
+h_apply_fci_pt2_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5141>`_
+h_apply_pt2_mono_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_pt2_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4504>`_
+h_apply_pt2_mono_delta_rho_diexc
   Undocumented
 
 
-`h_apply_pt2_mono_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4642>`_
+h_apply_pt2_mono_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4608>`_
+h_apply_pt2_mono_delta_rho_diexcp
   Undocumented
 
 
-`h_apply_pt2_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4951>`_
+h_apply_pt2_mono_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6906>`_
+h_apply_pt2_mono_di_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_pt2_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6271>`_
+h_apply_pt2_mono_di_delta_rho_diexc
   Undocumented
 
 
-`h_apply_pt2_mono_di_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6409>`_
+h_apply_pt2_mono_di_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_pt2_mono_di_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6375>`_
+h_apply_pt2_mono_di_delta_rho_diexcp
   Undocumented
 
 
-`h_apply_pt2_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6718>`_
+h_apply_pt2_mono_di_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4255>`_
+h_apply_select_mono_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_select_mono_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3582>`_
+h_apply_select_mono_delta_rho_diexc
   Undocumented
 
 
-`h_apply_select_mono_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3722>`_
+h_apply_select_mono_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L3687>`_
+h_apply_select_mono_delta_rho_diexcp
   Undocumented
 
 
-`h_apply_select_mono_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L4051>`_
+h_apply_select_mono_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_di_delta_rho <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L6022>`_
+h_apply_select_mono_di_delta_rho
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_select_mono_di_delta_rho_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5351>`_
+h_apply_select_mono_di_delta_rho_diexc
   Undocumented
 
 
-`h_apply_select_mono_di_delta_rho_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5491>`_
+h_apply_select_mono_di_delta_rho_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_select_mono_di_delta_rho_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5456>`_
+h_apply_select_mono_di_delta_rho_diexcp
   Undocumented
 
 
-`h_apply_select_mono_di_delta_rho_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/Full_CI/H_apply.irp.f_shell_43#L5820>`_
+h_apply_select_mono_di_delta_rho_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.

plugins/Full_CI/target_pt2.irp.f

@@ -27,6 +27,8 @@ program var_pt2_ratio_run
     call diagonalize_CI
     ratio = (CI_energy(1) - HF_energy) / (CI_energy(1)+pt2(1) - HF_energy)
     if (N_det > 20000) then
+      N_det = 20000
+      TOUCH N_det
       exit
     endif
   enddo

plugins/Hartree_Fock/diagonalize_fock.irp.f

@@ -73,19 +73,19 @@ BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
  END_DOC
  integer :: i,j
  double precision :: accu
- do i = 1,elec_alpha_num
+ do j = 1,elec_alpha_num
   accu = 0.d0
-  do j = 1, elec_alpha_num
+  do i = 1, elec_alpha_num
    accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
   enddo
-  diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i)
+  diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
  enddo
- do i = elec_alpha_num+1,mo_tot_num
+ do j = elec_alpha_num+1,mo_tot_num
   accu = 0.d0
-  do j = 1, elec_alpha_num
+  do i = 1, elec_alpha_num
    accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
   enddo
-  diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i)
+  diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
  enddo
 
 END_PROVIDER

plugins/Hartree_Fock/huckel.irp.f

@@ -4,7 +4,7 @@ subroutine huckel_guess
 ! Build the MOs using the extended Huckel model
   END_DOC
   integer                        :: i,j
-  double precision               :: tmp_matrix(ao_num_align,ao_num),accu
+  double precision               :: accu
   double precision               :: c
   character*(64)                 :: label
 
@@ -19,14 +19,12 @@ subroutine huckel_guess
   c = 0.5d0 * 1.75d0
 
   do j=1,ao_num
+    !DIR$ VECTOR ALIGNED
     do i=1,ao_num
-      if (i.ne.j) then
+      Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
-        Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral(i,i) + &
+                                                 ao_mono_elec_integral_diag(j))
-                                                 ao_mono_elec_integral(j,j))
-      else
-        Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j)
-      endif
     enddo
+    Fock_matrix_ao(j,j) = Fock_matrix_alpha_ao(j,j)
   enddo
   TOUCH Fock_matrix_ao
   mo_coef = eigenvectors_fock_matrix_mo

plugins/MP2/README.rst

@@ -40,3 +40,49 @@ Needed Modules
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
+* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
+* `SingleRefMethod <http://github.com/LCPQ/quantum_package/tree/master/plugins/SingleRefMethod>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+h_apply_mp2
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
+
+
+h_apply_mp2_diexc
+  Undocumented
+
+
+h_apply_mp2_diexcorg
+  Generate all double excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+h_apply_mp2_diexcp
+  Undocumented
+
+
+h_apply_mp2_monoexc
+  Generate all single excitations of key_in using the bit masks of holes and
+  particles.
+  Assume N_int is already provided.
+
+
+`mp2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/MP2/mp2.irp.f#L1>`_
+  Undocumented
+

plugins/MRCC_Utils/README.rst

@@ -283,26 +283,26 @@ Documentation
   Undocumented
 
 
-`h_apply_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L553>`_
+h_apply_mrcc
   Calls H_apply on the HF determinant and selects all connected single and double
   excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
 
 
-`h_apply_mrcc_diexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L3>`_
+h_apply_mrcc_diexc
   Undocumented
 
 
-`h_apply_mrcc_diexcorg <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L127>`_
+h_apply_mrcc_diexcorg
   Generate all double excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
 
 
-`h_apply_mrcc_diexcp <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L100>`_
+h_apply_mrcc_diexcp
   Undocumented
 
 
-`h_apply_mrcc_monoexc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/H_apply.irp.f_shell_27#L399>`_
+h_apply_mrcc_monoexc
   Generate all single excitations of key_in using the bit masks of holes and
   particles.
   Assume N_int is already provided.
@@ -312,15 +312,7 @@ Documentation
   Dressed H with Delta_ij
 
 
-`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L371>`_
+`h_u_0_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L374>`_
-  Computes v_0 = H|u_0>
-  .br
-  n : number of determinants
-  .br
-  H_jj : array of <j|H|j>
-
-
-`h_u_0_mrcc_org <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils/davidson.irp.f#L474>`_
   Computes v_0 = H|u_0>
   .br
   n : number of determinants

plugins/Molden/README.rst

@@ -38,3 +38,39 @@ Needed Modules
 * `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
 * `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`print_mos <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L1>`_
+  Undocumented
+
+
+`write_ao_basis <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L63>`_
+  Undocumented
+
+
+`write_geometry <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L45>`_
+  Undocumented
+
+
+`write_intro_gamess <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L26>`_
+  Undocumented
+
+
+`write_mo_basis <http://github.com/LCPQ/quantum_package/tree/master/plugins/Molden/print_mo.irp.f#L112>`_
+  Undocumented
+

plugins/Perturbation/README.rst

@@ -107,92 +107,92 @@ Documentation
   Undocumented
 
 
-`perturb_buffer_by_mono_delta_rho_one_point <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L896>`_
+perturb_buffer_by_mono_delta_rho_one_point
   Applly pertubration ``delta_rho_one_point`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L476>`_
+perturb_buffer_by_mono_dipole_moment_z
   Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L266>`_
+perturb_buffer_by_mono_epstein_nesbet
   Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L371>`_
+perturb_buffer_by_mono_epstein_nesbet_2x2
   Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L791>`_
+perturb_buffer_by_mono_epstein_nesbet_sc2
   Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L686>`_
+perturb_buffer_by_mono_epstein_nesbet_sc2_no_projected
   Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L581>`_
+perturb_buffer_by_mono_epstein_nesbet_sc2_projected
   Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L161>`_
+perturb_buffer_by_mono_h_core
   Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_by_mono_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L56>`_
+perturb_buffer_by_mono_moller_plesset
   Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_delta_rho_one_point <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L845>`_
+perturb_buffer_delta_rho_one_point
   Applly pertubration ``delta_rho_one_point`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_dipole_moment_z <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L425>`_
+perturb_buffer_dipole_moment_z
   Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L215>`_
+perturb_buffer_epstein_nesbet
   Applly pertubration ``epstein_nesbet`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L320>`_
+perturb_buffer_epstein_nesbet_2x2
   Applly pertubration ``epstein_nesbet_2x2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L740>`_
+perturb_buffer_epstein_nesbet_sc2
   Applly pertubration ``epstein_nesbet_sc2`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_sc2_no_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L635>`_
+perturb_buffer_epstein_nesbet_sc2_no_projected
   Applly pertubration ``epstein_nesbet_sc2_no_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L530>`_
+perturb_buffer_epstein_nesbet_sc2_projected
   Applly pertubration ``epstein_nesbet_sc2_projected`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_h_core <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L110>`_
+perturb_buffer_h_core
   Applly pertubration ``h_core`` to the buffer of determinants generated in the H_apply
   routine.
 
 
-`perturb_buffer_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation/perturbation.irp.f_shell_13#L5>`_
+perturb_buffer_moller_plesset
   Applly pertubration ``moller_plesset`` to the buffer of determinants generated in the H_apply
   routine.
 

plugins/QmcChem/.gitignore

@@ -20,4 +20,5 @@ ezfio_interface.irp.f
 irpf90.make
 irpf90_entities
 save_for_qmcchem
-tags
+tags
+target_pt2_qmc

plugins/QmcChem/README.rst

@@ -66,6 +66,14 @@ Documentation
   title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
 
 
+`compute_energy <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/target_pt2_qmc.irp.f#L80>`_
+  Compute an energy when a threshold is applied
+
+
+`e_curve <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/target_pt2_qmc.irp.f#L1>`_
+  Undocumented
+
+
 `mo_pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/plugins/QmcChem/pot_ao_pseudo_ints.irp.f#L56>`_
   Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
   .br

plugins/QmcChem/target_pt2_qmc.irp.f

@@ -0,0 +1,121 @@
+program e_curve
+  use bitmasks
+  implicit none
+  integer :: i,j,k, nab, m, l, n_up, n_dn, n
+  double precision :: norm, E, hij, num, ci, cj
+  integer, allocatable :: iorder(:)
+  double precision , allocatable :: norm_sort(:), psi_bilinear_matrix_values_save(:)
+  nab = n_det_alpha_unique+n_det_beta_unique
+
+  allocate ( norm_sort(0:nab), iorder(0:nab), psi_bilinear_matrix_values_save(N_det) )
+
+
+  norm_sort(0) = 0.d0
+  iorder(0) = 0
+  do i=1,n_det_alpha_unique
+   norm_sort(i) = det_alpha_norm(i)
+   iorder(i) = i
+  enddo
+  
+  do i=1,n_det_beta_unique
+   norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+   iorder(i+n_det_alpha_unique) = -i
+  enddo
+  
+  call dsort(norm_sort(1),iorder(1),nab)
+
+  if (.not.read_wf) then
+    stop 'Please set read_wf to true'
+  endif
+
+  psi_bilinear_matrix_values_save = psi_bilinear_matrix_values(:,1)
+  print *,  '=========================================================='
+  print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)',  'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
+  print *,  '=========================================================='
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+
+  double precision :: thresh, E_min, E_max, E_prev
+  thresh = 0.d0
+  call compute_energy(psi_bilinear_matrix_values_save,E_max,m,norm)
+  call i_h_j(psi_det_sorted(1,1,1), psi_det_sorted(1,1,1), N_int, E_min)
+  print *,  E_min, E_max
+
+  n_up = nab
+  n_dn = 0
+  do while (n_up > n_dn)
+    n = n_dn + (n_up-n_dn)/2
+    psi_bilinear_matrix_values_save = psi_bilinear_matrix_values(:,1)
+    do j=1,n
+      i = iorder(j)
+      if (i<0) then
+        do k=1,n_det
+          if (psi_bilinear_matrix_columns(k) == -i) then
+            psi_bilinear_matrix_values_save(k) = 0.d0
+          endif
+        enddo
+      else
+        do k=1,n_det
+          if (psi_bilinear_matrix_rows(k) == i) then
+            psi_bilinear_matrix_values_save(k) = 0.d0
+          endif
+        enddo
+      endif
+    enddo
+    call compute_energy(psi_bilinear_matrix_values_save,E,m,norm)
+    print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)',  norm_sort(n), m, &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * m ) / &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E
+    if (E < target_energy) then
+       n_dn = n+1
+    else
+       n_up = n
+    endif
+  enddo
+  print *,  '=========================================================='
+  print *,  norm_sort(n), target_energy
+
+  deallocate (iorder, norm_sort, psi_bilinear_matrix_values_save)
+end
+
+subroutine compute_energy(psi_bilinear_matrix_values_save, E, m, norm)
+  implicit none
+  BEGIN_DOC
+  ! Compute an energy when a threshold is applied
+  END_DOC
+  double precision, intent(in) :: psi_bilinear_matrix_values_save(n_det)
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  integer :: i,j, k, l, m
+  double precision :: num, norm, ci, cj, hij, E
+
+  
+    num = 0.d0
+    norm = 0.d0
+    m = 0
+    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
+    allocate( det_i(N_int,2), det_j(N_int,2))
+    !$OMP DO 
+    do k=1,n_det
+      if (psi_bilinear_matrix_values_save(k) == 0.d0) then
+        cycle
+      endif
+      ci = psi_bilinear_matrix_values_save(k)
+      det_i(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(k))
+      det_i(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(k))
+      do l=1,n_det
+        if (psi_bilinear_matrix_values_save(l) == 0.d0) then
+          cycle
+        endif
+        cj = psi_bilinear_matrix_values_save(l)
+        det_j(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(l))
+        det_j(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(l))
+        call i_h_j(det_i, det_j, N_int, hij)
+        num = num + ci*cj*hij
+      enddo
+     norm = norm + ci*ci
+     m = m+1
+    enddo
+    !$OMP END DO
+    deallocate (det_i,det_j)
+    !$OMP END PARALLEL
+    E = num / norm + nuclear_repulsion
+end

plugins/loc_cele/README.rst

@@ -20,3 +20,25 @@ Needed Modules
 
 * `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
 
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+.. image:: tree_dependency.png
+
+* `MO_Basis <http://github.com/LCPQ/quantum_package/tree/master/src/MO_Basis>`_
+
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+
+
+`loc_rasorb <http://github.com/LCPQ/quantum_package/tree/master/plugins/loc_cele/loc_cele.irp.f#L1>`_
+  This program performs a localization of the active orbitals
+  of a CASSCF wavefunction, reading the orbitals from a RASORB
+  file of molcas.
+  id1=max is the number of MO in a given symmetry.
+

plugins/loc_cele/loc_cele.irp.f

@@ -280,7 +280,7 @@
 
       enddo !big loop over symmetry
 
-      10 format (4E18.12)
+      10 format (4E20.12)
 
 
 !  Now we copyt the newcmo into the mo_coef

scripts/compilation/qp_create_ninja.py

@@ -707,7 +707,7 @@ def ninja_dot_tree_rule():
     l_string = [
         "rule build_dot_tree", "   command = {0}".format(" ; ".join(l_cmd)),
         "   generator = 1",
-        "   description = Generate Png representtion of the Tree Dependencies of $module_rel",
+        "   description = Generating Png representation of the Tree Dependencies of $module_rel",
         ""
     ]
 

scripts/generate_h_apply.py

@@ -99,7 +99,7 @@ class H_apply(object):
   deallocate(H_jj,iorder)
     """
 
-    s["size_max"] = str(1024*128) 
+    s["size_max"] = "256"
     s["copy_buffer"] = """call copy_H_apply_buffer_to_wf
   if (s2_eig) then
     call make_s2_eigenfunction
@@ -266,7 +266,7 @@ class H_apply(object):
       double precision, intent(inout) :: select_max_out"""
 
       self.data["params_post"] += ", select_max(min(i_generator,size(select_max,1)))"
-      self.data["size_max"] = str(1024*128) 
+      self.data["size_max"] = "256"
       self.data["copy_buffer"] = """
       call copy_H_apply_buffer_to_wf
       if (s2_eig) then

scripts/module/create_gitignore.sh

@@ -13,7 +13,6 @@ then
 fi
 source ${QP_ROOT}/scripts/qp_include.sh
 
-
 function do_gitingore()
 { 
   cat << EOF > .gitignore

scripts/module/module_handler.py

@@ -180,6 +180,11 @@ class ModuleHandler():
     def create_png(self, l_module):
         """Create the png of the dependency tree for a l_module"""
 
+        # Don't update if we are not in the main repository
+        from is_master_repository import is_master_repository
+        if not is_master_repository:
+            return
+
         basename = "tree_dependency"
         path = '{0}.png'.format(basename)
 
@@ -289,6 +294,12 @@ if __name__ == '__main__':
                         pass
 
             if arguments["create_git_ignore"]:
+
+                # Don't update if we are not in the main repository
+                from is_master_repository import is_master_repository
+                if not is_master_repository:
+                    sys.exit()
+
                 path = os.path.join(module_abs, ".gitignore")
 
                 with open(path, "w+") as f:

scripts/module/qp_update_readme.py

@@ -44,7 +44,7 @@ def get_url(path_module_rel):
     elif is_module(path_module_rel):
         url = "http://github.com/LCPQ/quantum_package/tree/master/src"
     else:
-        print "{0} Is not a valide module nor plugin".format(path_module_rel)
+        print "{0} Is not a valid module nor plugin".format(path_module_rel)
         sys.exit(1)
 
     return os.path.join(url, path_module_rel)
@@ -155,20 +155,32 @@ def update_documentation(d_readmen, root_module):
         l_doc = []
 
         for irp in d_info[path]:
-
             url = os.path.join(get_url(os.path.basename(path)), irp.file)
             doc = extract_doc(root_module, irp.provider)
 
-            l_doc += ["`{0} <{1}#L{2}>`_".format(irp.provider, url, irp.line),
+            if ".irp.f_shell_" in irp.file:
-                      doc,
+                l_doc += ["{0}".format(irp.provider),
-                      ""]
+                          doc,
+                          ""]
+            else:
+                l_doc += ["`{0} <{1}#L{2}>`_".format(irp.provider, url, irp.line),
+                          doc,
+                          ""]
 
         l_doc_section = [D_KEY["documentation"], '',
                          "\n".join(l_doc)]
 
         d_readme[path]["documentation"] = "\n".join(l_doc_section)
 
+
 if __name__ == '__main__':
+
+    # Update documentation only if the remote repository is
+    # the main repository
+    from is_master_repository import is_master_repository
+    if not is_master_repository:
+        sys.exit(0)
+
     arguments = docopt(__doc__)
 
     if arguments["--root_module"]:
@@ -188,8 +200,8 @@ if __name__ == '__main__':
         fetch_splitted_data(d_readme, l_module_readme)
     except IOError:
         print l_module_readme, "is not a module and/or",
-        print "have not a `README.rst` file inside"
+        print "has not a `README.rst` file inside"
-        print "Abort..."
+        print "Aborting..."
         sys.exit(1)
 
     update_needed(d_readme)

scripts/utility/is_master_repository.py

@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+
+import subprocess
+pipe = subprocess.Popen("git config --local --get remote.origin.url", \
+              shell=True, stdout=subprocess.PIPE)
+result = pipe.stdout.read()
+is_master_repository = "LCPQ/quantum_package" in result
+
+if __name__ == "__main__":
+    import sys
+    if is_master_repository:
+        sys.exit(0)
+    else:
+        sys.exit(-1)

src/Determinants/EZFIO.cfg

@@ -53,9 +53,10 @@ interface: ezfio,provider,ocaml
 default:        0.999
 
 [n_states_diag]
-type: integer
+type: States_number
 doc: n_states_diag
-interface: ezfio,provider
+default: 1
+interface: ezfio,provider,ocaml
 
 [n_int]
 interface: ezfio
@@ -89,24 +90,25 @@ doc: psi_det
 type: integer*8
 size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det)
 
-[det_num]
-interface: ezfio,provider
-doc: det_num
-type: integer
-
 [det_occ]
 interface: ezfio,provider
 doc: det_occ
 type: integer          
-size:  (electrons.elec_alpha_num,determinants.det_num,2)
+size:  (electrons.elec_alpha_num,determinants.n_det,2)
 
 [det_coef]
 interface: ezfio,provider
 doc: det_coef
 type: double precision 
-size:  (determinants.det_num)
+size:  (determinants.n_det)
 
 [expected_s2]
 interface: ezfio,provider
-doc: expcted_s2
+doc: Expected value of S^2
-type: double precision 
+type: double precision
+
+[target_energy]
+interface: ezfio,provider,ocaml
+doc: Energy that should be obtained when truncating the wave function (optional)
+type: Energy
+default: 0.

src/Determinants/README.rst

@@ -113,16 +113,16 @@ Documentation
   After calling this subroutine, N_det, psi_det and psi_coef need to be touched
 
 
-`create_wf_of_psi_bilinear_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L417>`_
+`create_wf_of_psi_bilinear_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L445>`_
   Generate a wave function containing all possible products
   of alpha and beta determinants
 
 
-`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L563>`_
+`davidson_converged <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L565>`_
   True if the Davidson algorithm is converged
 
 
-`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L553>`_
+`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L555>`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
 
@@ -145,7 +145,7 @@ Documentation
   Initial guess vectors are not necessarily orthonormal
 
 
-`davidson_diag_hjj <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L239>`_
+`davidson_diag_hjj <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L241>`_
   Davidson diagonalization with specific diagonal elements of the H matrix
   .br
   H_jj : specific diagonal H matrix elements to diagonalize de Davidson
@@ -174,7 +174,7 @@ Documentation
   Max number of Davidson sizes
 
 
-`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L554>`_
+`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L556>`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
 
@@ -186,22 +186,26 @@ Documentation
   degree : Degree of excitation
 
 
+`det_alpha_norm <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L353>`_
+  Norm of the alpha and beta spin determinants in the wave function:
+  .br
+  ||Da||_i \sum_j C_{ij}**2
+
+
+`det_beta_norm <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L354>`_
+  Norm of the alpha and beta spin determinants in the wave function:
+  .br
+  ||Da||_i \sum_j C_{ij}**2
+
+
 `det_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L138>`_
   det_coef
 
 
-`det_connections <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1380>`_
-  Build connection proxy between determinants
-
-
 `det_inf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L69>`_
   Undocumented
 
 
-`det_num <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L248>`_
-  det_num
-
-
 `det_occ <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L226>`_
   det_occ
 
@@ -278,26 +282,6 @@ Documentation
   idx(0) is the number of determinants that interact with key1
 
 
-`filter_connected_davidson <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L257>`_
-  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_davidson_shortcut <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L193>`_
-  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_davidson_warp <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L102>`_
   Filters out the determinants that are not connected by H
   returns the array idx which contains the index of the
@@ -308,7 +292,7 @@ Documentation
   key1 should come from psi_det_sorted_ab.
 
 
-`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L389>`_
+`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L390>`_
   returns the array idx which contains the index of the
   .br
   determinants in the array key1 that interact
@@ -318,7 +302,7 @@ Documentation
   idx(0) is the number of determinants that interact with key1
 
 
-`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L488>`_
+`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L489>`_
   standard filter_connected_i_H_psi but returns in addition
   .br
   the array of the index of the non connected determinants to key1
@@ -330,7 +314,7 @@ Documentation
   to repeat the excitations
 
 
-`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L472>`_
+`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L500>`_
   Create a wave function from all possible alpha x beta determinants
 
 
@@ -408,14 +392,6 @@ Documentation
   H_jj : array of <j|H|j>
 
 
-`h_u_0_org <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1303>`_
-  Computes v_0 = H|u_0>
-  .br
-  n : number of determinants
-  .br
-  H_jj : array of <j|H|j>
-
-
 `i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L356>`_
   Returns <i|H|j> where i and j are determinants
 
@@ -507,10 +483,6 @@ Documentation
   Energy of the reference bitmask used in Slater rules
 
 
-`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1372>`_
-  Number of integers to represent the connections between determinants
-
-
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L3>`_
   Number of determinants in the wave function
 
@@ -624,22 +596,22 @@ Documentation
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
 
-`psi_bilinear_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L400>`_
+`psi_bilinear_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L428>`_
   Coefficient matrix if the wave function is expressed in a bilinear form :
   D_a^t C D_b
 
 
-`psi_bilinear_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L362>`_
+`psi_bilinear_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L390>`_
   Sparse coefficient matrix if the wave function is expressed in a bilinear form :
   D_a^t C D_b
 
 
-`psi_bilinear_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L361>`_
+`psi_bilinear_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L389>`_
   Sparse coefficient matrix if the wave function is expressed in a bilinear form :
   D_a^t C D_b
 
 
-`psi_bilinear_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L360>`_
+`psi_bilinear_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L388>`_
   Sparse coefficient matrix if the wave function is expressed in a bilinear form :
   D_a^t C D_b
 
@@ -793,7 +765,7 @@ Documentation
   Reads the determinants from the EZFIO file
 
 
-`read_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L160>`_
+`read_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L116>`_
   If true, read the wave function from the EZFIO file
 
 
@@ -818,7 +790,7 @@ Documentation
   Undocumented
 
 
-`s2_eig <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L116>`_
+`s2_eig <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L248>`_
   Force the wave function to be an eigenfunction of S^2
 
 
@@ -872,11 +844,11 @@ Documentation
   for a given couple of hole/particle excitations i.
 
 
-`sort_dets_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L195>`_
+`sort_dets_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L197>`_
   Undocumented
 
 
-`sort_dets_ab_v <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L148>`_
+`sort_dets_ab_v <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L145>`_
   Undocumented
 
 
@@ -909,6 +881,10 @@ Documentation
   Undocumented
 
 
+`target_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ezfio_interface.irp.f#L160>`_
+  Energy that should be obtained when truncating the wave function (optional)
+
+
 `threshold_convergence_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L18>`_
   convergence of the correlation energy of SC2 iterations
 

src/Determinants/spindeterminants.irp.f

@@ -350,6 +350,34 @@ subroutine write_spindeterminants
   
 end
 
+ BEGIN_PROVIDER [ double precision, det_alpha_norm, (N_det_alpha_unique) ]
+&BEGIN_PROVIDER [ double precision, det_beta_norm, (N_det_beta_unique) ]
+ implicit none
+ BEGIN_DOC
+ ! Norm of the alpha and beta spin determinants in the wave function:
+ !
+ ! ||Da||_i \sum_j C_{ij}**2
+ END_DOC
+
+ integer :: i,j,k,l
+ double precision :: f
+
+ det_alpha_norm = 0.d0
+ det_beta_norm  = 0.d0
+ do k=1,N_det
+   i = psi_bilinear_matrix_rows(k)
+   j = psi_bilinear_matrix_columns(k)
+   do l=1,N_states
+    f = psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l)
+   enddo
+   det_alpha_norm(i) += f
+   det_beta_norm(j)  += f
+ enddo
+ det_alpha_norm = det_alpha_norm / dble(N_states)
+ det_beta_norm = det_beta_norm / dble(N_states)
+
+END_PROVIDER
+
 
 !==============================================================================!
 !                                                                              !

src/Ezfio_files/README.rst

@@ -195,103 +195,147 @@ Documentation
   .br
 
 
-`output_ao_basis <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L1>`_
+output_ao_basis
   Output file for AO_Basis
 
 
-`output_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L21>`_
+output_bitmask
   Output file for Bitmask
 
 
-`output_cas_sd <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L41>`_
+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 <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f#L2>`_
   Initial CPU and wall times when printing in the output files
 
 
-`output_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L61>`_
+output_ddci_selected
+  Output file for DDCI_selected
+
+
+output_determinants
   Output file for Determinants
 
 
-`output_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L81>`_
+output_electrons
   Output file for Electrons
 
 
-`output_ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L101>`_
+output_ezfio_files
   Output file for Ezfio_files
 
 
-`output_full_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L121>`_
+output_fcidump
+  Output file for FCIdump
+
+
+output_full_ci
   Output file for Full_CI
 
 
-`output_generators_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L141>`_
+output_generators_cas
   Output file for Generators_CAS
 
 
-`output_generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L161>`_
+output_generators_full
   Output file for Generators_full
 
 
-`output_hartree_fock <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L181>`_
+output_hartree_fock
   Output file for Hartree_Fock
 
 
-`output_integrals_bielec <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L201>`_
+output_integrals_bielec
   Output file for Integrals_Bielec
 
 
-`output_integrals_monoelec <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L221>`_
+output_integrals_monoelec
   Output file for Integrals_Monoelec
 
 
-`output_mo_basis <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L241>`_
+output_loc_cele
+  Output file for loc_cele
+
+
+output_mo_basis
   Output file for MO_Basis
 
 
-`output_moguess <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L261>`_
+output_moguess
   Output file for MOGuess
 
 
-`output_mrcc_cassd <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L281>`_
+output_molden
+  Output file for Molden
+
+
+output_mp2
+  Output file for MP2
+
+
+output_mrcc_cassd
   Output file for MRCC_CASSD
 
 
-`output_mrcc_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L301>`_
+output_mrcc_utils
   Output file for MRCC_Utils
 
 
-`output_nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L321>`_
+output_nuclei
   Output file for Nuclei
 
 
-`output_perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L341>`_
+output_perturbation
   Output file for Perturbation
 
 
-`output_properties <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L361>`_
+output_properties
   Output file for Properties
 
 
-`output_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L381>`_
+output_pseudo
   Output file for Pseudo
 
 
-`output_psiref_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L401>`_
+output_psiref_cas
   Output file for Psiref_CAS
 
 
-`output_psiref_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L421>`_
+output_psiref_utils
   Output file for Psiref_Utils
 
 
-`output_selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L441>`_
+output_qmcchem
+  Output file for QmcChem
+
+
+output_selectors_full
   Output file for Selectors_full
 
 
-`output_utils <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files/output.irp.f_shell_40#L461>`_
+output_singlerefmethod
+  Output file for SingleRefMethod
+
+
+output_utils
   Output file for Utils
 
 

src/Integrals_Bielec/ao_bi_integrals.irp.f

@@ -481,7 +481,7 @@ IRP_ENDIF COARRAY
   ao_bielec_integrals_in_map = .True.
   if (write_ao_integrals) then
     call dump_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin')
-    call ezfio_set_integrals_bielec_disk_access_ao_integrals(.True.)
+    call ezfio_set_integrals_bielec_disk_access_ao_integrals("Read")
   endif
   
 END_PROVIDER

src/Integrals_Bielec/mo_bi_integrals.irp.f

@@ -312,7 +312,7 @@ IRP_ENDIF
   
   if (write_mo_integrals) then
     call dump_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin')
-    call ezfio_set_integrals_bielec_disk_access_mo_integrals(.True.)
+    call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read")
   endif
   
   

src/Integrals_Monoelec/README.rst

@@ -93,6 +93,11 @@ Documentation
   : sum of the kinetic and nuclear electronic potential
 
 
+`ao_mono_elec_integral_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/ao_mono_ints.irp.f#L2>`_
+  array of the mono electronic hamiltonian on the AOs basis
+  : sum of the kinetic and nuclear electronic potential
+
+
 `ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_ints.irp.f#L1>`_
   interaction nuclear electron
 

src/Integrals_Monoelec/ao_mono_ints.irp.f

@@ -1,4 +1,5 @@
-BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
+ BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
+&BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_diag,(ao_num)]
   implicit none
   integer :: i,j,n,l
   BEGIN_DOC
@@ -6,9 +7,11 @@ BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
  ! : sum of the kinetic and nuclear electronic potential 
   END_DOC
   do j = 1, ao_num
+   !DIR$ VECTOR ALIGNED
    do i = 1, ao_num
     ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_pseudo_integral(i,j)
    enddo
+   ao_mono_elec_integral_diag(j) = ao_mono_elec_integral(j,j)
   enddo
 END_PROVIDER
 

src/Integrals_Monoelec/pseudopot.f90

@@ -214,13 +214,14 @@ integer :: l,k, nkl_max
 ! |_) | (_|   (_| |  | (_| \/ 
 !        _|                /  
 
-double precision, allocatable :: array_coefs_A(:,:,:)
+double precision, allocatable :: array_coefs_A(:,:)
-double precision, allocatable :: array_coefs_B(:,:,:)
+double precision, allocatable :: array_coefs_B(:,:)
 
 double precision, allocatable :: array_R(:,:,:,:,:)
 double precision, allocatable :: array_I_A(:,:,:,:,:)
 double precision, allocatable :: array_I_B(:,:,:,:,:)
 
+double precision :: f1, f2, f3
 
 !  _                
 ! /   _. |  _     | 
@@ -255,14 +256,14 @@ nkl_max=4
 ! A l l o c a t e !
 !=!=!=!=!=!=!=!=!=!
 
-allocate (array_coefs_A(0:ntot,0:ntot,0:ntot))
+allocate (array_coefs_A(0:ntot,3))
-allocate (array_coefs_B(0:ntot,0:ntot,0:ntot))
+allocate (array_coefs_B(0:ntot,3))
 
-allocate (array_R(0:ntot+nkl_max,kmax,0:lmax,0:lmax+ntot,0:lmax+ntot))
+allocate (array_R(kmax,0:ntot+nkl_max,0:lmax,0:lmax+ntot,0:lmax+ntot))
 
-allocate (array_I_A(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot))
+allocate (array_I_A(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot))
-
+                                                       
-allocate (array_I_B(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot))
+allocate (array_I_B(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot))
 
 if(ac.eq.0.d0.and.bc.eq.0.d0)then
 
@@ -310,108 +311,110 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
  phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc)
 
 
-
+ do lambdap=0,lmax+ntotB
-
- do ktot=0,ntotA+ntotB+nkl_max
    do lambda=0,lmax+ntotA
-   do lambdap=0,lmax+ntotB
+     do l=0,lmax
+       do ktot=0,ntotA+ntotB+nkl_max
         do k=1,kmax
-          do l=0,lmax
+           array_R(k,ktot,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
-          array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg)
+         enddo
-          enddo
+       enddo
-        enddo
+     enddo
-   enddo
    enddo
  enddo
-
+ 
  do k1=0,n_a(1)
+   array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)
+ enddo
  do k2=0,n_a(2)
+   array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)
+ enddo
  do k3=0,n_a(3)
-  array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
+   array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)
-                          *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)
- enddo
- enddo
  enddo
 
  do k1p=0,n_b(1)
+   array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p)
+ enddo
  do k2p=0,n_b(2)
+   array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p)
+ enddo
  do k3p=0,n_b(3)
-  array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)  &
+   array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p)
-                             *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
  enddo
- enddo
+ 
- enddo
-
  !=!=!=!=!=!=!=!
  ! c a l c u l !
  !=!=!=!=!=!=!=!
 
  accu=0.d0
  do l=0,lmax
-  do m=-l,l
+   do m=-l,l
-
+     
-       do lambda=0,l+ntotA
+     do k3=0,n_a(3)
-        do mu=-lambda,lambda
+       do k2=0,n_a(2)
-           do k1=0,n_a(1)
+         do k1=0,n_a(1)
-           do k2=0,n_a(2)
+           do lambda=0,l+ntotA
-           do k3=0,n_a(3)
+             do mu=-lambda,lambda
-                array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
+               array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
+             enddo
            enddo
-           enddo
+         enddo
-           enddo
-        enddo
        enddo
-
+     enddo
-       do lambdap=0,l+ntotB
+     
-        do mup=-lambdap,lambdap
+     do k3p=0,n_b(3)
-           do k1p=0,n_b(1)
+       do k2p=0,n_b(2)
-           do k2p=0,n_b(2)
+         do k1p=0,n_b(1)
-           do k3p=0,n_b(3)
+           do lambdap=0,l+ntotB
-                  array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
+             do mup=-lambdap,lambdap
-            enddo
+               array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
-            enddo
+             enddo
-            enddo
+           enddo
-        enddo
+         enddo
        enddo
-
+     enddo
-       do lambda=0,l+ntotA
+     
-        do mu=-lambda,lambda
+     do k3=0,n_a(3)
-
+       do k2=0,n_a(2)
-          do k1=0,n_a(1)
+         do k1=0,n_a(1)
-          do k2=0,n_a(2)
+           
-          do k3=0,n_a(3)
+           do lambda=0,l+ntotA
-
+             do mu=-lambda,lambda
-              prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3)
+               
-
+               prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*array_I_A(mu,lambda,k1,k2,k3)
-              do lambdap=0,l+ntotB
+               
-               do mup=-lambdap,lambdap
+               
-
+               do k3p=0,n_b(3)
-                    do k1p=0,n_b(1)
+                 do k2p=0,n_b(2)
-                    do k2p=0,n_b(2)
+                   do k1p=0,n_b(1)
-                    do k3p=0,n_b(3)
+                     do lambdap=0,l+ntotB
-
+                       do mup=-lambdap,lambdap
-                        prodp=ylm(lambdap,mup,theta_BC0,phi_BC0)*array_coefs_B(k1p,k2p,k3p)*array_I_B(lambdap,mup,k1p,k2p,k3p)
+                         
-
+                         prodp=prod*ylm(lambdap,mup,theta_BC0,phi_BC0)* &
-                        do k=1,kmax
+                            array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)* &
-                            ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l)
+                            array_I_B(mup,lambdap,k1p,k2p,k3p)
-                            accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,lambdap)
+                         
-                        enddo
+                         do k=1,kmax
-
+                           ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l)
-                    enddo
+                           accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,lambdap)
-                    enddo
+                         enddo
-                    enddo
+                         
-
+                       enddo
+                     enddo
+                   enddo
+                   
+                 enddo
                enddo
-              enddo
+               
-
+             enddo
-          enddo
+           enddo
-          enddo
+         enddo
-          enddo
+         
-
-        enddo
        enddo
-
+     enddo
-  enddo
+     
+   enddo
  enddo
 
  !=!=!=!=!
@@ -434,24 +437,24 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
  breal=2.d0*g_b*bc
  freal=dexp(-g_a*ac**2-g_b*bc**2)
 
- do ktot=0,ntotA+ntotB+nkl_max
+ do lambdap=0,lmax+ntotB
-  do lambdap=0,lmax+ntotB
+   do l=0,lmax
-    do k=1,kmax
+     do ktot=0,ntotA+ntotB+nkl_max
-      do l=0,lmax
+       do k=1,kmax
-        array_R(ktot,k,l,0,lambdap)=  int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
+         array_R(k,ktot,l,0,lambdap)=  int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg)
        enddo
-    enddo
+     enddo
-  enddo
+   enddo
  enddo
 
  do k1p=0,n_b(1)
+   array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p)
+ enddo
  do k2p=0,n_b(2)
+   array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p)
+ enddo
  do k3p=0,n_b(3)
-
+   array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p)
-  array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)  &
-                             *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
- enddo
- enddo
  enddo
 
  !=!=!=!=!=!=!=!
@@ -460,43 +463,43 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
 
  accu=0.d0
  do l=0,lmax
-  do m=-l,l
+   do m=-l,l
-
+     
-    do lambdap=0,l+ntotB
+     do k3p=0,n_b(3)
-      do mup=-lambdap,lambdap
+       do k2p=0,n_b(2)
-        do k1p=0,n_b(1)
-        do k2p=0,n_b(2)
-        do k3p=0,n_b(3)
-          array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
-        enddo
-        enddo
-        enddo
-      enddo
-     enddo
-    
-    prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
-    
-    do lambdap=0,l+ntotB
-      do mup=-lambdap,lambdap
          do k1p=0,n_b(1)
-         do k2p=0,n_b(2)
+           do lambdap=0,l+ntotB
-         do k3p=0,n_b(3)
+             do mup=-lambdap,lambdap
-    
+               array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p)
-          prodp=array_coefs_B(k1p,k2p,k3p)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(lambdap,mup,k1p,k2p,k3p)
+             enddo
-    
+           enddo
-          do k=1,kmax
-
-            ktot=ntotA+k1p+k2p+k3p+n_kl(k,l)
-            accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,0,lambdap)
-
-          enddo
-
-         enddo
-         enddo
          enddo
+       enddo
      enddo
-    enddo
+     
-  enddo
+     prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
+     
+     do k3p=0,n_b(3)
+       do k2p=0,n_b(2)
+         do k1p=0,n_b(1)
+           do lambdap=0,l+ntotB
+             do mup=-lambdap,lambdap
+               
+               prodp=prod*array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(mup,lambdap,k1p,k2p,k3p)
+               
+               do k=1,kmax
+                 
+                 ktot=ntotA+k1p+k2p+k3p+n_kl(k,l)
+                 accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,0,lambdap)
+                 
+               enddo
+               
+             enddo
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
  enddo
 
  !=!=!=!=!
@@ -519,26 +522,24 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
  breal=2.d0*g_b*bc
  freal=dexp(-g_a*ac**2-g_b*bc**2)
 
- do ktot=0,ntotA+ntotB+nkl_max
+ do lambda=0,lmax+ntotA
-  do lambda=0,lmax+ntotA
+   do l=0,lmax
-    do k=1,kmax
+     do ktot=0,ntotA+ntotB+nkl_max
-      do l=0,lmax
+       do k=1,kmax
-
+         array_R(k,ktot,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
-      array_R(ktot,k,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg)
+       enddo
-      enddo
+     enddo
-    enddo
+   enddo
-  enddo
  enddo
 
  do k1=0,n_a(1)
+   array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)
+ enddo
  do k2=0,n_a(2)
+   array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)
+ enddo
  do k3=0,n_a(3)
-
+   array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)
-  array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
-                          *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)
-
- enddo
- enddo
  enddo
 
  !=!=!=!=!=!=!=!
@@ -549,36 +550,36 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
  do l=0,lmax
   do m=-l,l
 
-    do lambda=0,l+ntotA
+    do k3=0,n_a(3)
-      do mu=-lambda,lambda
+      do k2=0,n_a(2)
         do k1=0,n_a(1)
-        do k2=0,n_a(2)
+          do lambda=0,l+ntotA
-        do k3=0,n_a(3)
+            do mu=-lambda,lambda
-          array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
+              array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
-        enddo
+            enddo
-        enddo
+          enddo
         enddo
       enddo
     enddo
 
-    do lambda=0,l+ntotA
+    do k3=0,n_a(3)
-      do mu=-lambda,lambda
+      do k2=0,n_a(2)
         do k1=0,n_a(1)
-        do k2=0,n_a(2)
+          do lambda=0,l+ntotA
-        do k3=0,n_a(3)
+            do mu=-lambda,lambda
-
+              
-          prod=array_coefs_A(k1,k2,k3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(lambda,mu,k1,k2,k3)
+              prod=array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(mu,lambda,k1,k2,k3)
-          prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
+              prodp=prod*bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))
-
+              
-          do k=1,kmax
+              do k=1,kmax
-            ktot=k1+k2+k3+ntotB+n_kl(k,l)
+                ktot=k1+k2+k3+ntotB+n_kl(k,l)
-            accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,0)
+                accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,0)
+              enddo
+              
+            enddo
           enddo
-
         enddo
-        enddo
+      enddo
-        enddo
-     enddo
     enddo
 
   enddo
@@ -850,22 +851,27 @@ implicit none
 integer lambda,mu,l,m,k1,k2,k3
 integer k,i,kp,ip
 double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm
+double precision sgn, sgnp
 pi=dacos(-1.d0)
 
 if(mu.gt.0.and.m.gt.0)then
 sum=0.d0
-factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
+factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
-factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
+factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
+sgn = 1.d0
 do k=0,mu/2
  do i=0,lambda-mu
+  sgnp = 1.d0
   do kp=0,m/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3)
    enddo
+   sgnp = -sgnp
   enddo
  enddo
+ sgn = -sgn
 enddo
 bigI=sum
 return
@@ -888,15 +894,17 @@ endif
 
 if(mu.eq.0.and.m.gt.0)then
 factor1=dsqrt((2*lambda+1)/(4.d0*pi))
-factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
+factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
 sum=0.d0
 do i=0,lambda
+ sgnp = 1.d0
  do kp=0,m/2
   do ip=0,l-m
    cylm=factor1*coef_pm(lambda,i)
-   cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+   cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
    sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3)
   enddo
+  sgnp = -sgnp
  enddo
 enddo
 bigI=sum
@@ -905,16 +913,18 @@ endif
 
 if(mu.gt.0.and.m.eq.0)then
 sum=0.d0
-factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
+factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
 factor2=dsqrt((2*l+1)/(4.d0*pi))
+sgn = 1.d0
 do k=0,mu/2
  do i=0,lambda-mu
   do ip=0,l
-   cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+   cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
    cylmp=factor2*coef_pm(l,ip)
    sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3)
   enddo
  enddo
+ sgn = -sgn
 enddo
 bigI=sum
 return
@@ -923,19 +933,23 @@ endif
 if(mu.lt.0.and.m.lt.0)then
 mu=-mu
 m=-m
-factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
+factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
-factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
+factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
 sum=0.d0
+sgn = 1.d0
 do k=0,(mu-1)/2
  do i=0,lambda-mu
+  sgnp = 1.d0
   do kp=0,(m-1)/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3)
    enddo
+   sgnp = -sgnp
   enddo
  enddo
+ sgn = -sgn
 enddo
 mu=-mu
 m=-m
@@ -946,15 +960,17 @@ endif
 if(mu.eq.0.and.m.lt.0)then
 m=-m
 factor1=dsqrt((2*lambda+1)/(4.d0*pi))
-factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
+factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
 sum=0.d0
 do i=0,lambda
+ sgnp = 1.d0
  do kp=0,(m-1)/2
   do ip=0,l-m
    cylm=factor1*coef_pm(lambda,i)
-   cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+   cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
    sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3)
   enddo
+  sgnp = -sgnp
  enddo
 enddo
 m=-m
@@ -965,16 +981,18 @@ endif
 if(mu.lt.0.and.m.eq.0)then
 sum=0.d0
 mu=-mu
-factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
+factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
 factor2=dsqrt((2*l+1)/(4.d0*pi))
+sgn = 1.d0
 do k=0,(mu-1)/2
  do i=0,lambda-mu
    do ip=0,l
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
     cylmp=factor2*coef_pm(l,ip)
     sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
    enddo
  enddo
+ sgn = -sgn
 enddo
 mu=-mu
 bigI=sum
@@ -983,19 +1001,23 @@ endif
 
 if(mu.gt.0.and.m.lt.0)then
 sum=0.d0
-factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
+factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
-factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
+factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
 m=-m
+sgn=1.d0
 do k=0,mu/2
  do i=0,lambda-mu
+  sgnp=1.d0
   do kp=0,(m-1)/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm =sgn *factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3)
    enddo
+   sgnp = -sgnp
   enddo
  enddo
+ sgn = -sgn
 enddo
 m=-m
 bigI=sum
@@ -1004,19 +1026,23 @@ endif
 
 if(mu.lt.0.and.m.gt.0)then
 mu=-mu
-factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu))))
+factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu))))
-factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m))))
+factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m))))
 sum=0.d0
+sgn = 1.d0
 do k=0,(mu-1)/2
  do i=0,lambda-mu
+  sgnp = 1.d0
   do kp=0,m/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=sgn*factor1  *binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3)
    enddo
+   sgnp = -sgnp
   enddo
  enddo
+ sgn = -sgn
 enddo
 bigI=sum
 mu=-mu
@@ -1128,28 +1154,49 @@ end
 !
         IMPLICIT DOUBLE PRECISION (P,X)
         DIMENSION PM(0:MM,0:(N+1))
-        DO 10 I=0,N
+        DOUBLE PRECISION, SAVE ::  INVERSE(100) = 0.D0
-        DO 10 J=0,M
+        DOUBLE PRECISION       :: LS, II, JJ
-10         PM(J,I)=0.0D0
+        IF (INVERSE(1) == 0.d0) THEN
+          DO I=1,100
+            INVERSE(I) = 1.D0/DBLE(I)
+          ENDDO
+        ENDIF
+        DO I=0,N
+          DO J=0,M
+            PM(J,I)=0.0D0
+          ENDDO
+        ENDDO
         PM(0,0)=1.0D0
         IF (DABS(X).EQ.1.0D0) THEN
-           DO 15 I=1,N
+           DO I=1,N
-15            PM(0,I)=X**I
+              PM(0,I)=X**I
+           ENDDO
            RETURN
         ENDIF
-        LS=1
+        LS=1.D0
-        IF (DABS(X).GT.1.0D0) LS=-1
+        IF (DABS(X).GT.1.0D0) LS=-1.D0
         XQ=DSQRT(LS*(1.0D0-X*X))
         XS=LS*(1.0D0-X*X)
-        DO 30 I=1,M
+        II = 1.D0
-30         PM(I,I)=-LS*(2.0D0*I-1.0D0)*XQ*PM(I-1,I-1)
+        DO I=1,M
-        DO 35 I=0,M
+          PM(I,I)=-LS*II*XQ*PM(I-1,I-1)
-35         PM(I,I+1)=(2.0D0*I+1.0D0)*X*PM(I,I)
+          II = II+2.D0
+        ENDDO
+        II = 1.D0
+        DO I=0,M
+          PM(I,I+1)=II*X*PM(I,I)
+          II = II+2.D0
+        ENDDO
 
-        DO 40 I=0,M
+        II = 0.D0
-        DO 40 J=I+2,N
+        DO I=0,M 
-           PM(I,J)=((2.0D0*J-1.0D0)*X*PM(I,J-1)- (I+J-1.0D0)*PM(I,J-2))/(J-I)
+          JJ = II+2.D0
-40      CONTINUE
+          DO J=I+2,N
+            PM(I,J)=((2.0D0*JJ-1.0D0)*X*PM(I,J-1)- (II+JJ-1.0D0)*PM(I,J-2))*INVERSE(J-I)
+            JJ = JJ+1.D0
+          ENDDO
+          II = II+1.D0
+        ENDDO
         END
 
 
@@ -1703,17 +1750,37 @@ end
 !c
 double precision function ylm(l,m,theta,phi)
 implicit none
-integer l,m
+integer l,m,i
-double precision theta,phi,pm,factor,pi,x,fact,sign
+double precision theta,phi,pm,factor,twopi,x,fact,sign
 DIMENSION PM(0:100,0:100)
-pi=dacos(-1.d0)
+twopi=2.d0*dacos(-1.d0)
 x=dcos(theta)
-sign=(-1.d0)**m
+if (iand(m,1) == 1) then
+  sign=-1.d0
+else
+  sign=1.d0
+endif
 CALL LPMN(100,l,l,X,PM)
-factor=dsqrt( (2*l+1)*fact(l-iabs(m)) /(4.d0*pi*fact(l+iabs(m))) )
+if (m > 0) then
-if(m.gt.0)ylm=sign*dsqrt(2.d0)*factor*pm(m,l)*dcos(dfloat(m)*phi)
+  factor=dsqrt((l+l+1)*fact(l-m) /(twopi*fact(l+m)) )
-if(m.eq.0)ylm=factor*pm(m,l)
+!  factor = dble(l+m)
-if(m.lt.0)ylm=sign*dsqrt(2.d0)*factor*pm(iabs(m),l)*dsin(dfloat(iabs(m))*phi)
+!  do i=-m,m-1
+!    factor = factor * (factor - 1.d0)
+!  enddo
+!  factor=dsqrt(dble(l+l+1)/(twopi*factor) )
+  ylm=sign*factor*pm(m,l)*dcos(dfloat(m)*phi)
+else if (m == 0) then
+  factor=dsqrt( 0.5d0*(l+l+1) /twopi )
+  ylm=factor*pm(m,l)
+else if (m < 0) then
+  factor=dsqrt( (l+l+1)*fact(l+m) /(twopi*fact(l-m)) )
+!  factor = dble(l-m)
+!  do i=m,-m-1
+!    factor = factor * (factor - 1.d0)
+!  enddo
+!  factor=dsqrt(dble(l+l+1)/(twopi*factor) )
+  ylm=sign*factor*pm(-m,l)*dsin(dfloat(-m)*phi)
+endif
 end
 
 !c Explicit representation of Legendre polynomials P_n(x) 
@@ -1829,11 +1896,12 @@ end
 double precision function binom_gen(alpha,n)
  implicit none
  integer :: n,k
- double precision :: fact,alpha,prod
+ double precision :: fact,alpha,prod, factn_inv
  prod=1.d0
+ factn_inv = 1.d0/fact(n)
  do k=0,n-1
    prod=prod*(alpha-k)
-   binom_gen = prod/(fact(n))
+   binom_gen = prod*factn_inv
  enddo
 end
 
@@ -1881,6 +1949,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
   double precision :: term_A, term_B, term_rap,  expo
   double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square
   double precision :: int_prod_bessel_loc
+  double precision :: inverses(0:300)
 
   logical done
 
@@ -1927,18 +1996,32 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
     ! Initialise the first recurence terme for the q loop
     s_q_0 = s_0_0
 
+
     ! Loop over q for the convergence of the sequence
     do while (.not.done)  
 
       ! Init
-      sum=0
       s_q_k=s_q_0
+      sum=s_q_0
 
-      ! Iteration of k
+      if (q>300) then
-      do k=0,q
+        stop 'pseudopot.f90 : q > 300'
+      endif
+
+      do k=0,q-1
+        s_q_k = ( dble(2*(q-k+m)+1)*dble(q-k)*inverses(k) ) * s_q_k
         sum=sum+s_q_k
-        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)/dble(2*(k+n)+3) ) * ( dble(q-k)/dble(k+1)) * s_q_k
       enddo
+      inverses(q) = a_over_b_square/(dble(2*(q+n)+3) * dble(q+1)) 
+!      do k=0,q
+!        sum=sum+s_q_k
+!        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k
+!      enddo
+      ! Iteration of k
+!      do k=0,q
+!        sum=sum+s_q_k
+!        s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k
+!      enddo
     
       int=int+sum
 
@@ -2120,15 +2203,15 @@ parameter (ntot_max=14)
 integer l,m
 double precision a(3),g_a,c(3)
 double precision prod,binom_func,accu,bigI,ylm,bessel_mod
-double precision theta_AC0,phi_AC0,ac,factor,fourpi,arg,r,areal
+double precision theta_AC0,phi_AC0,ac,ac2,factor,fourpi,arg,r,areal
 integer ntotA,mu,k1,k2,k3,lambda
 integer n_a(3)
-double precision &
+double precision y, f1, f2
-array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)
+double precision, allocatable :: array_coefs_A(:,:)
-double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max), y
 
-ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
+ac2=(a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2
-arg=g_a*(ac**2+r**2)
+ac=dsqrt(ac2)
+arg=g_a*(ac2+r*r)
 fourpi=4.d0*dacos(-1.d0)
 factor=fourpi*dexp(-arg)
 areal=2.d0*g_a*ac
@@ -2144,51 +2227,45 @@ else
  theta_AC0=dacos( (a(3)-c(3))/ac )
  phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
 
+ allocate (array_coefs_A(0:ntotA,3))
  do k1=0,n_a(1)
-  do k2=0,n_a(2)
+   array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)*r**(k1)
-   do k3=0,n_a(3)
-  array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
-  *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) &
-  *r**(k1+k2+k3)
-   enddo
-  enddo
  enddo
-
+ do k2=0,n_a(2)
- do lambda=0,l+ntotA
+   array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)*r**(k2)
-  do mu=-lambda,lambda
+ enddo
-   do k1=0,n_a(1)
+ do k3=0,n_a(3)
-   do k2=0,n_a(2)
+   array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)*r**(k3)
-   do k3=0,n_a(3)
-    array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
-   enddo
-   enddo
-   enddo
-  enddo
  enddo
 
  accu=0.d0
  do lambda=0,l+ntotA
-  do mu=-lambda,lambda
+   do mu=-lambda,lambda
-   y = ylm(lambda,mu,theta_AC0,phi_AC0)
+     y = ylm(lambda,mu,theta_AC0,phi_AC0)
-   if (y == 0.d0) then
+     if (y == 0.d0) then
-     cycle
+       cycle
-   endif
+     endif
-   do k1=0,n_a(1)
+     do k3=0,n_a(3)
-   do k2=0,n_a(2)
+       f1 = y*array_coefs_A(k3,3)
-   do k3=0,n_a(3)
+       if (f1 == 0.d0) cycle
-    prod=y*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3)
+       do k2=0,n_a(2)
-    if (prod == 0.d0) then
+         f2 = f1*array_coefs_A(k2,2)
-      cycle
+         if (f2 == 0.d0) cycle
-    endif
+         do k1=0,n_a(1)
-    if (areal*r < 100.d0) then ! overflow!
+           prod=f2*array_coefs_A(k1,1)*bigI(lambda,mu,l,m,k1,k2,k3)
-      accu=accu+prod*bessel_mod(areal*r,lambda)
+           if (prod == 0.d0) then
-    endif
+             cycle
+           endif
+           if (areal*r < 100.d0) then ! overflow!
+             accu=accu+prod*bessel_mod(areal*r,lambda)
+           endif
+         enddo
+       enddo
+     enddo
    enddo
-   enddo
-   enddo
-  enddo
  enddo
  ylm_orb=factor*accu
+ deallocate (array_coefs_A)
  return
 endif
 

src/Pseudo/EZFIO.cfg

@@ -51,7 +51,7 @@ size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax)
 
 [do_pseudo]
 type: logical
-doc:  Using pseudo potential integral of not
+doc:  Using pseudo potential integral or not
 interface: ezfio,provider,ocaml
 default: False
 

src/Pseudo/README.rst

@@ -29,7 +29,7 @@ Documentation
 
 
 `do_pseudo <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo/ezfio_interface.irp.f#L248>`_
-  Using pseudo potential integral of not
+  Using pseudo potential integral or not
 
 
 `pseudo_dz_k <http://github.com/LCPQ/quantum_package/tree/master/src/Pseudo/ezfio_interface.irp.f#L204>`_