QCaml/Nuclei/Nuclei.ml

open Util
open Xyz_ast

type t = (Element.t * Coordinate.t) array

let of_xyz_file ~filename =
  let lexbuf =
    let ic = open_in filename in
    Lexing.from_channel ic
  in
  let data = 
    Xyz_parser.input Nuclei_lexer.read_all lexbuf
  in

  let len = List.length data.nuclei in
  if len <> data.number_of_atoms then
    Printf.sprintf "Error: expected %d atoms but %d read" 
      data.number_of_atoms len
    |> failwith;

  List.map (fun nucleus ->
    nucleus.element, Coordinate.angstrom_to_bohr nucleus.coord
  ) data.nuclei
  |> Array.of_list



let of_zmt_file ~filename =
  let ic = open_in filename in
  let rec aux accu =
    try
      let line = input_line ic in
      aux (line::accu)
    with End_of_file ->
      close_in ic;
      List.rev accu
      |> String.concat "\n"
  in aux []
  |> Zmatrix.of_string
  |> Zmatrix.to_xyz
  |> Array.map (fun (e,x,y,z) -> (e, (Angstrom.make {Point.x ; y ; z} )))
    

let to_string atoms =
    "
                Nuclear Coordinates (Angstrom)
                ------------------------------

-----------------------------------------------------------------------
   Center   Atomic   Element          Coordinates (Angstroms)
            Number                 X             Y             Z
-----------------------------------------------------------------------
" ^
    (Array.mapi (fun i (e, coord) ->
         let coord =
           Coordinate.bohr_to_angstrom coord
         in
         Printf.sprintf " %5d    %5d    %5s   %12.6f  %12.6f  %12.6f"
            (i+1)  (Element.to_int e)  (Element.to_string e) 
            coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z
       ) atoms
     |> Array.to_list
     |> String.concat "\n" ) ^
"
-----------------------------------------------------------------------

"


let of_filename ~filename =
    of_xyz_file filename
  

let repulsion ~nuclei =
  let get_charge e = 
    Element.to_charge e
    |> Charge.to_float
  in
  Array.fold_left ( fun accu (e1, coord1) -> 
    accu +. 
    Array.fold_left (fun accu (e2, coord2) ->
      let r = Coordinate.(norm (coord1 |- coord2)) in
      if r > 0. then
          accu +. 0.5 *. (get_charge e2) *. (get_charge e1) /. r
      else accu
      ) 0. nuclei
    ) 0.  nuclei