Development of new coarse-grained models for chromonic liquid crystals: insights from top-down approaches

ABSTRACT Two top-down coarse-grained molecular simulation models for a chromonic liquid crystal, 3,6,7,10,11-hexa-(1,4,7-trioxa-octyl)-triphenylene, are tested. We use an extension of the well-known MARTINI model and develop a new coarse-grained model based on statistical associating fluid theory (SAFT)-γ perturbation theory. For both models, we demonstrate self-assembly in the isotropic phase of the chromonic and we test the effectiveness of both models in terms of the structures of the chromonic aggregates that are produced in solution and the thermodynamics of association. The latter is tested by calculations of the potential of mean force for dimers in solution, which measures the strength of molecular association. SAFT-γ provides valuable insights into the thermodynamics of assembly. Exploration of a range of interactions between unlike sites demonstrates that chromonic self-assembly only occurs in a small parameter space where the hydrophilic–lipophilic balance between aromatic core and ethylene oxide chains is optimal. Outside of this balance, chromonic self-assembly is replaced by the formation of conglomerates of molecules or short stacks. Graphical Abstract


Introduction
There has been considerable recent interest in chromonic liquid crystals. Chromonics are a non-conventional form of lyotropic liquid crystal [1][2][3][4]. Typically, chromonic mesogens are composed of a hydrophobic core (often a planar aromatic disc) with attached hydrophilic groups. This arrangement leads to molecules aggregating in aqueous solution to form stacks, which, at higher concentrations, self-organise to form lyotropic liquid crystal phases. Chromonic phases were initially found in groups of drug and dye molecules, with ionic functional groups providing the hydrophilic CONTACT Mark R. Wilson Mark.Wilson@durham.ac.uk Supplemental data for this article can be accessed here.
From a simulation perspective, chromonic systems are extremely challenging. Self-assembly into stacks occurs over periods of tens or, in most cases, hundreds of nanoseconds. Moreover, formation of liquid crystal mesophases, by self-organisation of stacks, takes place over time scales that are considerably longer. In two recent studies, progress has been made in understanding and validating the structures expected in short chromonic stacks via long atomistic simulations [15,16]. There have also been previous attempts to produce coarse-grained models of chromonics. Chromonic self-assembly and phases have been seen in simple disc models [17][18][19] and, at longer length scales, dissipative particle dynamics studies have recently demonstrated the first full simulations of chromonic liquid crystal phase diagrams using up to 5000 mesogens in water [20,21].
However, there is a real challenge to produce coarsegrained models that capture the key structural features responsible for the local molecular packing within stacks and yet are computationally tractable, allowing longer time scale events to be seen. Such a model would provide a bridge between atom the valistic models and simple phenomenological models and would potentially shed light on many of the intriguing questions posed by chromonic self-assembly. For example, can simulations rationalise the architecture of complex stacks with multiple molecule cross-sections and/or unusual self-assembly behaviour [22,23], or can simulations unravel the structure of large aggregate structures, such as hollow cylinders [24], which have been reported in some systems?
In this paper, we attempt to produce a new effective coarse-grained model for the chromonic mesogen, TP6EO2M, that bridges the gap between atomistic and highly coarse-grained simulation studies. Ideally, a new model should have chemical specificity so that results can be linked quantitatively to the system being studied. However, the model should also be sufficiently coarsegrained to allow for the simulation of large systems. In this study, we adapt the MARTINI force field, which was originally produced for lipids [25] and peptides [26] to act as a top-down coarse-grained force field for chromonic systems. Additionally, we attempt to develop a new top-down force field using statistical associating fluid theory (SAFT-γ) [27], in which the coarse-grained nonbonded terms are obtained by fitting of thermodynamic data to produce generalised Lennard-Jones parameters (Mie potentials) to be used in molecular simulation.

MARTINI model
The coarse-grained MARTINI force field was originally developed for the simulation of lipids [25] but has since been extended to other systems such as proteins [26] and polymers [28]. MARTINI consists of a number of distinct particle types, each of which represents a particular type of chemical group. These particles interact via Lennard-Jones 12:6 potentials, with the parameters used depending on the groups involved with A ij ¼ 4ε ij σ 12 and C ij ¼ 4ε ij σ 6 , for Lennard-Jones well depth ε ij and distance parameter, σ ij . MARTINI is parametrised to reproduce the water/ octanol partitioning free energies of a wide range of compounds. The mapping is usually four heavy atoms to one coarse-grained bead, although this has been altered for some particle types. Due to its relative ease of use and transferability, particularly compared to bottom-up coarse graining methods [29,30], MARTINI has been very widely used particularly in the simulation of biological systems. It has allowed simulations to reach higher time and length scales than fully atomistic models and usually provides qualitative or semi-quantitative accuracy for these systems.
Our initial coarse-grained model for TP6EO2M is based on the mapping shown in Figure 1, which uses four types of coarse-grained sites: CC, CO, AI, AO plus water beads P4 and BP4. Each water bead represents four water molecules within a MARTINI simulation. As MARTINI is a 12:6 potential, and hence rather more repulsive at close distances than many coarse-grained water interactions, it does not have the liquid stability range of real water and is apt to freeze in contact with organic systems. Hence, as previously employed in many biological simulations, we use anti-freeze particles, BP4 sites, which interact in the same way with the TP6EO2M, as normal water cluster sites but which have a larger effective radius when interacting with normal water, P4 sites. (Without the anti-freeze particles, the presence of a chromonic stack is found to nucleate the freezing of MARTINI water.) The nonbonded interaction parameters for the model were taken from the MARTINI benzene and polyethylene oxide (PEO) models.
In the original PEO model, each bead represents a single monomer. However, the hydrophilic arm in TP6EO2M is actually a methoxy group, and the mass of the end bead was altered to reflect this. In the benzene model, a bead represents a C 2 H 2 unit; the mass of the beads in the aromatic core was therefore altered according to the number of hydrogen atoms associated with that bead. Table 1 contains the Lennard-Jones parameters for the interactions in this system. Cross interactions were determined using the usual mixing rules: Bonds and angles were described by simple harmonic potentials, with bond lengths and angles taken from the energy minimised atomistic structure, mapped onto a coarse grained representation. Six improper dihedrals were required to keep the core of the molecule planar, but no other dihedral interactions were used in this model.
The GROMACS topology file for the system is provided in the Supporting information. Bond interactions are used in preference to LINCS constraint bonds, as the highly interconnected nature of the core makes satisfying all of the constraints extremely difficult. Intramolecular interaction parameters are given in Tables 3-4.

SAFT-γ model
The SAFT version used in this work is SAFT-γ Mie [31], which is based on a version of SAFT suitable for potentials of variable range (SAFT-VR [32]). The theory uses   Mie (generalised Lennard-Jones) potentials as developed by Lafitte et al. [33].
In the non-associating form of the theory, the Helmholtz free energy can be expressed as where A ideal accounts for the ideal kinetic and translational energy contributions, A monomer for the change in free energy due to attractive and repulsive interactions of the CG beads and A chain for the free energy contribution by forming chains. The molecular model in the SAFT-γ Mie equation of state (EoS) is a chain of tangentially connected beads interacting through Mie potentials where ε and σ are, respectively, the well depth and the segment diameter, λ r and λ a are the repulsive and attractive exponents, for which λ r ¼ 12 and λ a ¼ 6 is the usual form of the Lennard-Jones potential.
We follow the strategy used by Müller, Jackson and co-workers for using SAFT to parameterise the nonbonded parts of a coarse-grained model [34]. This has recently been successfully applied to lyotropic liquid crystal surfactants [35]. The SAFT model for TP6EO2M was built by combining Mie potential parameters that were fitted to experimental vapour-liquid equilibria over a wide range of temperatures. The 3bead benzene model by Lafitte et al. [36] was used for the core beads and the PEO bead by Lobanova [37] was used for the side chains. For the water beads, the coarse-grained biowater model by Lobanova et al. [38] was used, which is optimised for the saturatedliquid density and surface tension at ambient conditions. Ref. [37] also provides unlike interactions between the water and an EO bead, which were fitted to experimental excess enthalpies. However, in this work, we refitted all cross-interaction terms for our model, as described below, using the combining rules as a starting point to the fitting of Mie potentials to reproduce thermodynamic data.
The usual combining rules for contact distances give , for repulsive and attractive exponents [33], and well depth Here, the parameter k ij measures the extent of preferred interactions between species. For species that are chemically incompatible, for entropic or enthalpic reasons, k ij > 0, and k ij < 0 for species where one group is well solvated by another.
In the current work, the EO-aromatic (EO-Ar) site interaction was fitted to 1,2-dimethoxyethane/benzene ((EO) 2 -benzene) enthalpies [39,40] of mixing with k ¼ À0:0166 as the best fit. Aromatic-water (Ar-W) interactions were fitted to benzene-water liquid-liquid equilibrium compositions [41][42][43] with k ¼ þ0:392 as the best fit. EO-water (EO-W) interactions were fitted to triethylene glycol dimethyl ether/water ((EO) 4 -W) enthalpies of mixing [44,45] with k ¼ À0:167 as the best fit. Thus, interactions involving EO are more attractive (EO-Ar, EO-W) than the combining rules suggest, while the Ar-W interaction is considerably less attractive than the combining rule predicts. The solubility and excess enthalpy calculations of the SAFT models were performed using the SAFT-γ Mie EoS and numerical solvers provided by the gSAFT package (Process Systems Enterprise Ltd. [46]).
The Mie potentials (Table 5) were used as nonbonded interactions in the MD simulations, in combination with intramolecular parameters derived from the MARTINI force field. In our initial SAFT model, the σ values obtained from nonbonded potentials were used as equilibrium bond lengths (in the spirit of the SAFT theory which uses tangentially bonded spheres) in combination with the same force constants used in the MARTINI model. The angle potentials used were identical to the MARTINI model. The same improper dihedral potentials were also used in the SAFT model, but with a five times stronger (150 kJ mol −1 ) force constant to keep the core planar. In a later development of the SAFT model, we used the same bond lengths as used in MARTINI, in combination with SAFT nonbonded parameters.

Coarse-grained simulations
All molecular dynamics simulations were carried out using GROMACS 4.6.7. The equations of motion were solved using the leap-frog integrator with a time step of 2 fs (MARTINI) and 5 fs (SAFT). All simulations were carried out at a temperature of 280 K and a pressure of 1 bar, using the Nosé-Hoover thermostat and the Parrinello-Rahman barostat. The van der Waals, neighbour list and coulomb cut-offs were set to 1.2 nm (MARTINI) and 1.5 nm (SAFT). Simulations of a stacked MARTINI configuration were prepared by initially creating a stack of 10 TP6EO2M molecules, and solvating with 4709 water beads and 523 antifreeze beads. An energy minimisation was carried out, and the structure was equilibrated for 400 ps in the constant-NVT ensemble. A production run was then carried out in the constant-NPT ensemble for 100 ns.
MARTINI simulations were also carried out from a dispersed starting configuration. Ten TP6EO2M molecules were placed randomly in a simulation box and solvated with 3917 water beads and 435 antifreeze beads. The energy of the system was minimised, and a 10-ns equilibration run was carried out in the NVT ensemble. A 100-ns production was then carried out in the NPT ensemble.
For the SAFT-γ simulations, a randomly generated vapour configuration of 10 TP6EO2M molecules and 16,740 water beads was prepared. The system was compressed at 500 bar and 450 K for 1000 ps and equilibrated at 280 K and 1 bar for 500 ps using the Berendsen thermostat and Berendsen barostat in the NPT ensemble. A production run was carried out of at least 100 ns with the Nosé-Hoover thermostat and the Parrinello-Rahman barostat.

Potential of mean force calculations
The potential of mean force (PMF) for a TP6EO2M dimer was calculated by constraining the separation distance between the centres of mass (COMs) of the cores of two TP6EO2M molecules at a range of distances from 0.35 to 2.35 nm and calculating the average constraint force, hf C i. For each distance, V PMF was calculated by integrating over the separation distance, using: where k B is Boltzmann's constant and T is the temperature. First, a dimer was placed in a simulation box with 6914 water beads and 798 antifreeze beads, and the structure was energy minimised. The system was equilibrated in the constant-NPT ensemble for 1 ns with the Berendsen thermostat and barostat, followed by 5 ns with the Nosé-Hoover thermostat and the Parrinello-Rahman barostat. A further 5-ns constant-NVT equilibration was carried out at the equilibrium density from the constant-NPT simulations. An initial pull simulation was carried out, with the distance between the COMs of the aromatic cores of the two molecules constrained and increased from 0.35 nm at a rate of 0.01 nm ps −1 . Frames were then selected from this simulation at a range of distances, clustered more closely around the expected energy minimum. At each of these distances, a 25-ns simulation was carried out with COM-COM distance constrained, the last 20 ns of which was used to generate data for the PMF calculation. Figure 2 shows a series of snapshots taken from an initially dispersed TP6EO2M system in water. The system spontaneously self-assembles over time to form short chromonic stacks, with aggregates (Haggregates) growing through addition of a monomer or merging of two short stacks. If a stack is preassembled and solvated, it is stable for long simulation times (>100 ns). An equilibrated self-assembled stack shows the same structure as that seen atomistically for TP6EO2M, with the molecules exhibiting liquid crystalline, rather than crystalline, order.

MARTINI model
Two distinct measures of the distance between two molecules in a stack have been used in previous studies of TP6EO2M. The first, d com , is simply the distance between the COMs of two molecules, while the stacking distance, d S , is the distance between the molecules projected along the average of the vectors normal to the cores of the two molecules. Probability distributions for the two distances are shown in Figure 3. The maximum of the peak d S for adjacent molecules is 0.42 nm. The distribution for d com is shifted so that the maximum is at 0.46 nm. This reflects the fact that there is a slight offset between molecules in the stack. This stacking demonstrates H-aggregation, rather than J-aggregation. The latter would have a larger systematic offset. It is noted that the peaks do not have a long tail at higher distances, as there are no distinct bends in the stack. As expected, the stacking distances between molecules are close to the contact distance of the MARTINI nonbonded potentials for the CG sites. These separations are the minimum possible for this model due to the size of beads used but are (of course) slightly larger than the experimental and atomistic value of ,0:35 nm, which correspond to the distance of separation for single aromatic carbon atoms.
The PMF for the separation of a TP6EO2M dimer is given in Figure 4. The PMF has a similar shape to that reported for the atomistic system by Akinshina et al., with a short range repulsive and a long range attractive part. The free energy of association, ΔG, for the formation of a dimer can be obtained from the maximum well depth of the PMF. For this model, ΔG = −73.9 ± 1.5 kJ mol −1 at a concentration of 1:7 Â 10 À3 mol dm À3 , which is higher than the value of −40.7 ± 1.4 kJ mol −1 for a PMF obtained for the atomistic system. Therefore, the MARTINI model provides only semi-quantitative thermodynamic data.  We note that coarse graining, as desired, allows for a huge speed-up in observing chromonic assembly, with savings arising from a faster exploration of phase space in addition to a large reduction in the number of sites simulated. TP6EO2M reduces from 138 sites to 27 (a factor of ,5 Â and four water molecules (12 sites) reduce to 1 site.

SAFT-γ model
Using the initial SAFT-γ model developed above, which includes fitted cross-interactions, we see chromonic stacking, as seen in the MARTINI model. Molecular spacings in a stack are % σ CO corresponding to the smaller coarse-grained sites arising in SAFT-γ. It is interesting to note that if instead of fitted cross-interactions, we use instead interactions generated by standard combining rules for Lennard-Jones interactions, then SAFT-γ TP6EO2M molecules are perfectly soluble in water, with no association over % 190 ns. Association into chromonic stacks therefore arises when there are less favourable interactions between aromatic groups and water. Noting that part of this unfavourable interaction includes the hydrophobic interaction from release of water when aromatic groups interact, which is incorporated into the nonbonded interaction parameters within a coarse-grained model.
To investigate this further, we looked at a range of SAFT-γ models, as shown in Figure 5. For these models, we consider only combination of cases where k ij ¼ 0:2 or 0.0. The former value corresponds to unfavourable interactions relative to the combining rules. Here, as shown in Figure 5, a range of different aggregation behaviour is possible, from no aggregation where molecules are perfectly soluble, to conventional chromonic stacking, to effective phase separation with formation of non-equilibrium conglomerates of molecules or short stacks.
In SAFT, equilibrium bond lengths, r 0 , between coarse-grained sites are usually set to the same value as σ ij for a given pair of atom types. This means that while the bulk properties are correct, the internal structure is less faithful to an underlying atomistic framework. If SAFT is used to generate parameters for a simple polymer or oligomer chain or used for a simple linear surfactant, this is unlikely to be a problematic. However, for TP6EO2M, we noted the presence of some configurations where a dimer was able to sample the arrangement shown in Figure 6. We particularly noted this in PMF calculations, when the core-core distance was constrained to relatively small distances. The structure seen here is unphysical and can be attributed to the relatively large r 0 : σ ratio. To eliminate these configurations, we took the fitted SAFT model and further modified it to use the same bond lengths as used in the MARTINI model. After this modification, the interdigitated structure of Figure 6 was no longer observed (see Supplementary data).
The PMF generated with this final SAFT model is shown in Figure 4. As with MARTINI, a deep attractive potential well is seen, which is responsible for spontaneous aggregation in solution. The slight net repulsion seen at longer COM distances reflects the slight  preference for EO groups to be fully solvated in water rather than to be in contact with each other. In contrast to MARTINI, SAFT underestimates the free energy of association. It is possible that improvements to the SAFT well-depth may be possible if appropriate free energies were available for fitting of ethylene oxide cross-interaction terms (see Section 4), and also if we were able to move to the fitting of free energies for the interaction of extended aromatic systems instead of benzene.

Discussion
At first glance, the non-ionic chromonic system, TP6EO2M, appears to be quite a simple system to model, i.e. representable by two types of Lennard-Jones beads. However, the simplicity of the chemical structure hides a wealth of subtle complexities arising from modelling the right hydrophilic-lipophilic balance within different parts of the molecule. The initial clue to the complexity lies in the earlier atomistic simulation work of Akinshina et al. [16]. In this paper, it was shown that the contribution to association-free energy arises from both entropic and enthalpic contributions and these are similar in magnitude (ΔH ¼ À12:2 kJ mol À1 and À TΔS ¼ À16:3 to À 18:7 kJ mol À1 ). Favourable enthalpic interactions arise when discs self-assemble but there is also a strong hydrophobic interaction arising from the release of solvent. Moreover, the balance of interactions between chains and solvent, and between aromatic cores and solvent, is also important. For example, changing the atomistic force field from OPLS to GAFF leads to reduced hydration of ethylene oxide chains and consequently to association which is too strong to lead to the formation of fully hydrated chromonic stacks.
In this coarse-grained study, the same balance of interactions is found to be absolutely crucial. Moreover, within a coarse-grained representation, the coarse-graining procedure changes the balance between enthalpic and entropic interactions. Coarse-grained water has a lower gain in entropy, in comparison with real water, when freed from interaction with a chromonic molecule. Consequently, the free energy balance within a coarse-grained model must be represented by more favourable enthalpic interactions favouring association. One initial consequence of this observation is that in coarse graining, matching enthalpic interactions to an atomistic reference in a chromonic (or related) system is going to lead to poorer thermodynamics than if free-energies are used in the coarse-graining procedure. We note also the early work of Klein and co-workers who used transfer free energies as part of the development of early successful coarse-grained surfactant models [47].
The reason for the simple MARTINI model being successful, in terms of capturing chromonic assembly, is therefore very apparent. The nonbonded interactions within MARTINI are designed to model the partitioning (and hence the transfer free energy) from water to octanol. MARTINI is therefore able to do a reasonably good job of capturing both the overall hydrophiliclipophilic balance for a chromonic molecule in water, and the relative hydrophilic-lipophilic balance between the different parts of the molecule.
Coarse graining using SAFT-γ provides further insights. Figure 5 shows that across a small parameter range for the unlike interactions, one can achieve a range of aggregation behaviour, with chromonic selfassembly only occurring at the right balance of hydrophilic and lipophilic interactions for the aromatic and ethylene oxide components. If the balance of interactions is wrong, the self-assembly process may not be observed within a reasonable simulation time, either due to a lack of self-assembly, or the formation of other conglomerates.
The ability of the SAFT-γ procedure to fit to free energy, enthalpy, entropy and free-energy dependent data provides a powerful mechanism to match the correct hydrophilic-lipophilic balance. To ensure a correct balance, free energy (or free-energy dependent) experimental data have to be available for all interactions of interest. For the SAFT-γ model in this work, we have used free energy and free-energy dependent data for four out of six interactions (one out of three unlike interactions), while for the water-ethylene oxide and the aromatic-ethylene oxide interactions, we have used enthalpies of mixing due to lack of other experimental data. The shortcomings of the SAFT-γ are at least partially attributed to not ensuring the correct entropic contributions in those two ethylene oxide interactions. We have also highlighted another weakness to the SAFT-γ approach to coarse-graining. In order to obtain an accurate structure for the TP6EO2M molecule, we manually changed the bond lengths to match the atomistic structure. While successful for a dimer or stack system solvated with water, this approach would be likely to cause issues if dealing with a bulk system, because the model has been parametrised to match bulk quantities, such as densities, on the basis that r 0 and σ are equal. This suggests that the current ability to construct models for complex chemical systems with SAFT-γ is limited by its structural accuracy, and this is something that must be addressed in future work. Moreover, for complex chemical structures, such as TP6EO2M, a small contribution to the balance of hydrophilic/hydrophobic/lipophilic interactions arises from the parts of configurational space that are accessible (or inaccessible) in a dimer interaction. In this example, the accessible configurational space can be modified by the bond lengths used in representing complex structures.
Finally, it is worth commenting on a comparison of the methodologies adopted here and alternative bottom-up strategies, which rely on a fit to an underlying atomistic reference system. The two most commonly used bottom-up methods are the multi-scale coarsegrained (MS-CG) method of Voth and co-workers [48][49][50] (often referred to as force matching) and the iterative Boltzmann inversion (IBI) method [51]. IBI potentials are traditionally fitted to reproduce distribution functions in the atomistic reference system, so would be able to reproduce the structure of a chromonic stack given that structure within an atomistic system. However, good fitting is difficult to achieve in IBI when there are multiple partial radial distribution functions (RDFs) to fit arising from cross-interactions. Moreover, for systems such as this, the reproduction of thermodynamic quantities, particularly association free energies, tends to be poor, and it is unlikely that the same IBI model can reproduce both the correct structure of a chromonic stack and the correct selfassembly behaviour seen in dilute solution. The MS-CG method is, in principle, easier to fit to an atomistic reference system for a chromonic than multiple partial RDFs. However, here transferability between different concentrations and/or different aggregate states would be expected to be rather poor. In addition, as with IBI, there is no strong expectation that the thermodynamics of assembly, seen to be so crucial in the current study, could be easily reproduced.
In preliminary work, we have tested the MS-CG approach for the TP6EO2M/water system using an underlying atomistic reference composed of dispersed monomers, dimers and trimer stacks, solvated with water. This reference was chosen due to the tightly bound nature of the molecules in a chromonic stack, which led to poor sampling of the forces between atoms for many distances and orientations. The MS-CG scheme does not necessarily guarantee that the CG model will recover the reference atomistic structure, due to the assumptions made in trying to map all of the interactions in the atomistic system onto two-body coarse grained potentials. In practice, we found that using this MS-CG model, TP6EO2M tends to induce crystallisation of the CG water, so the correct aggregation behaviour is not seen. Running CG simulations without water, using the same MS-CG force field, led to the rapid collapse of a stacked starting configuration into a disordered conglomerate.
Hence, the modelling methods suggested in the present study currently seem to be the best realistic approach for designing coarse-grained models of chromonic systems when it is necessary to capture specific molecular interactions.

Conclusions
We have successfully produced two coarse-grained models to study molecular association in the non-ionic chromonic TP6EO2M/water system. The MARTINI model, extended to chromonics, is shown to work well in modelling self-assembly in solution, giving reasonably good PMF interactions for molecular association. In this model, considerable speed-ups are accessible in comparison to atomistic models. These are achieved via a considerable reduction in the number of interaction sites for the coarse-grained model: a reduction in sites of ,11:7 Â is achieved, depending on system size and concentrations, leading to a dramatic reduction in pair interactions , N at =N CG ð Þ 2 . Second, a speed-up in dynamics occurs as coarse-grained models move more quickly through phase space than their atomistic counterparts.
We show also that SAFT-γ perturbation theory provides a powerful methodology to carry out coarse-graining of the intermolecular interactions. We are able to use SAFT-γ fitting to experimental data to produce a coarsegrained model to describe TP6EO2M self-assembly in water. SAFT provides molecular insight into the nature of the interactions which are responsible for chromonic phase behaviour. Moreover, we see that the correct hydrophilic-lipophilic balance is absolutely crucial for chromonic association to take place in water, and this statement applies to the interactions of component parts of a chromonic molecule within water, in addition to the molecule as a whole.
In future, we suggest that the coarse-graining approaches we use here, which are both fast and effective, could be used for screening molecules for chromonic behaviour and designing new chromonic systems with desired self-assembly behaviour.