Added xyz types in symmetry (ocaml)

ocaml/atom.ml

@@ -5,6 +5,7 @@ exception AtomError of string
 module Charge : sig
   type t 
   val to_float : t -> float
+  val to_int   : t -> int
   val to_string: t -> string 
   val of_float : float -> t
   val of_int   : int -> t
@@ -12,6 +13,7 @@ module Charge : sig
 end = struct
   type t = float
   let to_float x = x
+  let to_int   x = Float.to_int x
   let to_string x = Float.to_string (to_float x)
   let of_float x = x
   let of_int   i = Float.of_int i

ocaml/basis.ml

@@ -0,0 +1,34 @@
+open Core.Std;;
+
+type t = Gto.t list;;
+
+(** Read all the basis functions of an element *)
+let read in_channel =
+  let rec read result =
+    try
+      let gto = Gto.read_one in_channel in
+      read (gto::result)
+    with
+    | Gto.End_Of_Basis -> List.rev result
+  in read []  
+;;
+
+(** Find an element in the basis set file *)
+let find in_channel element =
+  let element_read = ref Element.X in
+  while !element_read <> element
+  do
+    let buffer = input_line in_channel in
+    try
+      element_read := Element.of_string buffer
+    with
+    | Element.ElementError _ -> ()
+  done ;
+  !element_read
+;;
+
+let read_element in_channel element =
+  ignore (find in_channel element) ;
+  read in_channel ;
+;;
+

ocaml/create_ezfio_from_xyz.ml

@@ -0,0 +1,93 @@
+open Qputils;;
+open Qptypes;;
+open Core.Std;;
+
+let spec =
+  let open Command.Spec in
+  empty 
+  +> flag "o"  (optional string)
+     ~doc:"file Name of the created EZFIO file"
+  +> flag "b" (required string)
+     ~doc:"name Basis set."
+  +> flag "c" (optional_with_default 0 int)
+     ~doc:"int Total charge of the molecule. Default is 0."
+  +> flag "m"  (optional_with_default 1 int)
+     ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1"
+  +> anon ("xyz_file" %: string)
+;;
+
+let run ?o b c m xyz_file =
+
+  (* Open basis set channel *)
+  let basis_channel =
+    In_channel.create
+      (Qpackage.root / "data/basis" / (String.lowercase b))
+  in
+
+  (* Read molecule *)
+  let molecule =
+    Molecule.of_xyz_file xyz_file
+  in
+
+  (* Build EZFIO File name *)
+  let ezfio_file =
+    match o with
+    | Some x -> x
+    | None ->
+      begin
+        match String.rsplit2 ~on:'.' xyz_file with
+        | Some (x,"xyz") -> x^".ezfio"
+        | _ -> xyz_file^".ezfio"
+      end
+  in
+  if Sys.file_exists_exn ezfio_file then
+    failwith (ezfio_file^" already exists");
+
+  (* Create EZFIO *)
+  Ezfio.set_file ezfio_file;
+
+  (* Write Electrons *)
+  Ezfio.set_electrons_elec_alpha_num ( Positive_int.to_int
+    molecule.Molecule.elec_alpha ) ;
+  Ezfio.set_electrons_elec_beta_num  ( Positive_int.to_int
+    molecule.Molecule.elec_beta  ) ;
+
+  (* Write Nuclei *)
+  let nuclei = molecule.Molecule.nuclei in
+  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 Basis set *)
+;;
+
+let command = 
+    Command.basic 
+    ~summary: "Quantum Package command"
+    ~readme:(fun () ->
+      "Creates an EZFIO directory from a standard xyz file
+      ")
+    spec
+    (fun o b c m xyz_file () ->
+       run ?o b c m xyz_file )
+;;
+
+let () =
+    Command.run command
+;;
+
+
+

ocaml/gto.ml

@@ -65,35 +65,6 @@ let read_one in_channel =
   |> of_prim_coef_list
 ;;
 
-(** Read all the basis functions of an element *)
-let read_basis in_channel =
-  let rec read result =
-    try
-      let gto = read_one in_channel in
-      read (gto::result)
-    with
-    | End_Of_Basis -> List.rev result
-  in read []  
-;;
-
-(** Find an element in the basis set file *)
-let find in_channel element =
-  let element_read = ref Element.X in
-  while !element_read <> element
-  do
-    let buffer = input_line in_channel in
-    try
-      element_read := Element.of_string buffer
-    with
-    | Element.ElementError _ -> ()
-  done ;
-  !element_read
-;;
-
-let read_basis_of_element in_channel element =
-  ignore (find in_channel element) ;
-  read_basis in_channel ;
-;;
 
 (** Transform the gto to a string *)
 let to_string { sym = sym ; lc = lc } =

ocaml/molecule.ml

@@ -22,6 +22,10 @@ let get_multiplicity m =
   Multiplicity.of_alpha_beta m.elec_alpha m.elec_beta
 ;;
 
+let get_nucl_num m =
+  Strictly_positive_int.of_int (List.length m.nuclei)
+;;
+
 let name m = 
   let cm = Float.to_int (get_charge m) in
   let c = 

ocaml/qpackage.ml

@@ -0,0 +1,11 @@
+open Core.Std;;
+
+(** Variables related to the quantum package installation *)
+
+let root = 
+  match (Sys.getenv "QPACKAGE_ROOT") with
+  | None -> failwith "QPACKAGE_ROOT environment variable is not set.
+Please source the quantum_package.rc file."
+  | Some x -> x 
+;;
+

ocaml/qputils.ml

@@ -1 +1,3 @@
 let (/) (a:string) (b:string) = a^"/"^b;;
+
+

ocaml/symmetry.ml

@@ -1,3 +1,6 @@
+open Qptypes;;
+open Core.Std;;
+
 type t = S|P|D|F|G|H|I|J|K|L
 
 let to_string = function
@@ -38,3 +41,124 @@ let of_char = function
   | 'K' -> K
   | 'L' -> L
   | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L"))
+
+let to_l = function
+  | S -> Positive_int.of_int 0
+  | P -> Positive_int.of_int 1
+  | D -> Positive_int.of_int 2
+  | F -> Positive_int.of_int 3
+  | G -> Positive_int.of_int 4
+  | H -> Positive_int.of_int 5
+  | I -> Positive_int.of_int 6
+  | J -> Positive_int.of_int 7
+  | K -> Positive_int.of_int 8
+  | L -> Positive_int.of_int 9
+
+module Xyz : sig
+  type t
+  val  of_string : string -> t
+  val  to_string : t -> string 
+  val  get_l     : t -> Positive_int.t
+end = struct
+  type t = { x: Positive_int.t ;
+             y: Positive_int.t ;
+             z: Positive_int.t }
+
+  type state_type = Null | X | Y | Z
+  (* Input string is like "x2z3" *)
+
+  let of_string s =
+    let flush state accu number =
+      let n = 
+        if (number = "") then 0
+        else (Int.of_string number) 
+      in
+      match state with
+      | X -> { x= Positive_int.(of_int ( (to_int accu.x) +n));
+               y= accu.y ;
+               z= accu.z }
+      | Y -> { x= accu.x ; 
+               y= Positive_int.(of_int ( (to_int accu.y) +n));
+               z= accu.z }
+      | Z -> { x= accu.x ; 
+               y= accu.y ;
+               z= Positive_int.(of_int ( (to_int accu.z) +n))}
+      | Null -> accu
+    in
+    let rec do_work state accu number = function
+    | [] -> flush state accu number
+    | 'X'::rest | 'x'::rest -> 
+        let new_accu = flush state accu number in 
+        do_work X new_accu "" rest
+    | 'Y'::rest | 'y'::rest -> 
+        let new_accu = flush state accu number in 
+        do_work Y new_accu "" rest
+    | 'Z'::rest | 'z'::rest -> 
+        let new_accu = flush state accu number in 
+        do_work Z new_accu "" rest
+    | c::rest -> do_work state accu (number^(String.of_char c)) rest
+    in
+    String.to_list_rev s 
+    |> List.rev 
+    |> do_work Null 
+     { x=Positive_int.of_int 0 ; 
+       y=Positive_int.of_int 0 ;
+       z=Positive_int.of_int 0 } "" 
+  ;;
+
+  let to_string t = 
+    let x = match (Positive_int.to_int t.x) with
+    | 0 -> ""
+    | 1 -> "x"
+    | i -> Printf.sprintf "x%d" i
+    and y = match (Positive_int.to_int t.y) with
+    | 0 -> ""
+    | 1 -> "y"
+    | i -> Printf.sprintf "y%d" i
+    and z = match (Positive_int.to_int t.z) with
+    | 0 -> ""
+    | 1 -> "z"
+    | i -> Printf.sprintf "z%d" i
+    in
+    x^y^z
+  ;;
+
+  let get_l t =
+   let x = Positive_int.to_int t.x
+   and y = Positive_int.to_int t.y
+   and z = Positive_int.to_int t.z
+   in Positive_int.of_int (x+y+z)
+ ;;
+
+ let of_symmetry sym =
+   let l = Positive_int.to_int (to_l sym) in
+   let create_z xyz = 
+    { x=xyz.x ;
+      y=xyz.y ;
+      z=Positive_int.(of_int (l-((to_int xyz.x)+(to_int xyz.y))))
+    }
+   in
+   let rec create_y accu xyz =
+     let {x ; y ; z} = xyz in
+     match (Positive_int.to_int y) with
+     | 0 -> (create_z xyz)::accu
+     | i ->
+       let ynew = Positive_int.( (to_int y)-1 |> of_int) in
+       create_y ( (create_z xyz)::accu) { x ; y=ynew ; z}
+   in
+   let rec create_x accu xyz =
+     let {x ; y ; z} = xyz in
+     match (Positive_int.to_int x) with
+     | 0 -> (create_y [] xyz)@accu
+     | i ->
+       let xnew = Positive_int.( (to_int x)-1 |> of_int) in
+       let ynew = Positive_int.(l-(to_int xnew) |> of_int)
+       in
+       create_x ((create_y [] xyz)@accu) { x=xnew ; y=ynew ; z}
+   in
+   create_x [] { x=(to_l sym) ; y=Positive_int.of_int 0 ;
+     z=Positive_int.of_int 0 }
+   ;;
+
+end
+

ocaml/test_gto.ml

@@ -10,7 +10,7 @@ let test_prim () =
 ;;
   
 let test_gto_1 () =
-  let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pVDZ" in
+  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);
@@ -21,15 +21,15 @@ let test_gto_1 () =
 ;;
 
 let test_gto_2 () =
-  let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pVDZ" in
+  let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
   ignore (input_line in_channel);
-  let basis = Gto.read_basis in_channel in
+  let basis = Basis.read in_channel in
   List.iter basis ~f:(fun x-> Printf.printf "%s\n" (Gto.to_string x))
 ;;
 
 let test_gto () =
-  let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pVDZ" in
+  let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in
-  let basis = Gto.read_basis_of_element in_channel Element.C in
+  let basis = Basis.read_element in_channel Element.C in
   List.iter basis ~f:(fun x-> Printf.printf "%s\n" (Gto.to_string x))
 ;;
 

src/NEEDED_MODULES

@@ -0,0 +1 @@
+Utils

src/README.rst

@@ -116,3 +116,11 @@ Documentation
 
 
 
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+

src/Utils/README.rst

@@ -10,62 +10,6 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`apply_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L168>`_
-  Apply the rotation found by find_rotation
-
-`find_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L149>`_
-  Find A.C = B
-
-`get_pseudo_inverse <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L95>`_
-  Find C = A^-1
-
-`lapack_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L247>`_
-  Diagonalize matrix H
-  .br
-  H is untouched between input and ouptut
-  .br
-  eigevalues(i) = ith lowest eigenvalue of the H matrix
-  .br
-  eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
-  .br
-
-`lapack_diagd <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L180>`_
-  Diagonalize matrix H
-  .br
-  H is untouched between input and ouptut
-  .br
-  eigevalues(i) = ith lowest eigenvalue of the H matrix
-  .br
-  eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
-  .br
-
-`lapack_partial_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L310>`_
-  Diagonalize matrix H
-  .br
-  H is untouched between input and ouptut
-  .br
-  eigevalues(i) = ith lowest eigenvalue of the H matrix
-  .br
-  eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
-  .br
-
-`ortho_lowdin <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L1>`_
-  Compute C_new=C_old.S^-1/2 canonical orthogonalization.
-  .br
-  overlap : overlap matrix
-  .br
-  LDA : leftmost dimension of overlap array
-  .br
-  N : Overlap matrix is NxN (array is (LDA,N) )
-  .br
-  C : Coefficients of the vectors to orthogonalize. On exit,
-  orthogonal vectors
-  .br
-  LDC : leftmost dimension of C
-  .br
-  m : Coefficients matrix is MxN, ( array is (LDC,N) )
-  .br
-
 `abort_all <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/abort.irp.f#L1>`_
   If True, all the calculation is aborted