From b1ccb9dad6608ca4723892182c2fffe541cc2655 Mon Sep 17 00:00:00 2001 From: Ehsan Rahmatizad Date: Thu, 24 Jul 2025 15:20:58 +0200 Subject: [PATCH 1/4] ASCF Anlytical Derivative free boundary condition added --- dscribe/descriptors/acsf.py | 88 ++++++++++- dscribe/ext/acsf.cpp | 284 +++++++++++++++++++++++++++++++++++- dscribe/ext/acsf.h | 13 +- dscribe/ext/ext.cpp | 10 ++ setup.py | 2 +- 5 files changed, 392 insertions(+), 5 deletions(-) diff --git a/dscribe/descriptors/acsf.py b/dscribe/descriptors/acsf.py index 1eda8f80..36c14fc8 100644 --- a/dscribe/descriptors/acsf.py +++ b/dscribe/descriptors/acsf.py @@ -19,6 +19,7 @@ from dscribe.descriptors.descriptorlocal import DescriptorLocal from dscribe.ext import ACSFWrapper +from dscribe.ext import CellList class ACSF(DescriptorLocal): @@ -221,8 +222,17 @@ def validate_derivatives_method(self, method, attach): raise ValueError( "ACSF derivatives can only be calculated with attach=True." ) - return super().validate_derivatives_method(method, attach) - + if method == "auto": + method = "analytical" + return method + elif method == "analytical": + return method + elif method == "numerical": + return super().validate_derivatives_method(method, attach) + else: + raise ValueError( + "%s derivative method not implemented in ACSF derivatives"%method + ) @property def species(self): return self._species @@ -396,3 +406,77 @@ def g5_params(self): @g5_params.setter def g5_params(self, value): self.acsf_wrapper.g5_params = self.validate_g5_params(value) + + + def derivatives_analytical( + self, + d, + c, + system, + centers, + indices, + attach, + return_descriptor=True, + ): + """Return the analytical derivatives for the given system. + Args: + system (:class:`ase.Atoms`): Atomic structure. + indices (list): Indices of atoms for which the derivatives will be computed. + return_descriptor (bool): Whether to also calculate the descriptor + in the same function call. This is true by default as it + typically is faster to calculate both in one go. + Returns: + If return_descriptor is True, returns a tuple, where the first item + is the derivative array and the second is the descriptor array. + Otherwise only returns the derivatives array. The derivatives array + is a 3D numpy array. The dimensions are: [n_atoms, 3, n_features]. + The first dimension goes over the included atoms. The order is same + as the order of atoms in the given system. The second dimension + goes over the cartesian components, x, y and z. The last dimension + goes over the features in the default order. + """ + + # Validate and normalize system + positions = self.validate_positions(system.get_positions()) + atomic_numbers = self.validate_atomic_numbers(system.get_atomic_numbers()) + pbc = self.validate_pbc(system.get_pbc()) + cell = self.validate_cell(system.get_cell(), pbc) + + # Create C-compatible list of atomic indices for which the ACSF is + # calculated + if centers is None: + centers = np.arange(len(system), dtype=np.int32) + else: + centers = np.asarray(centers, dtype=np.int32) + """ + if desc_centers is None: + desc_centers = np.arange(len(system), dtype=np.int32) + else: + desc_centers = np.asarray(desc_centers, dtype=np.int32) + + if grad_centers is None: + grad_centers = np.arange(len(system), dtype=np.int32) + else: + grad_centers = np.asarray(grad_centers, dtype=np.int32) + """ + # --- CRUCIAL CHANGE HERE: Construct the CellList object in Python --- + # This assumes that CellList is exposed in your pybind11 wrapper as self.acsf_wrapper.CellList + cell_list_obj = CellList(positions, self.r_cut) + + # Call the C++ function with the correct arguments and order + self.acsf_wrapper.derivatives_analytical( + d, # 1st arg: py::array_t derivatives + c, # 2nd arg: py::array_t descriptor + atomic_numbers, # 3rd arg: py::array_t atomic_numbers + positions, # 4th arg: py::array_t atomic_positions + cell_list_obj, # 5th arg: CellList cell_list (the constructed object) + centers, # 6th arg: py::array_t desc_centers + centers, # 7th arg: py::array_t grad_centers + return_descriptor # 8th arg: const bool return_descriptor + ) + + if return_descriptor: + return d, c + else: + return d + diff --git a/dscribe/ext/acsf.cpp b/dscribe/ext/acsf.cpp index 50a9a375..02f5c845 100644 --- a/dscribe/ext/acsf.cpp +++ b/dscribe/ext/acsf.cpp @@ -109,6 +109,9 @@ void ACSF::create( CellList cell_list ) { + // This overload is likely for Python calls where centers are doubles, + // but the actual implementation expects int indices for centers. + // It's left empty as per the original acsf_new.cpp. } void ACSF::create( @@ -154,7 +157,7 @@ void ACSF::create( if (g4_params.size() != 0 || g5_params.size() != 0) { for (int k_neighbour = 0; k_neighbour < n_neighbours; ++k_neighbour) { int k = neighbours_i.indices[k_neighbour]; - if (k >= j) { + if (k >= j) { // Ensure k < j to avoid duplicate triplets and self-interaction continue; } @@ -259,3 +262,282 @@ inline void ACSF::compute_g5(py::detail::unchecked_mutable_reference offset++; } } + + +void ACSF::derivatives_analytical( + py::array_t derivatives, + py::array_t descriptor, + py::array_t atomic_numbers, + py::array_t atomic_positions, + CellList cell_list, + py::array_t desc_centers, + py::array_t grad_centers, // we want derivatives w.r.t. these atoms + const bool return_descriptor +) { + int n_desc_centers = desc_centers.shape(0); + + auto descriptor_mu = descriptor.mutable_unchecked<2>(); // [n_desc_centers, n_features] + auto derivatives_mu = derivatives.mutable_unchecked<4>(); // [n_desc_centers, n_grad_centers, 3, n_features] + + auto atomic_numbers_u = atomic_numbers.unchecked<1>(); + auto atomic_positions_u = atomic_positions.unchecked<2>(); + + auto desc_centers_u = desc_centers.unchecked<1>(); + auto grad_centers_u = grad_centers.unchecked<1>(); + + // d[idxi,idxj,c,l] = derivative of l-th ACSF of atom desc_centers_u(idxi), w.r.t. atom grad_centers_u(idxj) coordinate c + // descriptor_mu(idxi, l) = l-th ACSF of atom desc_centers_u(idxi) + + // Declare eta, zeta, lambda here to make them visible in the loops + double eta, zeta, lambda; + + for(int idxi=0; idxi= j_neighbour_idx) { // Ensure k index is less than j index to avoid duplicates (and self-k) + continue; + } + + // Calculate j-k distance: it is not contained in the cell lists. + double dx = atomic_positions_u(x, 0) - atomic_positions_u(y, 0); + double dy = atomic_positions_u(x, 1) - atomic_positions_u(y, 1); + double dz = atomic_positions_u(x, 2) - atomic_positions_u(y, 2); + double r_jk_square = dx*dx + dy*dy + dz*dz; + double r_jk = sqrt(r_jk_square); + + // Precompute values for G4 and G5 + double r_ik = neighbours_i.distances[k_neighbour_idx]; // Corrected access + double fc_ik = compute_cutoff(r_ik); + double r_ik_square = neighbours_i.distancesSquared[k_neighbour_idx]; // Corrected access + // r_ij_square is already defined from the outer loop + int index_k = atomic_number_to_index_map[atomic_numbers_u[y]]; + double costheta = 0.5*(r_ij_square+r_ik_square-r_jk_square)/(r_ij*r_ik); + double fc4 = fc_ij*fc_ik*compute_cutoff(r_jk); // G4 cutoff product, compute_cutoff for r_jk + double fc5 = fc_ij*fc_ik; // G5 cutoff product + + // Unit vectors for i-k and j-k + double e_ik[3], e_jk[3]; + for(int c=0; c<3; c++){ + e_ik[c] = (atomic_positions_u(i,c) - atomic_positions_u(y,c)) / r_ik; + e_jk[c] = (atomic_positions_u(x,c) - atomic_positions_u(y,c)) / r_jk; + } + + // Determine the location for this triplet of species + int its; + if (index_j >= index_k) its = (index_j*(index_j+1))/2 + index_k; + else its = (index_k*(index_k+1))/2 + index_j; + + // Offset for G4 and G5 features + offset = n_types * (1+n_g2+n_g3); // Skip G1, G2, G3 features for all types + offset += its * (n_g4+n_g5); // Skip G4 and G5 features for other type pairs + int o0_g4_g5 = offset; // Store initial offset for derivatives calculation + + // If return_descriptor is true, compute and add descriptor values + if(return_descriptor){ + int current_offset_for_descriptor = offset; // Use a temporary offset for descriptor calculation + compute_g4(descriptor_mu, idxi, current_offset_for_descriptor, costheta, r_jk, r_ij_square, r_ik_square, r_jk_square, fc_ij, fc_ik); + compute_g5(descriptor_mu, idxi, current_offset_for_descriptor, costheta, r_ij_square, r_ik_square, fc_ij, fc_ik); + } + + // Initialize variables to prevent uninitialized warnings + double dcostheta = 0.0; + double ang = 0.0; + double dang = 0.0; + double first_term = 0.0; + double secnd_term = 0.0; + double third_term = 0.0; + double final = 0.0; + + // Derivatives of cutoff functions w.r.t. r_ik and r_jk + double d_fc_ik = -(PI*sin((PI*r_ik)/r_cut))/(2.0*r_cut); + double d_fc_jk = -(PI*sin((PI*r_jk)/r_cut))/(2.0*r_cut); + + // G4 Derivatives + if (r_jk <= r_cut){ // Only compute if r_jk is within cutoff for G4 + for (auto params : g4_params){ // loop over G4 functions + eta = params[0]; // eta is declared outside the loop now + zeta = params[1]; // zeta is declared outside the loop now + lambda = params[2]; // lambda is declared outside the loop now + + ef = exp(-eta*(r_ij_square+r_ik_square+r_jk_square)); + ang = 2*pow(0.5*(1 + lambda*costheta), zeta); + dang = zeta*pow(0.5*(1 + lambda*costheta), zeta-1)*lambda; + + for(int idxj=0; idxj centers, CellList cell_list ); - + void create( py::array_t out, py::array_t positions, @@ -41,6 +41,17 @@ class ACSF : public DescriptorLocal { CellList cell_list ); + void derivatives_analytical( + py::array_t derivatives, + py::array_t descriptor, + py::array_t atomic_numbers, + py::array_t atomic_positions, + CellList cell_list, + py::array_t desc_centers, + py::array_t grad_centers, + const bool return_descriptor + ); + int get_number_of_features() const; void set_r_cut(double r_cut); diff --git a/dscribe/ext/ext.cpp b/dscribe/ext/ext.cpp index 95ceb474..5ca62592 100644 --- a/dscribe/ext/ext.cpp +++ b/dscribe/ext/ext.cpp @@ -75,6 +75,16 @@ PYBIND11_MODULE(ext, m) { py::class_(m, "ACSFWrapper") .def(py::init > , vector , vector > , vector >, vector , bool>()) .def("create", overload_cast_, py::array_t, py::array_t, py::array_t, py::array_t, py::array_t >()(&DescriptorLocal::create)) + .def("derivatives_analytical", &ACSF::derivatives_analytical, + py::arg("derivatives"), // py::array_t + py::arg("descriptor"), // py::array_t + py::arg("atomic_numbers"), // py::array_t + py::arg("atomic_positions"), // py::array_t + py::arg("cell_list"), // CellList + py::arg("desc_centers"), // py::array_t + py::arg("grad_centers"), // py::array_t + py::arg("return_descriptor") // const bool + ) .def("get_number_of_features", &ACSF::get_number_of_features) .def_readwrite("n_types", &ACSF::n_types) .def_readwrite("n_type_pairs", &ACSF::n_type_pairs) diff --git a/setup.py b/setup.py index d7f40e6f..8c99f585 100644 --- a/setup.py +++ b/setup.py @@ -77,7 +77,7 @@ def __str__(self): get_pybind_include(user=True), ], language='c++', - extra_compile_args=cpp_extra_compile_args + ["-fvisibility=hidden"], # the -fvisibility flag is needed by pybind11 + extra_compile_args=cpp_extra_compile_args + ["-fvisibility=hidden", "-ftemplate-depth=4096"], # the -fvisibility flag is needed by pybind11 extra_link_args=cpp_extra_link_args, ) ] From 1b5d6d2a1d383638e4bc26e85ff8dec8aa9f6318 Mon Sep 17 00:00:00 2001 From: Ehsan Rahmatizad Date: Fri, 25 Jul 2025 09:04:25 +0200 Subject: [PATCH 2/4] Starting implementation of periodic BC in acsf.cpp --- dscribe/descriptors/acsf.py | 13 ++++++------- dscribe/ext/acsf.cpp | 7 ++++++- dscribe/ext/acsf.h | 5 ++++- dscribe/ext/ext.cpp | 3 ++- 4 files changed, 18 insertions(+), 10 deletions(-) diff --git a/dscribe/descriptors/acsf.py b/dscribe/descriptors/acsf.py index 36c14fc8..87aed1f0 100644 --- a/dscribe/descriptors/acsf.py +++ b/dscribe/descriptors/acsf.py @@ -19,7 +19,6 @@ from dscribe.descriptors.descriptorlocal import DescriptorLocal from dscribe.ext import ACSFWrapper -from dscribe.ext import CellList class ACSF(DescriptorLocal): @@ -461,17 +460,17 @@ def derivatives_analytical( """ # --- CRUCIAL CHANGE HERE: Construct the CellList object in Python --- # This assumes that CellList is exposed in your pybind11 wrapper as self.acsf_wrapper.CellList - cell_list_obj = CellList(positions, self.r_cut) # Call the C++ function with the correct arguments and order self.acsf_wrapper.derivatives_analytical( d, # 1st arg: py::array_t derivatives c, # 2nd arg: py::array_t descriptor - atomic_numbers, # 3rd arg: py::array_t atomic_numbers - positions, # 4th arg: py::array_t atomic_positions - cell_list_obj, # 5th arg: CellList cell_list (the constructed object) - centers, # 6th arg: py::array_t desc_centers - centers, # 7th arg: py::array_t grad_centers + atomic_numbers, # 3rd arg: py::array_t desc_centers + cell, # 4th arg: py::array_t desc_centers + pbc, # 5th arg: py::array_t desc_centers + positions, # 6th arg: py::array_t atomic_positions + centers, # 7th arg: py::array_t desc_centers + centers, # 8th arg: py::array_t grad_centers return_descriptor # 8th arg: const bool return_descriptor ) diff --git a/dscribe/ext/acsf.cpp b/dscribe/ext/acsf.cpp index 02f5c845..49071350 100644 --- a/dscribe/ext/acsf.cpp +++ b/dscribe/ext/acsf.cpp @@ -268,12 +268,17 @@ void ACSF::derivatives_analytical( py::array_t derivatives, py::array_t descriptor, py::array_t atomic_numbers, + py::array_t cell, + py::array_t pbc, py::array_t atomic_positions, - CellList cell_list, py::array_t desc_centers, py::array_t grad_centers, // we want derivatives w.r.t. these atoms const bool return_descriptor ) { + + // Calculate neighbours with a cell list + CellList cell_list(atomic_positions, this->cutoff); + int n_desc_centers = desc_centers.shape(0); auto descriptor_mu = descriptor.mutable_unchecked<2>(); // [n_desc_centers, n_features] diff --git a/dscribe/ext/acsf.h b/dscribe/ext/acsf.h index be3f922a..170b1153 100644 --- a/dscribe/ext/acsf.h +++ b/dscribe/ext/acsf.h @@ -4,6 +4,8 @@ #include #include #include "descriptorlocal.h" +#include "celllist.h" +#include "geometry.h" #define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062 @@ -45,8 +47,9 @@ class ACSF : public DescriptorLocal { py::array_t derivatives, py::array_t descriptor, py::array_t atomic_numbers, + py::array_t cell, + py::array_t pbc, py::array_t atomic_positions, - CellList cell_list, py::array_t desc_centers, py::array_t grad_centers, const bool return_descriptor diff --git a/dscribe/ext/ext.cpp b/dscribe/ext/ext.cpp index 5ca62592..e47485e0 100644 --- a/dscribe/ext/ext.cpp +++ b/dscribe/ext/ext.cpp @@ -79,8 +79,9 @@ PYBIND11_MODULE(ext, m) { py::arg("derivatives"), // py::array_t py::arg("descriptor"), // py::array_t py::arg("atomic_numbers"), // py::array_t + py::arg("cell"), // py::array_t + py::arg("pbc"), // py::array_t py::arg("atomic_positions"), // py::array_t - py::arg("cell_list"), // CellList py::arg("desc_centers"), // py::array_t py::arg("grad_centers"), // py::array_t py::arg("return_descriptor") // const bool From c4f11e7a2ba3109b738aa4dee347449b19783c2c Mon Sep 17 00:00:00 2001 From: Ehsan Rahmatizad Date: Fri, 25 Jul 2025 11:13:16 +0200 Subject: [PATCH 3/4] Incomplete Periodic BC in acsf.cpp --- dscribe/ext/acsf.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/dscribe/ext/acsf.cpp b/dscribe/ext/acsf.cpp index 49071350..0a48cbd7 100644 --- a/dscribe/ext/acsf.cpp +++ b/dscribe/ext/acsf.cpp @@ -276,6 +276,15 @@ void ACSF::derivatives_analytical( const bool return_descriptor ) { + ///INCOMPLETE_PERIODIC_SYSTEMS///auto pbc_u = pbc.unchecked<1>(); + ///INCOMPLETE_PERIODIC_SYSTEMS///bool is_periodic = this->periodic && (pbc_u(0) || pbc_u(1) || pbc_u(2)); + ///INCOMPLETE_PERIODIC_SYSTEMS///if (is_periodic) { + ///INCOMPLETE_PERIODIC_SYSTEMS/// ExtendedSystem system_extended = extend_system(atomic_positions, atomic_numbers, cell, pbc, this->cutoff); + ///INCOMPLETE_PERIODIC_SYSTEMS/// atomic_positions = system_extended.positions; + ///INCOMPLETE_PERIODIC_SYSTEMS/// atomic_numbers = system_extended.atomic_numbers; + ///INCOMPLETE_PERIODIC_SYSTEMS///} + + // Calculate neighbours with a cell list CellList cell_list(atomic_positions, this->cutoff); From 0a1ae9b6b814a0c43d304819d1bfd530f4756445 Mon Sep 17 00:00:00 2001 From: Ehsan Rahmatizad Date: Fri, 25 Jul 2025 14:17:24 +0200 Subject: [PATCH 4/4] Minor fix in ACSF G4 analytical derivative --- dscribe/ext/acsf.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/dscribe/ext/acsf.cpp b/dscribe/ext/acsf.cpp index 0a48cbd7..3d93583d 100644 --- a/dscribe/ext/acsf.cpp +++ b/dscribe/ext/acsf.cpp @@ -469,7 +469,8 @@ void ACSF::derivatives_analytical( // G4 Derivatives if (r_jk <= r_cut){ // Only compute if r_jk is within cutoff for G4 - for (auto params : g4_params){ // loop over G4 functions + double fc_jk = compute_cutoff(r_jk); + for (auto params : g4_params) { eta = params[0]; // eta is declared outside the loop now zeta = params[1]; // zeta is declared outside the loop now lambda = params[2]; // lambda is declared outside the loop now @@ -488,26 +489,29 @@ void ACSF::derivatives_analytical( if(j == i){ // Derivative w.r.t. atom i dcostheta = 0.5*(r_ik*(r_ij_square-r_ik_square+r_jk_square)*e_ij[c] - r_ij*(r_ij_square-r_ik_square-r_jk_square)*e_ik[c])/(r_ij_square*r_ik_square); secnd_term = -2.0*eta * (r_ij*e_ij[c] + r_ik*e_ik[c])*ang*fc4; - third_term = (e_ij[c]*d_fc_ij*fc_ik*compute_cutoff(r_jk) + e_ik[c]*d_fc_ik*fc_ij*compute_cutoff(r_jk)) * ang + e_jk[c]*d_fc_jk*fc_ij*fc_ik * ang; + third_term = e_ij[c]*d_fc_ij*fc_ik*fc_jk + e_ik[c]*d_fc_ik*fc_ij*fc_jk; // + e_jk[c]*d_fc_jk*fc_ij*fc_ik ; } else if(j == x){ // Derivative w.r.t. atom x dcostheta = -0.5*(e_ik[c]*r_ij*r_ik - e_ij[c]*r_ik_square + e_jk[c]*r_ij*r_jk + e_ij[c]*r_jk_square)/(r_ij_square * r_ik); secnd_term = 2.0*eta * (r_ij*e_ij[c] - r_jk*e_jk[c])*ang*fc4; - third_term = (-e_ij[c]*d_fc_ij*fc_ik*compute_cutoff(r_jk) + e_jk[c]*d_fc_jk*fc_ij*fc_ik) * ang; + third_term = -e_ij[c]*d_fc_ij*fc_ik*fc_jk + e_jk[c]*d_fc_jk*fc_ij*fc_ik; } else if(j == y){ // Derivative w.r.t. atom y dcostheta = 0.5*(e_ik[c]*r_ij_square - r_ik*(e_ij[c]*r_ij - e_jk[c]*r_jk) - e_ik[c]*r_jk_square)/(r_ij*r_ik_square); secnd_term = 2.0*eta * (r_ik*e_ik[c] + r_jk*e_jk[c])*ang*fc4; - third_term = (-e_ik[c]*d_fc_ik*fc_ij*compute_cutoff(r_jk) - e_jk[c]*d_fc_jk*fc_ij*fc_ik) * ang; + third_term = -e_ik[c]*d_fc_ik*fc_ij*fc_jk - e_jk[c]*d_fc_jk*fc_ij*fc_ik; } first_term = dang*dcostheta*fc4; - final = (first_term + secnd_term + third_term) * ef; + third_term *= ang; + + final = first_term + secnd_term + third_term; + final *= ef; derivatives_mu(idxi,idxj,c,o0_g4_g5) += final; + } } + o0_g4_g5++; // Move to next G4 feature } - o0_g4_g5++; // Move to next G4 feature - } } else { o0_g4_g5 += g4_params.size(); // Skip G4 features if r_jk is outside cutoff }