Simulation setup. Molecular Dynamics (MD) simulations were performed in explicit solvent. Water has been represented with the TIP3P model and ions with the OPLS-AA force field . Our analysis is carried out, apart from water, in N,N-dimethylformamide (DMF) using force field parameters compatible with OPLS-AA 60 obtained from the virtualchemistry.org database . The MOF half-SBUs were modelled using force field parameters obtained from a previous work . Force field parameters are also included as Supporting Information (SI), in a text file (SI_FF.txt). A brief discussion of the force field is provided in the SI, section I, p. 2. The leapfrog integrator was used to propagate dynamics of the system with a time step of 2 fs. The LINCS algorithm preserved the distances of bonds involving hydrogen atoms. The cut-off for non-bonded interactions is 10 Å. Long range electrostatics were treated using the Particle-Mesh Ewald (PME) scheme. The Bussi-Donadio-Parrinello thermostat and the Berendsen barostat 66 preserved the temperature and pressure at 493 K and 3,500 bar respectively for an equilibration period of 5 ns. Production molecular dynamics simulations followed using the Parrinello-Rahman Barostat 67 with a relaxation time of 2 ps in water and 4 ps in DMF. Production simulations were carried out for 100 ns. The system consists of 132 half-SBUs in explicit solvent. The SBUs are introduced in varying proportions of three Metal-Ligand (ML) isomers each with different orientations of the ligands about the Cr3O3 core denoted as, MLA, MLB, and MLC . These isomers can be depicted in Figure . Systems consist of more than 780,000 atoms in water and 500,000 atoms in DMF randomly dispersed in a cubic box with an edge of 200 Å. Two initial compositions of half-SBUs are considered. In the first one (purely MLA), there are 132 MLA half-SBUs in solution. The second one (equiprobable MLA, MLB, and MLC), consists of 44 MLA, 44 MLB and 44 MLC in solution. Simulations were performed with periodic boundary conditions in three dimensions, using GROMACS 2018 . Chemical structures were visualized with VMD 69 . Distance and adjacency matrices. At first, a (n x n) distance matrix is constructed, where n is the number of half-SBUs. The generic 𝑖𝑗 element of the matrix corresponds to the Euclidean distance between the central oxygen atoms of the i th and j th half-SBUs units. This distance was used as the argument of a step function to define the adjacency of two half-SBUs. A cutoff distance of 15 Å is chosen as this value lies between the first and second coordination shells emerging from the calculation of the pair radial distribution function between central oxygen atoms of half-SBUs. This is provided in the SI, section II, p. 3, Figure . In the adjacency matrix, the element 𝑖𝑗 is equal to the unity if the distance between 𝑖 and 𝑗 is below the cutoff otherwise it is set to zero. Examples of distance and adjacency matrices are discussed in the SI, section II, p. 3. Graph-based clustering. The clustering analysis of half-SBUs is performed by converting molecular structures into lower dimension graph representations. We consider the central oxygen atoms of the building units as the nodes of the graph. Then we connect these nodes advised by the adjacency matrix discussed in the previous paragraph. Then, we identify clusters as connected components in the graph, where the smallest cluster consists of two half-SBUs (dimer). Departing form depth first search (DFS) that can be used to identify clusters, we analyze the properties of clusters as components of a graph; hence enrich the information we have for each cluster. Also, we calculate the mass of each cluster as the sum of the masses of its constituent particles. Furthermore, the local environment of each half-SBU is characterized by enumerating the neighbors of each node. This is a measure of the coordination of half-SBUs in the cluster. Also, we calculate the number of "free" half-SBUs as isolated nodes, without any edges attached. This allows us to further calculate certain properties of the graph, such as the number of connections (average neighborhood degree of each node, named graph degree for simplicity), the extent of interconnected triplets (transitivity) and the number of nodes which connect with similar ones based on their degree (assortativity coefficient). The graph representation was constructed using the NetworkX Python library . Spherical radius and radius of gyration. The spherical radius and the radius of gyration are calculated for each cluster. The first follows from the ideal process of including each cluster in a sphere. The radius of that sphere would then be half the maximum distance between any two metal centers in this cluster. The radius of gyration provides insight into the distribution of mass in complex structures and it is used to calculate the fractal dimension of each cluster. Periodic boundary conditions in three dimensions are appropriately considered in all these calculations. Fractal Dimension. The fractal dimension D ! was computed using a power law approach in the form,