diff --git a/Basis/AtomicShell.ml b/Basis/AtomicShell.ml index 438caa8..a0d9c7c 100644 --- a/Basis/AtomicShell.ml +++ b/Basis/AtomicShell.ml @@ -4,12 +4,12 @@ open Constants type t = { expo : float array array; coef : float array array; - center : Coordinate.t; - ang_mom : AngularMomentum.t; norm_coef : float array array; norm_coef_scale : float array; - index : int; contr : ContractedShell.t array; + index : int; + center : Coordinate.t; + ang_mom : AngularMomentum.t; } module Am = AngularMomentum diff --git a/Basis/AtomicShellPair.ml b/Basis/AtomicShellPair.ml index 63db58b..db7b242 100644 --- a/Basis/AtomicShellPair.ml +++ b/Basis/AtomicShellPair.ml @@ -5,9 +5,9 @@ exception Null_contribution type t = { + contracted_shell_pairs : ContractedShellPair.t list; atomic_shell_a : AtomicShell.t; atomic_shell_b : AtomicShell.t; - contracted_shell_pairs : ContractedShellPair.t list; } diff --git a/Basis/AtomicShellPairCouple.ml b/Basis/AtomicShellPairCouple.ml index e43078a..c993a1a 100644 --- a/Basis/AtomicShellPairCouple.ml +++ b/Basis/AtomicShellPairCouple.ml @@ -2,14 +2,14 @@ open Util type t = { + contracted_shell_pair_couples : ContractedShellPairCouple.t list ; atomic_shell_pair_p: AtomicShellPair.t ; atomic_shell_pair_q: AtomicShellPair.t ; atomic_shell_a : AtomicShell.t ; atomic_shell_b : AtomicShell.t ; atomic_shell_c : AtomicShell.t ; atomic_shell_d : AtomicShell.t ; - ang_mom : AngularMomentum.t ; - contracted_shell_pair_couples : ContractedShellPairCouple.t list ; + ang_mom : AngularMomentum.t ; } module Am = AngularMomentum diff --git a/Basis/ContractedShell.ml b/Basis/ContractedShell.ml index f8de68d..96725a6 100644 --- a/Basis/ContractedShell.ml +++ b/Basis/ContractedShell.ml @@ -4,12 +4,12 @@ open Constants type t = { expo : float array; coef : float array; - center : Coordinate.t; - ang_mom : AngularMomentum.t; norm_coef : float array; norm_coef_scale : float array; - index : int; prim : PrimitiveShell.t array; + center : Coordinate.t; + ang_mom : AngularMomentum.t; + index : int; } module Am = AngularMomentum diff --git a/Basis/ContractedShellPair.ml b/Basis/ContractedShellPair.ml index 16392da..85c1f82 100644 --- a/Basis/ContractedShellPair.ml +++ b/Basis/ContractedShellPair.ml @@ -5,9 +5,9 @@ exception Null_contribution type t = { + coefs_and_shell_pairs : (float * PrimitiveShellPair.t) list; shell_a : ContractedShell.t; shell_b : ContractedShell.t; - coefs_and_shell_pairs : (float * PrimitiveShellPair.t) list; } diff --git a/Basis/ContractedShellPairCouple.ml b/Basis/ContractedShellPairCouple.ml index 2b47302..f371e5b 100644 --- a/Basis/ContractedShellPairCouple.ml +++ b/Basis/ContractedShellPairCouple.ml @@ -2,6 +2,7 @@ open Util type t = { + coefs_and_shell_pair_couples : (float * PrimitiveShellPairCouple.t) list ; shell_pair_p: ContractedShellPair.t ; shell_pair_q: ContractedShellPair.t ; shell_a : ContractedShell.t ; @@ -9,7 +10,6 @@ type t = shell_c : ContractedShell.t ; shell_d : ContractedShell.t ; ang_mom : AngularMomentum.t ; - coefs_and_shell_pair_couples : (float * PrimitiveShellPairCouple.t) list ; } module Am = AngularMomentum diff --git a/Basis/PrimitiveShell.ml b/Basis/PrimitiveShell.ml index d709149..b560679 100644 --- a/Basis/PrimitiveShell.ml +++ b/Basis/PrimitiveShell.ml @@ -3,9 +3,9 @@ open Constants open Coordinate type t = { + norm_scales : float array lazy_t; exponent : float; normalization : float; - norm_scales : float array lazy_t; center : Coordinate.t; ang_mom : AngularMomentum.t; } diff --git a/Basis/PrimitiveShellPair.ml b/Basis/PrimitiveShellPair.ml index 3c07783..1084705 100644 --- a/Basis/PrimitiveShellPair.ml +++ b/Basis/PrimitiveShellPair.ml @@ -3,15 +3,15 @@ open Constants type t = { + norm_scales : float array lazy_t; exponent : float; (* {% $\alpha + \beta$ %} *) exponent_inv : float; (* {% $1/(\alpha + \beta)$ %} *) + a_minus_b_sq : float; (* {% $|A-B|^2$ %} *) + normalization : float; (* [norm_coef_a * norm_coef_b * g], with + {% $g = (\pi/(\alpha+\beta))^(3/2) \exp (-|A-B|^2 \alpha\beta/(\alpha+\beta))$ %} *) center : Coordinate.t; (* {% $P = (\alpha A + \beta B)/(\alpha+\beta)$ %} *) center_minus_a : Coordinate.t; (* {% $P - A$ %} *) a_minus_b : Coordinate.t; (* {% $A - B$ %} *) - a_minus_b_sq : float; (* {% $|A-B|^2$ %} *) - norm_scales : float array lazy_t; - normalization : float; (* [norm_coef_a * norm_coef_b * g], with - {% $g = (\pi/(\alpha+\beta))^(3/2) \exp (-|A-B|^2 \alpha\beta/(\alpha+\beta))$ %} *) ang_mom : AngularMomentum.t; shell_a : PrimitiveShell.t; shell_b : PrimitiveShell.t; diff --git a/Basis/PrimitiveShellPairCouple.ml b/Basis/PrimitiveShellPairCouple.ml index cf07d73..fc8c21e 100644 --- a/Basis/PrimitiveShellPairCouple.ml +++ b/Basis/PrimitiveShellPairCouple.ml @@ -2,6 +2,7 @@ exception NullContribution type t = { + zkey_array : Zkey.t array lazy_t; shell_pair_p: PrimitiveShellPair.t ; shell_pair_q: PrimitiveShellPair.t ; shell_a : PrimitiveShell.t ; @@ -9,7 +10,6 @@ type t = shell_c : PrimitiveShell.t ; shell_d : PrimitiveShell.t ; ang_mom : AngularMomentum.t ; - zkey_array : Zkey.t array lazy_t; } module Am = AngularMomentum diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index b976945..2b9ed04 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -157,52 +157,25 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = let of_rhf hf = - let mo_num = Vec.dim hf.HF.eigenvalues in let mo_coef = hf.HF.eigenvectors in let simulation = hf.HF.simulation in let mo_type = RHF in - let nocc = hf.HF.nocc in - let mo_occupation = - Array.init mo_num (fun i -> - if i < nocc then 2. else 0.) - |> Vec.of_array - in - let result = - make ~simulation ~mo_type ~mo_occupation ~mo_coef () - in + let mo_occupation = hf.HF.occupation in + make ~simulation ~mo_type ~mo_occupation ~mo_coef () -(* - let () = - let e = ref 0. in - let t = KinInt.matrix (Lazy.force result.kin_ints) in - let v = NucInt.matrix (Lazy.force result.eN_ints) in - let g = Lazy.force result.ee_ints in - for i = 1 to 5 do - e := !e +. 2. *. (t.{i,i} +. v.{i,i}) - done; - Printf.printf "Energy one-e = %20.15f\n" !e; - let e2 = ref 0. in - for i = 1 to 5 do - for j = i+1 to 5 do - e2 := !e2 +. 2. *. (ERI.get_phys g i j i j -. - ERI.get_phys g i j j i) - done; - done; - for i = 1 to 5 do - for j = 1 to 5 do - e2 := !e2 +. ERI.get_phys g i j i j - done; - done; - Printf.printf "Energy two-e = %20.15f\n" !e2; - Printf.printf "Energy = %20.15f\n" (Si.nuclear_repulsion simulation +. !e +. !e2) - in - *) - result + +let of_rohf hf = + let mo_coef = hf.HF.eigenvectors in + let simulation = hf.HF.simulation in + let mo_type = ROHF in + let mo_occupation = hf.HF.occupation in + make ~simulation ~mo_type ~mo_occupation ~mo_coef () let of_hartree_fock = function - | HF.RHF hf -> of_rhf hf - | _ -> assert false + | HF.RHF hf -> of_rhf hf + | HF.ROHF hf -> of_rohf hf + | HF.UHF _ -> assert false diff --git a/SCF/HartreeFock_type.ml b/SCF/HartreeFock_type.ml index ca19393..47a7802 100644 --- a/SCF/HartreeFock_type.ml +++ b/SCF/HartreeFock_type.ml @@ -5,16 +5,17 @@ type s = { simulation : Simulation.t; guess : Guess.t; - eigenvectors : Lacaml.D.Mat.t ; - eigenvalues : Lacaml.D.Vec.t ; - nocc : int; + eigenvectors : Mat.t ; + eigenvalues : Vec.t ; + occupation : Vec.t ; + iterations : (float * float * float) array; energy : float ; nuclear_repulsion : float ; kin_energy : float ; eN_energy : float ; coulomb_energy : float ; exchange_energy : float ; - iterations : (float * float * float) array; + nocc : int; } diff --git a/SCF/HartreeFock_type.mli b/SCF/HartreeFock_type.mli index a79d046..285da5f 100644 --- a/SCF/HartreeFock_type.mli +++ b/SCF/HartreeFock_type.mli @@ -1,19 +1,22 @@ +open Lacaml.D + (** Data structure representing the output of a Hartree-Fock caculation *) type s = { simulation : Simulation.t; (** Simulation which was used for HF calculation *) guess : Guess.t; (** Initial guess *) - eigenvectors : Lacaml.D.Mat.t ; (** Final eigenvectors *) - eigenvalues : Lacaml.D.Vec.t ; (** Final eigenvalues *) - nocc : int ; (** Number of occupied MOs *) + eigenvectors : Mat.t ; (** Final eigenvectors *) + eigenvalues : Vec.t ; (** Final eigenvalues *) + occupation : Vec.t ; (** Diagonal of the density matrix *) + iterations : (float * float * float) array; energy : float ; (** Final energy *) nuclear_repulsion : float ; (** Nucleus-Nucleus potential energy *) kin_energy : float ; (** Kinetic energy *) eN_energy : float ; (** Electron-nucleus potential energy *) coulomb_energy : float ; (** Electron-Electron potential energy *) exchange_energy : float ; (** Exchange energy *) - iterations : (float * float * float) array; + nocc : int ; (** Number of occupied MOs *) (** Energy, convergence and HOMO-LUMO gap of all iterations *) } diff --git a/SCF/RHF.ml b/SCF/RHF.ml index 3545a9a..b4fbbe4 100644 --- a/SCF/RHF.ml +++ b/SCF/RHF.ml @@ -174,7 +174,8 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= eN_energy = Mat.gemm_trace m_P m_V; coulomb_energy = 0.5 *. Mat.gemm_trace m_P m_J; exchange_energy = 0.5 *. Mat.gemm_trace m_P m_K; - }) + occupation = Mat.copy_diag m_P; + }) in diff --git a/SCF/ROHF.ml b/SCF/ROHF.ml index c4aafdd..33c44d2 100644 --- a/SCF/ROHF.ml +++ b/SCF/ROHF.ml @@ -265,6 +265,7 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= 0.5 *. (Mat.gemm_trace m_P_b m_J_b); exchange_energy = 0.5 *. (Mat.gemm_trace m_P_a m_K_a) +. 0.5 *. (Mat.gemm_trace m_P_b m_K_b); + occupation = Mat.copy_diag m_P; }) in