From bb8654b17c4ae375125cfc5505604d8dd1c67070 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 20 Nov 2019 16:42:28 +0100 Subject: [PATCH] Debugging --- Basis/GeneralBasis.ml | 2 +- CI/CI.ml | 4 +- CI/F12CI.ml | 109 ++++++++++++++++++++++++++++++++++++++-- MOBasis/MOBasis.ml | 2 + SCF/Guess.ml | 2 - SCF/HartreeFock.ml | 4 ++ Utils/FourIdxStorage.ml | 8 +-- Utils/Util.ml | 9 ++++ Utils/Util.mli | 3 ++ 9 files changed, 129 insertions(+), 14 deletions(-) diff --git a/Basis/GeneralBasis.ml b/Basis/GeneralBasis.ml index c5e186f..a6f3990 100644 --- a/Basis/GeneralBasis.ml +++ b/Basis/GeneralBasis.ml @@ -83,7 +83,7 @@ let combine basis_list = List.concat basis_list |> List.iter (fun (k,v) -> let l = Hashtbl.find_all h k in - Hashtbl.replace h k (Array.concat (v::l) ) + Hashtbl.replace h k (Array.concat (l@[v]) ) ) ; Hashtbl.fold (fun k v accu -> (k, v)::accu) h [] diff --git a/CI/CI.ml b/CI/CI.ml index bb0e626..ab1ef95 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 - (Conventions.rephase eigenvectors), eigenvalues + (Conventions.rephase @@ Util.remove_epsilons 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 - (Conventions.rephase eigenvectors), eigenvalues + (Conventions.rephase @@ Util.remove_epsilons eigenvectors), eigenvalues in match algo with diff --git a/CI/F12CI.ml b/CI/F12CI.ml index f5104b8..5ee8586 100644 --- a/CI/F12CI.ml +++ b/CI/F12CI.ml @@ -70,6 +70,7 @@ let single_matrices hf12_integrals density = done; *) +(*--- let hf_ij_non_zero mo_basis hf12_integrals deg_a deg_b ki kj = let integrals = [ @@ -185,6 +186,104 @@ let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) +--- *) + + +let hf_ij_non_zero hf12_integrals deg_a deg_b ki kj = + let integrals = [ + + let { HF12. + simulation ; aux_basis ; + hf12 ; hf12_anti ; + hf12_single ; hf12_single_anti } = hf12_integrals + in + + let kia = De.alfa ki and kib = De.beta ki in + let kja = De.alfa kj and kjb = De.beta kj in + + let mo_a = + Bitstring.logand (Sp.bitstring kia) (Sp.bitstring kja) + |> Bitstring.to_list + and mo_b = + Bitstring.logand (Sp.bitstring kib) (Sp.bitstring kjb) + |> Bitstring.to_list + in + + let one_e _ _ _ = 0. in + + let two_e i j k l s s' = + if s = s' then + hf12_anti.{i,j,k,l} -. ( + (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,i,j,k,l}) 0. mo_a) +. + (List.fold_left (fun accu m -> accu +. hf12_single_anti.{m,j,i,l,k}) 0. mo_b) + ) + else + hf12.{i,j,k,l} -. ( + (List.fold_left (fun accu m -> accu +. hf12_single.{m,i,j,k,l}) 0. mo_a) +. + (List.fold_left (fun accu m -> accu +. hf12_single.{m,j,i,l,k}) 0. mo_b) + ) + in + + let h = MOBasis.ee_ints aux_basis in + let two_e_h i j k l s s' = + if s' <> s then + ERI.get_phys h l k j i + else + (ERI.get_phys h l k j i) -. (ERI.get_phys h k l j i) + in + + let f = MOBasis.f12_ints aux_basis in + let two_e_f i j k l s s' = + if s' <> s then + F12.get_phys f i j k l + else + (F12.get_phys f i j k l) -. (F12.get_phys f i j l k) + in + + let mo_of_s = function + | Spin.Alfa -> mo_a + | Spin.Beta -> mo_b + in + + let three_e i j k l m n s s' s'' = + List.fold_left (fun accu a -> accu +. two_e_h i j l a s s' *. two_e_f a k m n s' s'') 0. (mo_of_s s' ) + -. List.fold_left (fun accu a -> accu +. two_e_h j i m a s' s *. two_e_f a k l n s s'') 0. (mo_of_s s ) + +. List.fold_left (fun accu a -> accu +. two_e_h j k m a s' s'' *. two_e_f a i n l s'' s ) 0. (mo_of_s s'') + -. List.fold_left (fun accu a -> accu +. two_e_h k j n a s'' s' *. two_e_f a i m l s' s ) 0. (mo_of_s s' ) + +. List.fold_left (fun accu a -> accu +. two_e_h k i n a s'' s *. two_e_f a j l m s s' ) 0. (mo_of_s s ) + -. List.fold_left (fun accu a -> accu +. two_e_h i k l a s s'' *. two_e_f a j n m s'' s' ) 0. (mo_of_s s'') + in + + (one_e, two_e, Some three_e) + ] + in + CIMatrixElement.non_zero integrals deg_a deg_b ki kj + |> List.hd + + +let dressing_vector ~frozen_core hf12_integrals f12_amplitudes ci = + + if Parallel.master then + Printf.printf "Building matrix\n%!"; + let det_space = + ci.CI.det_space + in + + let m_HF = + + let f = + match Ds.determinants det_space with + | Ds.Arbitrary _ -> CI.create_matrix_arbitrary + | Ds.Spin _ -> CI.create_matrix_spin_computed ~nmax:3 + in + f (fun deg_a deg_b ki kj -> + hf_ij_non_zero hf12_integrals deg_a deg_b ki kj + ) det_space + + in + + Matrix.mm (Lazy.force m_HF) (Matrix.dense_of_mat f12_amplitudes) + let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename ?(state=1) () = @@ -210,10 +309,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; - *) +(* +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 = @@ -284,7 +383,7 @@ let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filen in let eigenvectors = - Conventions.rephase eigenvectors + Conventions.rephase @@ Util.remove_epsilons eigenvectors in diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index eb5274c..42431aa 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -49,6 +49,7 @@ let mo_energies t = let m_P = x_o_xt m_N m_C in match t.mo_type with | RHF -> Fock.make_rhf ~density:m_P (ao_basis t) + | Projected | ROHF -> (Mat.scal 0.5 m_P; Fock.make_uhf ~density_same:m_P ~density_other:m_P (ao_basis t)) | _ -> failwith "Not implemented" @@ -139,6 +140,7 @@ let of_mo_basis simulation other = in (* Gram-Schmidt Orthonormalization *) gemm m_X @@ (Util.qr_ortho @@ gemm ~transa:`T m_X result) + |> Util.remove_epsilons |> Conventions.rephase in diff --git a/SCF/Guess.ml b/SCF/Guess.ml index d0ef751..80d332f 100644 --- a/SCF/Guess.ml +++ b/SCF/Guess.ml @@ -17,7 +17,6 @@ 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 = @@ -47,7 +46,6 @@ let huckel_guess ao_basis = else diag.{i} ) - |> Conventions.rephase let make ?(nocc=0) ~guess ao_basis = diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index b0304fe..7d8c126 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -346,6 +346,7 @@ let make (* MOs in AO basis *) let m_C = gemm m_X m_C' + |> Util.remove_epsilons |> Conventions.rephase in @@ -541,6 +542,7 @@ let make (* MOs in AO basis *) let m_C = gemm m_X m_C' + |> Util.remove_epsilons |> Conventions.rephase in @@ -691,6 +693,8 @@ let pp_summary ppf t = ; print "Hartree-Fock energy" (energy t) ; Format.fprintf ppf "@]" ; Format.fprintf ppf "@[%s@]@;" (Printing.line ~c:'=' linewidth) + ; Util.debug_matrix "MOs" (eigenvectors t) + diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index d224c98..d2188a8 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -344,7 +344,7 @@ let to_stream d = (** Write all integrals to a file with the convention *) -let to_file ?(cutoff=Constants.epsilon) ~filename data = +let to_file ?(cutoff=Constants.integrals_cutoff) ~filename data = let oc = open_out filename in to_stream data |> Stream.iter (fun {i_r1 ; j_r2 ; k_r1 ; l_r2 ; value} -> @@ -496,7 +496,7 @@ let four_index_transform_dense_sparse ds coef source = List.iter (fun beta -> List.iter (fun alpha -> let x = u.{alpha,beta,gamma} in - if x <> 0. then + if abs_float x > Constants.integrals_cutoff then result := (alpha, beta, gamma, delta, x) :: !result; ) (list_range 1 beta) ) range_mo @@ -541,7 +541,7 @@ let svd_of_dense t = scal (1. /. d.{1}) w; let rec aux i = if i = Vec.dim w then i else - if w.{i} < 1.e-15 then i else + if w.{i} < epsilon_float then i else aux (i+1) in aux 1 in @@ -668,7 +668,7 @@ let four_index_transform_svd coef source = for c=1 to mo_num do for d=c to mo_num do let x = result.{a,b,c,d} in - if abs_float x > epsilon then + if abs_float x > Constants.integrals_cutoff then set_chem destination a b c d x; done done diff --git a/Utils/Util.ml b/Utils/Util.ml index 1d6faff..3b33011 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -302,6 +302,15 @@ let x_o_xt ~o ~x = gemm o x ~transb:`T |> gemm x +let remove_epsilons m = + let vecs = + Mat.to_col_vecs m + in + Array.map (fun v -> + let m = abs_float (2. *. amax v) in + Vec.map (fun x -> if abs_float x < m *. epsilon_float then 0. else x) v + ) vecs + |> Mat.of_col_vecs let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = let d, u, _ = gesvd (lacpy overlap) in diff --git a/Utils/Util.mli b/Utils/Util.mli index 7f083ad..8aedfa5 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -130,6 +130,9 @@ val xt_o_x : o:Mat.t -> x:Mat.t -> Mat.t val x_o_xt : o:Mat.t -> x:Mat.t -> Mat.t (** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *) +val remove_epsilons : Mat.t -> Mat.t +(** Remove values below machine precision *) + val qr_ortho : Mat.t -> Mat.t (** QR orthogonalization of the input matrix *)