diff --git a/Basis/ERI.ml b/Basis/ERI.ml index 0f6563a..44337fa 100644 --- a/Basis/ERI.ml +++ b/Basis/ERI.ml @@ -38,6 +38,10 @@ let cutoff = integrals_cutoff *) let zero_m ~maxm ~expo_pq_inv ~norm_pq_sq = + assert (expo_pq_inv <> 0.); + let norm_pq_sq = + if norm_pq_sq > integrals_cutoff then norm_pq_sq else 0. + in let exp_pq = 1. /. expo_pq_inv in let t = norm_pq_sq *. exp_pq in let f = two_over_sq_pi *. (sqrt exp_pq) in @@ -124,7 +128,7 @@ let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls let xl = to_powers powers_l in let key = Zkey.of_powers_twelve xi xj xk xl in let value = Zmap.find cls key in - set_chem data i_c j_c k_c l_c value; + set_chem data i_c j_c k_c l_c value ) (Cs.zkey_array (Csp.shell_b shell_q)) ) (Cs.zkey_array (Csp.shell_a shell_q)) ) (Cs.zkey_array (Csp.shell_b shell_p)) @@ -135,7 +139,7 @@ let store_class ?(cutoff=integrals_cutoff) data contracted_shell_pair_couple cls -let of_basis_serial basis = +let of_basis basis = let n = Bs.size basis and shell = Bs.contracted_shells basis @@ -188,8 +192,11 @@ let of_basis_serial basis = in match cspc with - | Some cspc -> let cls = class_of_contracted_shell_pair_couple cspc in - store_class ~cutoff eri_array cspc cls + | Some cspc -> + let cls = + class_of_contracted_shell_pair_couple cspc + in + store_class ~cutoff eri_array cspc cls | None -> () ) shell_pairs with Exit -> () @@ -199,7 +206,7 @@ let of_basis_serial basis = - +(* let of_basis_parallel basis = let store_class ?(cutoff=integrals_cutoff) push_socket contracted_shell_pair_couple cls = @@ -256,7 +263,7 @@ let of_basis_parallel basis = |> filter_contracted_shell_pairs ~cutoff in - Printf.printf "%d significant shell pairs computed in %f seconds\n" + Printf.printf "%d significant shell pairs computed in %f seconds\n%!" (List.length shell_pairs) (Unix.gettimeofday () -. t0); @@ -265,20 +272,29 @@ let of_basis_parallel basis = - let zmq_addr = Printf.sprintf "ipc://%d" (Unix.getpid ()) in - let () = + let zmq_port = 12345 in + begin + match Unix.fork () with + | 0 -> Printf.printf "pouet\n%!" + | pid -> Printf.printf "coucou\n%!" + end; + begin + match Unix.fork () with | 0 -> begin + let zmq_addr = Printf.sprintf "tcp://localhost:%d" zmq_port in let zmq = ref None in +Printf.printf "PID %d OK\n%!" 0; Parmap.pariter ~chunksize:1 - ~init:(fun _ -> + ~init:(fun rank -> let zmq_context = Zmq.Context.create () in let push_socket = Zmq.Socket.create zmq_context Zmq.Socket.push in + Printf.printf "Init %d OK\n%!" rank; Zmq.Socket.connect push_socket zmq_addr; zmq := Some (zmq_context, push_socket) ) @@ -341,6 +357,8 @@ let of_basis_parallel basis = end | pid -> begin +Printf.printf "PID %d OK\n%!" pid; + let zmq_addr = Printf.sprintf "tcp://*:%d" zmq_port in let zmq_context = Zmq.Context.create () in @@ -364,10 +382,8 @@ let of_basis_parallel basis = Zmq.Context.terminate zmq_context; ignore (Unix.wait ()) end - in + end; Printf.printf "Computed ERIs in %f seconds\n%!" (Unix.gettimeofday () -. t0); eri_array +*) - - -let of_basis = of_basis_parallel diff --git a/Basis/TwoElectronRRVectorized.ml b/Basis/TwoElectronRRVectorized.ml index d8abcf1..af043a0 100644 --- a/Basis/TwoElectronRRVectorized.ml +++ b/Basis/TwoElectronRRVectorized.ml @@ -631,8 +631,9 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q in zero_m_array.(0) with NullQuartet -> 0. - ) in - Mat.gemm_trace zm_array coef + ) + in + Mat.gemm_trace zm_array coef with (Invalid_argument _) -> 0. end | _ -> diff --git a/INSTALL.md b/INSTALL.md index 351ef28..b8afc44 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -21,15 +21,6 @@ export LACAML_LIBS="-L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_rt -lpthr opam install lacaml ``` -# Parmap - -Multicore library. - -```bash -opam install parmap -``` - - # odoc-ltxhtml This plugin allows to embed equations in the documentation generated by Ocamldoc. diff --git a/SCF/Fock.ml b/SCF/Fock.ml index 0f51f20..35ca589 100644 --- a/SCF/Fock.ml +++ b/SCF/Fock.ml @@ -1,6 +1,7 @@ open Lacaml.D open Simulation open Constants +open Util type t = @@ -61,6 +62,15 @@ let make ~density ao_basis = let m_J = Mat.of_array m_J and m_K = Mat.of_array m_K in - { fock = Mat.add m_Hc @@ Mat.add m_J m_K ; + { fock = Mat.add m_Hc (Mat.add m_J m_K) ; core = m_Hc ; coulomb = m_J ; exchange = m_K } + +let pp_fock ppf a = + Format.fprintf ppf "@[<2>"; + Format.fprintf ppf "@[ Fock matrix:@[<2>@[%a@]@.]@]" pp_matrix a.fock; + Format.fprintf ppf "@[ Core Hamiltonian:@[<2>@[%a@]@.]@]" pp_matrix a.core; + Format.fprintf ppf "@[ Coulomb matrix:@[<2>@[%a@]@.]@]" pp_matrix a.coulomb; + Format.fprintf ppf "@[ Exchange matrix:@[<2>@[%a@]@.]@]" pp_matrix a.exchange; + Format.fprintf ppf "@]" + diff --git a/SCF/Guess.ml b/SCF/Guess.ml index e347ccf..a3cb000 100644 --- a/SCF/Guess.ml +++ b/SCF/Guess.ml @@ -33,7 +33,9 @@ let huckel_guess ao_basis = function | 0 -> invalid_arg "Huckel guess needs a non-zero number of occupied MOs." | nocc -> - let m_F = (Fock.make ~density:(gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X) ao_basis).Fock.fock in + let density = gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X in + let fock = Fock.make ~density ao_basis in + let m_F = fock.Fock.fock in for j=1 to ao_num do for i=1 to ao_num do if (i <> j) then diff --git a/Utils/DIIS.mli b/Utils/DIIS.mli index bdf057c..cfef058 100644 --- a/Utils/DIIS.mli +++ b/Utils/DIIS.mli @@ -19,7 +19,7 @@ with a Langrange multiplier {% $\lambda$ %}. {% \begin{align*} \mathcal{L} & = ||\sum_i c_i \mathbf{e}_i||^2 - \lambda \left(\sum_i c_i - 1\right) \\ - & = \sum_{ij} c_i c_j B_{ij} - \lambda \left(\sum_i c_i - 1\right) + & = \sum_{ij} c_i c_j B_{ij} - \lambda \left(\sum_{i=1}^m c_i - 1\right) \end{align*} %} with diff --git a/Utils/FourIdxStorage.ml b/Utils/FourIdxStorage.ml index 2c8b57e..7a4f4e3 100644 --- a/Utils/FourIdxStorage.ml +++ b/Utils/FourIdxStorage.ml @@ -104,11 +104,21 @@ let increment_four_index ~r1 ~r2 ~value t = in Hashtbl.replace a key (old_value +. value) -let get ~r1 ~r2 = - get_four_index ~r1 ~r2 +let get ~r1 ~r2 a = + get_four_index ~r1 ~r2 a -let set ~r1 ~r2 = - set_four_index ~r1 ~r2 +let set ~r1 ~r2 ~value = + match classify_float value with + | FP_normal -> set_four_index ~r1 ~r2 ~value + | FP_zero + | FP_subnormal -> fun _ -> () + | FP_infinite + | FP_nan -> + let msg = + Printf.sprintf "FourIdxStorage.ml : set : r1 = (%d,%d) ; r2 = (%d,%d)" + r1.first r1.second r2.first r2.second + in + raise (Invalid_argument msg) let increment ~r1 ~r2 = increment_four_index ~r1 ~r2 diff --git a/Utils/Util.ml b/Utils/Util.ml index be4f8fd..bcbf32f 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -133,6 +133,9 @@ let boys_function ~maxm t = let result = Array.init (maxm+1) (fun m -> 1. /. float_of_int (2*m+1)) in + (* + assert (abs_float t > 1.e-10); + *) if t <> 0. then begin let fmax = @@ -140,7 +143,12 @@ let boys_function ~maxm t = let n = float_of_int maxm in let dm = 0.5 +. n in let f = (pow t_inv (maxm+maxm+1) ) in - (incomplete_gamma dm t) *. 0.5 *. f + match classify_float f with + | FP_zero + | FP_subnormal + | FP_normal -> + (incomplete_gamma dm t) *. 0.5 *. f + | _ -> invalid_arg "zero_m overflow" in let emt = exp (-. t) in result.(maxm) <- fmax; @@ -199,34 +207,6 @@ let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = gemm c u -let string_of_matrix m = - let open Lacaml.Io in - let rows = Mat.dim1 m - and cols = Mat.dim2 m - in - let rec aux accu first last = - - if (first > last) then String.concat "\n" (List.rev accu) - else - let nw = - Format.asprintf "\n\n %a\n" (Lacaml.Io.pp_lfmat - ~row_labels: - (Array.init rows (fun i -> Printf.sprintf "%d " (i + 1))) - ~col_labels: - (Array.init (min 5 (cols-first+1)) (fun i -> Printf.sprintf "-- %d --" (i + first) )) - ~print_right:false - ~print_foot:false - () ) (lacpy ~ac:first ~n:(min 5 (cols-first+1)) m) - in - aux (nw :: accu) (first+5) last - in - aux [] 1 cols - - -let debug_matrix name a = - Printf.printf "%s =\n%s\n" name (string_of_matrix a) - - (** {2 Printers} *) @@ -250,3 +230,33 @@ let pp_float_2darray_size ppf a = Array.iter (fun f -> Format.fprintf ppf "@[%a@]@ " pp_float_array_size f) a; Format.fprintf ppf "]@]@]" +let pp_matrix ppf m = + let open Lacaml.Io in + let rows = Mat.dim1 m + and cols = Mat.dim2 m + in + let rec aux first last = + if (first <= last) then begin + Format.fprintf ppf "@[\n\n %a@]@ " (Lacaml.Io.pp_lfmat + ~row_labels: + (Array.init rows (fun i -> Printf.sprintf "%d " (i + 1))) + ~col_labels: + (Array.init (min 5 (cols-first+1)) (fun i -> Printf.sprintf "-- %d --" (i + first) )) + ~print_right:false + ~print_foot:false + () ) (lacpy ~ac:first ~n:(min 5 (cols-first+1)) m); + aux (first+5) last + end + in + aux 1 cols + + + + +let string_of_matrix m = + Format.asprintf "%a" pp_matrix m + +let debug_matrix name a = + Format.printf "@[%s =\n@[%a@]@]@ " name pp_matrix a + + diff --git a/Utils/Util.mli b/Utils/Util.mli index b94d358..21c3cfd 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -122,3 +122,5 @@ val pp_float_2darray : Format.formatter -> float array array -> unit ]} *) +val pp_matrix : Format.formatter -> Lacaml.D.mat -> unit + diff --git a/_tags b/_tags index 9568e5b..dca8492 100644 --- a/_tags +++ b/_tags @@ -1,4 +1,4 @@ -true: package(str,unix,bigarray,lacaml,parmap,zmq) +true: package(str,unix,bigarray,lacaml) <*.byte> : linkdep(Utils/math_functions.o), custom <*.native>: linkdep(Utils/math_functions.o) : not_hygienic