QCaml/run_cas.ml

let () =
  let open Command_line in
  begin
    set_header_doc (Sys.argv.(0) ^ " - QCaml command");
    set_description_doc "Runs a CAS-CI 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='f' ; long="frozen-core" ; opt=Optional;
          arg=Without_arg ;
          doc="Freeze core MOs. Default is false."; } ;

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

        { short='n' ; long="n_elec" ; opt=Mandatory;
          arg=With_arg "<int>";
          doc="Number of electrons in the active space."; } ;

        { short='o' ; long="n_orb" ; opt=Mandatory;
          arg=With_arg "<int>";
          doc="Number of MOs in the active space."; } ;
      ]
  end;

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

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

  let n = Util.of_some @@ Command_line.get "n_elec" |> int_of_string in
  let m = Util.of_some @@ Command_line.get "n_orb" |> int_of_string 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 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 s =
    Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
  in

  let hf = HartreeFock.make s in
  if Parallel.master then
    Format.fprintf ppf "@[%a@]@." HartreeFock.pp hf;

  let mo_basis =
    MOBasis.of_hartree_fock hf
  in

  let frozen_core =
     Command_line.get_bool "frozen-core"
  in

  let space =
    DeterminantSpace.cas_of_mo_basis ~frozen_core mo_basis n m
  in
  let ci = CI.make space ~n_states:state in
  if Parallel.master then
    Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);

  let pt2 = CI.pt2_mp ci in
  Format.fprintf ppf "CAS-MP2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);


  let pt2 = CI.pt2_en ci in
  Format.fprintf ppf "CAS-EN2 energy             : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);

  let pt2 = CI.pt2_en_reference ci in
  Format.fprintf ppf "CAS-EN2 energy (reference) : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2);

(*
  let variance = CI.variance ci in
  Format.fprintf ppf "CAS variance   : %20.16f@." variance
  *)

  (*
  let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in
  Util.list_range 1 (DeterminantSpace.size space)
  |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i});
  *)

  Format.fprintf ppf "@.";