10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-07-16 08:00:22 +02:00
QCaml/run_Localisation.ml

152 lines
4.4 KiB
OCaml

(* ./run_Localisation.native -b basis -x xyz -c charge -o name.mos -i EZFIOdirectory > name.out *)
(*
#.mos contains the localised orbitales
#.out contains the localization convergence and the analysis of the spatial extent of the orbitales
*)
open Lacaml.D
let read_qp_mo dirname =
let ic = Unix.open_process_in ("zcat "^dirname^"/mo_basis/mo_coef.gz") in
let check = String.trim (input_line ic) in
assert (check = "2");
let int_list =
input_line ic
|> String.split_on_char ' '
|> List.filter (fun x -> x <> "")
|> List.map int_of_string
in
let n_ao, n_mo =
match int_list with
| [ x ; y ] -> x, y
| _ -> assert false
in
let result =
Mat.init_cols n_ao n_mo (fun i j ->
let s = input_line ic in
Scanf.sscanf s " %f " (fun x -> x)
)
in
let exit_code = Unix.close_process_in ic in
assert (exit_code = Unix.WEXITED 0);
result
let () =
let open Command_line in
begin
set_header_doc (Sys.argv.(0) ^ " - QCaml command");
set_description_doc "Localizes MOs";
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='o' ; long="output" ; opt=Mandatory;
arg=With_arg "<string>";
doc="Name of the file containing the localized MOs"; } ;
{ short='i' ; long="import" ; opt=Optional;
arg=With_arg "<string>";
doc="Name of the EZFIO directory containing MOs"; } ;
{ short='s' ; long="separation" ; opt=Optional;
arg=With_arg "<string>";
doc="Type of separation: [ov: occ/vir | sp: sigma/pi | custom] "; } ;
]
end;
let basis_file = Util.of_some @@ Command_line.get "basis" in
let nuclei_file = Util.of_some @@ Command_line.get "xyz" in
let ezfio_file = Command_line.get "import" 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 output_filename = Util.of_some @@ Command_line.get "output" in
let separation =
match Command_line.get "separation" with
| None
| Some "ov" -> Localisation.Occ_vir
| Some "sp" -> Localisation.Sigma_pi
| Some "custom" ->
let list_of_lists =
print_endline "Enter ranges, one per line:";
let rec loop accu =
try
let r =
input_line stdin
|> Range.of_string
|> Range.to_int_list
in
loop (r::accu)
with End_of_file -> List.rev accu
in loop []
in
Localisation.Custom list_of_lists
| _ -> invalid_arg "separation"
in
(* II : Hartree-Fock *)
(* 1. Def pour HF *)
let simulation =
Simulation.of_filenames
~charge ~multiplicity ~nuclei:nuclei_file basis_file
in
(* 2. Calcul de Hartree-Fock*)
let mo_basis =
match ezfio_file with
| Some ezfio_file ->
let mo_coef = read_qp_mo ezfio_file in
let mo_type = MOBasis.Localized "Boys" in
let elec = Simulation.electrons simulation in
let n_mo = Mat.dim2 mo_coef in
let mo_occupation = Vec.init n_mo (fun i ->
if i <= Electrons.n_beta elec then 2.
else if i <= Electrons.n_alfa elec then 1.
else 0.) in
MOBasis.make ~simulation ~mo_type ~mo_occupation ~mo_coef ()
| None -> HartreeFock.make simulation
|> MOBasis.of_hartree_fock
in
let local_mos = Localisation.localize mo_basis separation in
let oc = open_out output_filename in
Printf.fprintf oc "[\n";
Mat.as_vec local_mos
|> Vec.iter (fun x -> Printf.fprintf oc "%20.15e,\n" x);
Printf.fprintf oc "]\n";
close_out oc