Improved AO basis printing

ocaml/Basis.ml

@@ -36,8 +36,30 @@ let read_element in_channel at_number element =
 ;;
 
 let to_string b =
+  let new_nucleus n = 
+    Printf.sprintf "Atom %d:" n
+  in
+  
+  let rec do_work accu current_nucleus = function
+  | [] -> List.rev accu
+  | (g,n)::tail -> 
+    let n = Nucl_number.to_int n
+    in
+    let accu = 
+       if (n <> current_nucleus) then
+         (new_nucleus n)::""::accu
+       else
+         accu
+    in
+    do_work ((Gto.to_string g)::accu) n tail
+  in
+  do_work [new_nucleus 1] 1 b
+  |> String.concat ~sep:"\n"
+;;
+(*
   List.map ~f:(fun (g,n) ->
      let n = Nucl_number.to_int n in
      (Int.to_string n)^":"^(Gto.to_string g)) b
-  |> String.concat ~sep:"\n"
 ;;
+*)
+

ocaml/Gto.ml

@@ -74,3 +74,23 @@ let to_string { sym = sym ; lc = lc } =
   Printf.sprintf "( %s, %s )" (Symmetry.to_string sym)
     (String.concat (List.map ~f:f lc) ~sep:", ")
 ;;
+
+let to_string { sym = sym ; lc = lc } =
+  let result = 
+    Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc)
+  in
+  let rec do_work accu i = function
+  | [] -> List.rev accu
+  | (p,c)::tail -> 
+    let p = AO_expo.to_float p.Primitive.expo
+    and c = AO_coef.to_float c
+    in
+    let result = 
+      Printf.sprintf "%3d %16f  %16f" i p c
+    in
+    do_work (result::accu) (i+1) tail
+  in
+  (do_work [result] 1 lc)
+  |> String.concat ~sep:"\n"
+;;
+

ocaml/Input_ao_basis.ml

@@ -134,27 +134,25 @@ end = struct
         (Array.to_list gto_array)
         (Array.to_list b.ao_nucl)
     in
-    Printf.sprintf "
-# AO Basis
-# ========
-%s" (Long_basis.to_string long_basis)
+    let short_basis =
+      Long_basis.to_basis long_basis
+    in
+    Printf.sprintf "Basis name : %s\n\n%s" b.ao_basis
+    (Basis.to_string short_basis)
   ;;
     
       
   let debug b =
     Printf.sprintf "
-# AO Basis
-# ========
-#
-# ao_basis        = %s
-# ao_num          = %s
-# ao_prim_num     = %s
-# ao_prim_num_max = %s
-# ao_nucl         = %s
-# ao_power        = %s
-# ao_coef         = %s
-# ao_expo         = %s
-#"
+ao_basis        = %s
+ao_num          = %s
+ao_prim_num     = %s
+ao_prim_num_max = %s
+ao_nucl         = %s
+ao_power        = %s
+ao_coef         = %s
+ao_expo         = %s
+"
     b.ao_basis
     (AO_number.to_string b.ao_num)
     (b.ao_prim_num |> Array.to_list |> List.map

ocaml/Input_bi_integrals.ml

@@ -15,6 +15,7 @@ module Bielec_integrals : sig
   ;;
   val read : unit -> t
   val to_string : t -> string
+  val of_string : string -> t
 end = struct
   type t = 
     { read_ao_integrals  : bool;
@@ -114,8 +115,7 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "
-read_ao_integrals  = %s
+    Printf.sprintf "read_ao_integrals  = %s
 read_mo_integrals  = %s
 write_ao_integrals = %s
 write_mo_integrals = %s
@@ -130,6 +130,13 @@ direct             = %s
         (Threshold.to_string b.threshold_ao)
         (Threshold.to_string b.threshold_mo)
         (Bool.to_string b.direct)
+  ;;
+
+  let of_string s =
+    input_to_sexp s
+    |> t_of_sexp
+  ;;
+
 end
 
 

ocaml/Long_basis.ml

@@ -19,11 +19,26 @@ let of_basis b =
   |> List.rev
 ;;
 
+let to_basis b =
+  let rec do_work accu = function
+  | [] -> List.rev accu
+  | (s,g,n)::tail -> 
+    let first_sym = 
+      Symmetry.Xyz.of_symmetry g.Gto.sym
+      |> List.hd_exn
+    in
+    let new_accu = 
+      if ( s = first_sym ) then
+        (g,n)::accu
+      else
+        accu
+    in
+    do_work new_accu tail
+  in
+  do_work [] b 
+;;
 
 let to_string b =
-  Sexp.to_string (sexp_of_t b)
-;;
-(*
   let middle = List.map ~f:(fun (x,y,z) -> 
      "( "^((Int.to_string (Nucl_number.to_int z)))^", "^
      (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y)
@@ -32,4 +47,3 @@ let to_string b =
   |> String.concat ~sep:",\n"
   in "("^middle^")"
 ;;
-*)

ocaml/Long_basis.mli

@@ -5,12 +5,16 @@ open Qptypes;;
   * all the D orbitals are converted to xx, xy, xz, yy, yx
   * etc
 *)
-type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list
+type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list with sexp
 
 (** Transform a basis to a long basis *)
 val of_basis :
   (Gto.t * Nucl_number.t) list -> (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list
 
+(** Transform a long basis to a basis *)
+val to_basis :
+  (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list 
+
 (** Convert the basis into its string representation *)
 val to_string :
   (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> string

ocaml/Qputils.ml

@@ -1,3 +1,6 @@
+open Core.Std;;
+
+(*
 let rec transpose = function
 | []          -> []
 | []::tail    -> transpose tail
@@ -7,5 +10,20 @@ let rec transpose = function
     in 
     new_head @ new_tail
 ;;
+*)
 
 
+let input_to_sexp s =
+    let result = 
+      String.split_lines s 
+      |> List.filter ~f:(fun x->
+         (String.strip x) <> "")
+      |> List.map ~f:(fun x->
+         "("^(String.tr '=' ' ' x)^")")
+      |> String.concat 
+    in
+    print_endline ("("^result^")");
+    "("^result^")"
+    |> Sexp.of_string
+;;
+

ocaml/qp_edit.ml

@@ -0,0 +1,95 @@
+open Qputils;;
+open Qptypes;;
+open Core.Std;;
+
+let instructions filename = Printf.sprintf
+"# ===============
+# Quantum Package
+# ===============
+# 
+# File : %s
+#
+# Lines starting with a '#' sign are commented.
+#" filename
+  
+type keyword = 
+| Ao_basis
+| Bielec_integrals
+;;
+
+let keyword_to_string = function
+| Ao_basis -> "AO basis"
+| Bielec_integrals -> "Two electron integrals"
+;;
+
+let make_header kw =
+  let s = keyword_to_string kw in
+  let l = String.length s in
+  "\n\n# "^s^"\n"^"# "^(String.init l ~f:(fun _ -> '='))^"\n\n"
+;;
+
+let get_bielec () = 
+  (make_header Bielec_integrals)^
+  (Input.Bielec_integrals.(to_string (read ())))
+;;
+
+let get_ao_basis () = 
+  (make_header Ao_basis)^
+  (Input.Ao_basis.(to_string (read ())))
+;;
+
+let run ezfio_filename =
+
+  (* Open EZFIO *)
+  if (not (Sys.file_exists_exn ezfio_filename)) then
+    failwith (ezfio_filename^" does not exists");
+
+  Ezfio.set_file ezfio_filename;
+
+  let output = [
+    (instructions ezfio_filename) ;
+    (get_ao_basis ()) ;
+    (get_bielec ()) 
+  ] in
+  String.concat output
+  |> print_string
+;;
+
+
+
+let spec =
+  let open Command.Spec in
+  empty 
+(*
+  +> flag "i"  (optional string)
+     ~doc:"Prints input data"
+  +> flag "o" (optional string)
+     ~doc:"Prints output data"
+*)
+  +> anon ("ezfio_file" %: string)
+;;
+
+let command = 
+    Command.basic 
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+      "
+Edit input data
+      ")
+    spec
+(*    (fun i o ezfio_file () -> *)
+    (*fun ezfio_file () -> 
+       try 
+           run ezfio_file
+       with
+       | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n")
+    *)
+    (fun ezfio_file () -> run ezfio_file)
+;;
+
+let () =
+    Command.run command
+;;
+
+
+

ocaml/test_basis.ml

@@ -21,9 +21,20 @@ let test_module () =
       (Basis.read_element basis_channel (Nucl_number.of_int 2) Element.F) 
     in
 
-    Long_basis.of_basis basis 
-    |> Long_basis.to_string 
-    |> print_endline
+    print_string "Long basis\n==========\n";
+    let long_basis = 
+      Long_basis.of_basis basis 
+    in
+    print_endline (Long_basis.to_string long_basis);
+
+    let short_basis = 
+      Long_basis.to_basis long_basis
+    in
+    if (short_basis <> basis) then
+      print_endline "(short_basis <> basis)"
+    ;
+    print_string "Short basis\n===========\n";
+    print_endline (Basis.to_string basis);
 ;;
 
 test_module ();

ocaml/test_gto.ml

@@ -13,11 +13,20 @@ let test_gto_1 () =
   let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
   ignore (input_line in_channel);
   let gto = Gto.read_one in_channel in
-  print_string (Gto.to_string gto);
-  let gto = Gto.read_one in_channel in
-  print_string (Gto.to_string gto);
-  let gto = Gto.read_one in_channel in
-  print_string (Gto.to_string gto);
+  print_endline (Gto.to_string gto);
+  In_channel.seek in_channel 0L;
+  ignore (input_line in_channel);
+  let gto2 = Gto.read_one in_channel in
+  print_endline (Gto.to_string gto2);
+  let gto3 = Gto.read_one in_channel in
+  print_endline (Gto.to_string gto3);
+  if (gto2 = gto) then
+    print_endline "gto2 = gto";
+  if (gto3 = gto) then
+    print_endline "gto3 = gto";
+  if (gto3 = gto3) then
+    print_endline "gto3 = gto3";
+  
 ;;
 
 let test_gto_2 () =
@@ -34,7 +43,7 @@ let test_gto () =
 ;;
 
 let test_module () =
-  test_gto()
+  test_gto_1()
 ;;
 
 test_module ();;

ocaml/test_input.ml

@@ -10,7 +10,15 @@ let test_bielec_intergals () =
   Ezfio.set_file "F2.ezfio" ;
   let b = Input.Bielec_integrals.read ()
   in
-  print_endline (Input.Bielec_integrals.to_string b);
+  let output = Input.Bielec_integrals.to_string b
+  in
+  print_endline output;
+  let b2 = Input.Bielec_integrals.of_string output in
+  if (b = b2) then
+   print_endline "OK"
+  else
+   print_endline "b <> b2"
+  ;
 ;;
 
 let test_bitmasks () =
@@ -77,15 +85,15 @@ let test_nucl () =
 ;;
 
 (*
-test_hf ();;
 test_ao ();;
 test_bielec_intergals ();;
 test_bitmasks ();
 test_cis ();
-test_dets ();
 test_cisd_sc2 ();
+test_dets ();
+test_hf ();;
 test_mo ();;
 test_nucl ();
 *)
-test_ao ();;
+test_bielec_intergals ();;