QCaml/run_mp2.ml

open Lacaml.D 

let () =
  let open Command_line in
  begin
    set_header_doc (Sys.argv.(0) ^ " - QCaml command");
    set_description_doc "Runs a Hartree-Fock calculation";
    set_specs
      [ { short='b' ; long="basis" ; opt=Mandatory;
          arg=With_arg "<string>";
          doc="Name of the file containing the 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='g' ; long="guess" ; opt=Optional;
          arg=With_arg "<string>";
          doc="Guess with another basis set."; } ;
      ]
  end;

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

  let basis_file  = Util.of_some @@ Command_line.get "basis" in

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

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

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

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

  let guess = 
    match Command_line.get "guess" with
    | None -> `Huckel
    | Some filename ->
      let s_guess = Simulation.of_filenames
          ~charge ~multiplicity ~nuclei:nuclei_file filename in
      let hf = HartreeFock.make s_guess in
      Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf;
      let guess_mos =
        MOBasis.of_hartree_fock hf
        |> MOBasis.of_mo_basis s
      in
      `Matrix (MOBasis.mo_coef guess_mos)
  in

  let hf = HartreeFock.make ~guess s in

  Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf;

  let e_hf = HartreeFock.energy hf in

  let mp2 = MP2.make ~frozen_core:true hf in
  Format.fprintf ppf "@[MP2   = %15.10f@]@." mp2;
  Format.fprintf ppf "@[E+MP2 = %15.10f@]@." (mp2 +. e_hf)