Modelling structural properties of cyanine dye nanotubes at coarse-grained level

Self-assembly is a ubiquitous process spanning from biomolecular aggregates to nanomaterials. Even though the resulting aggregates can be studied through experimental techniques, the dynamic pathways of the process and the molecular details of the final structures are not necessarily easy to resolve. Consequently, rational design of self-assembling aggregates and their properties remains extremely challenging. At the same time, modelling the self-assembly with computational methods is not trivial, because its spatio-temporal scales are usually beyond the limits of all-atom based simulations. The use of coarse-grained (CG) models can alleviate this limitation, but usually suffers from the lack of optimised parameters for the molecular constituents. In this work, we describe the procedure of parametrizing a CG Martini model for a cyanine dye (C8S3) that self-assembles into hollow double-walled nanotubes. First, we optimised the model based on quantum mechanics calculations and all-atom reference simulations, in combination with available experimental data. Then, we conducted random self-assembly simulations, and the performance of our model was tested on preformed assemblies. Our simulations provide information on the time-dependent local arrangement of this cyanine dye, when aggregates are being formed. Furthermore, we provide guidelines for designing and optimising parameters for similar self-assembling nanomaterials.


Parameter optimisation
The bonded terms for the polymethine bridge and the partial charges of the aromatic core (C0C0) were optimised based on QM calculations. Specifically, the torsional potentials to be added to the dihedral between atoms C2-C3-C5-N9 and C4-C2-C3-C5 have been determined based on the energy difference between the QM and the MD dihedral profiles. For this purpose, models for the C0C0 (Figure 1b of the main text, where R represents hydrogen atoms) and C1C1 (Figure 1b, where R represents methyl groups) molecules were used. The atoms of these aromatic cores are the same as the C8S3 core, but the simple side chains allow us to focus solely on the polymethine bridge and the conjugated aromatic core. QM dihedral profiles were obtained by performing a relaxed scan with angle increments of 5 degrees using the Gaussian16 software ? . For these scans, ωB97xD functional with 6-311G(d,p) basis set was used. MD dihedral profiles were obtained by performing similar relaxed dihedral scans, where the existing dihedral potentials on the dihedrals of interest were removed. Then, the difference between the QM and MD profiles was fitted to a Ryckaert-Bellemans (RB) type of function, Eq. ??. The point charges for each system were generated after optimising the C0C0 and C1C1 structures with the Hartree-Fock (HF) method and 6-31G* basis set in two different ways: Dipole Preserving Analysis (DPA) ? using the GAMESS-UK software ? , and Restricted Electrostatic Potential (RESP) ? method using the Gaussian16 software ? . The force constants for optimised MD parameters are reported in Table ??.
n is the number of added potentials, C n is the RB coefficient for each potential and θ is the dihedral angle.  An overlap of the energy profiles between the QM and MD methods for C0C0 and C1C1 is shown in Figure ??. Since the values for the RB coefficients of the dihedral angles are similar between C0C0 and C1C1, we can assume that the contribution of the substituents does not significantly affect the potentials for the rotation around the bonds of the polymethine bridge. To maintain simplicity, the respective C0C0 force constant coefficients were used for the simulated cyanine dyes.  Figure S1: Potential energy profiles for dihedrals of C0C0 for the HF-RESP method. a) C5-C3-C2-C4. b) N9-C5-C3-C2. c) Atom index of C8S3 atomistic model. The atom numbers of the aromatic core among the generated models are identical.
After obtaining an optimal set of parameters for the dihedrals of the polymethine bridge at QM level, we compared different methods for calculating partial charges based on their ability to reproduce structural features from experimental data. Specifically, atomistic simulations of the crystal structures of two cyanine dye, C2C2 ? and C8O3 ? , were used to evaluate the performance of three different models: • Model 1: DPA charges from QM calculations (B3LYP/6-31G*) on C1C1 and cosinoid dihedral definitions ?
The structural features of the crystal structures that were evaluated were the dimensions of the supercell (simulation box) and its density. Additionally, the root-mean-square deviation (RMSD) from the initial conformation and the position of the first peak in Radial Distribution Function (RDF) calculations are reported. The results from the crystal simulations are reported in Table ??.
In general, even though the differences are quite small, the models with HF-RESP charges performed better in terms of maintaining the structural properties of the cyanine crystal structures. Model 1 performed worse than the others, especially in the C8O3 crystal. The results for Model 2 and Model 3 are quite similar, but Model 3 reproduced slightly better the position of the first peak in the RDF calculation. The parameters of Model 3 were used for the atomistic simulations of the C8S3 monomer in water. The partial charges for C8S3 are reported in Table ??.

C8Scoarse-grained model
In Martini 3, the mapping of aromatic moieties is based on the centre of geometry of all atoms (heavy and hydrogen atoms) that participate in each chemical group. However, this rule is slightly flexible, since parametrization aims at optimising the surface area of the molecules, and different definitions can be used to maintain the overall shape. In the C8S3 model, constraints hold together the aromatic rings, whereas bonds connect the polymethine bridge with the benzimidazoles, Figures S2 and S4. Each benzimidazole ring is constituted by four normal beads, one virtual site at their centre of geometry and one dummy particle. Virtual sites have no mass, so the mass of the respective bead has been evenly distributed to the four normal beads. In contrast to virtual sites, dummy particles do not interact via non-bonded interactions with any bead, and only act as supportive particles to allow specific conformations. The angles and dihedrals that control the orientation of the aromatic core are presented in Figures S3, S5 and S6. There were no angles between the side chains and the central bead to allow free rotation around the defined axis. No dihedral definitions were used between the core and side chains, since the angle definitions were sufficient to describe the movement of the side chains. Finally, the aromatic core of these cyanine dyes is positively charged. The extra charge was assigned to the beads that represent atoms of the polymethine chain. The solvent accessible surface area (SASA) was calculated for the atomistic and the CG model, in order to compare the final size of the models. The initial mapping underestimated the volume of the aromatic core. Overlapping the volume of the atomistic and CG models revealed that the area around the chlorines was slightly smaller in the initial CG models (Model 1). Consequently, the mapping was modified to increase the overlap of the two surfaces. Instead of mapping the SX3 bead on the centre of geometry of the CHCCl group (Model 1), the bead was placed at the centre of the C-Cl bond (Model 2). The surface of the aromatic core with the new mapping is almost identical to the atomistic value, Figure ?  C1A-C2A-C3A Figure S5: Comparison of atomistic and CG angle distributions for the tails of C8S3. Figure S6: Comparison of atomistic and CG dihedral distributions for C2C2. Figure S7: Solvent accessible surface area for the aromatic core of the C8S3 molecule. The histograms represent the SASA for one monomer in solution for a trajectory of 100 ns. The number of dots for SASA in the histogram analysis was set to 10000, whereas for the surface representation the number of dots were 50.
4 Self-assembly of C8S3      Figure S10: Snapshots from a C8S3 nanotube simulation with the herringbone arrangement after 500 ns in the production phase (left panel). Initial arrangement of C8S3 molecules (right panel).

C8S3 bundle preparation
The C8S3 bundle construction can be summarised in following steps: i) preparation of a single C8S3 nanotube (System 1), ii) replication of the single nanotube and translation of the nanotubes' position until the optimal thickness (2.5 nm) is achieved, and iii) removal of the overlapping C8S3 molecules of the outer wall. Figure S11: Schematic representation of the C8S3 bundle preparation.