mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-09-01 17:23:41 +02:00
merge localization
This commit is contained in:
parent
7b9bb6aca6
commit
666fd04369
@ -24,6 +24,15 @@ let matrix_z2 t = t.(5)
|
||||
let matrix_xy t = t.(6)
|
||||
let matrix_yz t = t.(7)
|
||||
let matrix_zx t = t.(8)
|
||||
let matrix_x3 t = t.(9)
|
||||
let matrix_y3 t = t.(10)
|
||||
let matrix_z3 t = t.(11)
|
||||
let matrix_x4 t = t.(12)
|
||||
let matrix_y4 t = t.(13)
|
||||
let matrix_z4 t = t.(14)
|
||||
|
||||
|
||||
|
||||
|
||||
let cutoff = integrals_cutoff
|
||||
|
||||
@ -37,7 +46,7 @@ let to_powers x =
|
||||
let contracted_class shell_a shell_b : float Zmap.t array =
|
||||
|
||||
match Csp.make shell_a shell_b with
|
||||
| None -> Array.init 9 (fun _ -> Zmap.create 0)
|
||||
| None -> Array.init 15 (fun _ -> Zmap.create 0)
|
||||
| Some shell_p ->
|
||||
begin
|
||||
|
||||
@ -45,7 +54,7 @@ let contracted_class shell_a shell_b : float Zmap.t array =
|
||||
let class_indices = Csp.zkey_array shell_p in
|
||||
|
||||
let contracted_class =
|
||||
Array.init 9 (fun _ -> Array.make (Array.length class_indices) 0.)
|
||||
Array.init 15 (fun _ -> Array.make (Array.length class_indices) 0.)
|
||||
in
|
||||
|
||||
let a_minus_b =
|
||||
@ -97,10 +106,26 @@ let contracted_class shell_a shell_b : float Zmap.t array =
|
||||
Overlap_primitives.hvrr (Po.get xyz angMomA + 2, Po.get xyz angMomB)
|
||||
expo_inv
|
||||
(Co.get xyz a_minus_b, Co.get xyz center_pa)
|
||||
|
||||
in
|
||||
(* 1D <i|(x-Xa)^3|j> *)
|
||||
let j k =
|
||||
let xyz = xyz_of_int k in
|
||||
Overlap_primitives.hvrr (Po.get xyz angMomA + 3, Po.get xyz angMomB)
|
||||
expo_inv
|
||||
(Co.get xyz a_minus_b, Co.get xyz center_pa)
|
||||
in
|
||||
(* 1D <i|(x-Xa)^4|j> *)
|
||||
let l k =
|
||||
let xyz = xyz_of_int k in
|
||||
Overlap_primitives.hvrr (Po.get xyz angMomA + 4, Po.get xyz angMomB)
|
||||
expo_inv
|
||||
(Co.get xyz a_minus_b, Co.get xyz center_pa)
|
||||
|
||||
in
|
||||
let norm = norm_coef_scales.(i) in
|
||||
let f0, f1, f2, g0, g1, g2, h0, h1, h2 =
|
||||
f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2
|
||||
let f0, f1, f2, g0, g1, g2, h0, h1, h2, j0, j1, j2 , l0, l1, l2 =
|
||||
f 0, f 1, f 2, g 0, g 1, g 2, h 0, h 1, h 2, j 0, j 1, j 2, l 0, l 1, l 2
|
||||
in
|
||||
let x = g0 +. f0 *. xa in
|
||||
let y = g1 +. f1 *. ya in
|
||||
@ -108,6 +133,13 @@ let contracted_class shell_a shell_b : float Zmap.t array =
|
||||
let x2 = h0 +. xa *. (2. *. x -. xa *. f0) in
|
||||
let y2 = h1 +. ya *. (2. *. y -. ya *. f1) in
|
||||
let z2 = h2 +. za *. (2. *. z -. za *. f2) in
|
||||
let x3 = j0 +. xa *. f0 *. (3. *. x2 -. 3. *. x *. xa +. xa *. xa) in
|
||||
let y3 = j1 +. ya *. f1 *. (3. *. y2 -. 3. *. y *. ya +. ya *. ya) in
|
||||
let z3 = j2 +. za *. f2 *. (3. *. z2 -. 3. *. z *. za +. za *. za) in
|
||||
let x4 = l0 +. xa *. f0 *. ( 4. *. x3 -. 6. *. x2 *. xa +. 4. *. x *. xa *. xa -. xa *. xa *. xa) in
|
||||
let y4 = l1 +. ya *. f1 *. ( 4. *. y3 -. 6. *. y2 *. ya +. 4. *. y *. ya *. ya -. ya *. ya *. ya) in
|
||||
let z4 = l2 +. za *. f2 *. ( 4. *. z3 -. 6. *. z2 *. za +. 4. *. z *. za *. za -. za *. za *. za) in
|
||||
|
||||
let c = contracted_class in
|
||||
let d = coef_prod *. norm in
|
||||
c.(0).(i) <- c.(0).(i) +. d *. x *. f1 *. f2;
|
||||
@ -119,6 +151,13 @@ let contracted_class shell_a shell_b : float Zmap.t array =
|
||||
c.(6).(i) <- c.(6).(i) +. d *. x *. y *. f2;
|
||||
c.(7).(i) <- c.(7).(i) +. d *. f0 *. y *. z;
|
||||
c.(8).(i) <- c.(8).(i) +. d *. x *. f1 *. z;
|
||||
c.(9).(i) <- c.(9).(i) +. d *. x3 *. f1 *. f2;
|
||||
c.(10).(i) <- c.(10).(i) +. d *. f0 *. y3 *. f2;
|
||||
c.(11).(i) <- c.(11).(i) +. d *. f0 *. f1 *. z3;
|
||||
c.(12).(i) <- c.(12).(i) +. d *. x4 *. f1 *. f2;
|
||||
c.(13).(i) <- c.(13).(i) +. d *. f0 *. y4 *. f2;
|
||||
c.(14).(i) <- c.(14).(i) +. d *. f0 *. f1 *. z4;
|
||||
|
||||
) class_indices
|
||||
end
|
||||
) (Csp.coefs_and_shell_pairs shell_p);
|
||||
@ -147,7 +186,7 @@ let of_basis basis =
|
||||
and shell = Bs.contracted_shells basis
|
||||
in
|
||||
|
||||
let result = Array.init 9 (fun _ -> Mat.create n n) in
|
||||
let result = Array.init 15 (fun _ -> Mat.create n n) in
|
||||
for j=0 to (Array.length shell) - 1 do
|
||||
for i=0 to j do
|
||||
(* Compute all the integrals of the class *)
|
||||
@ -155,7 +194,7 @@ let of_basis basis =
|
||||
contracted_class shell.(i) shell.(j)
|
||||
in
|
||||
|
||||
for k=0 to 8 do
|
||||
for k=0 to 14 do
|
||||
Array.iteri (fun j_c powers_j ->
|
||||
let j_c = Cs.index shell.(j) + j_c + 1 in
|
||||
let xj = to_powers powers_j in
|
||||
@ -176,7 +215,7 @@ let of_basis basis =
|
||||
done;
|
||||
done;
|
||||
done;
|
||||
for k=0 to 8 do
|
||||
for k=0 to 14 do
|
||||
Mat.detri result.(k);
|
||||
done;
|
||||
result
|
||||
|
@ -31,4 +31,22 @@ val matrix_y2 : t -> Mat.t
|
||||
val matrix_z2 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | z^2 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val matrix_x3 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | x^3 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val matrix_y3 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | y^3 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val matrix_z3 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | z^3 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val matrix_x4 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | x^4 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val matrix_y4 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | y^4 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val matrix_z4 : t -> Mat.t
|
||||
(** {% $$ \langle \chi_i | z^4 | \chi_j \rangle $$ %} *)
|
||||
|
||||
val of_basis : Basis.t -> t
|
||||
|
1376
MOBasis/Localisation.ml
Normal file
1376
MOBasis/Localisation.ml
Normal file
File diff suppressed because it is too large
Load Diff
124
run_Localisation.ml
Normal file
124
run_Localisation.ml
Normal file
@ -0,0 +1,124 @@
|
||||
(* ./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"; } ;
|
||||
|
||||
]
|
||||
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
|
||||
|
||||
|
||||
(* 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 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
|
Loading…
Reference in New Issue
Block a user