quantum_package/ocaml/Molecule.ml

open Core.Std ;;
open Qptypes ;;

exception MultiplicityError of string;;

type t = {
  nuclei     : Atom.t list ;
  elec_alpha : Elec_alpha_number.t ;
  elec_beta  : Elec_beta_number.t ;
} with sexp

let get_charge { nuclei  ; elec_alpha ; elec_beta } =
  let result =
     (Elec_alpha_number.to_int elec_alpha) +
     (Elec_beta_number.to_int elec_beta)
  in
  let rec nucl_charge = function
  | a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest
  | [] -> 0.
  in
  Charge.of_float (nucl_charge nuclei  -. (Float.of_int result))


let get_multiplicity m = 
  let elec_alpha =
     m.elec_alpha
  in
  Multiplicity.of_alpha_beta elec_alpha m.elec_beta


let get_nucl_num m =
  let nmax =
    List.length m.nuclei 
  in
  Nucl_number.of_int nmax ~max:nmax


let name m = 
  let cm =
    get_charge m
    |> Charge.to_int
  in
  let c = 
     match cm with
     | 0 -> ""
     | 1 -> " (+)"
     | (-1) -> " (-)"
     | i when i>1 -> Printf.sprintf " (%d+)" i
     | i -> Printf.sprintf " (%d-)" (-i)
  in
  let mult =
    get_multiplicity m
    |> Multiplicity.to_string
  in
  let { nuclei  ; elec_alpha ; elec_beta } = m
  in
  let rec build_list accu = function
  | a::rest ->
      begin
        let e = a.Atom.element in
        match (List.Assoc.find accu e) with
        | None   -> build_list (List.Assoc.add accu e 1) rest
        | Some i -> build_list (List.Assoc.add accu e (i+1)) rest
      end
  | [] -> accu
  in
  let rec build_name accu = function
  | (a, n)::rest ->
    let a =
      Element.to_string a
    in
    begin
      match n with 
      | 1 -> build_name (a::accu) rest
      | i when i>1 -> 
        let tmp = Printf.sprintf "%s%d" a i
        in build_name (tmp::accu) rest
      | _ -> assert false
    end
  | [] -> accu
  in
  let result =
     build_list [] nuclei |> build_name [c ; ", " ; mult]
  in
  String.concat (result)


let to_string_general ~f m = 
  let { nuclei  ; elec_alpha ; elec_beta } = m
  in
  let n =
    List.length nuclei
  in
  let title =
     name m
  in
  [ Int.to_string n ; title ] @  (List.map ~f nuclei)
  |> String.concat ~sep:"\n"

let to_string =
  to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x)

let to_xyz = 
  to_string_general ~f:Atom.to_xyz 


let of_xyz_string
    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
    ?(units=Units.Angstrom)
    s =
  let l = String.split s ~on:'\n'
       |> List.filter ~f:(fun x -> x <> "")
       |> List.map ~f:(fun x -> Atom.of_string units x)
  in
  let ne = ( get_charge { 
        nuclei=l ;
        elec_alpha=(Elec_alpha_number.of_int 1) ;
        elec_beta=(Elec_beta_number.of_int 0) } 
      |> Charge.to_int 
      ) + 1 - (Charge.to_int charge)
      |> Elec_number.of_int 
  in
  let (na,nb) =
     Multiplicity.to_alpha_beta ne multiplicity
  in
  let result = 
  { nuclei = l ;
    elec_alpha = na ;
    elec_beta  = nb }
  in
  if ((get_multiplicity result) <> multiplicity) then
     let msg = Printf.sprintf
      "With %d electrons multiplicity %d is impossible"
      (Elec_number.to_int ne)
      (Multiplicity.to_int multiplicity)
     in
     raise (MultiplicityError msg);
  else () ;
  result



let of_xyz_file
    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
    ?(units=Units.Angstrom)
    filename =
  let (_,buffer) = In_channel.read_all filename 
  |> String.lsplit2_exn ~on:'\n' in
  let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
  of_xyz_string ~charge ~multiplicity ~units buffer


let of_zmt_file  
    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
    ?(units=Units.Angstrom)
    filename = 
  In_channel.read_all filename
  |> Zmatrix.of_string
  |> Zmatrix.to_xyz_string
  |> of_xyz_string ~charge ~multiplicity ~units


let of_file 
    ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
    ?(units=Units.Angstrom)
    filename = 
  try
    of_xyz_file ~charge ~multiplicity ~units filename
  with _ ->
    of_zmt_file ~charge ~multiplicity ~units filename


let distance_matrix molecule =
  let coord =
    molecule.nuclei
    |> List.map ~f:(fun x -> x.Atom.coord)
    |> Array.of_list
  in
  let n = 
    Array.length coord
  in
  let result =
    Array.make_matrix ~dimx:n  ~dimy:n 0.
  in
  for i = 0 to (n-1)
  do
    for j = 0 to (n-1)
    do
      result.(i).(j) <- Point3d.distance coord.(i) coord.(j)
    done;
  done;
  result
  
  
  

include To_md5
let to_md5 = to_md5 sexp_of_t