To initialize a first swarm particle at the very first step of the PSO process, bonded parameters (BP) are estimated via BI to obtain initial values for bond, angle and dihedral equilibrium force constants. Bond lengths, values and dihedral torsions are set to their average values calculated from the AA-mapped MD trajectory. When using execution mode 2, force constants guessed via BI are applied together with angles and dihedrals equilibrium values defined in the preliminary CG model file. Accordingly, execution mode 2 produces more flexible CG models with respect to the AA-mapped trajectories (i.e. CG distributions are typically broader than the AA ones with large Rg fluctuations) and models are corrected during optimization to reach appropriate average Rg value. Potential functions used for BI and optimization are all user-specified via the preliminary model file. Histograms are generated using 50 bins on the domains of bonds and angles, and 100 bins for dihedrals domains. For bonds and angles, three steps of exponential smoothing are applied using an alpha parameter of 0.55 and a window of 5 bins, while no transformation is performed for dihedrals. These histograms are then transformed using Boltzmann inversion as 𝐵𝐼(𝑣) = -𝐾 𝐵 𝑇 log(𝑣), with 𝐾 𝐵 the Boltzmann constant, 𝑇 the temperature and 𝑣 as bins values. GROMACS potential functions are then fit using SciPy 3 to obtain guesses of the force constants. To ensure these initial guesses of bonded parameters will produce a stable enough model to initialize the optimization procedure, force constants are capped to 17000 kJ.mol -1 for bonds using function 1, 1200 kJ.mol -1 for angles function 2, 3.5 kJ.mol -1 for dihedrals functions 1, 4 and 9 and 250 kJ.mol -1 .rad -2 for dihedrals function 2 (see GROMACS manual 4 sections 4.8 "Bonded interactions" and 5.2 "File formats").