diff --git a/CI/CI.ml b/CI/CI.ml index 171fa6d..bb0e626 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -676,7 +676,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = Parallel.broadcast (lazy result) in let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in - eigenvectors, eigenvalues + (Conventions.rephase eigenvectors), eigenvalues in let eigensystem_direct () = @@ -702,7 +702,7 @@ let make ?(n_states=1) ?(algo=`Direct) det_space = Parallel.broadcast (lazy result) in let eigenvalues = Vec.map (fun x -> x +. e_shift) eigenvalues in - eigenvectors, eigenvalues + (Conventions.rephase eigenvectors), eigenvalues in match algo with diff --git a/CI/F12CI.ml b/CI/F12CI.ml index afa9506..f5104b8 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -210,6 +210,10 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen let rec iteration ~state psi = + (* + Format.printf "%a@." DeterminantSpace.pp_det_space @@ CI.det_space ci; + Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat psi; + *) let column_idx = iamax (Mat.to_col_vecs psi).(state-1) in let delta = @@ -217,6 +221,10 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen dressing_vector ~frozen_core hf12_integrals psi ci |> Matrix.to_mat in + (* + Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat delta; + *) + Printf.printf "Cmax : %e\n" psi.{column_idx,state}; Printf.printf "Norm : %e\n" (sqrt (gemm ~transa:`T delta delta).{state,state}); @@ -271,10 +279,13 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen Parallel.broadcast (lazy ( - Davidson.make ~threshold:1.e-9 ~guess:psi ~n_states:state diagonal matrix_prod + Davidson.make ~threshold:1.e-10 ~guess:psi ~n_states:state diagonal matrix_prod )) in + let eigenvectors = + Conventions.rephase eigenvectors + in Vec.iter (fun energy -> Printf.printf "%g\t" energy) eigenvalues; diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index 800a58b..eb5274c 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -41,6 +41,7 @@ let two_e_ints t = Lazy.force t.ee_ints let f12_ints t = Lazy.force t.f12_ints let one_e_ints t = Lazy.force t.one_e_ints + let mo_energies t = let m_C = mo_coef t in let f = @@ -138,6 +139,7 @@ let of_mo_basis simulation other = in (* Gram-Schmidt Orthonormalization *) gemm m_X @@ (Util.qr_ortho @@ gemm ~transa:`T m_X result) + |> Conventions.rephase in let mo_occupation = @@ -152,8 +154,6 @@ let of_mo_basis simulation other = - - let pp_mo ?(start=1) ?finish ppf t = let open Lacaml.Io in let rows = Mat.dim1 t.mo_coef diff --git a/MOBasis/MOBasis.mli b/MOBasis/MOBasis.mli index 0c6e055..f6e6091 100644 --- a/MOBasis/MOBasis.mli +++ b/MOBasis/MOBasis.mli @@ -55,7 +55,6 @@ val size : t -> int val mo_energies : t -> Vec.t (** Fock MO energies *) - (** {1 Creators} *) val make : simulation:Simulation.t -> diff --git a/SCF/Guess.ml b/SCF/Guess.ml index 5e96d6c..d0ef751 100644 --- a/SCF/Guess.ml +++ b/SCF/Guess.ml @@ -17,6 +17,7 @@ let hcore_guess ao_basis = and kin_ints = Ao.kin_ints ao_basis |> KinInt.matrix in Mat.add eN_ints kin_ints + |> Conventions.rephase let huckel_guess ao_basis = @@ -46,6 +47,7 @@ let huckel_guess ao_basis = else diag.{i} ) + |> Conventions.rephase let make ?(nocc=0) ~guess ao_basis = @@ -66,6 +68,7 @@ let test_case ao_basis = Lacaml.D.Mat.add (AOBasis.eN_ints ao_basis |> NucInt.matrix) (AOBasis.kin_ints ao_basis |> KinInt.matrix) + |> Conventions.rephase |> Lacaml.D.Mat.to_array in Array.iteri (fun i x -> diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 7ea890f..b0304fe 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -188,8 +188,6 @@ let exchange_energy t = - - let make ?kind ?guess:(guess=`Huckel) @@ -348,6 +346,7 @@ let make (* MOs in AO basis *) let m_C = gemm m_X m_C' + |> Conventions.rephase in (* Hartree-Fock energy *) @@ -542,6 +541,7 @@ let make (* MOs in AO basis *) let m_C = gemm m_X m_C' + |> Conventions.rephase in (* Hartree-Fock energy *) @@ -672,25 +672,26 @@ let pp_summary ppf t = let ha_to_ev = Constants.ha_to_ev in let e = eigenvalues t in - Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth); - Format.fprintf ppf "@["; - print "One-electron energy" (kin_energy t +. eN_energy t) ; - print "Kinetic" (kin_energy t) ; - print "Potential" (eN_energy t) ; - line () ; - print "Two-electron energy" (coulomb_energy t +. exchange_energy t) ; - print "Coulomb" (coulomb_energy t) ; - print "Exchange" (exchange_energy t) ; - line (); - print "HF HOMO" (ha_to_ev *. e.{nocc t}); - print "HF LUMO" (ha_to_ev *. e.{nocc t + 1}); - print "HF LUMO-LUMO" (ha_to_ev *. (e.{nocc t + 1} -. e.{nocc t })); - line (); - print "Electronic energy" (energy t -. nuclear_repulsion t) ; - print "Nuclear repulsion" (nuclear_repulsion t) ; - print "Hartree-Fock energy" (energy t) ; - Format.fprintf ppf "@]"; Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth) + ; Format.fprintf ppf "@[" + ; print "One-electron energy" (kin_energy t +. eN_energy t) + ; print "Kinetic" (kin_energy t) + ; print "Potential" (eN_energy t) + ; line () + ; print "Two-electron energy" (coulomb_energy t +. exchange_energy t) + ; print "Coulomb" (coulomb_energy t) + ; print "Exchange" (exchange_energy t) + ; line () + ; print "HF HOMO" (ha_to_ev *. e.{nocc t}) + ; print "HF LUMO" (ha_to_ev *. e.{nocc t + 1}) + ; print "HF LUMO-LUMO" (ha_to_ev *. (e.{nocc t + 1} -. e.{nocc t })) + ; line () + ; print "Electronic energy" (energy t -. nuclear_repulsion t) + ; print "Nuclear repulsion" (nuclear_repulsion t) + ; print "Hartree-Fock energy" (energy t) + ; Format.fprintf ppf "@]" + ; Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth) + let pp_hf ppf t = diff --git a/Utils/Davidson.ml b/Utils/Davidson.ml index df3634b..78c2caf 100644 --- a/Utils/Davidson.ml +++ b/Utils/Davidson.ml @@ -136,11 +136,14 @@ let make Mat.of_col_vecs_list u_proposed in let maxu = lange u_proposed ~norm:`M in - let thr = maxu *. 0.001 in + let thr = maxu *. 0.00001 in let u_proposed = Mat.map (fun x -> if abs_float x < thr then 0. else x) u_proposed |> Mat.to_col_vecs_list in + (* + Format.printf "%a@." Matrix.pp_matrix @@ Matrix.dense_of_mat m_new_U; + *) if residual_norm > threshold && macro > 0 then diff --git a/Utils/Util.ml b/Utils/Util.ml index 03203d7..1d6faff 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -20,15 +20,15 @@ external popcnt : int64 -> int32 = "popcnt_bytecode" "popcnt" [@@unboxed] [@@noalloc] (** popcnt instruction *) -let popcnt i = popcnt i |> Int32.to_int +let popcnt i = (popcnt [@inlined] ) i |> Int32.to_int -external trailz : int64 -> int32 = "trailz_bytecode" "trailz" +external trailz : int64 -> int32 = "trailz_bytecode" "trailz" "int" [@@unboxed] [@@noalloc] (** ctz instruction *) let trailz i = trailz i |> Int32.to_int -external leadz : int64 -> int32 = "leadz_bytecode" "leadz" +external leadz : int64 -> int32 = "leadz_bytecode" "leadz" "int" [@@unboxed] [@@noalloc] (** bsf instruction *) diff --git a/run_hartree_fock.ml b/run_hartree_fock.ml index c985076..bb67a14 100644 --- a/run_hartree_fock.ml +++ b/run_hartree_fock.ml @@ -71,7 +71,7 @@ let () = Format.fprintf ppf "@[%a@]@." HartreeFock.pp_hf hf (* - let mos = MOBasis.of_hartree_fock hf in + ; let mos = MOBasis.of_hartree_fock hf in Format.fprintf ppf "@[%a@]@." (fun ppf x -> MOBasis.pp_mo ppf x) mos *)