From bdece8ea6029b8949a2ed9f9bb1f7d837f25785a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Feb 2021 10:47:34 +0100 Subject: [PATCH] Moving all MOs is not optimal --- docs/mo.html | 90 +++++++++---------- docs/top.html | 10 +-- examples/dune | 1 + examples/ex_localization.ml | 170 ++++++++++++++++++++++++++++++++++++ mo/lib/localization.ml | 32 ++++--- mo/localization.org | 32 ++++--- 6 files changed, 257 insertions(+), 78 deletions(-) create mode 100644 examples/ex_localization.ml diff --git a/docs/mo.html b/docs/mo.html index 8efb47d..64999fe 100644 --- a/docs/mo.html +++ b/docs/mo.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Molecular orbitals @@ -272,46 +272,46 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summmary

+
+

1 Summmary

-
-

2 Frozen core

+
+

2 Frozen core

Defines how the core electrons are frozen, for each atom.

-
-

2.1 Type

+
+

2.1 Type

-
type kind =
+
type kind =
   | All_electron
   | Small
   | Large
@@ -325,8 +325,8 @@ Defines how the core electrons are frozen, for each atom.
-
-

2.2 Creation

+
+

2.2 Creation

val make : kind -> Particles.Nuclei.t -> t
@@ -362,7 +362,7 @@ Defines how the core electrons are frozen, for each atom.
 
 
 
-
+
 let f = Frozen_core.(make Small nuclei) ;;
 val f : Frozen_core.t = [|0; 2; 2; 0|]
 
@@ -372,8 +372,8 @@ val f : Frozen_core.t = [|0; 2; 2; 0|]
-
-

2.3 Access

+
+

2.3 Access

val num_elec : t -> int
@@ -402,7 +402,7 @@ val f : Frozen_core.t = [|0; 2; 2; 0|]
 
 
 
-
+
 Frozen_core.num_elec f ;;
 - : int = 4
 
@@ -412,8 +412,8 @@ Frozen_core.num_mos f ;;
-
-

2.4 Printers

+
+

2.4 Printers

val pp : Format.formatter -> t -> unit
@@ -423,8 +423,8 @@ Frozen_core.num_mos f ;;
-
-

3 Orbital localization

+
+

3 Orbital localization

Molecular orbital localization function. @@ -440,11 +440,11 @@ Edmiston-Rudenberg:

-
-

3.1 Type

+
+

3.1 Type

-
open Linear_algebra
+
open Linear_algebra
 
 type localization_kind =
   | Edmiston
@@ -464,16 +464,16 @@ Edmiston-Rudenberg:
-
-

3.2 Edmiston-Rudenberg

+
+

3.2 Edmiston-Rudenberg

-
-

3.3 Boys

+
+

3.3 Boys

-
-

3.4 Access

+
+

3.4 Access

val kind         : t -> localization_kind
@@ -522,8 +522,8 @@ Edmiston-Rudenberg:
-
-

3.5 Printers

+
+

3.5 Printers

val pp : Format.formatter -> t -> unit
@@ -532,14 +532,14 @@ Edmiston-Rudenberg:
-
-

3.6 Tests

+
+

3.6 Tests

Author: Anthony Scemama

-

Created: 2021-02-03 Wed 09:30

+

Created: 2021-02-03 Wed 10:47

Validate

diff --git a/docs/top.html b/docs/top.html index 6de9436..94b59bb 100644 --- a/docs/top.html +++ b/docs/top.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Top-level @@ -250,18 +250,18 @@ org_html_manager.setup(); // activate after the parameters are set

Table of Contents

-
-

1 Summmary

+
+

1 Summmary

Author: Anthony Scemama

-

Created: 2021-02-03 Wed 09:30

+

Created: 2021-02-03 Wed 10:47

Validate

diff --git a/examples/dune b/examples/dune index 3e5d340..17d6727 100644 --- a/examples/dune +++ b/examples/dune @@ -2,6 +2,7 @@ (names ex_integrals ex_hartree_fock + ex_localization ) (libraries qcaml)) diff --git a/examples/ex_localization.ml b/examples/ex_localization.ml new file mode 100644 index 0000000..4d258bc --- /dev/null +++ b/examples/ex_localization.ml @@ -0,0 +1,170 @@ +(* ./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 Qcaml + +module Command_line = Common.Command_line +module Util = Common.Util +module Range = Common.Range +module Electrons = Particles.Electrons + +open Linear_algebra + +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 = + Matrix.init_cols n_ao n_mo (fun _ _ -> + 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 ""; + 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='o' ; long="output" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the localized MOs"; } ; + + { short='i' ; long="import" ; opt=Optional; + arg=With_arg ""; + doc="Name of the EZFIO directory containing MOs"; } ; + + { short='r' ; long="range" ; opt=Mandatory; + arg=With_arg ""; + doc="Range of orbitals to include in localization"; } ; + + ] + end; + + + + (* II : Hartree-Fock *) + + (* 1. Def pour HF *) + +let basis_file = Util.of_some @@ Command_line.get "basis" in +let nuclei_file = Util.of_some @@ Command_line.get "xyz" in + +let charge = + match Command_line.get "charge" with + | Some x -> ( if x.[0] = 'm' then + ~- (int_of_string (String.sub x 1 (String.length x - 1))) + else + 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 range = + Command_line.get "range" + |> Util.of_some + |> Range.of_string + |> Range.to_int_list + in + + let ezfio_file = Command_line.get "import" in + + let output_filename = Util.of_some @@ Command_line.get "output" in +(* Interpretation:1 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:1]] *) +let nuclei = + Qcaml.Particles.Nuclei.of_xyz_file nuclei_file +in +(* Computation:1 ends here *) + +(* [[file:~/QCaml/examples/ex_hartree_fock.org::*Computation][Computation:2]] *) +let ao_basis = + Qcaml.Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file +in +(* Computation:2 ends here *) + + + let simulation = + Simulation.make ~charge ~multiplicity ~nuclei ao_basis + 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 = Mo.Basis.Localized "Boys" in + let elec = Simulation.electrons simulation in + let n_mo = Matrix.dim2 mo_coef in + let mo_occupation = Vector.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 + Mo.Basis.make ~simulation ~mo_type ~mo_occupation ~mo_coef () + | None -> Mo.Hartree_fock.make simulation + |> Mo.Basis.of_hartree_fock + in + + + let localization = + Mo.Localization.make ~kind:Mo.Localization.Boys mo_basis range + in + + Format.printf "@[%a@]" (Mo.Localization.pp) localization; + + let local_mo_basis = + Mo.Localization.to_basis localization + in + + let local_mos = + Mo.Basis.mo_coef local_mo_basis + in + + let oc = open_out output_filename in + Printf.fprintf oc "[\n"; + Matrix.as_vec local_mos + |> Vector.iter (fun x -> Printf.fprintf oc "%20.15e,\n" x); + Printf.fprintf oc "]\n"; + close_out oc diff --git a/mo/lib/localization.ml b/mo/lib/localization.ml index 61416f4..31850ce 100644 --- a/mo/lib/localization.ml +++ b/mo/lib/localization.ml @@ -253,7 +253,7 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos = let data_initial = let iteration = 0 - and scaling = 0.5 + and scaling = 0.1 in let kappa, loc_value = kappa_loc x in let convergence = abs_float (Matrix.amax kappa) in @@ -263,20 +263,24 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos = in let iteration data = - let iteration = data.iteration + 1 - and x = data.coefficients + let iteration = data.iteration+1 in + let new_kappa, new_loc_value = kappa_loc data.coefficients in + let new_convergence = abs_float (Matrix.amax new_kappa) in + let accept = + new_loc_value >= data.loc_value*.0.98 in - let kappa, loc_value = kappa_loc x in - let convergence = abs_float (Matrix.amax kappa) in - let scaling = - if convergence <= data.convergence then - min 1. (data.scaling *. 1.1) - else - data.scaling *. 0.75 - in - let kappa = Matrix.scale scaling kappa in - let coefficients = next_coef kappa x in - { coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration } + if accept then + let coefficients = next_coef new_kappa data.coefficients in + let scaling = min 1. (data.scaling *. 1.1) in + let kappa = Matrix.scale scaling new_kappa in + let convergence = new_convergence in + let loc_value = new_loc_value in + { coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration } + else + let scaling = + data.scaling *. 0.5 + in + { data with scaling ; iteration } in let array_data = diff --git a/mo/localization.org b/mo/localization.org index c51b93f..c8bd442 100644 --- a/mo/localization.org +++ b/mo/localization.org @@ -313,7 +313,7 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos = let data_initial = let iteration = 0 - and scaling = 0.5 + and scaling = 0.1 in let kappa, loc_value = kappa_loc x in let convergence = abs_float (Matrix.amax kappa) in @@ -323,20 +323,24 @@ let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos = in let iteration data = - let iteration = data.iteration + 1 - and x = data.coefficients + let iteration = data.iteration+1 in + let new_kappa, new_loc_value = kappa_loc data.coefficients in + let new_convergence = abs_float (Matrix.amax new_kappa) in + let accept = + new_loc_value >= data.loc_value*.0.98 in - let kappa, loc_value = kappa_loc x in - let convergence = abs_float (Matrix.amax kappa) in - let scaling = - if convergence <= data.convergence then - min 1. (data.scaling *. 1.1) - else - data.scaling *. 0.75 - in - let kappa = Matrix.scale scaling kappa in - let coefficients = next_coef kappa x in - { coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration } + if accept then + let coefficients = next_coef new_kappa data.coefficients in + let scaling = min 1. (data.scaling *. 1.1) in + let kappa = Matrix.scale scaling new_kappa in + let convergence = new_convergence in + let loc_value = new_loc_value in + { coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration } + else + let scaling = + data.scaling *. 0.5 + in + { data with scaling ; iteration } in let array_data =