diff --git a/run_cas.ml b/run_cas.ml new file mode 100644 index 0000000..05b7f90 --- /dev/null +++ b/run_cas.ml @@ -0,0 +1,82 @@ +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 ""; + doc="Name of the file containing the basis set"; } ; + + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; + + { short='m' ; long="multiplicity" ; opt=Optional; + arg=With_arg ""; + doc="Spin multiplicity (2S+1). Default is singlet"; } ; + + { short='c' ; long="charge" ; opt=Optional; + arg=With_arg ""; + doc="Total charge of the molecule. Default is 0"; } ; + + { short='n' ; long="n_elec" ; opt=Mandatory; + arg=With_arg ""; + doc="Number of electrons in the active space."; } ; + + { short='o' ; long="n_orb" ; opt=Mandatory; + arg=With_arg ""; + 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 -> 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 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 + Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf; + + let mos = + MOBasis.of_hartree_fock hf + in + + let space = + DeterminantSpace.cas_of_mo_basis mos n m + in + let ci = CI.make space in + Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); + + let pt2 = CI.pt2 ci in + Format.fprintf ppf "CAS-EN2 energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s +. pt2) + (* + 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}); + *) +