diff --git a/.gitignore b/.gitignore index 1a33ea8..81626aa 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ _build .merlin *.install +qcaml.opam diff --git a/ao_basis/test/ao_basis_gaussian.ml b/ao_basis/test/ao_basis_gaussian.ml new file mode 100644 index 0000000..9fa0b9f --- /dev/null +++ b/ao_basis/test/ao_basis_gaussian.ml @@ -0,0 +1,125 @@ +open Alcotest +open Qcaml_common +open Qcaml_gaussian_integrals +open Qcaml_gaussian_basis +open Qcaml_ao_basis +open Qcaml_particles +open Qcaml_operators +open Qcaml_linear_algebra +open Ao_basis_gaussian + +let wd = Qcaml.root ^ Filename.dir_sep ^ "test" + +let make_tests name t = + + let check_matrix title a r = + let ax = + Matrix.as_vec_inplace a + |> Vector.to_bigarray_inplace + in + Matrix.as_vec_inplace r + |> Vector.iteri (fun i x -> + let message = + Printf.sprintf "%s line %d" title i + in + check (float 1.e-10) message ax.{i} x + ) + in + + let check_eri a r = + let f { Eri.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } = + (i_r1, (j_r2, (k_r1, (l_r2, value)))) + in + let a = Eri.to_list a |> List.rev_map f in + let r = Eri.to_list r |> List.rev_map f in + check (list (pair int (pair int (pair int (pair int (float 1.e-10)))))) "ERI" a r + in + + let check_eri_lr a r = + let f { Eri_long_range.i_r1 ; j_r2 ; k_r1 ; l_r2 ; value } = + (i_r1, (j_r2, (k_r1, (l_r2, value)))) + in + let a = Eri_long_range.to_list a |> List.rev_map f in + let r = Eri_long_range.to_list r |> List.rev_map f in + check (list (pair int (pair int (pair int (pair int (float 1.e-10)))))) "ERI_lr" a r + in + + let test_overlap () = + let reference = + Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_overlap.ref") + in + let overlap = + Lazy.force t.overlap |> Overlap.matrix + in + check_matrix "Overlap" overlap reference + in + + let test_eN_ints () = + let reference = + Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_nuc.ref") + in + let eN_ints = + Lazy.force t.eN_ints |> Electron_nucleus.matrix + in + check_matrix "eN_ints" eN_ints reference + in + + let test_kin_ints () = + let reference = + Matrix.sym_matrix_of_file (wd ^ Filename.dir_sep ^ name ^ "_kin.ref") + in + let kin_ints = + Lazy.force t.kin_ints |> Kinetic.matrix + in + check_matrix "kin_ints" kin_ints reference + in + + let test_ee_ints () = + let reference = + Eri.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri.ref") ~sparsity:`Dense ~size:(Basis.size t.basis) + in + let ee_ints = + Lazy.force t.ee_ints + in + check_eri ee_ints reference + ; + in + + let test_ee_lr_ints () = + let reference = + Eri_long_range.of_file (wd ^ Filename.dir_sep ^ name ^ "_eri_lr.ref") ~sparsity:`Dense + ~size:(Basis.size t.basis) + in + let ee_lr_ints = + Lazy.force t.ee_lr_ints + in + check_eri_lr ee_lr_ints reference + in + + [ + "Overlap", `Quick, test_overlap; + "eN_ints", `Quick, test_eN_ints; + "kin_ints", `Quick, test_kin_ints; + "ee_ints", `Quick, test_ee_ints; + "ee_lr_ints", `Quick, test_ee_lr_ints; + ] + +let tests = + List.concat [ + let nuclei = + wd ^ Filename.dir_sep ^ "water.xyz" + |> Nuclei.of_xyz_file + in + let basis = + wd ^ Filename.dir_sep ^ "cc-pvdz" + |> Basis.of_nuclei_and_basis_filename ~nuclei + in + let rs = 0.5 in + let operators = [ Operator.of_range_separation rs ] in + let ao_basis = + Ao_basis_gaussian.make ~basis ~operators nuclei + in + make_tests "water" ao_basis ; + + ] + diff --git a/ao_basis/test/dune b/ao_basis/test/dune new file mode 100644 index 0000000..1af4abe --- /dev/null +++ b/ao_basis/test/dune @@ -0,0 +1,10 @@ +(library + (name test_ao_basis) + (libraries + alcotest + qcaml.common + qcaml.linear_algebra + qcaml.gaussian_integrals + qcaml.ao_basis + ) + (synopsis "Tests for the AO basis")) diff --git a/linear_algebra/lib/conventions.ml b/linear_algebra/lib/conventions.ml new file mode 100644 index 0000000..234394d --- /dev/null +++ b/linear_algebra/lib/conventions.ml @@ -0,0 +1,28 @@ +(** Contains the conventions relative to the program. + + The phase convention is given by: + {% $\sum_i c_i > 0$ %} or {% $\min_i c_i \ge 0$ %} + +*) + +let in_phase vec = + let s = Vector.sum vec in + if s = 0. then + let rec first_non_zero k = + if k > Vector.dim vec then + k-1 + else if Vector.at vec k = 0. then + first_non_zero (k+1) + else k + in + let k = first_non_zero 1 in + Vector.at vec k >= 0. + else + s > 0. + + +let rephase mat = + Matrix.to_col_vecs mat + |> Array.map (fun v -> if in_phase v then v else Vector.neg v) + |> Matrix.of_col_vecs + diff --git a/linear_algebra/lib/conventions.mli b/linear_algebra/lib/conventions.mli new file mode 100644 index 0000000..60da2d3 --- /dev/null +++ b/linear_algebra/lib/conventions.mli @@ -0,0 +1,14 @@ +(** Contains the conventions relative to the program. + + The phase convention is given by: + {% $\sum_i c_i > 0$ %} or {% $\min_i c_i \ge 0$ %} + +*) + + +val in_phase : 'a Vector.t -> bool +(** Checks if one MO respects the phase convention *) + +val rephase : ('a,'b) Matrix.t -> ('a,'b) Matrix.t +(** Apply the phase convention to the MOs. *) +