From a03e31e293ef9eecb319fa1df503b3f17b61b9e7 Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:02:50 -0400 Subject: [PATCH 01/19] solve conflict --- src/ace_binding.cpp | 26 ++++++++++++++++++- src/flare_pp/descriptors/b1.cpp | 21 ++++++++++++++-- src/flare_pp/descriptors/b1.h | 6 +++++ src/flare_pp/descriptors/b3.cpp | 44 +++++++++++++++++++++++++++++++++ src/flare_pp/descriptors/b3.h | 6 +++++ 5 files changed, 100 insertions(+), 3 deletions(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index edd2255a0..39ad08c86 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -110,7 +110,17 @@ PYBIND11_MODULE(_C_flare, m) { py::class_(m, "B1") .def(py::init &, const std::vector &, - const std::vector &>()); + const std::vector &>()) + .def(py::init &, const std::vector &, + const std::vector &, + const Eigen::MatrixXd &>()) + .def_readonly("radial_basis", &B1::radial_basis) + .def_readonly("cutoff_function", &B1::cutoff_function) + .def_readonly("radial_hyps", &B1::radial_hyps) + .def_readonly("cutoff_hyps", &B1::cutoff_hyps) + .def_readonly("cutoffs", &B1::cutoffs) + .def_readonly("descriptor_settings", &B1::descriptor_settings); py::class_(m, "B2") .def(py::init(m, "B3") .def(py::init &, const std::vector &, +<<<<<<< HEAD const std::vector &>()); +======= + const std::vector &>()) + .def(py::init &, const std::vector &, + const std::vector &, + const Eigen::MatrixXd &>()) + .def_readonly("radial_basis", &B3::radial_basis) + .def_readonly("cutoff_function", &B3::cutoff_function) + .def_readonly("radial_hyps", &B3::radial_hyps) + .def_readonly("cutoff_hyps", &B3::cutoff_hyps) + .def_readonly("cutoffs", &B3::cutoffs) + .def_readonly("descriptor_settings", &B3::descriptor_settings); +>>>>>>> f251f1db (add multiple cutoffs input to b1 and b3) // Kernel functions py::class_(m, "Kernel"); diff --git a/src/flare_pp/descriptors/b1.cpp b/src/flare_pp/descriptors/b1.cpp index 7568d5ed4..d3ee61ede 100644 --- a/src/flare_pp/descriptors/b1.cpp +++ b/src/flare_pp/descriptors/b1.cpp @@ -31,9 +31,26 @@ B1 ::B1(const std::string &radial_basis, const std::string &cutoff_function, cutoffs = Eigen::MatrixXd::Constant(n_species, n_species, cutoff_val); } -void B1 ::write_to_file(std::ofstream &coeff_file, int coeff_size) { - coeff_file << "B1" << "\n"; +B1 ::B1(const std::string &radial_basis, const std::string &cutoff_function, + const std::vector &radial_hyps, + const std::vector &cutoff_hyps, + const std::vector &descriptor_settings, + const Eigen::MatrixXd &cutoffs) { + + this->radial_basis = radial_basis; + this->cutoff_function = cutoff_function; + this->radial_hyps = radial_hyps; + this->cutoff_hyps = cutoff_hyps; + this->descriptor_settings = descriptor_settings; + + set_radial_basis(radial_basis, this->radial_pointer); + set_cutoff(cutoff_function, this->cutoff_pointer); + + // Assign cutoff matrix. + this->cutoffs = cutoffs; +} +void B1 ::write_to_file(std::ofstream &coeff_file, int coeff_size) { // Report radial basis set. coeff_file << radial_basis << "\n"; diff --git a/src/flare_pp/descriptors/b1.h b/src/flare_pp/descriptors/b1.h index df9ef42f6..69f3f31c2 100644 --- a/src/flare_pp/descriptors/b1.h +++ b/src/flare_pp/descriptors/b1.h @@ -29,6 +29,12 @@ class B1 : public Descriptor { const std::vector &cutoff_hyps, const std::vector &descriptor_settings); + B1(const std::string &radial_basis, const std::string &cutoff_function, + const std::vector &radial_hyps, + const std::vector &cutoff_hyps, + const std::vector &descriptor_settings, + const Eigen::MatrixXd &cutoffs); + DescriptorValues compute_struc(Structure &structure); void write_to_file(std::ofstream &coeff_file, int coeff_size); diff --git a/src/flare_pp/descriptors/b3.cpp b/src/flare_pp/descriptors/b3.cpp index d2325823e..d0b5b5b51 100644 --- a/src/flare_pp/descriptors/b3.cpp +++ b/src/flare_pp/descriptors/b3.cpp @@ -26,6 +26,50 @@ B3 ::B3(const std::string &radial_basis, const std::string &cutoff_function, set_cutoff(cutoff_function, this->cutoff_pointer); } +B3 ::B3(const std::string &radial_basis, const std::string &cutoff_function, + const std::vector &radial_hyps, + const std::vector &cutoff_hyps, + const std::vector &descriptor_settings, + const Eigen::MatrixXd &cutoffs) { + + this->radial_basis = radial_basis; + this->cutoff_function = cutoff_function; + this->radial_hyps = radial_hyps; + this->cutoff_hyps = cutoff_hyps; + this->descriptor_settings = descriptor_settings; + + set_radial_basis(radial_basis, this->radial_pointer); + set_cutoff(cutoff_function, this->cutoff_pointer); + + // Assign cutoff matrix. + this->cutoffs = cutoffs; +} + +void B3 ::write_to_file(std::ofstream &coeff_file, int coeff_size) { + // Report radial basis set. + coeff_file << radial_basis << "\n"; + + // Record number of species, nmax, lmax, and the cutoff. + int n_species = descriptor_settings[0]; + int n_max = descriptor_settings[1]; + int l_max = descriptor_settings[2]; + double cutoff = radial_hyps[1]; + + coeff_file << n_species << " " << n_max << " " << l_max << " "; + coeff_file << coeff_size << "\n"; + coeff_file << cutoff_function << "\n"; + + // Report cutoffs to 2 decimal places. + coeff_file << std::fixed << std::setprecision(2); + for (int i = 0; i < n_species; i ++){ + for (int j = 0; j < n_species; j ++){ + coeff_file << cutoffs(i, j) << " "; + } + } + coeff_file << "\n"; +} + + DescriptorValues B3 ::compute_struc(Structure &structure) { // Initialize descriptor values. diff --git a/src/flare_pp/descriptors/b3.h b/src/flare_pp/descriptors/b3.h index b6da7ee8f..461c35b6d 100644 --- a/src/flare_pp/descriptors/b3.h +++ b/src/flare_pp/descriptors/b3.h @@ -29,6 +29,12 @@ class B3 : public Descriptor { const std::vector &cutoff_hyps, const std::vector &descriptor_settings); + B3(const std::string &radial_basis, const std::string &cutoff_function, + const std::vector &radial_hyps, + const std::vector &cutoff_hyps, + const std::vector &descriptor_settings, + const Eigen::MatrixXd &cutoffs); + DescriptorValues compute_struc(Structure &structure); nlohmann::json return_json(); From 57e50a5d4b67df035d6a3ba755ce035765d411ea Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:14:19 -0400 Subject: [PATCH 02/19] b3 complex multiple cutoffs --- src/flare_pp/descriptors/b3.cpp | 24 +++++++++++++++++------- src/flare_pp/descriptors/b3.h | 5 +++-- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/flare_pp/descriptors/b3.cpp b/src/flare_pp/descriptors/b3.cpp index d0b5b5b51..72a112350 100644 --- a/src/flare_pp/descriptors/b3.cpp +++ b/src/flare_pp/descriptors/b3.cpp @@ -88,7 +88,7 @@ DescriptorValues B3 ::compute_struc(Structure &structure) { complex_single_bond(single_bond_vals, force_dervs, neighbor_coords, unique_neighbor_count, cumulative_neighbor_count, descriptor_indices, radial_pointer, cutoff_pointer, nos, - N, lmax, radial_hyps, cutoff_hyps, structure); + N, lmax, radial_hyps, cutoff_hyps, structure, cutoffs); // Compute descriptor values. Eigen::MatrixXd B3_vals, B3_force_dervs; @@ -281,7 +281,7 @@ void compute_B3(Eigen::MatrixXd &B3_vals, Eigen::MatrixXd &B3_force_dervs, } } -void complex_single_bond( +void complex_single_bond_multiple_cutoffs( Eigen::MatrixXcd &single_bond_vals, Eigen::MatrixXcd &force_dervs, Eigen::MatrixXd &neighbor_coordinates, Eigen::VectorXi &neighbor_count, Eigen::VectorXi &cumulative_neighbor_count, @@ -293,14 +293,12 @@ void complex_single_bond( std::vector)> cutoff_function, int nos, int N, int lmax, const std::vector &radial_hyps, - const std::vector &cutoff_hyps, const Structure &structure) { + const std::vector &cutoff_hyps, const Structure &structure, + const Eigen::MatrixXd &cutoffs) { int n_atoms = structure.noa; int n_neighbors = structure.n_neighbors; - // TODO: Make rcut an attribute of the descriptor calculator. - double rcut = radial_hyps[1]; - // Count atoms inside the descriptor cutoff. neighbor_count = Eigen::VectorXi::Zero(n_atoms); Eigen::VectorXi store_neighbors = Eigen::VectorXi::Zero(n_neighbors); @@ -308,9 +306,12 @@ void complex_single_bond( for (int i = 0; i < n_atoms; i++) { int i_neighbors = structure.neighbor_count(i); int rel_index = structure.cumulative_neighbor_count(i); + int central_species = structure.species[i]; for (int j = 0; j < i_neighbors; j++) { int current_count = neighbor_count(i); int neigh_index = rel_index + j; + int neighbor_species = structure.neighbor_species(neigh_index); + double rcut = cutoffs(central_species, neighbor_species); double r = structure.relative_positions(neigh_index, 0); // Check that atom is within descriptor cutoff. if (r <= rcut) { @@ -356,6 +357,10 @@ void complex_single_bond( int i_neighbors = structure.neighbor_count(i); int rel_index = structure.cumulative_neighbor_count(i); int neighbor_index = cumulative_neighbor_count(i); + int central_species = structure.species[i]; + + // Initialize radial hyperparameters. + std::vector new_radial_hyps = radial_hyps; // Initialize radial and spherical harmonic vectors. std::vector g = std::vector(N, 0); @@ -370,6 +375,8 @@ void complex_single_bond( int s, neigh_index, descriptor_counter, unique_ind; for (int j = 0; j < i_neighbors; j++) { neigh_index = rel_index + j; + int neighbor_species = structure.neighbor_species(neigh_index); + double rcut = cutoffs(central_species, neighbor_species); r = structure.relative_positions(neigh_index, 0); if (r > rcut) continue; // Skip if outside cutoff. @@ -378,6 +385,9 @@ void complex_single_bond( z = structure.relative_positions(neigh_index, 3); s = structure.neighbor_species(neigh_index); + // Reset the endpoint of the radial basis set. + new_radial_hyps[1] = rcut; + // Store neighbor coordinates. neighbor_coordinates(neighbor_index, 0) = x; neighbor_coordinates(neighbor_index, 1) = y; @@ -385,7 +395,7 @@ void complex_single_bond( // Compute radial basis values and spherical harmonics. calculate_radial(g, gx, gy, gz, radial_function, cutoff_function, x, y, z, - r, rcut, N, radial_hyps, cutoff_hyps); + r, rcut, N, new_radial_hyps, cutoff_hyps); get_complex_Y(h, hx, hy, hz, x, y, z, lmax); // Store the products and their derivatives. diff --git a/src/flare_pp/descriptors/b3.h b/src/flare_pp/descriptors/b3.h index 461c35b6d..f1e2924d6 100644 --- a/src/flare_pp/descriptors/b3.h +++ b/src/flare_pp/descriptors/b3.h @@ -49,7 +49,7 @@ void compute_B3(Eigen::MatrixXd &B3_vals, Eigen::MatrixXd &B3_force_dervs, const Eigen::VectorXi &descriptor_indices, int nos, int N, int lmax, const Eigen::VectorXd &wigner3j_coeffs); -void complex_single_bond( +void complex_single_bond_multiple_cutoffs( Eigen::MatrixXcd &single_bond_vals, Eigen::MatrixXcd &force_dervs, Eigen::MatrixXd &neighbor_coordinates, Eigen::VectorXi &neighbor_count, Eigen::VectorXi &cumulative_neighbor_count, @@ -61,6 +61,7 @@ void complex_single_bond( std::vector)> cutoff_function, int nos, int N, int lmax, const std::vector &radial_hyps, - const std::vector &cutoff_hyps, const Structure &structure); + const std::vector &cutoff_hyps, const Structure &structure, + const Eigen::MatrixXd &cutoffs); #endif From e0c4d51301e4c5c228e90b648739bf10d301e1af Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:17:27 -0400 Subject: [PATCH 03/19] fix syntax --- src/flare_pp/descriptors/b3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/flare_pp/descriptors/b3.cpp b/src/flare_pp/descriptors/b3.cpp index 72a112350..3e49db1fc 100644 --- a/src/flare_pp/descriptors/b3.cpp +++ b/src/flare_pp/descriptors/b3.cpp @@ -85,7 +85,7 @@ DescriptorValues B3 ::compute_struc(Structure &structure) { int N = descriptor_settings[1]; int lmax = descriptor_settings[2]; - complex_single_bond(single_bond_vals, force_dervs, neighbor_coords, + complex_single_bond_multiple_cutoffs(single_bond_vals, force_dervs, neighbor_coords, unique_neighbor_count, cumulative_neighbor_count, descriptor_indices, radial_pointer, cutoff_pointer, nos, N, lmax, radial_hyps, cutoff_hyps, structure, cutoffs); From 0dcd1c40bc40aac5bf12704ea08aae1d2a17a085 Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:33:35 -0400 Subject: [PATCH 04/19] solve conflict --- src/ace_binding.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/ace_binding.cpp b/src/ace_binding.cpp index 39ad08c86..16e671020 100644 --- a/src/ace_binding.cpp +++ b/src/ace_binding.cpp @@ -150,9 +150,6 @@ PYBIND11_MODULE(_C_flare, m) { py::class_(m, "B3") .def(py::init &, const std::vector &, -<<<<<<< HEAD - const std::vector &>()); -======= const std::vector &>()) .def(py::init &, const std::vector &, @@ -164,7 +161,6 @@ PYBIND11_MODULE(_C_flare, m) { .def_readonly("cutoff_hyps", &B3::cutoff_hyps) .def_readonly("cutoffs", &B3::cutoffs) .def_readonly("descriptor_settings", &B3::descriptor_settings); ->>>>>>> f251f1db (add multiple cutoffs input to b1 and b3) // Kernel functions py::class_(m, "Kernel"); From 5361d59f065f4307d6e5a8fa8a9031e52f41a01a Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:41:41 -0400 Subject: [PATCH 05/19] add cutoffs to b3 --- src/flare_pp/descriptors/b3.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/flare_pp/descriptors/b3.h b/src/flare_pp/descriptors/b3.h index f1e2924d6..8a79b5d03 100644 --- a/src/flare_pp/descriptors/b3.h +++ b/src/flare_pp/descriptors/b3.h @@ -22,6 +22,8 @@ class B3 : public Descriptor { std::string descriptor_name = "B3"; + Eigen::MatrixXd cutoffs; + B3(); B3(const std::string &radial_basis, const std::string &cutoff_function, From 07368bf0d56289ff6ec8b3ec1cd4cc6996f276f7 Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:45:12 -0400 Subject: [PATCH 06/19] add write_to_file --- src/flare_pp/descriptors/b3.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/flare_pp/descriptors/b3.h b/src/flare_pp/descriptors/b3.h index 8a79b5d03..32a144e83 100644 --- a/src/flare_pp/descriptors/b3.h +++ b/src/flare_pp/descriptors/b3.h @@ -39,6 +39,8 @@ class B3 : public Descriptor { DescriptorValues compute_struc(Structure &structure); + void write_to_file(std::ofstream &coeff_file, int coeff_size); + nlohmann::json return_json(); }; From dbc59f77dab7d3ea250fc6434a6ff6f5a72e1560 Mon Sep 17 00:00:00 2001 From: YuuuXie Date: Tue, 16 May 2023 14:49:16 -0400 Subject: [PATCH 07/19] add iostream to b3.cpp --- src/flare_pp/descriptors/b3.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/flare_pp/descriptors/b3.cpp b/src/flare_pp/descriptors/b3.cpp index 3e49db1fc..0b05f26fb 100644 --- a/src/flare_pp/descriptors/b3.cpp +++ b/src/flare_pp/descriptors/b3.cpp @@ -5,6 +5,8 @@ #include "structure.h" #include "wigner3j.h" #include "y_grad.h" +#include // File operations +#include // setprecision #include B3 ::B3() {} From 8f07592aa4bf50d1a6d3318424287ce94a9be8fa Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 14:58:58 -0400 Subject: [PATCH 08/19] fix cutoffs of b3 --- src/flare_pp/descriptors/b3.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/flare_pp/descriptors/b3.cpp b/src/flare_pp/descriptors/b3.cpp index 0b05f26fb..d89225e51 100644 --- a/src/flare_pp/descriptors/b3.cpp +++ b/src/flare_pp/descriptors/b3.cpp @@ -26,6 +26,11 @@ B3 ::B3(const std::string &radial_basis, const std::string &cutoff_function, set_radial_basis(radial_basis, this->radial_pointer); set_cutoff(cutoff_function, this->cutoff_pointer); + + // Create cutoff matrix. + int n_species = descriptor_settings[0]; + double cutoff_val = radial_hyps[1]; + cutoffs = Eigen::MatrixXd::Constant(n_species, n_species, cutoff_val); } B3 ::B3(const std::string &radial_basis, const std::string &cutoff_function, From 8d1bea7ab4699bb061f4e1afd17addeea026ae67 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 15:08:18 -0400 Subject: [PATCH 09/19] add randomness in unit tests --- tests/get_sgp.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/tests/get_sgp.py b/tests/get_sgp.py index d42d1a147..f1f869ef1 100644 --- a/tests/get_sgp.py +++ b/tests/get_sgp.py @@ -8,21 +8,16 @@ from ase.build import make_supercell # Define kernel. -sigma = 2.0 +sigma = np.random.rand() + 1.0 power = 1.0 dotprod_kernel = DotProduct(sigma, power) normdotprod_kernel = NormalizedDotProduct(sigma, power) # Define remaining parameters for the SGP wrapper. -sigma_e = 0.3 -sigma_f = 0.2 -sigma_s = 0.1 -species_map = {6: 0, 8: 1} -single_atom_energies = {0: -5, 1: -6} variance_type = "local" max_iterations = 20 opt_method = "L-BFGS-B" -bounds = [(None, None), (sigma_e, None), (None, None), (None, None)] +bounds = [(None, None), (None, None), (None, None), (None, None)] def get_random_atoms(a=2.0, sc_size=2, numbers=[6, 8], set_seed: int = None): @@ -89,6 +84,18 @@ def get_empty_sgp(n_types=2, power=2, multiple_cutoff=False, kernel_type="Normal cutoff_matrix, ) + sigma_e = np.random.rand() + sigma_f = np.random.rand() + sigma_s = np.random.rand() + print("Hyps", kernel.sigma, sigma_e, sigma_f, sigma_s) + + if n_types == 1: + species_map = {6: 0} + single_atom_energies = {0: np.random.rand()} + elif n_types == 2: + species_map = {6: 0, 8: 1} + single_atom_energies = {0: np.random.rand(), 1: np.random.rand()} + empty_sgp = SGP_Wrapper( [kernel], [b2_calc], @@ -123,16 +130,18 @@ def get_updated_sgp(n_types=2, power=2, multiple_cutoff=False, kernel_type="Norm energy = training_structure.get_potential_energy() stress = training_structure.get_stress() + size = max(1, np.random.randint(len(training_structure))) + custom_range = np.random.choice(len(training_structure), size=size, replace=False).tolist() sgp.update_db( training_structure, forces, - custom_range=(1, 2, 3, 4, 5), + custom_range=custom_range, energy=energy, stress=stress, mode="specific", - rel_e_noise=0.1, - rel_f_noise=0.2, - rel_s_noise=0.1, + rel_e_noise=np.random.rand(), + rel_f_noise=np.random.rand(), + rel_s_noise=np.random.rand(), ) # add an isolated atom to the training data From 34ae88dc287d3786de7c806f9bb5b7ac5b089797 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 15:10:13 -0400 Subject: [PATCH 10/19] change varmap in norm dotproduct to use qr --- .../kernels/normalized_dot_product.cpp | 58 +++++-------------- 1 file changed, 14 insertions(+), 44 deletions(-) diff --git a/src/flare_pp/kernels/normalized_dot_product.cpp b/src/flare_pp/kernels/normalized_dot_product.cpp index caecef892..587e93171 100644 --- a/src/flare_pp/kernels/normalized_dot_product.cpp +++ b/src/flare_pp/kernels/normalized_dot_product.cpp @@ -784,8 +784,8 @@ Eigen::MatrixXd NormalizedDotProduct ::compute_varmap_coefficients( int beta_size = p_size * (p_size + 1) / 2; int n_species = gp_model.sparse_descriptors[kernel_index].n_types; int n_sparse = gp_model.sparse_descriptors[kernel_index].n_clusters; - //mapping_coeffs = Eigen::MatrixXd::Zero(n_species * n_species, p_size * p_size); // can be reduced by symmetry mapping_coeffs = Eigen::MatrixXd::Zero(n_species, p_size * p_size); // can be reduced by symmetry + double jitter_root = pow(gp_model.Kuu_jitter, 0.5); // Get alpha index. @@ -797,58 +797,28 @@ Eigen::MatrixXd NormalizedDotProduct ::compute_varmap_coefficients( // Loop over types. for (int s = 0; s < n_species; s++){ int n_types = gp_model.sparse_descriptors[kernel_index].n_clusters_by_type[s]; - int c_types = - gp_model.sparse_descriptors[kernel_index].cumulative_type_count[s]; - int K_ind = alpha_ind + c_types; - + Eigen::MatrixXd pi = gp_model.sparse_descriptors[kernel_index].descriptors[s]; // Loop over clusters within each type. for (int i = 0; i < n_types; i++){ - Eigen::VectorXd pi_current = - gp_model.sparse_descriptors[kernel_index].descriptors[s].row(i); double pi_norm = gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](i); // Skip empty environments. if (pi_norm < empty_thresh) continue; + pi.row(i) /= pi_norm; // normalize descriptor by norm + } - // TODO: include symmetry of i & j - // Loop over clusters within each type. - for (int j = 0; j < n_types; j++){ - Eigen::VectorXd pj_current = - gp_model.sparse_descriptors[kernel_index].descriptors[s].row(j); - double pj_norm = - gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](j); - - // Skip empty environments. - if (pj_norm < empty_thresh) - continue; - - double Kuu_inv_ij = gp_model.Kuu_inverse(K_ind + i, K_ind + j); - double Kuu_inv_ij_normed = Kuu_inv_ij / pi_norm / pj_norm; -// double Sigma_ij = gp_model.Sigma(K_ind + i, K_ind + j); -// double Sigma_ij_normed = Sigma_ij / pi_norm / pj_norm; - int beta_count = 0; - - // First loop over descriptor values. - for (int k = 0; k < p_size; k++) { - double p_ik = pi_current(k); - - // Second loop over descriptor values. - for (int l = 0; l < p_size; l++){ - double p_jl = pj_current(l); - - // Update beta vector. - double beta_val = sig2 * sig2 * p_ik * p_jl * (- Kuu_inv_ij_normed); // + Sigma_ij_normed); // To match the compute_cluster_uncertainty function - mapping_coeffs(s, beta_count) += beta_val; - - if (k == l && i == 0 && j == 0) { - mapping_coeffs(s, beta_count) += sig2; // the self kernel term - } - - beta_count++; - } - } + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(p_size + n_types, n_types); + A.block(0, 0, p_size, n_types) = pi.transpose(); + A.block(p_size, 0, n_types, n_types) = jitter_root * Eigen::MatrixXd::Identity(n_types, n_types) / sigma; + Eigen::HouseholderQR qr(A); + Eigen::MatrixXd Q = qr.householderQ(); // (p_size + n_types, n_types) + int beta_count = 0; + for (int i = 0; i < p_size; i++) { + for (int j = 0; j < p_size; j++) { + if (j < n_types) mapping_coeffs(s, beta_count) += Q(i, j); + beta_count++; } } } From 9bfb46a61034a49abd6af4642a5f4a2e8bd5cb9f Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 15:14:32 -0400 Subject: [PATCH 11/19] change varmap in dotproduct to use qr --- src/flare_pp/kernels/dot_product.cpp | 65 +++++----------------------- 1 file changed, 12 insertions(+), 53 deletions(-) diff --git a/src/flare_pp/kernels/dot_product.cpp b/src/flare_pp/kernels/dot_product.cpp index f61482ef3..040b20783 100644 --- a/src/flare_pp/kernels/dot_product.cpp +++ b/src/flare_pp/kernels/dot_product.cpp @@ -786,8 +786,8 @@ Eigen::MatrixXd DotProduct ::compute_varmap_coefficients( int beta_size = p_size * (p_size + 1) / 2; int n_species = gp_model.sparse_descriptors[kernel_index].n_types; int n_sparse = gp_model.sparse_descriptors[kernel_index].n_clusters; - //mapping_coeffs = Eigen::MatrixXd::Zero(n_species * n_species, p_size * p_size); // can be reduced by symmetry mapping_coeffs = Eigen::MatrixXd::Zero(n_species, p_size * p_size); // can be reduced by symmetry + double jitter_root = pow(gp_model.Kuu_jitter, 0.5); // Get alpha index. @@ -799,58 +799,17 @@ Eigen::MatrixXd DotProduct ::compute_varmap_coefficients( // Loop over types. for (int s = 0; s < n_species; s++){ int n_types = gp_model.sparse_descriptors[kernel_index].n_clusters_by_type[s]; - int c_types = - gp_model.sparse_descriptors[kernel_index].cumulative_type_count[s]; - int K_ind = alpha_ind + c_types; - - // Loop over clusters within each type. - for (int i = 0; i < n_types; i++){ - Eigen::VectorXd pi_current = - gp_model.sparse_descriptors[kernel_index].descriptors[s].row(i); - double pi_norm = - gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](i); - - // Skip empty environments. - if (pi_norm < empty_thresh) - continue; - - // TODO: include symmetry of i & j - // Loop over clusters within each type. - for (int j = 0; j < n_types; j++){ - Eigen::VectorXd pj_current = - gp_model.sparse_descriptors[kernel_index].descriptors[s].row(j); - double pj_norm = - gp_model.sparse_descriptors[kernel_index].descriptor_norms[s](j); - - // Skip empty environments. - if (pj_norm < empty_thresh) - continue; - - double Kuu_inv_ij = gp_model.Kuu_inverse(K_ind + i, K_ind + j); - double Kuu_inv_ij_normed = Kuu_inv_ij; // / pi_norm / pj_norm; -// double Sigma_ij = gp_model.Sigma(K_ind + i, K_ind + j); -// double Sigma_ij_normed = Sigma_ij / pi_norm / pj_norm; - int beta_count = 0; - - // First loop over descriptor values. - for (int k = 0; k < p_size; k++) { - double p_ik = pi_current(k); - - // Second loop over descriptor values. - for (int l = 0; l < p_size; l++){ - double p_jl = pj_current(l); - - // Update beta vector. - double beta_val = sig2 * sig2 * p_ik * p_jl * (- Kuu_inv_ij_normed); // + Sigma_ij_normed); // To match the compute_cluster_uncertainty function - mapping_coeffs(s, beta_count) += beta_val; - - if (k == l && i == 0 && j == 0) { - mapping_coeffs(s, beta_count) += sig2; // the self kernel term - } - - beta_count++; - } - } + Eigen::MatrixXd pi = gp_model.sparse_descriptors[kernel_index].descriptors[s]; + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(p_size + n_types, n_types); + A.block(0, 0, p_size, n_types) = pi.transpose(); + A.block(p_size, 0, n_types, n_types) = jitter_root * Eigen::MatrixXd::Identity(n_types, n_types) / sigma; + Eigen::HouseholderQR qr(A); + Eigen::MatrixXd Q = qr.householderQ(); // (p_size + n_types, n_types) + int beta_count = 0; + for (int i = 0; i < p_size; i++) { + for (int j = 0; j < p_size; j++) { + if (j < n_types) mapping_coeffs(s, beta_count) += Q(i, j); + beta_count++; } } } From 889da121b6276ee0f7a1c012805c25de8a4da5ce Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 15:18:22 -0400 Subject: [PATCH 12/19] change lmp mapped var to use qr --- lammps_plugins/compute_flare_std_atom.cpp | 33 ++++++++++++----------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/lammps_plugins/compute_flare_std_atom.cpp b/lammps_plugins/compute_flare_std_atom.cpp index 16f01c091..1427f0e28 100644 --- a/lammps_plugins/compute_flare_std_atom.cpp +++ b/lammps_plugins/compute_flare_std_atom.cpp @@ -180,10 +180,17 @@ void ComputeFlareStdAtom::compute_peratom() { continue; if (use_map) { - int power = 2; - compute_energy_and_u(B2_vals, B2_norm_squared, single_bond_vals, power, - n_species, n_max, l_max, beta_matrices[itype - 1], u, &variance, normalized); - variance /= sig2; + Eigen::VectorXd Q_desc; + double K_self; + if (normalized) { + K_self = 1.0; + double B2_norm = pow(B2_norm_squared, 0.5); + Q_desc = beta_matrices[itype - 1].transpose() * B2_vals / B2_norm; + } else { + K_self = B2_norm_squared / size_norm; // only power 1 is supported + Q_desc = beta_matrices[itype - 1].transpose() * B2_vals; + } + variance = K_self - Q_desc.dot(Q_desc); } else { Eigen::VectorXd kernel_vec = Eigen::VectorXd::Zero(n_clusters); double K_self; @@ -462,18 +469,14 @@ void ComputeFlareStdAtom::read_file(char *filename) { int beta_count = 0; double beta_val; for (int k = 0; k < n_species; k++) { -// for (int l = 0; l < n_species; l++) { - - beta_matrix = Eigen::MatrixXd::Zero(n_descriptors, n_descriptors); - for (int i = 0; i < n_descriptors; i++) { - for (int j = 0; j < n_descriptors; j++) { - //beta_matrix(k * n_descriptors + i, l * n_descriptors + j) = beta[beta_count]; - beta_matrix(i, j) = beta[beta_count]; - beta_count++; - } + beta_matrix = Eigen::MatrixXd::Zero(n_descriptors, n_descriptors); + for (int i = 0; i < n_descriptors; i++) { + for (int j = 0; j < n_descriptors; j++) { + beta_matrix(i, j) = beta[beta_count]; + beta_count++; } - beta_matrices.push_back(beta_matrix); -// } + } + beta_matrices.push_back(beta_matrix); } } From 66b42149933b4b3942ca99d33184a4db574680b0 Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 16:03:41 -0400 Subject: [PATCH 13/19] add kernel type to lammps unit tests --- tests/test_lammps.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_lammps.py b/tests/test_lammps.py index fb7729678..c822596b6 100644 --- a/tests/test_lammps.py +++ b/tests/test_lammps.py @@ -28,7 +28,8 @@ @pytest.mark.parametrize("struc", struc_list) @pytest.mark.parametrize("multicut", [False, True]) @pytest.mark.parametrize("n_cpus", n_cpus_list) -def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus): +@pytest.mark.parametrize("kernel_type", ["NormalizedDotProduct", "DotProduct"]) +def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus, kernel_type): """Test the flare pair style.""" if n_species > n_types: @@ -100,11 +101,11 @@ def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus): energy_lmp += sgp_model.gp_model.single_atom_energies[coded_spec] thresh = 1e-6 + r_thresh = 1e-3 print(energy, energy_lmp) - assert np.allclose(energy, energy_lmp, atol=thresh) - # print(forces, forces_lmp) - assert np.allclose(forces, forces_lmp, atol=thresh) - assert np.allclose(stress, stress_lmp, atol=thresh) + assert np.allclose(energy, energy_lmp, atol=thresh, rtol=r_thresh) + assert np.allclose(forces, forces_lmp, atol=thresh, rtol=r_thresh) + assert np.allclose(stress, stress_lmp, atol=thresh, rtol=r_thresh) # Remove files. os.remove(potential_name) @@ -119,7 +120,6 @@ def test_write_potential(n_species, n_types, power, struc, multicut, n_cpus): from ase.calculators.lammpsrun import LAMMPS from flare.md.lammps import LAMMPS_MOD, LAMMPS_MD, get_kinetic_stress - @pytest.mark.skipif( not os.environ.get("lmp", False), reason=( From 7c347d9912b5d09f56457ea444fdf25a91e319eb Mon Sep 17 00:00:00 2001 From: Yu Xie Date: Tue, 16 May 2023 16:08:54 -0400 Subject: [PATCH 14/19] rm norm squared for dot product in lmp --- lammps_plugins/compute_flare_std_atom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lammps_plugins/compute_flare_std_atom.cpp b/lammps_plugins/compute_flare_std_atom.cpp index 1427f0e28..d63cf09c7 100644 --- a/lammps_plugins/compute_flare_std_atom.cpp +++ b/lammps_plugins/compute_flare_std_atom.cpp @@ -187,7 +187,7 @@ void ComputeFlareStdAtom::compute_peratom() { double B2_norm = pow(B2_norm_squared, 0.5); Q_desc = beta_matrices[itype - 1].transpose() * B2_vals / B2_norm; } else { - K_self = B2_norm_squared / size_norm; // only power 1 is supported + K_self = B2_norm_squared; // only power 1 is supported Q_desc = beta_matrices[itype - 1].transpose() * B2_vals; } variance = K_self - Q_desc.dot(Q_desc); From 01754f57bec7ac8c5817d6aef7b9ed93f51f500f Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 20 Jul 2023 11:31:34 +0000 Subject: [PATCH 15/19] Bump pygments from 2.11.2 to 2.15.0 Bumps [pygments](https://github.com/pygments/pygments) from 2.11.2 to 2.15.0. - [Release notes](https://github.com/pygments/pygments/releases) - [Changelog](https://github.com/pygments/pygments/blob/master/CHANGES) - [Commits](https://github.com/pygments/pygments/compare/2.11.2...2.15.0) --- updated-dependencies: - dependency-name: pygments dependency-type: direct:production ... Signed-off-by: dependabot[bot] --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 4f6cf6233..c05b69b03 100644 --- a/requirements.txt +++ b/requirements.txt @@ -15,6 +15,6 @@ wandb docutils==0.17.1 alabaster==0.7.12 babel -pygments==2.11.2 +pygments==2.15.0 nptyping nbsphinx From c249459a06363e4d5c27d04ebc3de177ccb69a74 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 15 Sep 2024 14:24:49 +0000 Subject: [PATCH 16/19] Bump mattnotmitt/doxygen-action from 1.1.0 to 1.9.8 Bumps [mattnotmitt/doxygen-action](https://github.com/mattnotmitt/doxygen-action) from 1.1.0 to 1.9.8. - [Release notes](https://github.com/mattnotmitt/doxygen-action/releases) - [Commits](https://github.com/mattnotmitt/doxygen-action/compare/v1.1.0...v1.9.8) --- updated-dependencies: - dependency-name: mattnotmitt/doxygen-action dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/flare.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flare.yml b/.github/workflows/flare.yml index cf09ba3df..eaeb03f8a 100644 --- a/.github/workflows/flare.yml +++ b/.github/workflows/flare.yml @@ -118,7 +118,7 @@ jobs: sudo apt-get install python3-sphinx python3-sphinx-rtd-theme python3-breathe python3-nbsphinx - name: Run Doxygen - uses: mattnotmitt/doxygen-action@v1.1.0 + uses: mattnotmitt/doxygen-action@v1.9.8 with: # Path to Doxyfile doxyfile-path: "./Doxyfile" # default is ./Doxyfile From 67d15e79abddda20b266b51ff8933f46c955aa81 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 15 Sep 2024 14:24:51 +0000 Subject: [PATCH 17/19] Bump actions/checkout from 2 to 4 Bumps [actions/checkout](https://github.com/actions/checkout) from 2 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v2...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/flare.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flare.yml b/.github/workflows/flare.yml index cf09ba3df..bed8086a0 100644 --- a/.github/workflows/flare.yml +++ b/.github/workflows/flare.yml @@ -30,7 +30,7 @@ jobs: # Steps represent a sequence of tasks that will be executed as part of the job steps: # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v4 with: From 8a741a2d06e588cd6213627b4456626916a95fb8 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 15 Sep 2024 14:24:53 +0000 Subject: [PATCH 18/19] Bump actions/setup-python from 4 to 5 Bumps [actions/setup-python](https://github.com/actions/setup-python) from 4 to 5. - [Release notes](https://github.com/actions/setup-python/releases) - [Commits](https://github.com/actions/setup-python/compare/v4...v5) --- updated-dependencies: - dependency-name: actions/setup-python dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/flare.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flare.yml b/.github/workflows/flare.yml index cf09ba3df..47fb50aa8 100644 --- a/.github/workflows/flare.yml +++ b/.github/workflows/flare.yml @@ -32,7 +32,7 @@ jobs: # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it - uses: actions/checkout@v2 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} From 413a2341ded6b0fa9c300154a18b8c988db19864 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 15 Sep 2024 14:24:55 +0000 Subject: [PATCH 19/19] Bump peaceiris/actions-gh-pages from 3 to 4 Bumps [peaceiris/actions-gh-pages](https://github.com/peaceiris/actions-gh-pages) from 3 to 4. - [Release notes](https://github.com/peaceiris/actions-gh-pages/releases) - [Changelog](https://github.com/peaceiris/actions-gh-pages/blob/main/CHANGELOG.md) - [Commits](https://github.com/peaceiris/actions-gh-pages/compare/v3...v4) --- updated-dependencies: - dependency-name: peaceiris/actions-gh-pages dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/flare.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flare.yml b/.github/workflows/flare.yml index cf09ba3df..08564da84 100644 --- a/.github/workflows/flare.yml +++ b/.github/workflows/flare.yml @@ -134,7 +134,7 @@ jobs: make html - name: Publish the docs - uses: peaceiris/actions-gh-pages@v3 + uses: peaceiris/actions-gh-pages@v4 with: github_token: ${{ secrets.GITHUB_TOKEN }} # Default Doxyfile build documentation to html directory.