Working on qp_edit

ocaml/Basis.ml

@@ -37,7 +37,7 @@ let read_element in_channel at_number element =
 
 let to_string b =
   let new_nucleus n = 
-    Printf.sprintf "Atom %d:" n
+    Printf.sprintf "Atom %d" n
   in
   
   let rec do_work accu current_nucleus = function

ocaml/Determinant.ml

@@ -0,0 +1,54 @@
+open Core.Std;;
+open Qptypes;;
+
+type t = int64 array with sexp
+
+let to_int64_array (x:t) = (x:int64 array)
+;;
+
+let of_int64_array n_int x =
+   assert ((Array.length x) = (N_int_number.to_int n_int)*2) ;
+   x
+;;
+
+let to_alpha_beta x = 
+  let x = to_int64_array x in
+  let n_int = (Array.length x)/2 in
+  ( Array.init n_int ~f:(fun i -> x.(i)) ,
+    Array.init n_int ~f:(fun i -> x.(i+n_int)) )
+;;
+
+let to_bitlist_couple x =
+  let (xa,xb) = to_alpha_beta x in
+  let xa = to_int64_array xa
+  |> Array.to_list
+  |> Bitlist.of_int64_list
+  and xb = to_int64_array xb
+  |> Array.to_list
+  |> Bitlist.of_int64_list
+  in (xa,xb)
+;;
+
+let of_bitlist_couple (xa,xb) =
+  let ba = Bitlist.to_int64_list xa in
+  let bb = Bitlist.to_int64_list xb in
+  let n_int = Bitlist.n_int_of_mo_tot_num (List.length ba) in
+  of_int64_array n_int (Array.of_list (ba@bb))
+;;
+
+let bitlist_to_string ~mo_tot_num x =
+  List.map x ~f:(fun i -> match i with
+    | Bit.Zero -> "-"
+    | Bit.One  -> "+" )
+  |> String.concat
+  |> String.sub ~pos:0 ~len:(MO_number.to_int mo_tot_num)
+;;
+
+let to_string ~mo_tot_num x = 
+  let (xa,xb) = to_bitlist_couple x in
+  [ bitlist_to_string ~mo_tot_num:mo_tot_num xa ;
+    bitlist_to_string ~mo_tot_num:mo_tot_num xb ]
+  |> String.concat ~sep:"\n"
+;;
+
+

ocaml/Determinant.mli

@@ -0,0 +1,26 @@
+(** Determinants are stored as follows :
+  *    <-------- N_int ---------->
+  * [|  i1_alpha ; i2_alpha ; ... ;
+  *     i1_beta  ; i2_beta  ; ... ; |]
+  * where each int64 is a list of 64 MOs. When the bit is set
+  * to 1, the MO is occupied.
+  *)
+type t = int64 array with sexp
+
+(** Transform to an int64 array *)
+val to_int64_array : t -> int64 array
+
+(** Create from an int64 array *)
+val of_int64_array : Qptypes.N_int_number.t -> int64 array -> t
+
+(** Split into an alpha-only and a beta-only determinant *)
+val to_alpha_beta : t -> (int64 array)*(int64 array)
+
+(** Transform to a bit list *)
+val to_bitlist_couple : t -> Bitlist.t * Bitlist.t
+
+(** Create from a bit list *)
+val of_bitlist_couple : Bitlist.t * Bitlist.t -> t
+
+(** String representation *)
+val to_string : mo_tot_num:Qptypes.MO_number.t -> t -> string

ocaml/Input_ao_basis.ml

@@ -142,8 +142,21 @@ end = struct
       
   let to_string b =
     let short_basis = to_basis b in
-    Printf.sprintf "Basis name : %s\n\n%s\n" b.ao_basis
-    (Basis.to_string short_basis)
+    Printf.sprintf "
+Name of the AO basis ::
+
+  ao_basis = %s
+
+Basis set ::
+  
+%s
+
+" b.ao_basis
+    (Basis.to_string short_basis 
+     |> String.split ~on:'\n'
+     |> List.map ~f:(fun x-> "  "^x)
+     |> String.concat ~sep:"\n"
+    )
   ;;
 
   let to_md5 b =

ocaml/Input_bi_integrals.ml

@@ -115,13 +115,26 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "read_ao_integrals  = %s
-read_mo_integrals  = %s
-write_ao_integrals = %s
-write_mo_integrals = %s
-threshold_ao       = %s
-threshold_mo       = %s
-direct             = %s
+    Printf.sprintf "
+Read AO/MO integrals from disk ::
+
+  read_ao_integrals = %s
+  read_mo_integrals = %s
+
+Write AO/MO integrals to disk ::
+
+  write_ao_integrals = %s
+  write_mo_integrals = %s
+
+Thresholds on integrals ::
+
+  threshold_ao = %s
+  threshold_mo = %s
+
+Direct calculation of integrals ::
+
+  direct = %s
+
 "
         (Bool.to_string b.read_ao_integrals)
         (Bool.to_string b.read_mo_integrals)

ocaml/Input_cis.ml

@@ -95,12 +95,27 @@ end = struct
 
   let to_string b =
     Printf.sprintf "
-n_state_cis        = %s
-n_core_cis         = %s
-n_act_cis          = %s
-mp2_dressing       = %s
-standard_doubles   = %s
-en_2_2             = %s
+Number of states ::
+
+  n_state_cis = %s
+
+Core and active MOs ::
+
+  n_core_cis = %s
+  n_act_cis  = %s
+
+Dress with MP2 perturbation ::
+
+  mp2_dressing       = %s
+
+Use standard double-excitations ::
+
+  standard_doubles   = %s
+
+Epstein-Nesbet 2x2 diagonalization ::
+
+  en_2_2             = %s
+
 "
         (States_number.to_string b.n_state_cis)
         (Positive_int.to_string b.n_core_cis)

ocaml/Input_cisd_sc2.ml

@@ -61,9 +61,19 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "n_det_max_cisd_sc2 = %s
-pt2_max            = %s
-do_pt2_end         = %s
+    Printf.sprintf "
+Stop when the `n_det` > `n_det_max_cisd_sc2` ::
+
+  n_det_max_cisd_sc2 = %s
+
+Stop when -E(PT2) < `pt2_max` ::
+
+  pt2_max = %s
+
+Compute E(PT2) at the end ::
+
+  do_pt2_end = %s
+
 "
         (Det_number_max.to_string b.n_det_max_cisd_sc2)
         (PT2_energy.to_string b.pt2_max)

ocaml/Input_determinants.ml

@@ -22,6 +22,7 @@ module Determinants : sig
   ;;
   val read : unit -> t
   val to_string : t -> string
+  val debug : t -> string
 end = struct
   type t = 
     { n_int                  : N_int_number.t;
@@ -191,12 +192,12 @@ end = struct
         Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; 1 |] ~data:data
         |> Ezfio.set_determinants_psi_det 
       end  ;
-    let n_int = N_int_number.to_int n_int in
+      (*
     let rec transform accu1 accu2 n_rest = function 
       | [] ->
           let accu1 = List.rev accu1
           |> Array.of_list 
-          |> Determinant.of_int64_array
+          |> Determinant.of_int64_array n_int
           in
           List.rev (accu1::accu2) |> Array.of_list
       | i::rest ->
@@ -205,14 +206,25 @@ end = struct
           else
             let accu1 = List.rev accu1
             |> Array.of_list 
-            |> Determinant.of_int64_array
+            |> Determinant.of_int64_array n_int
             in
+            let n_int = N_int_number.to_int n_int in
             transform [] (accu1::accu2) (2*n_int) rest
     in
-    (Ezfio.get_determinants_psi_det ()).Ezfio.data
-    |> Ezfio.flattened_ezfio_data
-    |> Array.to_list
-    |> transform [] [] (2*n_int)
+    *)
+    let n_int = N_int_number.to_int n_int in
+    let psi_det_array = Ezfio.get_determinants_psi_det () in
+    let dim = psi_det_array.Ezfio.dim
+    and data =  Ezfio.flattened_ezfio_data psi_det_array.Ezfio.data 
+    in
+    assert (n_int = dim.(0));
+    assert (dim.(1) = 2);
+    assert (dim.(2) = (Det_number.to_int (read_n_det ())));
+    List.init dim.(2) ~f:(fun i ->
+      Array.sub ~pos:(2*n_int*i) ~len:(2*n_int) data)
+    |> List.map ~f:(Determinant.of_int64_array
+      (N_int_number.of_int n_int))
+    |> Array.of_list
   ;;
 
   let read () =
@@ -234,30 +246,74 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "read_wf                = %s
-n_det                  = %s
-n_states               = %s
-threshold_generators   = %s
-threshold_selectors    = %s
-n_det_max_jacobi       = %s
-n_states_diag          = %s
-s2_eig                 = %s
-expected_s2            = %s
-mo_label               = \"%s\"
+    let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () 
+    |> MO_number.of_int in
+    let det_text = 
+      List.map2_exn ~f:(fun coef det ->
+        Printf.sprintf "  %f\n%s\n"
+        (Det_coef.to_float coef)
+        (Determinant.to_string ~mo_tot_num:mo_tot_num det 
+         |> String.split ~on:'\n'
+         |> List.map ~f:(fun x -> "  "^x)
+         |> String.concat ~sep:"\n"
+        )
+      ) (Array.to_list b.psi_coef) (Array.to_list b.psi_det)
+      |> String.concat ~sep:"\n"
+    in
+    Printf.sprintf "
+Read the current wave function ::
+
+  read_wf = %s
+
+Label of the MOs on which the determinants were computed ::
+
+  mo_label = %s
+
+Force the selected wave function to be an eigenfunction of S^2.
+If true, input the expected value of S^2 ::
+
+  s2_eig                 = %s
+  expected_s2            = %s
+
+Thresholds on generators and selectors (fraction of the norm) ::
+
+  threshold_generators   = %s
+  threshold_selectors    = %s
+
+Number of requested states, and number of states used for the
+Davidson diagonalization ::
+
+  n_states               = %s
+  n_states_diag          = %s
+
+Maximum size of the Hamiltonian matrix that will be fully diagonalized ::
+
+  n_det_max_jacobi       = %s
+
+Number of determinants ::
+
+  n_det                  = %s
+
+Determinants ::
+
+%s
 "
      (b.read_wf       |> Bool.to_string)
-     (b.n_det         |> Det_number.to_string)
-     (b.n_states      |> States_number.to_string)
-     (b.threshold_generators |> Threshold.to_string)
-     (b.threshold_selectors |> Threshold.to_string)
-     (b.n_det_max_jacobi |> Strictly_positive_int.to_string)
-     (b.n_states_diag |> States_number.to_string)
+     (b.mo_label      |> Non_empty_string.to_string)
      (b.s2_eig        |> Bool.to_string)
      (b.expected_s2   |> Positive_float.to_string)
-     (b.mo_label      |> Non_empty_string.to_string)
+     (b.threshold_generators |> Threshold.to_string)
+     (b.threshold_selectors |> Threshold.to_string)
+     (b.n_states      |> States_number.to_string)
+     (b.n_states_diag |> States_number.to_string)
+     (b.n_det_max_jacobi |> Strictly_positive_int.to_string)
+     (b.n_det         |> Det_number.to_string)
+     det_text
   ;;
 
   let debug b =
+    let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () 
+    |> MO_number.of_int in
     Printf.sprintf "
 n_int                  = %s
 bit_kind               = %s
@@ -288,11 +344,8 @@ psi_det                = %s
      (b.s2_eig        |> Bool.to_string)
      (b.psi_coef  |> Array.to_list |> List.map ~f:Det_coef.to_string
       |> String.concat ~sep:", ")
-     (b.psi_det   |> Array.map ~f:(fun x -> Determinant.to_int64_array x
-      |> Array.map ~f:(fun x-> 
-          Int64.to_string x )|> Array.to_list |>
-       String.concat ~sep:", ") |> Array.to_list
-      |> String.concat ~sep:" | ")
+     (b.psi_det   |> Array.to_list |> List.map ~f:(Determinant.to_string
+       ~mo_tot_num:mo_tot_num) |> String.concat ~sep:"\n\n")
   ;;
 
 end

ocaml/Input_electrons.ml

@@ -47,9 +47,17 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "elec_alpha_num     = %s
-elec_beta_num      = %s
+    Printf.sprintf "
+Spin multiplicity is %s.
+
+Number of alpha and beta electrons ::
+
+  elec_alpha_num = %s
+  elec_beta_num  = %s
+
 "
+        (Multiplicity.of_alpha_beta b.elec_alpha_num b.elec_beta_num
+         |> Multiplicity.to_string)
         (Elec_alpha_number.to_string b.elec_alpha_num)
         (Elec_beta_number.to_string b.elec_beta_num)
   ;;

ocaml/Input_full_ci.ml

@@ -59,9 +59,19 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "n_det_max_fci      = %s
-pt2_max            = %s
-do_pt2_end         = %s
+    Printf.sprintf "
+Stop when the `n_det` > `n_det_max_fci` ::
+
+  n_det_max_fci = %s
+
+Stop when -E(PT2) < `pt2_max` ::
+
+  pt2_max = %s
+
+Compute E(PT2) at the end ::
+
+  do_pt2_end = %s
+
 "
         (Det_number_max.to_string b.n_det_max_fci)
         (PT2_energy.to_string b.pt2_max)

ocaml/Input_hartree_fock.ml

@@ -46,8 +46,15 @@ end = struct
   ;;
 
   let to_string b =
-    Printf.sprintf "n_it_scf_max       = %s
-thresh_scf         = %s
+    Printf.sprintf "
+Max number of SCF iterations ::
+
+  n_it_scf_max = %s
+
+SCF convergence criterion (on energy) ::
+
+  thresh_scf = %s
+
 "
         (Strictly_positive_int.to_string b.n_it_scf_max)
         (Threshold.to_string b.thresh_scf)

ocaml/Input_mo_basis.ml

@@ -71,8 +71,14 @@ end = struct
 
   let to_string b =
     Printf.sprintf "
-mo_label        = %s
-mo_tot_num      = %s
+Label of the molecular orbitals ::
+
+  mo_label = %s
+
+Total number of MOs ::
+
+  mo_tot_num = %s
+
 "
     (Non_empty_string.to_string b.mo_label)
     (MO_number.to_string b.mo_tot_num)

ocaml/Input_nuclei.ml

@@ -11,7 +11,7 @@ module Nuclei : sig
     } with sexp
   ;;
   val read : unit -> t
-  val to_string : t -> string
+  val debug : t -> string
 end = struct
   type t = 
     { nucl_num        : Nucl_number.t ;
@@ -65,7 +65,7 @@ end = struct
     }
   ;;
 
-  let to_string b =
+  let debug b =
     Printf.sprintf "
 nucl_num         = %s
 nucl_label       = %s

ocaml/qp_edit.ml

@@ -3,14 +3,14 @@ open Qptypes;;
 open Core.Std;;
 
 let instructions filename = Printf.sprintf
-"# ===============
-# Quantum Package
-# ===============
-# 
-# File : %s
-#
-# Lines starting with a '#' sign are commented.
-#" filename
+"
+==================================================================
+                       Quantum Package
+==================================================================
+
+Editing file `%s`
+
+" filename
   
 type keyword = 
 | Ao_basis
@@ -37,7 +37,7 @@ let keyword_to_string = function
 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"
+  "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n"
 ;;
 
 let get s = 

ocaml/qptypes_generator.ml

@@ -132,6 +132,7 @@ let input_data = "
 ;;
 
 let untouched = "
+(*
 module Determinant : sig
   type t with sexp
   val to_int64_array : t -> int64 array
@@ -151,6 +152,7 @@ end = struct
     |> List.map ~f:Int64.to_string
     |> String.concat ~sep:\", \"
 end
+*)
 
 "
 

ocaml/test_input.ml

@@ -39,6 +39,7 @@ let test_dets () =
   Ezfio.set_file "F2.ezfio" ;
   let b = Input.Determinants.read ()
   in
+  (*print_endline (Input.Determinants.debug b);*)
   print_endline (Input.Determinants.to_string b);
 ;;
 
@@ -95,5 +96,5 @@ test_hf ();;
 test_mo ();;
 test_nucl ();
 *)
-test_ao ();;
+test_dets();;