diff --git a/dscribe/descriptors/acsf.py b/dscribe/descriptors/acsf.py
index 1eda8f80..87aed1f0 100644
--- a/dscribe/descriptors/acsf.py
+++ b/dscribe/descriptors/acsf.py
@@ -221,8 +221,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 +405,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
+        
+        # Call the C++ function with the correct arguments and order
+        self.acsf_wrapper.derivatives_analytical(
+            d,                 # 1st arg: py::array_t<double> derivatives
+            c,                 # 2nd arg: py::array_t<double> descriptor
+            atomic_numbers,    # 3rd arg: py::array_t<int> desc_centers
+            cell,              # 4th arg: py::array_t<double> desc_centers
+            pbc,               # 5th arg: py::array_t<int> desc_centers
+            positions,         # 6th arg: py::array_t<double> atomic_positions
+            centers,           # 7th arg: py::array_t<int> desc_centers
+            centers,           # 8th arg: py::array_t<int> 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..3d93583d 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,300 @@ inline void ACSF::compute_g5(py::detail::unchecked_mutable_reference<double, 2>
 		offset++;
 	}
 }
+
+
+void ACSF::derivatives_analytical(
+    py::array_t<double> derivatives,
+    py::array_t<double> descriptor,
+    py::array_t<int> atomic_numbers,
+    py::array_t<double> cell,
+    py::array_t<int> pbc,
+    py::array_t<double> atomic_positions,
+    py::array_t<int> desc_centers,
+    py::array_t<int> grad_centers, // we want derivatives w.r.t. these atoms
+    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);
+    
+    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<n_desc_centers; idxi++) {
+        
+        int i = desc_centers_u(idxi);
+        // Get neighbours and distances for atom i using CellList
+        CellListResult neighbours_i = cell_list.getNeighboursForIndex(i);
+        int n_neighbours = neighbours_i.indices.size();
+
+        // Loop over neighbours of atom i (j)
+        for (int j_neighbour_idx = 0; j_neighbour_idx < n_neighbours; ++j_neighbour_idx) {
+            int x = neighbours_i.indices[j_neighbour_idx]; // x is the j atom
+            if(i == x) continue; // Skip self-interaction
+
+            // Precompute some values for G1, G2, G3
+            double r_ij = neighbours_i.distances[j_neighbour_idx]; // Corrected access
+            double r_ij_square = neighbours_i.distancesSquared[j_neighbour_idx]; // Add this for G4/G5 later
+            double fc_ij = compute_cutoff(r_ij);
+            int index_j = atomic_number_to_index_map[atomic_numbers_u[x]];
+            // Initial offset for G1, G2, G3 for the current atomic species (index_j)
+            int offset = index_j * (1+n_g2+n_g3);  
+            int o0_g1_g3 = offset; // Store initial offset for derivatives calculation
+
+            // Compute unit vector e_ij = (pos_i - pos_j) / r_ij
+            double e_ij[3];
+            for(int c=0; c<3; c++)
+                e_ij[c] = (atomic_positions_u(i,c) - atomic_positions_u(x,c)) / r_ij;
+            
+            // If return_descriptor is true, compute and add descriptor values
+            if(return_descriptor){
+                // Note: The compute_gX functions in acsf_new.cpp take out_mu, index, offset.
+                // Here, out_mu is descriptor_mu, index is idxi (the index in the descriptor array),
+                // and offset is the current feature offset.
+                int current_offset_for_descriptor = offset; // Use a temporary offset for descriptor calculation
+                compute_g1(descriptor_mu, idxi, current_offset_for_descriptor, fc_ij);
+                compute_g2(descriptor_mu, idxi, current_offset_for_descriptor, r_ij, fc_ij);
+                compute_g3(descriptor_mu, idxi, current_offset_for_descriptor, r_ij, fc_ij);
+            }
+
+            // Derivative of cutoff function w.r.t. r_ij
+            double d_fc_ij = -(PI*sin((PI*r_ij)/r_cut))/(2.0*r_cut);
+
+            // G1 Derivatives
+            for(int idxj=0; idxj<grad_centers.shape(0); idxj++){
+                int j = grad_centers_u(idxj);
+                if(j != i && j != x) continue; // Only i and x contribute to G1 derivative
+                if(j == i){
+                    for(int c=0;c<3;c++)
+                        derivatives_mu(idxi,idxj,c,o0_g1_g3) += e_ij[c]*d_fc_ij;
+                }
+                else if(j == x){
+                    for(int c=0;c<3;c++)
+                        derivatives_mu(idxi,idxj,c,o0_g1_g3) -= e_ij[c]*d_fc_ij;
+                }
+            }
+            o0_g1_g3++; // Move to next feature for G2
+
+            // G2 Derivatives
+            double Rs, ef, der, tmp;
+            for (auto params : g2_params) {
+                eta = params[0]; // eta is declared outside the loop now
+                Rs = params[1];
+                tmp = (r_ij - Rs);
+                ef = exp(-eta * tmp*tmp);
+                der = ef*(-eta*2*tmp*fc_ij + d_fc_ij); // Derivative w.r.t r_ij
+
+                for(int idxj=0; idxj<grad_centers.shape(0); idxj++){
+                    int j = grad_centers_u(idxj);
+                    if(j != i && j != x) continue;
+                    if(j == i){
+                        for(int c=0;c<3;c++)
+                            derivatives_mu(idxi,idxj,c,o0_g1_g3) += e_ij[c]*der;
+                    }
+                    else if(j==x){
+                        for(int c=0;c<3;c++)
+                            derivatives_mu(idxi,idxj,c,o0_g1_g3) -= e_ij[c]*der;
+                    }
+                }
+                o0_g1_g3++; // Move to next feature for G3
+            }
+
+            // G3 Derivatives
+            for (auto param : g3_params) {
+                der = param * (-sin(r_ij*param)*fc_ij) + cos(r_ij*param)*d_fc_ij; // Derivative w.r.t r_ij
+
+                for(int idxj=0; idxj<grad_centers.shape(0); idxj++){
+                    int j = grad_centers_u(idxj);
+                    if(j != i && j != x) continue;
+                    if(j == i){
+                        for(int c=0;c<3;c++)
+                            derivatives_mu(idxi,idxj,c,o0_g1_g3) += e_ij[c]*der;
+                    }
+                    else if(j==x){
+                        for(int c=0;c<3;c++)
+                            derivatives_mu(idxi,idxj,c,o0_g1_g3) -= e_ij[c]*der;
+                    }
+                }
+                o0_g1_g3++; // Move to next feature (if any, for G4/G5)
+            }
+            
+            // If no G4 or G5 parameters, continue to next j_neighbour
+            if (g4_params.size()==0 && g5_params.size()==0) continue;
+            
+            // Loop over second neighbours of atom i (k or y) for 3-body terms (G4, G5)
+            for (int k_neighbour_idx = 0; k_neighbour_idx < n_neighbours; ++k_neighbour_idx) {
+                int y = neighbours_i.indices[k_neighbour_idx]; // y is the k atom
+                if (y == i || k_neighbour_idx >= 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
+                    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
+
+                        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<grad_centers.shape(0); idxj++){
+                            int j = grad_centers_u(idxj);
+                            // Only i, x, y contribute to the gradient of G4/G5 for this triplet
+                            if(j!=i && j!=x && j!=y) continue; 
+
+                            for(int c=0;c<3;c++) {
+                                // Calculate dcostheta and other terms based on which atom (i, x, or y) is the gradient center
+                                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*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*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*fc_jk - e_jk[c]*d_fc_jk*fc_ij*fc_ik;
+                                }
+
+                                first_term = dang*dcostheta*fc4;
+                                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
+                        }
+                } else {
+                    o0_g4_g5 += g4_params.size(); // Skip G4 features if r_jk is outside cutoff
+                }
+
+                // G5 Derivatives
+                for (auto params : g5_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
+
+                    ef = exp(-eta*(r_ij_square+r_ik_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<grad_centers.shape(0); idxj++){
+                        int j = grad_centers_u(idxj);
+                        if(j!=i && j!=x && j!=y) continue; 
+
+                        for(int c=0;c<3;c++) {
+                            
+                            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*fc5;
+                                third_term = (e_ij[c]*d_fc_ij*fc_ik + e_ik[c]*d_fc_ik*fc_ij) * ang;
+
+                            } 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])*ang*fc5;
+                                third_term = (-e_ij[c]*d_fc_ij*fc_ik) * ang;
+
+                            } 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])*ang*fc5;
+                                third_term = (-e_ik[c]*d_fc_ik*fc_ij) * ang;
+                            }
+
+                            first_term = dang*dcostheta*fc5;
+                            final = (first_term + secnd_term + third_term) * ef;
+                            derivatives_mu(idxi,idxj,c,o0_g4_g5) += final;
+                        }
+                    }
+                    o0_g4_g5++; // Move to next G5 feature
+                }
+            }
+        }
+    }
+}
diff --git a/dscribe/ext/acsf.h b/dscribe/ext/acsf.h
index fd6e0e3b..170b1153 100644
--- a/dscribe/ext/acsf.h
+++ b/dscribe/ext/acsf.h
@@ -4,6 +4,8 @@
 #include <unordered_map>
 #include <vector>
 #include "descriptorlocal.h"
+#include "celllist.h"
+#include "geometry.h"
 
 #define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062
 
@@ -32,7 +34,7 @@ class ACSF : public DescriptorLocal {
             py::array_t<int> centers,
             CellList cell_list
         );
-
+        
         void create(
             py::array_t<double> out, 
             py::array_t<double> positions,
@@ -41,6 +43,18 @@ class ACSF : public DescriptorLocal {
             CellList cell_list
         );
 
+        void derivatives_analytical(
+            py::array_t<double> derivatives,
+            py::array_t<double> descriptor,
+            py::array_t<int> atomic_numbers,
+            py::array_t<double> cell,
+            py::array_t<int> pbc,
+            py::array_t<double> atomic_positions,
+            py::array_t<int> desc_centers,
+            py::array_t<int> 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..e47485e0 100644
--- a/dscribe/ext/ext.cpp
+++ b/dscribe/ext/ext.cpp
@@ -75,6 +75,17 @@ PYBIND11_MODULE(ext, m) {
     py::class_<ACSF>(m, "ACSFWrapper")
         .def(py::init<double, vector<vector<double> > , vector<double> , vector<vector<double> > , vector<vector<double> >, vector<int> , bool>())
         .def("create", overload_cast_<py::array_t<double>, py::array_t<double>, py::array_t<int>, py::array_t<double>, py::array_t<bool>, py::array_t<int> >()(&DescriptorLocal::create))
+        .def("derivatives_analytical", &ACSF::derivatives_analytical,
+             py::arg("derivatives"),          // py::array_t<double>
+             py::arg("descriptor"),           // py::array_t<double>
+             py::arg("atomic_numbers"),       // py::array_t<int>
+             py::arg("cell"),                 // py::array_t<double>
+             py::arg("pbc"),                  // py::array_t<int>
+             py::arg("atomic_positions"),     // py::array_t<double>
+             py::arg("desc_centers"),         // py::array_t<int>
+             py::arg("grad_centers"),         // py::array_t<int>
+             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,
     )
 ]