The solution-phase conformations and dynamics of flexible polymers are generally well understood. The conformation of a single flexible polymer chain in solution shows predictable scaling with chain length, with a scaling exponent determined by the solvent quality, i.e. the relative strength of the polymer-polymer, solvent-solvent, and polymer-solvent interactions. In a good solvent, in which polymer-solvent attractions dominate, the chain is in an extended conformation, whereas a more compact, collapsed structure is formed under poor solvent conditions in order to minimise unfavourable interactions with the solvent or to maximise favourable intramolecular polymer interactions. Metrics such as the radius of gyration (R g ) are therefore expected to decrease as solvent quality decreases and chains become more compact. This behaviour has been observed in a number of simulation studies of sufficiently flexible single chains, with a transition from an extended coil to collapsed globule structure observed with decreasing solvent quality. However, semiconducting polymers typically have stiffer semiflexible backbones, which, coupled with the more anisotropic shape and interactions imparted by the conjugated backbone, means that the aggregate structure may deviate from that predicted for flexible chains. (See ref. 33 for a more comprehensive review of the behaviour of semiflexbile polymers in dilute solutions.) Indeed, for both single- and multi-chain systems, both the backbone stiffness and solvent quality have been shown to be important for predicting the types of structures formed by semiflexible polymers. For single chains (i.e. in very dilute solution), as chain stiffness increases, different structures, such as hairpins or toroids, take the place of disordered globules in poorer solvents, with the exact structure depending on both the solvent quality and chain flexibility. At higher concentrations where multiple chains are able to interact, bundles of fully overlapping chains, rather than collapsed globules, are expected for semiflexible polymers in a poor solvent due to the unfavourable bending energy. Indeed, Monte Carlo simulations have shown the equilibrium structure of multichain aggregates in a poor solvent to shift from a disordered globular morphology to fully overlapped twisted or folded rodlike bundles as the chain stiffness is increased. However, although these fully overlapped bundles are expected to be the equilibrium structure under these conditions, due to maximising favourable polymer-polymer interactions while minimising unfavourable polymer-solvent interactions and bending energy, such structures would not lead to rod-like aggregates significantly longer than a single chain. Thus, the known equilibrium solutionphase behaviour of semiflexible polymers is not consistent with experimental observations on P(NDI2OD-T2). While a multitude of studies have examined single-chain behaviour using simulations or theory, those examining multi-chain systems, which are more relevant for the behaviour of realistic OSC systems in which chains are rarely so isolated, are less common. Although solution-phase molecular simulations of multi-chain systems of semiconducting polymers are relatively rare, owing to the need for often prohibitively large systems to explicitly account for solvent, especially for long polymer chains, studies examining OSC solubility using all-atom (AA) molecular simulation methods can be found. Some have used mean-field solution theories such as the Flory-Huggins theory, in which simulations of short oligomers were used to estimate the Flory-Huggins parameter, which is used as a measure of solvent quality and thus the propensity for aggregation. Others have examined the aggregation mechanism and effect of solvent and polymer properties, again using short chains or implicit solvent models. While these studies provide valuable insights into some of the many factors affecting the solution-phase morphology, the Flory-Huggins theory provides a fairly simplistic model of the effects of polymer chain length and the relative strength of the solvent-solvent, solvent-polymer, and polymer-polymer interactions on solubility; it does not capture the roles of chain stiffness and conformation and so cannot account for extended aggregates expected for P(NDI2OD-T2). Other simulation studies that have more accurately calculated solubility through free energy perturbation methods, and examined aggregation mechanisms and the effect of various molecular properties on the solution-phase behaviour, have not been able to reach experimental chain lengths and have generally considered only the equilibrium behaviour as the time scales relevant to the kinetic processes are generally not accessible to detailed AA models. In this work we have developed a systematically coarse-grained (CG) model of P(NDI2OD-T2) in order to investigate the reported formation of large extended aggregates in poor solvents, to reconcile discrepancies between experimental findings on the solution-phase morphology of this polymer in such solvents, and, more broadly, to clarify the general factors that control the solution-phase morphology of semiconducting polymers. By combining atoms with correlated motion into a single CG site, and accounting implicitly for the solvent, the number of degrees of freedom of the system can be greatly reduced. This allows access to polymer length and time scales on the order of those studied experimentally. A similar model has previously been developed for the commonly studied semiconducting polymer P3HT and used to accurately predict the experimental solution-phase conformation of this polymer, giving results consistent with a more computationally expensive AA model and experiment.