open Lacaml.D

let () =
  let open Command_line in
  begin
    set_header_doc (Sys.argv.(0) ^ " - QCaml command");
    set_description_doc "Runs a F12-Full CI calculation";
    set_specs
      [ { short='b' ; long="basis" ; opt=Mandatory;
          arg=With_arg "<string>";
          doc="Name of the file containing the basis set"; } ;

        { short='a' ; long="aux_basis" ; opt=Mandatory;
          arg=With_arg "<string>";
          doc="Name of the file containing the auxiliary basis set"; } ;

        { short='x' ; long="xyz" ; opt=Mandatory;
          arg=With_arg "<string>";
          doc="Name of the file containing the nuclear coordinates in xyz format"; } ;

        { short='m' ; long="multiplicity" ; opt=Optional;
          arg=With_arg "<int>";
          doc="Spin multiplicity (2S+1). Default is singlet"; } ;

        { short='c' ; long="charge" ; opt=Optional;
          arg=With_arg "<int>";
          doc="Total charge of the molecule. Default is 0"; } ;

(*
        { short='e' ; long="expo" ; opt=Optional;
          arg=With_arg "<float>";
          doc="Exponent of the Gaussian geminal"; } ;
*)

        { short='s' ; long="state" ; opt=Optional;
          arg=With_arg "<int>";
          doc="Requested state. Default is 1."; } ;
      ]
  end;

  (* Handle options *)
  let basis_file  = Util.of_some @@ Command_line.get "basis" in

  let aux_basis_filename  = Util.of_some @@ Command_line.get "aux_basis" in

  let nuclei_file = Util.of_some @@ Command_line.get "xyz" in

  let charge = 
    match Command_line.get "charge" with
    | Some x -> ( if x.[0] = 'm' then
                    ~- (int_of_string (String.sub x 1 (String.length x - 1)))
                  else
                    int_of_string x )
    | None   -> 0
  in

(*
  let expo = 
    match Command_line.get "expo" with
    | Some x -> float_of_string x
    | None   -> 1.0
  in
*)

  let state =
    match Command_line.get "state" with
    | Some x -> int_of_string x
    | None -> 1
  in

  let multiplicity =
    match Command_line.get "multiplicity" with
    | Some x -> int_of_string x
    | None -> 1
  in

  let ppf = 
    if Parallel.master then Format.std_formatter
    else Printing.ppf_dev_null
  in

(*
  let f12 =
    F12factor.gaussian_geminal expo
  in
*)
  let f12 =
    F12factor.gaussian_geminal 1.0
  in

  let simulation =
    Simulation.of_filenames ~f12 ~charge ~multiplicity ~nuclei:nuclei_file basis_file
  in

  let hf = HartreeFock.make simulation in
  Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf;

  let mo_basis = 
    MOBasis.of_hartree_fock hf
  in


  let fcif12 =
    F12CI.make ~simulation ~frozen_core:false ~mo_basis ~aux_basis_filename ~state ()
  in

  let ci = F12CI.ci fcif12 in
  Format.fprintf ppf "FCI energy     : ";
  Vec.iteri (fun i x -> if i <= state then
    Format.fprintf ppf "%20.16f@;  " (x +. Simulation.nuclear_repulsion simulation) )
    (CI.eigenvalues ci);
  Format.fprintf ppf "@.";

  let _, e_cif12 = F12CI.eigensystem fcif12 in
  Format.fprintf ppf "FCI-F12 energy : ";
  Vec.iteri (fun i x -> if i <= state then
    Format.fprintf ppf "%20.16f@;  " (x +. Simulation.nuclear_repulsion simulation) )
    e_cif12;
  Format.fprintf ppf "@."