diff --git a/mlcolvar/core/transform/continuous_hist.py b/mlcolvar/core/transform/continuous_hist.py
deleted file mode 100644
index 4977825e..00000000
--- a/mlcolvar/core/transform/continuous_hist.py
+++ /dev/null
@@ -1,53 +0,0 @@
-import torch
-
-from mlcolvar.core.transform import Transform
-from mlcolvar.core.transform.utils import easy_KDE
-
-from typing import Union
-
-__all__ = ["ContHist"]
-
-class ContHist(Transform):
-    """
-    Compute continuous histogram with KDE-like method
-    """
-
-    def __init__(self,
-                 in_features : int,
-                 min : float,
-                 max : float,
-                 bins : int,
-                 sigma_to_center : float) -> torch.Tensor :
-       
-        super().__init__(in_features=in_features, out_features=bins)
-
-        self.in_features = in_features
-        self.min = min
-        self.max = max
-        self.bins = bins
-        self.sigma_to_center = sigma_to_center
-    
-    def compute_hist(self, x):
-        hist = easy_KDE(x=x,
-                        n_input=self.in_features, 
-                        min_max=[self.min, self.max], 
-                        n=self.bins, 
-                        sigma_to_center=self.sigma_to_center)
-        # if self.normalize = 
-        return hist
-
-    def forward(self, x: torch.Tensor):
-        x = self.compute_hist(x)
-        return x
-    
-def test_continuous_histogram():
-    x = torch.randn((5,100))
-    hist = ContHist(in_features=100,
-                    min=-1,
-                    max=1,
-                    bins=10,
-                    sigma_to_center=1)
-    out = hist(x)
-    
-if __name__ == "__main__":
-    test_continuous_histogram()
\ No newline at end of file
diff --git a/mlcolvar/core/transform/eigs_adjacency_matrix.py b/mlcolvar/core/transform/eigs_adjacency_matrix.py
deleted file mode 100644
index 5dfc1b76..00000000
--- a/mlcolvar/core/transform/eigs_adjacency_matrix.py
+++ /dev/null
@@ -1,115 +0,0 @@
-import torch
-
-from mlcolvar.core.transform import Transform
-from mlcolvar.core.transform.utils import compute_adjacency_matrix
-
-from typing import Union
-
-__all__ = ["EigsAdjMat"]
-
-class EigsAdjMat(Transform):
-    """
-    Eigenvalues of adjacency matrix transform, compute the eigenvalues of the adjacency matrix for a set of atoms from their positions
-    """
-
-    def __init__(self,
-                 mode : str,
-                 cutoff : float, 
-                 n_atoms : int,
-                 PBC: bool,
-                 real_cell: Union[float, list],
-                 scaled_coords : bool,
-                 switching_function = None) -> torch.Tensor:
-        """Initialize an eigenvalues of an adjacency matrix object.
-
-        Parameters
-        ----------
-        mode : str
-            Mode for cutoff application, either:
-            - 'continuous': applies a switching function to the distances which can be specified with switching_function keyword, has stable derivatives
-            - 'discontinuous': set at zero everything above the cutoff and one below, derivatives may be be incorrect
-        cutoff : float
-            Cutoff for the adjacency criterion 
-        n_atoms : int
-            Number of atoms in the system
-        PBC : bool
-            Switch for Periodic Boundary Conditions use
-        real_cell : Union[float, list]
-            Dimensions of the real cell, orthorombic-like cells only
-        scaled_coords : bool
-            Switch for coordinates scaled on cell's vectors use
-        switching_function : _type_, optional
-            Switching function to be applied for the cutoff, can be either initialized as a switching_functions/SwitchingFunctions class or a simple function, by default None
-
-        Returns
-        -------
-        torch.Tensor
-            Adjacency matrix of all the n_atoms according to cutoff
-        """
-        super().__init__(in_features=int(n_atoms*3), out_features=n_atoms)
-
-        # parse args
-        self.mode = mode
-        self.cutoff = cutoff 
-        self.n_atoms = n_atoms
-        self.PBC = PBC
-        self.real_cell = real_cell
-        self.scaled_coords = scaled_coords
-        self.switching_function = switching_function
-
-    def compute_adjacency_matrix(self, x):
-        x = compute_adjacency_matrix(pos=x,
-                                     mode = self.mode,
-                                     cutoff = self.cutoff, 
-                                     n_atoms = self.n_atoms,
-                                     PBC = self.PBC,
-                                     real_cell = self.real_cell,
-                                     scaled_coords = self.scaled_coords,
-                                     switching_function=self.switching_function)
-        # x = torch.squeeze(x)
-        return x
-    
-    def get_eigenvalues(self, x):
-        eigs = torch.linalg.eigvalsh(x)
-        return eigs
-
-    def forward(self, x: torch.Tensor):
-        x = self.compute_adjacency_matrix(x)
-        eigs = self.get_eigenvalues(x)
-        return eigs
-
-def test_eigs_of_adj_matrix():
-    from mlcolvar.core.transform.switching_functions import SwitchingFunctions
-    
-    n_atoms=2
-    pos = torch.Tensor([ [ [0., 0., 0.],
-                           [1., 1., 1.] ],
-                         [ [0., 0., 0.],
-                           [1., 1.1, 1.] ] ]
-                      )
-    real_cell = torch.Tensor([1., 2., 1.])
-    cutoff = 1.8
-    switching_function=SwitchingFunctions(in_features=n_atoms*3, name='Fermi', cutoff=cutoff, options={'q':0.05})
-  
-    model = EigsAdjMat(mode = 'continuous',
-                       cutoff = cutoff, 
-                       n_atoms = n_atoms,
-                       PBC = True,
-                       real_cell = real_cell,
-                       scaled_coords = False,
-                       switching_function=switching_function)
-    out = model(pos)
-
-    pos = torch.einsum('bij,j->bij', pos, 1/real_cell)
-    model = EigsAdjMat(mode = 'continuous',
-                       cutoff = cutoff, 
-                       n_atoms = n_atoms,
-                       PBC = True,
-                       real_cell = real_cell,
-                       scaled_coords = True,
-                       switching_function=switching_function)
-    out = model(pos)
-    assert(out.shape[-1] == model.out_features)
-
-if __name__ == "__main__":
-    test_eigs_of_adj_matrix()
\ No newline at end of file
diff --git a/mlcolvar/core/transform/pairwise_distances.py b/mlcolvar/core/transform/pairwise_distances.py
deleted file mode 100644
index 345239ed..00000000
--- a/mlcolvar/core/transform/pairwise_distances.py
+++ /dev/null
@@ -1,86 +0,0 @@
-import torch
-
-from mlcolvar.core.transform import Transform
-from mlcolvar.core.transform.utils import compute_distances_matrix
-
-from typing import Union
-
-__all__ = ["PairwiseDistances"]
-
-class PairwiseDistances(Transform):
-    """
-    Pairwise distances transform, compute all the non duplicated pairwise distances for a set of atoms from their positions
-    """
-
-    def __init__(self, 
-                 n_atoms : int,
-                 PBC: bool,
-                 real_cell: Union[float, list],
-                 scaled_coords : bool) -> torch.Tensor:
-        """Initialize a pairwise distances matrix object.
-
-        Parameters
-        ----------
-        n_atoms : int
-            Number of atoms in the system
-        PBC : bool
-            Switch for Periodic Boundary Conditions use
-        real_cell : Union[float, list]
-            Dimensions of the real cell, orthorombic-like cells only
-        scaled_coords : bool
-            Switch for coordinates scaled on cell's vectors use
-
-        Returns
-        -------
-        torch.Tensor
-            Non duplicated pairwise distances between all the atoms
-        """
-        super().__init__(in_features=int(n_atoms*3), out_features=int(n_atoms*(n_atoms-1) / 2))
-
-        # parse args
-        self.n_atoms = n_atoms
-        self.PBC = PBC
-        self.real_cell = real_cell
-        self.scaled_coords = scaled_coords
-
-    def compute_pairwise_distances(self, pos):
-        dist = compute_distances_matrix(pos=pos,
-                                        n_atoms=self.n_atoms,
-                                        PBC=self.PBC,
-                                        real_cell=self.real_cell,
-                                        scaled_coords=self.scaled_coords)
-        batch_size = dist.shape[0]
-        device = pos.device
-        # mask out diagonal elements
-        aux_mask = torch.ones_like(dist, device=device) - torch.eye(dist.shape[-1], device=device)
-        # keep upper triangular part to avoid duplicates
-        unique = aux_mask.triu().nonzero(as_tuple=True)
-        pairwise_distances = dist[unique].reshape((batch_size, -1)) 
-        return pairwise_distances
-
-    def forward(self, x: torch.Tensor):
-        x = self.compute_pairwise_distances(x)
-        return x
-
-def test_pairwise_distances():
-
-    pos = torch.Tensor([ [ [0., 0., 0.],
-                           [1., 1., 1.],
-                           [1., 1., 1.1] ],
-                         [ [0., 0., 0.],
-                           [1., 1.1, 1.],
-                           [1., 1., 1.] ] ]
-                      )
-    
-    real_cell = torch.Tensor([1., 2., 1.])
-    cutoff = 1.8
-  
-    model = PairwiseDistances(n_atoms = 3,
-                              PBC = True,
-                              real_cell = real_cell,
-                              scaled_coords = False)
-    out = model(pos)
-    assert(out.reshape(pos.shape[0], -1).shape[-1] == model.out_features)
-
-if __name__ == "__main__":
-    test_pairwise_distances()
\ No newline at end of file
diff --git a/mlcolvar/core/transform/radial_distribution_function.py b/mlcolvar/core/transform/radial_distribution_function.py
deleted file mode 100644
index ed6b7c4f..00000000
--- a/mlcolvar/core/transform/radial_distribution_function.py
+++ /dev/null
@@ -1,116 +0,0 @@
-import torch
-
-from mlcolvar.core.transform import Transform,PairwiseDistances
-from mlcolvar.core.transform.utils import easy_KDE
-
-from typing import Union
-
-__all__ = ["RDF"]
-
-class RDF(Transform):
-    """
-    Compute radial distribution function. Optionally returns the integral over some specific peaks.
-    """
-
-    def __init__(self,
-                 n_atoms : int,
-                 PBC : bool,
-                 real_cell : torch.Tensor,
-                 scaled_coords : bool,
-                 min_max : list,
-                 n_bins : int,
-                 sigma_to_center : float,
-                 peak_integral : list = None,
-                 integral_baseline : bool = False) -> torch.Tensor :
-        # check output dimension
-        if peak_integral:
-            super().__init__(in_features=int(n_atoms*(n_atoms-1) / 2), out_features=1)
-        else:
-            super().__init__(in_features=int(n_atoms*(n_atoms-1) / 2), out_features=n_bins)
-        
-        # parse arguments
-        self.in_features = int(n_atoms*(n_atoms-1) / 2)
-        self.n_atoms = n_atoms
-        self.PBC = PBC
-        self.scaled_coords = scaled_coords
-
-        if real_cell.shape[0] != 1:
-           if real_cell[0] == real_cell[1] and real_cell[1] == real_cell[2]:
-                pass
-           else:
-                raise ValueError('Radial distribution function implemented only for cubic cells!') 
-        self.real_cell = real_cell
-
-        self.min_max = min_max
-        self.n_bins = n_bins
-        self.sigma_to_center = sigma_to_center
-
-        self.peak_integral = peak_integral
-        self.integral_baseline = integral_baseline
-
-        # initialize pairwise distances transform
-        self.pairDist = PairwiseDistances(n_atoms=self.n_atoms,
-                                 PBC=self.PBC,
-                                 real_cell = self.real_cell,
-                                 scaled_coords=self.scaled_coords)
-        
-    def compute_distances(self, x):
-        dist = self.pairDist(x)
-        return dist
-    
-    def compute_hist(self, x):
-        hist, bins = easy_KDE(x=x,
-                              n_input=self.in_features, 
-                              min_max=self.min_max, 
-                              n=self.n_bins, 
-                              sigma_to_center=self.sigma_to_center,
-                              return_bins=True)
-        
-        # store the bins, shell and integration attributes if not present
-        if not hasattr(self, 'bins'):
-            self.bins = bins
-            self.bins_size = bins[1] - bins[0]
-
-            bins_ext = torch.linspace(self.min_max[0], self.min_max[1] + self.bins_size.item(), self.n_bins + 1, device=x.device)
-            self.shell = 4/3 * torch.pi * (self.n_atoms / self.real_cell**3) * (bins_ext[1:]**3 - bins_ext[:-1]**3)
-
-            if self.peak_integral is not None:
-                self.peak_lower_bound =  torch.where(self.bins == torch.max(self.bins[(self.bins < self.peak_integral[0])]))[0].item()
-                self.peak_upper_bound =  torch.where(self.bins == torch.min(self.bins[(self.bins > self.peak_integral[1])]))[0].item()
-
-        return hist
-    
-    def compute_rdf(self, x):
-        x = self.compute_distances(x)
-        x = self.compute_hist(x)
-        x = x / (self.shell*self.n_atoms)
-        return x
-        
-
-    def forward(self, x: torch.Tensor):
-        x = self.compute_rdf(x)
-        if self.peak_integral is not None:
-            if self.integral_baseline:
-                integral = torch.sum((x[:, self.peak_lower_bound:self.peak_upper_bound] - 1) * self.bins_size, dim = 1, keepdim=True)
-            else:
-                integral = torch.sum((x[:, self.peak_lower_bound:self.peak_upper_bound]) * self.bins_size, dim = 1, keepdim=True)
-            return integral
-        else:
-            return x
-
-def test_radial_distribution_function():
-    x = torch.randn((5,300))
-    cell = 1.0
-    rdf = RDF(n_atoms=100,
-              PBC=True,
-              real_cell=cell,
-              scaled_coords=False,
-              min_max=[0,1],
-              n_bins=10,
-              sigma_to_center=1,
-              peak_integral=[0.2,0.8],
-              integral_baseline = True)
-    out = rdf(x)
-    
-if __name__ == "__main__":
-    test_radial_distribution_function()
\ No newline at end of file
diff --git a/mlcolvar/core/transform/radius_graph.py b/mlcolvar/core/transform/radius_graph.py
deleted file mode 100644
index 3c42af7f..00000000
--- a/mlcolvar/core/transform/radius_graph.py
+++ /dev/null
@@ -1,128 +0,0 @@
-import torch
-
-from mlcolvar.core.transform import Transform
-from mlcolvar.core.transform.utils import compute_distances_matrix,apply_cutoff
-
-from typing import Union
-
-__all__ = ["RadiusGraph"]
-
-class RadiusGraph(Transform):
-    """
-    Radius Graph transform, compute the elements to build a distance-cutoff based graph for a set of atoms from their positions
-    """
-
-    def __init__(self,
-                 mode : str,
-                 cutoff : float, 
-                 n_atoms : int,
-                 PBC: bool,
-                 real_cell: Union[float, list],
-                 scaled_coords : bool,
-                 switching_function = None,
-                 zero_threshold : float = 1e-4) -> torch.Tensor:
-        """Initialize a radius graph object
-
-        Parameters
-        ----------
-        mode : str
-            Mode for cutoff application, either:
-            - 'continuous': applies a switching function to the distances which can be specified with switching_function keyword and keep only what is >= zero_threshold, has stable derivatives
-            - 'discontinuous': keep what is >= cutoff and discard the rest, derivatives may be be incorrect
-        cutoff : float
-            Radial cutoff for the connectivity criterion
-        n_atoms : int
-            Number of atoms in the system
-        PBC : bool
-            Switch for Periodic Boundary Conditions use
-        real_cell : Union[float, list]
-            Dimensions of the real cell, orthorombic-like cells only
-        scaled_coords : bool
-            Switch for coordinates scaled on cell's vectors use
-        switching_function : _type_, optional
-            Switching function to be applied for the cutoff, can be either initialized as a switching_functions/SwitchingFunctions class or a simple function, by default None
-        zero_threshold : float, optional
-            Threshold to be considered zero when in continuous mode, it is ignored if in discontinuous mode, by default 1e-4
-
-        Returns
-        -------
-        torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor
-            Elements to build a distance-cutoff based graph:
-            - Distances for each edge
-            - Batch indeces
-            - Edge source indeces
-            - Edge destination indeces
-
-
-        """
-        super().__init__(in_features=int(n_atoms*3), out_features=None)
-
-        # parse args
-        self.mode = mode
-        self.cutoff = cutoff 
-        self.n_atoms = n_atoms
-        self.PBC = PBC
-        self.real_cell = real_cell
-        self.scaled_coords = scaled_coords
-        self.switching_function = switching_function
-        self.zero_threshold = zero_threshold
-        
-    def compute_radius_graph(self, pos, mode):
-        dist = compute_distances_matrix(pos=pos,
-                                        n_atoms=self.n_atoms,
-                                        PBC=self.PBC,
-                                        real_cell=self.real_cell,
-                                        scaled_coords=self.scaled_coords)
-        
-        # we want to smooth dist at cutoff, we first get a switch and then apply that to dist
-        # we need a clone! Otherwise dist would be modified!
-        aux_switch = torch.clone(dist)
-        aux_switch = apply_cutoff(x=aux_switch,
-                                  cutoff=self.cutoff, 
-                                  mode=mode, 
-                                  switching_function = self.switching_function)          
-
-        if mode == 'continuous':
-            # smooth dist
-            dist = dist * aux_switch 
-            # discard what is almost zero --> use self.zero_threshold
-            unique = torch.nonzero(torch.ge(dist.triu(), self.zero_threshold), as_tuple=True)
-
-        elif mode == 'discontinuous': 
-            # discard zeros entries
-            unique = dist.triu().nonzero(as_tuple=True) 
-        
-        distances = dist[unique]
-        batch_indeces, edge_src, edge_dst = unique
-        
-        return distances,batch_indeces,edge_src,edge_dst
-
-    def forward(self, x : torch.Tensor):
-        x = self.compute_radius_graph(x, mode=self.mode)
-        return x
-
-def test_radiusgraph():
-    from mlcolvar.core.transform.switching_functions import SwitchingFunctions
-
-    n_atoms=3
-    pos = torch.Tensor([ [ [0., 0., 0.],
-                           [1., 1., 1.] ],
-                         [ [0., 0., 0.],
-                           [1., 1.1, 1.] ] ]
-                      )
-    
-    real_cell = torch.Tensor([1., 2., 1.])
-    cutoff = 1.8
-    switching_function=SwitchingFunctions(in_features=n_atoms*3, name='Fermi', cutoff=cutoff, options={'q':0.01})
-    model = RadiusGraph(mode = 'continuous',
-                        cutoff = cutoff, 
-                        n_atoms = 2,
-                        PBC = True,
-                        real_cell = real_cell,
-                        scaled_coords = False,
-                        switching_function=switching_function)
-    distances,batch_indeces,edge_src,edge_dst = model(pos)
-
-if __name__ == "__main__":
-    test_radiusgraph()
-              
\ No newline at end of file
diff --git a/mlcolvar/core/transform/switching_functions.py b/mlcolvar/core/transform/switching_functions.py
deleted file mode 100644
index 473b87bf..00000000
--- a/mlcolvar/core/transform/switching_functions.py
+++ /dev/null
@@ -1,80 +0,0 @@
-import torch
-
-from mlcolvar.core.transform import Transform
-
-
-__all__ = ["SwitchingFunctions"]
-
-class SwitchingFunctions(Transform):
-    """
-    Transform class with some common switching functions
-    """
-    SWITCH_FUNCS = ['Fermi', 'Rational']
-
-    def __init__(self,
-                 in_features : int,
-                 name : str, 
-                 cutoff : float, 
-                 options : dict = None):
-        f"""Initialize switching function object
-
-        Parameters
-        ----------
-        name : str
-            Name of the switching function to be used, available {",".join(self.SWITCH_FUNCS)}
-        cutoff : float
-            Cutoff for the swtiching functions
-        options : dict, optional
-            Dictionary with all the arguments of the switching function, by default None
-        """
-        super().__init__(in_features=in_features, out_features=in_features)
-
-        self.name = name
-        self.cutoff = cutoff
-        if options is None:
-            options = {}
-        self.options = options
-       
-        if name not in self.SWITCH_FUNCS:
-            raise NotImplementedError(f'''The switching function {name} is not implemented in this class. The available options are: {",".join(self.SWITCH_FUNCS)}.
-                    You can initialize it as a method of the SwitchingFunctions class and tell us on Github, contributions are welcome!''')  
-
-    def forward(self, x : torch.Tensor):
-        switch_function = self.__getattribute__(f'{self.name}_switch')
-        y = switch_function(x, self.cutoff, **self.options)
-        return y
-    
-    # ========================== define here switching functions ==========================
-    def Fermi_switch(self,
-                     x : torch.Tensor, 
-                     cutoff : float, 
-                     q : float = 0.01, 
-                     prefactor_cutoff : float = 1.0):
-        y = torch.div( 1, ( 1 + torch.exp( torch.div((x - prefactor_cutoff*cutoff) , q ))))
-        return y
-
-    def Rational_switch(self,
-                        x : torch.Tensor, 
-                        cutoff : float, 
-                        n : int = 6, 
-                        m : int = 12, 
-                        eps : float = 1e-8, 
-                        prefactor_cutoff : float = 1.0):
-        y = torch.div((1 - torch.pow(x/(prefactor_cutoff*cutoff), n) + eps) , (1 - torch.pow(x/(prefactor_cutoff*cutoff), m)  + 2*eps) )
-        return y
-
-
-def test_switchingfunctions():
-    x = torch.Tensor([1., 2., 3.])
-    cutoff = 2
-    switch = SwitchingFunctions(in_features=len(x), name='Fermi', cutoff=cutoff)
-    out = switch(x)
-
-    switch = SwitchingFunctions(in_features=len(x), name='Fermi', cutoff=cutoff, options = {'q' : 0.5})
-    out = switch(x)
-
-    switch = SwitchingFunctions(in_features=len(x), name='Rational', cutoff=cutoff, options = {'n' : 6, 'm' : 12})
-    out = switch(x)
-
-if __name__ == "__main__":
-    test_switchingfunctions()
\ No newline at end of file