Merge pull request #132 from scemama/master

Corrected bug when generating pseudopotentials

configure

@@ -533,7 +533,7 @@ def recommendation():
     print "   source {0}".format(path)
     print ""
     print "Then, install the modules you want to install using :"
-    print "    qp_install_module.py "
+    print "    qp_module.py"
     print ""
     print "Finally :"
     print "   ninja"

ocaml/Qputils.ml

@@ -1,4 +1,4 @@
-open Core.Std;;
+open Core.Std
 
 (*
 let rec transpose = function
@@ -24,5 +24,21 @@ let input_to_sexp s =
     print_endline ("("^result^")");
     "("^result^")"
     |> Sexp.of_string
-;;
+
+let rmdir dirname =
+  let rec remove_one dir =
+    Sys.chdir dir;
+    Sys.readdir "."
+    |> Array.iter ~f:(fun x ->
+      match (Sys.is_directory x, Sys.is_file x) with
+      | (`Yes, _) -> remove_one x
+      | (_, `Yes) -> Sys.remove x
+      | _ -> failwith ("Unable to remove file "^x^".")
+    );
+    Sys.chdir "..";
+    Unix.rmdir dir
+  in
+  remove_one dirname
+  
+    
 

ocaml/qp_create_ezfio_from_xyz.ml

@@ -1,6 +1,6 @@
-open Qputils;;
+open Qputils
-open Qptypes;;
+open Qptypes
-open Core.Std;;
+open Core.Std
 
 let spec =
   let open Command.Spec in
@@ -316,289 +316,316 @@ let run ?o b c d m p cart xyz_file =
   if Sys.file_exists_exn ezfio_file then
     failwith (ezfio_file^" already exists");
 
-  (* Create EZFIO *)
+  let write_file () = 
-  Ezfio.set_file ezfio_file;
+      (* Create EZFIO *)
+      Ezfio.set_file ezfio_file;
 
-  (* Write Pseudo *)
+      (* Write Pseudo *)
-  let pseudo =
+      let pseudo =
-    List.map nuclei ~f:(fun x ->
+        List.map nuclei ~f:(fun x ->
-         match pseudo_channel x.Atom.element with
+            match pseudo_channel x.Atom.element with
-         | Some channel -> Pseudo.read_element channel x.Atom.element
+            | Some channel -> Pseudo.read_element channel x.Atom.element
-         | None -> Pseudo.empty x.Atom.element
+            | None -> Pseudo.empty x.Atom.element
-       ) 
+          ) 
-  in
+      in
 
-  let molecule =     
+      let molecule =     
-    let n_elec_to_remove =
+        let n_elec_to_remove =
-      List.fold pseudo ~init:0 ~f:(fun accu x ->
+          List.fold pseudo ~init:0 ~f:(fun accu x ->
-        accu + (Positive_int.to_int x.Pseudo.n_elec))
+            accu + (Positive_int.to_int x.Pseudo.n_elec))
-    in
-    { Molecule.elec_alpha = 
-        (Elec_alpha_number.to_int molecule.Molecule.elec_alpha)
-         - n_elec_to_remove/2
-        |> Elec_alpha_number.of_int;
-      Molecule.elec_beta = 
-        (Elec_beta_number.to_int molecule.Molecule.elec_beta)
-         - (n_elec_to_remove - n_elec_to_remove/2)
-        |> Elec_beta_number.of_int;
-      Molecule.nuclei = 
-        let charges =
-          List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec
-            |> Float.of_int)
-          |> Array.of_list
         in
-        List.mapi molecule.Molecule.nuclei ~f:(fun i x ->
+        { Molecule.elec_alpha = 
-          { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i)
+            (Elec_alpha_number.to_int molecule.Molecule.elec_alpha)
-            |> Charge.of_float }
+            - n_elec_to_remove/2
-        )
+            |> Elec_alpha_number.of_int;
-    }
+          Molecule.elec_beta = 
-  in
+            (Elec_beta_number.to_int molecule.Molecule.elec_beta)
-  let nuclei =
+            - (n_elec_to_remove - n_elec_to_remove/2)
-    molecule.Molecule.nuclei @ dummy
+            |> Elec_beta_number.of_int;
-  in
+          Molecule.nuclei = 
-    
+            let charges =
-
+              List.map pseudo ~f:(fun x -> Positive_int.to_int x.Pseudo.n_elec
-  (* Write Electrons *)
+                |> Float.of_int)
-  Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int
+              |> Array.of_list
-    molecule.Molecule.elec_alpha ) ;
+            in
-  Ezfio.set_electrons_elec_beta_num  ( Elec_beta_number.to_int
+            List.mapi molecule.Molecule.nuclei ~f:(fun i x ->
-    molecule.Molecule.elec_beta  ) ;
+              { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i)
-
+                |> Charge.of_float }
-  (* Write Nuclei *)
+            )
-  let labels =
+        }
-     List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei
-  and charges = 
-     List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei 
-  and coords  = 
-    (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @
-    (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @
-    (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in
-  let nucl_num = (List.length labels) in
-  Ezfio.set_nuclei_nucl_num nucl_num ;
-  Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list 
-    ~rank:1 ~dim:[| nucl_num |] ~data:labels);
-  Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list 
-    ~rank:1 ~dim:[| nucl_num |] ~data:charges);
-  Ezfio.set_nuclei_nucl_coord  (Ezfio.ezfio_array_of_list 
-    ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords);
-
-
-  (* Write pseudopotential *)
-  let () = 
-    match p with
-    | None -> Ezfio.set_pseudo_do_pseudo false
-    | _    -> Ezfio.set_pseudo_do_pseudo true
-  in
-
-  let klocmax = 
-    List.fold pseudo ~init:0 ~f:(fun accu x -> 
-      let x = 
-        List.length x.Pseudo.local
       in
-      if (x > accu) then x
+      let nuclei =
-      else accu
+        molecule.Molecule.nuclei @ dummy
-    )
-  and kmax = 
-    List.fold pseudo ~init:0 ~f:(fun accu x -> 
-      let x = 
-        List.length x.Pseudo.non_local
       in
-      if (x > accu) then x
+        
-      else accu
+
-    )
+      (* Write Electrons *)
-  and lmax = 
+      Ezfio.set_electrons_elec_alpha_num ( Elec_alpha_number.to_int
-    List.fold pseudo ~init:0 ~f:(fun accu x -> 
+        molecule.Molecule.elec_alpha ) ;
-      let x = 
+      Ezfio.set_electrons_elec_beta_num  ( Elec_beta_number.to_int
-        List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) ->
+        molecule.Molecule.elec_beta  ) ;
+
+      (* Write Nuclei *)
+      let labels =
+        List.map ~f:(fun x->Element.to_string x.Atom.element) nuclei
+      and charges = 
+        List.map ~f:(fun x-> Atom.(Charge.to_float x.charge)) nuclei 
+      and coords  = 
+        (List.map ~f:(fun x-> x.Atom.coord.Point3d.x) nuclei) @
+        (List.map ~f:(fun x-> x.Atom.coord.Point3d.y) nuclei) @
+        (List.map ~f:(fun x-> x.Atom.coord.Point3d.z) nuclei) in
+      let nucl_num = (List.length labels) in
+      Ezfio.set_nuclei_nucl_num nucl_num ;
+      Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list 
+        ~rank:1 ~dim:[| nucl_num |] ~data:labels);
+      Ezfio.set_nuclei_nucl_charge (Ezfio.ezfio_array_of_list 
+        ~rank:1 ~dim:[| nucl_num |] ~data:charges);
+      Ezfio.set_nuclei_nucl_coord  (Ezfio.ezfio_array_of_list 
+        ~rank:2 ~dim:[| nucl_num ; 3 |] ~data:coords);
+
+
+      (* Write pseudopotential *)
+      let () = 
+        match p with
+        | None -> Ezfio.set_pseudo_do_pseudo false
+        | _    -> Ezfio.set_pseudo_do_pseudo true
+      in
+
+      let klocmax = 
+        List.fold pseudo ~init:0 ~f:(fun accu x -> 
           let x = 
-             Positive_int.to_int x.Pseudo.Primitive_non_local.proj
+            List.length x.Pseudo.local
+          in
+          if (x > accu) then x
+          else accu
+        )
+      and lmax = 
+        List.fold pseudo ~init:0 ~f:(fun accu x -> 
+          let x = 
+            List.fold x.Pseudo.non_local ~init:0 ~f:(fun accu (x,_) ->
+              let x = 
+                Positive_int.to_int x.Pseudo.Primitive_non_local.proj
+              in
+              if (x > accu) then x
+              else accu
+            )
           in
           if (x > accu) then x
           else accu
         )
       in
-      if (x > accu) then x
-      else accu
-    )
-  in
 
-  let () = 
+     let kmax = 
-    Ezfio.set_pseudo_pseudo_klocmax klocmax;
+         Array.init (lmax+1) ~f:(fun i->
-    Ezfio.set_pseudo_pseudo_kmax kmax;
+           List.map pseudo ~f:(fun x -> 
-    Ezfio.set_pseudo_pseudo_lmax lmax;
+             List.filter x.Pseudo.non_local ~f:(fun (y,_) -> 
-    let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k =
+               (Positive_int.to_int y.Pseudo.Primitive_non_local.proj) = i)
-      Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
+             |> List.length )
-      Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
+           |> List.fold ~init:0 ~f:(fun accu x ->
-      Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0
+             if accu > x then accu else x)
-    in
+           )
-    List.iteri pseudo ~f:(fun j x -> 
+         |> Array.fold ~init:0 ~f:(fun accu i ->
-      List.iteri x.Pseudo.local ~f:(fun i (y,c) ->
+            if i > accu then i else accu)
-        tmp_array_v_k.(i).(j)  <- AO_coef.to_float c;
+      in
-        let y, z =
+ 
-          AO_expo.to_float y.Pseudo.Primitive_local.expo, 
+              
-          R_power.to_int y.Pseudo.Primitive_local.r_power 
+      let () = 
+        Ezfio.set_pseudo_pseudo_klocmax klocmax;
+        Ezfio.set_pseudo_pseudo_kmax kmax;
+        Ezfio.set_pseudo_pseudo_lmax lmax;
+        let tmp_array_v_k, tmp_array_dz_k, tmp_array_n_k =
+          Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
+          Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0. ,
+          Array.make_matrix ~dimx:klocmax ~dimy:nucl_num 0
         in
-        tmp_array_dz_k.(i).(j) <- y;
+        List.iteri pseudo ~f:(fun j x -> 
-        tmp_array_n_k.(i).(j)  <- z;
+          List.iteri x.Pseudo.local ~f:(fun i (y,c) ->
-      )
+            tmp_array_v_k.(i).(j)  <- AO_coef.to_float c;
-    );
+            let y, z =
-    let concat_2d tmp_array = 
+              AO_expo.to_float y.Pseudo.Primitive_local.expo, 
-      let data = 
+              R_power.to_int y.Pseudo.Primitive_local.r_power 
-        Array.map tmp_array ~f:Array.to_list
-        |> Array.to_list
-        |> List.concat
-      in
-      Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data
-    in  
-    concat_2d tmp_array_v_k 
-    |> Ezfio.set_pseudo_pseudo_v_k ;
-    concat_2d tmp_array_dz_k 
-    |> Ezfio.set_pseudo_pseudo_dz_k;
-    concat_2d tmp_array_n_k 
-    |> Ezfio.set_pseudo_pseudo_n_k;
-
-    let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = 
-      Array.create ~len:(lmax+1) 
-       (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ),
-      Array.create ~len:(lmax+1) 
-       (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. ),
-      Array.create ~len:(lmax+1) 
-       (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 )
-    in
-    List.iteri pseudo ~f:(fun j x -> 
-      List.iteri x.Pseudo.non_local ~f:(fun i (y,c) ->
-        let k, y, z =
-          Positive_int.to_int y.Pseudo.Primitive_non_local.proj,
-          AO_expo.to_float y.Pseudo.Primitive_non_local.expo,
-          R_power.to_int y.Pseudo.Primitive_non_local.r_power
-        in
-        tmp_array_v_kl.(k).(i).(j)  <- AO_coef.to_float c;
-        tmp_array_dz_kl.(k).(i).(j) <- y;
-        tmp_array_n_kl.(k).(i).(j)  <- z;
-      ) 
-    );    
-    let concat_3d tmp_array = 
-      let data = 
-        Array.map tmp_array ~f:(fun x -> 
-          Array.map x ~f:Array.to_list
-          |> Array.to_list
-          |> List.concat)
-        |> Array.to_list
-        |> List.concat
-      in
-      Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data
-    in
-    concat_3d tmp_array_v_kl
-    |> Ezfio.set_pseudo_pseudo_v_kl ;
-    concat_3d tmp_array_dz_kl
-    |> Ezfio.set_pseudo_pseudo_dz_kl ;
-    concat_3d tmp_array_n_kl
-    |> Ezfio.set_pseudo_pseudo_n_kl ;
-  in
-    
-
-
-
-  (* Write Basis set *)
-  let basis =
-
-    let nmax =
-       Nucl_number.get_max ()
-    in
-    let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
-    | [] -> accu
-    | e::tail ->
-      let new_accu =
-        (e,(Nucl_number.of_int ~max:nmax n))::accu
-      in
-      do_work new_accu (n+1) tail
-    in
-    let result = do_work [] 1  nuclei
-    |> List.rev
-    |> List.map ~f:(fun (x,i) ->
-       try 
-         let e =
-           match x.Atom.element with
-           | Element.X -> Element.H
-           | e -> e
-         in
-         Basis.read_element (basis_channel x.Atom.element) i e
-       with
-       | End_of_file -> failwith
-             ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
-       ) 
-    |> List.concat
-    in
-    (* close all in_channels *)
-    result
-  in
-  let long_basis = Long_basis.of_basis basis in
-  let ao_num = List.length long_basis in
-  Ezfio.set_ao_basis_ao_num ao_num;
-  Ezfio.set_ao_basis_ao_basis b;
-  let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) 
-  and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n)
-  and ao_power= 
-    let l = List.map long_basis ~f:(fun (x,_,_) -> x) in
-    (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@
-    (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@
-    (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) 
-  in
-  let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x ->
-    if x > s then x
-    else s) ao_prim_num
-  in
-  let gtos = 
-     List.map long_basis ~f:(fun (_,x,_) -> x)
-  in
-
-  let create_expo_coef ec = 
-      let coefs = 
-        begin match ec with
-        | `Coefs -> List.map gtos ~f:(fun x->
-          List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) )
-        | `Expos -> List.map gtos ~f:(fun x->
-          List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float
-          prim.Primitive.expo) )
-        end
-      in
-      let rec get_n n accu = function
-        | [] -> List.rev accu  
-        | h::tail ->
-            let y =
-            begin match List.nth h n with
-            | Some x -> x
-            | None -> 0.
-            end
             in
-            get_n n (y::accu) tail
+            tmp_array_dz_k.(i).(j) <- y;
-      in
+            tmp_array_n_k.(i).(j)  <- z;
-      let rec build accu = function
+          )
-        | n when n=ao_prim_num_max -> accu
+        );
-        | n -> build ( accu @ (get_n n [] coefs) ) (n+1)
+        let concat_2d tmp_array = 
-      in
+          let data = 
-      build [] 0
+            Array.map tmp_array ~f:Array.to_list
-  in
+            |> Array.to_list
+            |> List.concat
+          in
+          Ezfio.ezfio_array_of_list ~rank:2 ~dim:[|nucl_num ; klocmax|] ~data
+        in  
+        concat_2d tmp_array_v_k 
+        |> Ezfio.set_pseudo_pseudo_v_k ;
+        concat_2d tmp_array_dz_k 
+        |> Ezfio.set_pseudo_pseudo_dz_k;
+        concat_2d tmp_array_n_k 
+        |> Ezfio.set_pseudo_pseudo_n_k;
 
-  let ao_coef = create_expo_coef `Coefs
+        let tmp_array_v_kl, tmp_array_dz_kl, tmp_array_n_kl = 
-  and ao_expo = create_expo_coef `Expos
+          Array.init (lmax+1) ~f:(fun _ ->
+          (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )),
+          Array.init (lmax+1) ~f:(fun _ ->
+          (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0. )),
+          Array.init (lmax+1) ~f:(fun _ ->
+          (Array.make_matrix ~dimx:kmax ~dimy:nucl_num 0 ))
+        in
+        List.iteri pseudo ~f:(fun j x -> 
+          let last_idx =
+            Array.create ~len:(lmax+1) 0
+          in
+          List.iter x.Pseudo.non_local ~f:(fun (y,c) -> 
+            let k, y, z =
+              Positive_int.to_int y.Pseudo.Primitive_non_local.proj,
+              AO_expo.to_float y.Pseudo.Primitive_non_local.expo,
+              R_power.to_int y.Pseudo.Primitive_non_local.r_power
+            in
+            let i = 
+              last_idx.(k) 
+            in
+            tmp_array_v_kl.(k).(i).(j)  <- AO_coef.to_float c;
+            tmp_array_dz_kl.(k).(i).(j) <- y;
+            tmp_array_n_kl.(k).(i).(j)  <- z;
+            last_idx.(k) <- i+1;
+           )
+        );    
+        let concat_3d tmp_array = 
+          let data = 
+            Array.map tmp_array ~f:(fun x -> 
+              Array.map x ~f:Array.to_list
+              |> Array.to_list
+              |> List.concat)
+            |> Array.to_list
+            |> List.concat
+          in
+          Ezfio.ezfio_array_of_list ~rank:3 ~dim:[|nucl_num ; kmax ; lmax+1|] ~data
+        in
+        concat_3d tmp_array_v_kl
+        |> Ezfio.set_pseudo_pseudo_v_kl ;
+        concat_3d tmp_array_dz_kl
+        |> Ezfio.set_pseudo_pseudo_dz_kl ;
+        concat_3d tmp_array_n_kl
+        |> Ezfio.set_pseudo_pseudo_n_kl ;
+      in
+        
+
+
+
+      (* Write Basis set *)
+      let basis =
+
+        let nmax =
+          Nucl_number.get_max ()
+        in
+        let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function
+        | [] -> accu
+        | e::tail ->
+          let new_accu =
+            (e,(Nucl_number.of_int ~max:nmax n))::accu
+          in
+          do_work new_accu (n+1) tail
+        in
+        let result = do_work [] 1  nuclei
+        |> List.rev
+        |> List.map ~f:(fun (x,i) ->
+          try 
+            let e =
+              match x.Atom.element with
+              | Element.X -> Element.H
+              | e -> e
+            in
+            Basis.read_element (basis_channel x.Atom.element) i e
+          with
+          | End_of_file -> failwith
+                ("Element "^(Element.to_string x.Atom.element)^" not found in basis set.")
+          ) 
+        |> List.concat
+        in
+        (* close all in_channels *)
+        result
+      in
+      let long_basis = Long_basis.of_basis basis in
+      let ao_num = List.length long_basis in
+      Ezfio.set_ao_basis_ao_num ao_num;
+      Ezfio.set_ao_basis_ao_basis b;
+      let ao_prim_num = List.map long_basis ~f:(fun (_,g,_) -> List.length g.Gto.lc) 
+      and ao_nucl = List.map long_basis ~f:(fun (_,_,n) -> Nucl_number.to_int n)
+      and ao_power= 
+        let l = List.map long_basis ~f:(fun (x,_,_) -> x) in
+        (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) )@
+        (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) )@
+        (List.map l ~f:(fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) ) 
+      in
+      let ao_prim_num_max = List.fold ~init:0 ~f:(fun s x ->
+        if x > s then x
+        else s) ao_prim_num
+      in
+      let gtos = 
+        List.map long_basis ~f:(fun (_,x,_) -> x)
+      in
+
+      let create_expo_coef ec = 
+          let coefs = 
+            begin match ec with
+            | `Coefs -> List.map gtos ~f:(fun x->
+              List.map x.Gto.lc ~f:(fun (_,coef) -> AO_coef.to_float coef) )
+            | `Expos -> List.map gtos ~f:(fun x->
+              List.map x.Gto.lc ~f:(fun (prim,_) -> AO_expo.to_float
+              prim.Primitive.expo) )
+            end
+          in
+          let rec get_n n accu = function
+            | [] -> List.rev accu  
+            | h::tail ->
+                let y =
+                begin match List.nth h n with
+                | Some x -> x
+                | None -> 0.
+                end
+                in 
+                get_n n (y::accu) tail
+          in
+          let rec build accu = function
+            | n when n=ao_prim_num_max -> accu
+            | n -> build ( accu @ (get_n n [] coefs) ) (n+1)
+          in
+          build [] 0
+      in
+
+      let ao_coef = create_expo_coef `Coefs
+      and ao_expo = create_expo_coef `Expos
+      in
+      let () =
+        Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
+        Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
+        Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list
+        ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
+        Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list
+        ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ;
+        Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list
+        ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ;
+        Ezfio.set_ao_basis_ao_cartesian(cart);
+      in
+      match Input.Ao_basis.read () with
+      | None -> failwith "Error in basis"
+      | Some x -> Input.Ao_basis.write x
   in
   let () =
-    Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
+    try write_file () with
-      ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
+    | ex -> 
-    Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list
+      begin
-      ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
+        begin
-    Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list
+          match Sys.is_directory ezfio_file with
-    ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
+          | `Yes ->  rmdir ezfio_file
-    Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list
+          | _ -> ()
-    ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ;
+        end; 
-    Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list
+        raise ex;
-    ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ;
+      end
-    Ezfio.set_ao_basis_ao_cartesian(cart);
+  in ()
-  in
-  match Input.Ao_basis.read () with
-  | None -> failwith "Error in basis"
-  | Some x -> Input.Ao_basis.write x
  
 
 

plugins/MRCC_Utils/H_apply.irp.f

@@ -24,5 +24,27 @@ s.data["size_max"] = "3072"
 print s
 
 
+s = H_apply("mrcepa")
+s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_ref, Ndet_non_ref"
+s.data["declarations"] += """
+    integer, intent(in) :: Ndet_ref,Ndet_non_ref
+    double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
+    double precision, intent(in) :: delta_ii_(Ndet_ref,*)
+"""
+s.data["keys_work"] = "call mrcepa_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
+s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
+s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
+s.data["decls_main"] += """
+    integer, intent(in) :: Ndet_ref,Ndet_non_ref
+    double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
+    double precision, intent(in) :: delta_ii_(Ndet_ref,*)
+"""
+s.data["finalization"] = ""
+s.data["copy_buffer"] = ""
+s.data["generate_psi_guess"] = ""
+s.data["size_max"] = "3072"
+# print s
+
+
 END_SHELL
 

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -26,7 +26,11 @@
 ! TODO --- Test perturbatif  ------
    do k=1,N_states
      lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
+     ! TODO : i_h_psi peut sortir de la boucle?
      call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
+     if (ihpsi_current(k) == 0.d0) then
+       ihpsi_current(k) = 1.d-32
+     endif
      tmp = psi_non_ref_coef(i,k)/ihpsi_current(k)
      i_pert = 0
      ! Perturbation only if 1st order < 0.5 x second order

plugins/MRCC_Utils/mrcepa_dress.irp.f

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

plugins/MRCC_Utils/mrcepa_general.irp.f

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

scripts/module/qp_module.py

@@ -175,7 +175,11 @@ def main(arguments):
         d_child = d_local.copy()
         d_child.update(d_plugin)
 
-        l_name = arguments["<name>"]
+        normalize_case = {}
+        for name in d_local.keys() + d_plugin.keys():
+            normalize_case [ name.lower() ] = name
+
+        l_name = [ normalize_case[name.lower()] for name in arguments["<name>"] ]
 
         for name in l_name:
             if name in d_local:
@@ -206,6 +210,7 @@ def main(arguments):
                 raise
 
             print "[ OK ]"
+            print ""
             print "You can now compile as usual"
             print "`cd {0} ; ninja` for exemple".format(QP_ROOT)
             print " or --in developement mode-- you can cd in a directory and compile here"