Molecular Dynamics Simulations. The LAMMPS package was used for carrying out all MD simulations of this NP-polymer system. All simulations were carried out in the canonical (NVT) ensemble at a temperature of 𝜀/𝑘 * and a polymer density of 0.85 beads/𝜎 % . Under these conditions, the polymer matrix exhibits a melt-like state. Periodic boundary conditions were employed in all three directions of the simulation box. A velocity-Verlet algorithm with a time step of 0.002 (𝑚𝜎 ! 𝜀 ⁄ ) 4/! and a Nosé-Hoover thermostat of time constant 1.0 (𝑚𝜎 ! 𝜀 ⁄ ) 4/! was used for integrating the equations of motion and keeping the system temperature fixed. The simulations were initialized by placing the grafted NPs and the polymer chains in a cubic simulation box 50 times larger than the required dimensions to prevent overlap among the chain beads and NP cores. The box was then gradually shrunk in each direction over a period of 0.8 million timesteps until the targeted polymer density was reached. Simultaneously, the NPs were slowly brought into a configuration used for initiating the PMF calculations or NP assembly simulations. After initialization, we continued the MD simulations for an additional 1 million timesteps to further equilibrate the system while keeping the NPs fixed. During the production period, the net linear and angular momentums of the entire system were zeroed at every step to avoid the "unreal" drifting or rotation induced by the periodic boundary conditions.