Merge branch 'master' into 'master'

Localization

See merge request scemama/QCaml!1

Basis/Multipole.ml

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

Basis/Multipole.mli

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

run_Localisation.ml

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