From 57f8c0710c88ae53557ec2b1021ce526cbb0d5c8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 28 May 2018 19:58:40 +0200 Subject: [PATCH] Canonical orthogonalization OK in cart. TODO: spherical. --- Basis/ERI.ml | 1 + Basis/Orthonormalization.ml | 23 +++++++++++++++++++++-- Basis/Overlap.ml | 2 -- HartreeFock/RHF.ml | 1 + Simulation.ml | 15 +++++++++++---- Utils/SphericalToCartesian.ml | 2 +- Utils/SphericalToCartesian.mli | 6 ++++-- Utils/Util.ml | 14 ++++++++++++++ Utils/Util.mli | 4 ++++ 9 files changed, 57 insertions(+), 11 deletions(-) diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 577b85c..2282246 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -77,6 +77,7 @@ let filter_contracted_shell_pairs ?(cutoff=integrals_cutoff) shell_pairs = | Some cspc -> let cls = class_of_contracted_shell_pair_couple cspc in (pair, Zmap.fold (fun key value accu -> max (abs_float value) accu) cls 0. ) + (* TODO \sum_k |coef_k * integral_k| *) | None -> (pair, -1.) ) shell_pairs |> List.filter (fun (_, schwartz_p_max) -> schwartz_p_max >= cutoff) diff --git a/Basis/Orthonormalization.ml b/Basis/Orthonormalization.ml index 9d0df20..6fb3e25 100644 --- a/Basis/Orthonormalization.ml +++ b/Basis/Orthonormalization.ml @@ -7,7 +7,21 @@ type t = | Svd of Mat.t -let make_lowdin ?(thresh=1.e-12) ~overlap = + + +let make_canonical ~cartesian ~thresh ~overlap = + let result = + if cartesian then + canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap) + else + (* TODO *) + canonical_ortho ~thresh ~overlap (Mat.identity @@ Mat.dim1 overlap) + in + Canonical result + + + +let make_lowdin ~thresh ~overlap = let u_vec, u_val = diagonalize_symm overlap in @@ -26,4 +40,9 @@ let make_lowdin ?(thresh=1.e-12) ~overlap = Lowdin result -let make = make_lowdin + +let make ~cartesian ?(thresh=1.e-12) overlap = + (* + make_lowdin ~thresh ~overlap + *) + make_canonical ~cartesian ~thresh ~overlap diff --git a/Basis/Overlap.ml b/Basis/Overlap.ml index 2e98a8d..3ab474b 100644 --- a/Basis/Overlap.ml +++ b/Basis/Overlap.ml @@ -127,7 +127,6 @@ let of_basis basis = - (** Write all overlap integrals to a file *) let to_file ~filename overlap = @@ -146,4 +145,3 @@ let to_file ~filename overlap = - diff --git a/HartreeFock/RHF.ml b/HartreeFock/RHF.ml index 64df9e3..dfa5d67 100644 --- a/HartreeFock/RHF.ml +++ b/HartreeFock/RHF.ml @@ -23,6 +23,7 @@ let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64) | Svd x -> x in + (* Core Hamiltonian *) let m_Hc = let m_T = Lazy.force simulation.kin_ints diff --git a/Simulation.ml b/Simulation.ml index 193c0fa..66a3d72 100644 --- a/Simulation.ml +++ b/Simulation.ml @@ -9,9 +9,15 @@ type t = { kin_ints : KinInt.t lazy_t; ee_ints : ERI.t lazy_t; nuclear_repulsion : float; + cartesian : bool; } -let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis = +let make ?cartesian:(cartesian=true) + ?multiplicity:(multiplicity=1) + ?charge:(charge=0) + ~nuclei + basis + = (* Tune Garbage Collector *) let gc = Gc.get () in @@ -29,20 +35,21 @@ let make ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis = { charge ; basis ; nuclei ; electrons ; overlap ; - overlap_ortho = lazy (Orthonormalization.make (Lazy.force overlap)); + overlap_ortho = lazy (Orthonormalization.make ~cartesian (Lazy.force overlap)); eN_ints = lazy (NucInt.of_basis_nuclei basis nuclei); kin_ints = lazy (KinInt.of_basis basis); ee_ints = lazy (ERI.of_basis basis); nuclear_repulsion = Nuclei.repulsion nuclei; + cartesian ; } -let of_filenames ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis_filename = +let of_filenames ?cartesian:(cartesian=false) ?multiplicity:(multiplicity=1) ?charge:(charge=0) ~nuclei basis_filename = let nuclei = Nuclei.of_filename nuclei in let basis = Basis.of_nuclei_and_basis_filename ~nuclei basis_filename in - make ~charge ~multiplicity ~nuclei basis + make ~cartesian ~charge ~multiplicity ~nuclei basis diff --git a/Utils/SphericalToCartesian.ml b/Utils/SphericalToCartesian.ml index 37c0d0f..c2c9142 100644 --- a/Utils/SphericalToCartesian.ml +++ b/Utils/SphericalToCartesian.ml @@ -647,7 +647,7 @@ let matrix = function | _ -> invalid_arg "Not implemented" -let index nx ny nz = +let index ~nx ~ny ~nz = let l = nx + ny + nz in ((l-nx)*(l-nx+1))/2 + nz + 1 diff --git a/Utils/SphericalToCartesian.mli b/Utils/SphericalToCartesian.mli index 0d0e67d..d8f555f 100644 --- a/Utils/SphericalToCartesian.mli +++ b/Utils/SphericalToCartesian.mli @@ -1,3 +1,5 @@ +(** Conversion from spherical coordinate to cartesian corrdinates. *) + val matrix : AngularMomentum.t -> float array (** Returns a transformation matrix to rotate between the basis of atom-centered spherical coordinates to x,y,z coordinates. @@ -12,12 +14,12 @@ val matrix : AngularMomentum.t -> float array ]} *) -val index : int -> int -> int -> int +val index : nx:int -> ny:int -> nz:int -> int (** Unique index given to a triplet of powers. Used to identify a cartesian shell. Example: {[ let nx, ny, nz = 3, 2, 1 in - SphericalToCartesian.index nx ny nz -> 8 + SphericalToCartesian.index ~nx ~ny ~nz -> 8 ]} *) diff --git a/Utils/Util.ml b/Utils/Util.ml index 6641098..ed58607 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -182,6 +182,20 @@ let xt_o_x ~o ~x = |> gemm ~transa:`T x +let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = + let d, u, _ = gesvd (lacpy overlap) in + let d0 = Vec.sqrt d in + let n = Vec.fold (fun accu x -> if x > thresh then accu + 1 else accu) 0 d in + let d = Vec.map (fun x -> + if x >= thresh then 1. /. x + else 0. ) ~y:d d0 + in + if n < Vec.dim d0 then + Printf.printf "Removed linear dependencies below %f\n" (1. /. d.{n}) + ; + Mat.scal_cols u d ; + gemm c u + (** {2 Printers} *) diff --git a/Utils/Util.mli b/Utils/Util.mli index 4d99554..6a3228d 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -69,6 +69,10 @@ val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat (** Computes X{^T}.O.X *) +val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat +(** Canonical orthogonalization. [overlap] is the overlap matrix, and the last argument + contains the vectors to orthogonalize. *) + (** {2 Printers} *) val pp_float_array_size : Format.formatter -> float array -> unit