Moving all MOs is not optimal

This commit is contained in:
Anthony Scemama 2021-02-03 10:47:34 +01:00
parent 254e8b5dc8
commit bdece8ea60
6 changed files with 257 additions and 78 deletions

View File

@ -3,7 +3,7 @@
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2021-02-03 Wed 09:30 -->
<!-- 2021-02-03 Wed 10:47 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Molecular orbitals</title>
@ -272,46 +272,46 @@ org_html_manager.setup(); // activate after the parameters are set
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org0347fac">1. Summmary</a></li>
<li><a href="#org3e52ccb">2. Frozen core</a>
<li><a href="#org16f6c6b">1. Summmary</a></li>
<li><a href="#orgebbcfe0">2. Frozen core</a>
<ul>
<li><a href="#org314a379">2.1. Type</a></li>
<li><a href="#org3f41ebf">2.2. Creation</a></li>
<li><a href="#org30d9e55">2.3. Access</a></li>
<li><a href="#orgbfbb067">2.4. Printers</a></li>
<li><a href="#org61f3b2e">2.1. Type</a></li>
<li><a href="#orgfff2d7e">2.2. Creation</a></li>
<li><a href="#org2519af7">2.3. Access</a></li>
<li><a href="#orga5fe616">2.4. Printers</a></li>
</ul>
</li>
<li><a href="#org71f8a17">3. Orbital localization</a>
<li><a href="#org03637a8">3. Orbital localization</a>
<ul>
<li><a href="#orgd6cd073">3.1. Type</a></li>
<li><a href="#org5abbc4c">3.2. Edmiston-Rudenberg</a></li>
<li><a href="#org1360ab2">3.3. Boys</a></li>
<li><a href="#org44850fb">3.4. Access</a></li>
<li><a href="#orgb1a2fc0">3.5. Printers</a></li>
<li><a href="#orgc751f7e">3.6. Tests</a></li>
<li><a href="#orge723841">3.1. Type</a></li>
<li><a href="#org6933918">3.2. Edmiston-Rudenberg</a></li>
<li><a href="#org44c9408">3.3. Boys</a></li>
<li><a href="#org64dada8">3.4. Access</a></li>
<li><a href="#org5b86753">3.5. Printers</a></li>
<li><a href="#orgcfa3655">3.6. Tests</a></li>
</ul>
</li>
</ul>
</div>
</div>
<div id="outline-container-org0347fac" class="outline-2">
<h2 id="org0347fac"><span class="section-number-2">1</span> Summmary</h2>
<div id="outline-container-org16f6c6b" class="outline-2">
<h2 id="org16f6c6b"><span class="section-number-2">1</span> Summmary</h2>
</div>
<div id="outline-container-org3e52ccb" class="outline-2">
<h2 id="org3e52ccb"><span class="section-number-2">2</span> Frozen core</h2>
<div id="outline-container-orgebbcfe0" class="outline-2">
<h2 id="orgebbcfe0"><span class="section-number-2">2</span> Frozen core</h2>
<div class="outline-text-2" id="text-2">
<p>
Defines how the core electrons are frozen, for each atom.
</p>
</div>
<div id="outline-container-org314a379" class="outline-3">
<h3 id="org314a379"><span class="section-number-3">2.1</span> Type</h3>
<div id="outline-container-org61f3b2e" class="outline-3">
<h3 id="org61f3b2e"><span class="section-number-3">2.1</span> Type</h3>
<div class="outline-text-3" id="text-2-1">
<div class="org-src-container">
<pre class="src src-ocaml" id="org46bfdf8"><span class="org-tuareg-font-lock-governing">type</span> <span class="org-type">kind</span> <span class="org-tuareg-font-lock-operator">=</span>
<pre class="src src-ocaml" id="org786e5f6"><span class="org-tuareg-font-lock-governing">type</span> <span class="org-type">kind</span> <span class="org-tuareg-font-lock-operator">=</span>
<span class="org-tuareg-font-lock-operator">|</span> <span class="org-tuareg-font-lock-constructor">All_electron</span>
<span class="org-tuareg-font-lock-operator">|</span> <span class="org-tuareg-font-lock-constructor">Small</span>
<span class="org-tuareg-font-lock-operator">|</span> Large
@ -325,8 +325,8 @@ Defines how the core electrons are frozen, for each atom.
</div>
</div>
<div id="outline-container-org3f41ebf" class="outline-3">
<h3 id="org3f41ebf"><span class="section-number-3">2.2</span> Creation</h3>
<div id="outline-container-orgfff2d7e" class="outline-3">
<h3 id="orgfff2d7e"><span class="section-number-3">2.2</span> Creation</h3>
<div class="outline-text-3" id="text-2-2">
<div class="org-src-container">
<pre class="src src-ocaml"><span class="org-tuareg-font-lock-governing">val</span> <span class="org-function-name">make</span> <span class="org-tuareg-font-lock-operator">:</span> kind <span class="org-tuareg-font-lock-operator">-&gt;</span> <span class="org-tuareg-font-lock-module">Particles.Nuclei.</span>t <span class="org-tuareg-font-lock-operator">-&gt;</span> t
@ -362,7 +362,7 @@ Defines how the core electrons are frozen, for each atom.
</tbody>
</table>
<pre class="example" id="org41e22ed">
<pre class="example" id="org863654b">
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|]
</div>
</div>
<div id="outline-container-org30d9e55" class="outline-3">
<h3 id="org30d9e55"><span class="section-number-3">2.3</span> Access</h3>
<div id="outline-container-org2519af7" class="outline-3">
<h3 id="org2519af7"><span class="section-number-3">2.3</span> Access</h3>
<div class="outline-text-3" id="text-2-3">
<div class="org-src-container">
<pre class="src src-ocaml"><span class="org-tuareg-font-lock-governing">val</span> <span class="org-function-name">num_elec</span> <span class="org-tuareg-font-lock-operator">:</span> t <span class="org-tuareg-font-lock-operator">-&gt;</span> int
@ -402,7 +402,7 @@ val f : Frozen_core.t = [|0; 2; 2; 0|]
</tbody>
</table>
<pre class="example" id="org64fff1a">
<pre class="example" id="org517ff86">
Frozen_core.num_elec f ;;
- : int = 4
@ -412,8 +412,8 @@ Frozen_core.num_mos f ;;
</div>
</div>
<div id="outline-container-orgbfbb067" class="outline-3">
<h3 id="orgbfbb067"><span class="section-number-3">2.4</span> Printers</h3>
<div id="outline-container-orga5fe616" class="outline-3">
<h3 id="orga5fe616"><span class="section-number-3">2.4</span> Printers</h3>
<div class="outline-text-3" id="text-2-4">
<div class="org-src-container">
<pre class="src src-ocaml"><span class="org-tuareg-font-lock-governing">val</span> <span class="org-function-name">pp</span> <span class="org-tuareg-font-lock-operator">:</span> <span class="org-tuareg-font-lock-module">Format.</span>formatter <span class="org-tuareg-font-lock-operator">-&gt;</span> t <span class="org-tuareg-font-lock-operator">-&gt;</span> unit
@ -423,8 +423,8 @@ Frozen_core.num_mos f ;;
</div>
</div>
<div id="outline-container-org71f8a17" class="outline-2">
<h2 id="org71f8a17"><span class="section-number-2">3</span> Orbital localization</h2>
<div id="outline-container-org03637a8" class="outline-2">
<h2 id="org03637a8"><span class="section-number-2">3</span> Orbital localization</h2>
<div class="outline-text-2" id="text-3">
<p>
Molecular orbital localization function.
@ -440,11 +440,11 @@ Edmiston-Rudenberg:
</div>
<div id="outline-container-orgd6cd073" class="outline-3">
<h3 id="orgd6cd073"><span class="section-number-3">3.1</span> Type</h3>
<div id="outline-container-orge723841" class="outline-3">
<h3 id="orge723841"><span class="section-number-3">3.1</span> Type</h3>
<div class="outline-text-3" id="text-3-1">
<div class="org-src-container">
<pre class="src src-ocaml" id="orgbd40984"><span class="org-tuareg-font-lock-governing">open </span><span class="org-tuareg-font-lock-module">Linear_algebra</span>
<pre class="src src-ocaml" id="org1b6afe0"><span class="org-tuareg-font-lock-governing">open </span><span class="org-tuareg-font-lock-module">Linear_algebra</span>
<span class="org-tuareg-font-lock-governing">type</span> <span class="org-type">localization_kind</span> <span class="org-tuareg-font-lock-operator">=</span>
<span class="org-tuareg-font-lock-operator">|</span> <span class="org-tuareg-font-lock-constructor">Edmiston</span>
@ -464,16 +464,16 @@ Edmiston-Rudenberg:
</div>
</div>
<div id="outline-container-org5abbc4c" class="outline-3">
<h3 id="org5abbc4c"><span class="section-number-3">3.2</span> Edmiston-Rudenberg</h3>
<div id="outline-container-org6933918" class="outline-3">
<h3 id="org6933918"><span class="section-number-3">3.2</span> Edmiston-Rudenberg</h3>
</div>
<div id="outline-container-org1360ab2" class="outline-3">
<h3 id="org1360ab2"><span class="section-number-3">3.3</span> Boys</h3>
<div id="outline-container-org44c9408" class="outline-3">
<h3 id="org44c9408"><span class="section-number-3">3.3</span> Boys</h3>
</div>
<div id="outline-container-org44850fb" class="outline-3">
<h3 id="org44850fb"><span class="section-number-3">3.4</span> Access</h3>
<div id="outline-container-org64dada8" class="outline-3">
<h3 id="org64dada8"><span class="section-number-3">3.4</span> Access</h3>
<div class="outline-text-3" id="text-3-4">
<div class="org-src-container">
<pre class="src src-ocaml"><span class="org-tuareg-font-lock-governing">val</span> <span class="org-function-name">kind</span> <span class="org-tuareg-font-lock-operator">:</span> t <span class="org-tuareg-font-lock-operator">-&gt;</span> localization_kind
@ -522,8 +522,8 @@ Edmiston-Rudenberg:
</div>
</div>
<div id="outline-container-orgb1a2fc0" class="outline-3">
<h3 id="orgb1a2fc0"><span class="section-number-3">3.5</span> Printers</h3>
<div id="outline-container-org5b86753" class="outline-3">
<h3 id="org5b86753"><span class="section-number-3">3.5</span> Printers</h3>
<div class="outline-text-3" id="text-3-5">
<div class="org-src-container">
<pre class="src src-ocaml"><span class="org-tuareg-font-lock-governing">val</span> <span class="org-function-name">pp</span> <span class="org-tuareg-font-lock-operator">:</span> <span class="org-tuareg-font-lock-module">Format.</span>formatter <span class="org-tuareg-font-lock-operator">-&gt;</span> t <span class="org-tuareg-font-lock-operator">-&gt;</span> unit
@ -532,14 +532,14 @@ Edmiston-Rudenberg:
</div>
</div>
<div id="outline-container-orgc751f7e" class="outline-3">
<h3 id="orgc751f7e"><span class="section-number-3">3.6</span> Tests</h3>
<div id="outline-container-orgcfa3655" class="outline-3">
<h3 id="orgcfa3655"><span class="section-number-3">3.6</span> Tests</h3>
</div>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Anthony Scemama</p>
<p class="date">Created: 2021-02-03 Wed 09:30</p>
<p class="date">Created: 2021-02-03 Wed 10:47</p>
<p class="validation"><a href="https://validator.w3.org/check?uri=referer">Validate</a></p>
</div>
</body>

View File

@ -3,7 +3,7 @@
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2021-02-03 Wed 09:30 -->
<!-- 2021-02-03 Wed 10:47 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Top-level</title>
@ -250,18 +250,18 @@ org_html_manager.setup(); // activate after the parameters are set
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org30505a3">1. Summmary</a></li>
<li><a href="#orgbfde666">1. Summmary</a></li>
</ul>
</div>
</div>
<div id="outline-container-org30505a3" class="outline-2">
<h2 id="org30505a3"><span class="section-number-2">1</span> Summmary</h2>
<div id="outline-container-orgbfde666" class="outline-2">
<h2 id="orgbfde666"><span class="section-number-2">1</span> Summmary</h2>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Anthony Scemama</p>
<p class="date">Created: 2021-02-03 Wed 09:30</p>
<p class="date">Created: 2021-02-03 Wed 10:47</p>
<p class="validation"><a href="https://validator.w3.org/check?uri=referer">Validate</a></p>
</div>
</body>

View File

@ -2,6 +2,7 @@
(names
ex_integrals
ex_hartree_fock
ex_localization
)
(libraries qcaml))

170
examples/ex_localization.ml Normal file
View File

@ -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 "<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='r' ; long="range" ; opt=Mandatory;
arg=With_arg "<range>";
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

View File

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

View File

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