LCPQ
/
QCaml
mirror of https://gitlab.com/scemama/QCaml.git
10
1
Fork 0
Code Releases Activity

Merge branch 'master' of gitlab.com:scemama/QCaml

This commit is contained in:
Anthony Scemama 2020-07-09 22:56:34 +02:00
parent a5c4851814 9df479a02f
commit 0f86ece4f3
4 changed files with 1564 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
Basis/Multipole.ml
View File

@ -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

18
Basis/Multipole.mli
View File

@ -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
View File

File diff suppressed because it is too large Load Diff

124
run_Localisation.ml Normal file
View 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