Modeling of excitonic properties in tubular molecular aggregates Bondarenko, Anna

Modeling of excitonic properties in tubular molecular aggregates Bondarenko, Anna

DOI:

10.33612/diss.98528598

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date:

2019

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Bondarenko, A. (2019). Modeling of excitonic properties in tubular molecular aggregates. University of Groningen. https://doi.org/10.33612/diss.98528598

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

Download date: 28-06-2021

4

Multiscale modeling of complex molecular aggregates

λ (nm)

A(λ)

Experiment Modelling

Published as A.S. Bondarenko, I. Patmanidis, R. Alessandri, P.C.T. Souza, T.L.C. Jansen, A.H. de Vries, S.J. Marrink, and J. Knoester, 2019, submitted.

55

4

Supramolecular aggregates of synthetic dye molecules offer great per- spectives to prepare biomimetic functional materials for light-harvest- ing and energy transport. The design is complicated by the fact that structure—property relationships are hard to establish, because the molecular packing results from a delicate balance of interactions and the excitonic properties that dictate the optics and excited state dynam- ics, in turn sensitively depend on this packing. Here we show how an iterative multiscale approach combining molecular dynamics and quantum mechanical exciton modeling can be used to obtain accurate insight into the packing of thousands of dye molecules in a complex double-walled tubular aggregate. The thus obtained optical spectra compare very well with experiment. The microscopic model devel- oped opens the route to accurate predictions of energy transport and other dynamic properties.

4.1. Introduction

S upramolecular structures may self-assemble from a variety of building blocks, resulting in a wide range of advanced materials with attractive biomimetic, sens- ing, catalytic, optoelectronic and photonic functionalities.

The close-packed nanoscale organization of the individual molecules within a supramolecular system, held together via noncovalent interactions, gives rise to the aggregate’s (collective) properties. In case of natural or synthetic dyes, the assemblies often exhibit unique collective optical properties and are of interest for opto-electronic applications as well as artificial light-harvesting complexes that mimic natural antennae systems of photo- synthetic bacteria and plants.

For example, chlorosomal antennae complexes of photosynthetic green sulfur bacteria are self-assembled into multilayer tubular structures having bacteriochlorophyll pigments as building blocks.

The struc- ture of these antennae complexes and the underlying molecular arrangement ensures that the process of light-harvesting and excitation energy transport is very efficient, even under extremely low light conditions.

The quest to recreate such efficiency under laboratory conditions has sparked numerous studies of synthetic self-assembled systems mimicking natural chlorosomes, e.g. using porphyrins,

zinc chlorin,

and cyanine dyes.

Of particular interest are the tubular aggregates of 3,3’-bis(2- sulfopropyl)-5,5’,6,6’-tetrachloro-1,1’-dioctylbenzimidacarbocyanine (C8S3).

Cryo-TEM reveals a hierarchy of supramolecular architectures, including double-

walled nanotubes; under certain conditions, bundles of nanotubes arise.

Thus, this

4

system allows for the occurrence of electronic excitation energy transport at various levels: within one wall, between walls of one tube, and between different tubes, similar to the situation in natural systems.

To understand how such supramolecular systems work, as well as propose design rules for new materials, it is essential to determine the relationship between molecular structure and optical properties. Current experimental techniques, however, are unable to resolve the structure at the molecular level. This, in combination with the sensitivity of spectral properties to the details of the molecular packing, leads to a crucial role for theoretical modeling.

For example, molecular dynamics (MD) simulations have been used to predict the molecular packing within a variety of supramolecular assemblies.

However, synthetic amphiphiles with aromatic groups, such as cyanine dyes—often used to prepare aggregates with optical functionality—tend to fall into kinetic traps during spontaneous self-assembly simulations and the packing of the aromatic chromophores remains highly disordered on the accessible time scale, leading to predicted (optical) spectra that are not consistent with experimental data.

This problem can be overcome by building assemblies based upon proposed architectures and assessing their stability in relatively short MD simulations.

The drawback of this approach is the requirement of a thorough understanding of what to use as a starting point and how to validate the structure. In any case, proper validation requires the modeling of the optical spectra of the obtained structure, and finally, comparing it to the experiment. The demanding character of such methods explains why an important role is played by phenomenological modeling, in which a molecular packing is guessed and the optics is obtained from parametrizing an exciton model that describes the collective excited states of the assembly with interactions dictated by the guessed packing. By comparing the calculated spectra to experimental ones, the structure and exciton model may be fine-tuned. While this method has been successful in describing spectra,

it is limited in its predictive power and therefore not very useful for design of new materials.

In this work, we use an advanced multiscale approach to determine structure—

optical property relationships for the C8S3 double-walled nanotubes, guided by

comparison to experiments. The optical spectrum of these aggregates suggests a rather

complex underlying molecular packing, making it an exceptionally challenging system

to resolve. To do so, we iteratively combine the use of MD simulations to capture the

details of molecular packing and structural disorder, an exciton Hamiltonian approach

to calculate optical signatures, and explicit microelectrostatic calculations to estimate

energetic disorder and solvent shifts. Previous attempts to reveal the structure of

cyanine-based nanotubes were limited to small-scale system sizes,

modeling

optical features phenomenologically rather than using atomistic information

or

featuring simpler, single-walled systems.

Our current approach not only allows us

to unravel the molecular packing details of the C8S3 double-walled nanotubes, but also

4

to understand how local structure subtly influences collective optical functionality. In a broader sense, our study opens the way to explain and predict at an unprecedented level of detail the functional properties of highly complex molecular materials.

4.2. Results and discussions

Multiscale approach to probe molecular packing. Our multiscale approach to reveal the molecular packing of the complex double-walled tubular aggregate includes the following steps (Figure 4.1; see Section 4.3 and 4.4 for details). 1) We construct an initial set of ideal—perfectly ordered—single-walled tubular structures with different structural parameters taking into account a number of restraints (also see below) from available experimental and modeling data, namely inner wall (IW) and outer wall (OW) radii from cryo-TEM,

wall-to-wall thickness from MD simulations on bilayers,

a starting unit cell of a similar molecule,

and the herringbone lattice from phenomenological modeling

(Figure 4.1a-c). 2) Based on the comparison of the calculated absorption spectra of this set of ideal single-walled tubes to the experimental spectrum, we select candidates for the IW and the OW and preassemble them into double-walled aggregates (up to 100 nm long,

7,000 C8S3 molecules) that are used as starting point for MD simulations in an explicit solvent environment. During the MD simulations (20-100 ns), the modeled systems “acquire” realistic structural disorder, i.e., disorder in molecular packing due to the interaction with the solvent and within the aggregate itself (Figure 4.1d). 3) From the obtained MD snapshots, we estimate how the environment of each chromophore changes the gas-phase monomer excitation energy via atomistic microelectrostatic calculations, resulting in an average energy shift and so-called static energy (diagonal) disorder. 4) Next, we model the absorption spectrum of the obtained aggregate. For each snapshot, taken every 100 ps along the MD trajectory, we translate the structural information into a Frenkel exciton Hamiltonian that describes the assembly’s collective excited states responsible for the optical response; this Hamiltonian accounts for interaction (off-diagonal) disorder arising from orientational and positional disorder. For calculating the intermolecular excitonic couplings, we use the extended dipole model,

where the transition dipole is mapped on the polymethine bridge of each C8S3 molecule (Figure 4.1a).

Furthermore, for each snapshot a random energy shift taken from a Gaussian with mean and width obtained from the above microelectrostatic calculations is added to the gas phase monomer excitation energy for each molecule. By averaging the absorption spectra calculated for each snapshot, we obtain the final spectrum. By comparing this spectrum to the experimental one, structural parameters for the double-walled tube are adjusted and tested by repeating one or more of the above steps.

Figure 4.1e,f give the calculated spectrum for the “optimal” structure whose

spectrum best compares to the experimental one taken from Ref. 103. The agreement

4

<latexit sha1_base64="k/AN3iXODITolNbwP576Vc7DU5I=">AAAB7XicbVDLSgNBEOz1GeMr6tHLYBA8hd0o6DHoxWME84BkCbOT2WTM7Mwy0yuEkH/w4kERr/6PN//GyQPRxIKGoqqb7q4olcKi7395K6tr6xubua389s7u3n7h4LBudWYYrzEttWlG1HIpFK+hQMmbqeE0iSRvRIObid945MYKre5xmPIwoT0lYsEoOqnexj5H2ikU/ZI/BfkhwSIpwhzVTuGz3dUsS7hCJqm1rcBPMRxRg4JJPs63M8tTyga0x1uOKppwG46m147JqVO6JNbGlUIyVX9PjGhi7TCJXGdCsW8XvYn4n9fKML4KR0KlGXLFZoviTBLUZPI66QrDGcqhI5QZ4W4lrE8NZegCyrsQll5eJvVyKTgvle8uipXreRw5OIYTOIMALqECt1CFGjB4gCd4gVdPe8/em/c+a13x5jNH8Afexzems48t</latexit>

540 560 580 600 620

Wavelength (nm) 0

0.2 0.4 0.6 0.8 1

Absorption

IW and OW components

Experiment IW modelling OW modelling

(f)

540 560 580 600 620

Wavelength (nm) 0

0.2 0.4 0.6 0.8 1

Absorption

Total spectrum

Experiment Modelling

(e)

(a) (d)

(c) (b)

5 nm

Figure 4.1 | The model system and comparison of modeled and experimental absorption spectra. (a) A C8S3 molecule overlaid with its transition dipole moment (green arrow) and a basic unit cell including two molecules with a herringbone arrangement. (b) Two-dimensional herringbone lattice and rolling vector (blue arrow) with corresponding rolling angleµ. (c) Single cylinder obtained from the rolling of the two-dimensional lattice along the rolling vector; in the inset, the arrows represent the transition dipole moments used in the Frenkel exciton Hamiltonian—their positions and orientations are imposed by the molecular structure. (d) Snapshot of the double-walled C8S3 nanotube for the optimal (defined as best reproducing the measured absorption spectrum) structure after 100 ns MD simulation in water; sulfonate groups are in red, carbon atoms in grey (light (dark) grey for inner (outer) wall molecules), water in cyan, and Na⁺in blue. (e) Absorption spectrum calculated for the modeled system with the optimal structure, compared to the experimental one.¹⁰³(f) Contributions to the absorption spectrum from the inner and outer walls for the optimal structure.

4

is very good, not only in the overal spectral shape, but also in essential details: (i) The contributions from the inner and outer walls are correctly assigned, i.e., the lowest- energy peak from the inner wall is lower in energy than the one from the outer wall, correctly reproducing the assignment found experimentally.

(ii) The relative oscillator strength of these two peaks is reproduced correctly, as are the polarization directions of the optical transitions: the two lowest-energy absorption bands are polarized parallel to the tube’s axis, while higher bands also show regions of perpendicular polarization.

(iii) The experimental linewidths are reproduced well, where we find that energy and interaction disorder are comparable, resulting in a disorder strength of 220 cm

and 370 cm

, respectively. This compares well with expectations from previous nonlinear optical experiments.

Lifetime broadening was introduced following the procedure of previous modeling,

in accordance with expected lifetimes based on experimental observations.

The spectrum’s absolute position was fixed by using the gas-phase transition energy of a single C8S3 molecule to be 19,498 cm

; this quantity is unknown experimentally, while quantum theoretical calculations on its absolute value have large error bars.

Following the multiscale modeling procedure outlined above, we obtained insights into the subtle interplay between molecular packing motives on the one hand, and collective properties reflected in the absorption spectrum on the other hand, including testing hypotheses regarding the splitting of the IW and OW peaks, whose origin has hitherto been elusive. Our detailed account of solvent shifts and energy and structural disorder, allows us to conclude that the origin of this splitting lies in the IW’s packing density being slightly larger than the OW’s density. Below, more detailed information is given about the rationale and sensitivity of the different steps in the procedure, as well as some of the essential numbers.

Structural restraints. It is well-known that an ideal single-walled tubular aggregate with one molecule per unit cell possesses two excitonic absorption bands, polarized parallel and perpendicular to the tube’s axis, respectively.

Previous optical studies of C8S3 aggregates

revealed more than 4 bands for the double-walled system, suggesting a structure with two molecules per unit cell, leading to a Davydov splitting of the optical bands.

A phenomenological herringbone lattice was able to reproduce the complex optical spectrum of the C8S3 aggregate

and its changes upon halogen substitution of the dye,

and therefore is used as starting point for our ideal structures. Other packings, such as brickwork or staircase lattices cannot explain the multitude of experimentally observed peaks in the absorption spectrum.

For the wall radii we use values of 3.2

0.5 nm and 6.5

0.5 nm, as defined by the positions of the sulfonate groups of the IW and OW, respectively, obtained from cryo-TEM imaging.

For all the tubes, we use the optimal wall-to-wall thickness of 1.7 nm measured from the cores of the molecules.

Probing the effect of the rolling angle. The optical spectra of the tubular aggregate

4

540 560 580 600 620

Wavelength (nm) 0

0.2 0.4 0.6 0.8 1

Experiment IW MD OW MD

540 560 580 600 620

Wavelength (nm) 0

0.2 0.4 0.6 0.8 1

Experiment MD

540 560 580 600 620

Wavelength (nm) 0

0.2 0.4 0.6 0.8

1 Experiment

IW MD OW MD

540 560 580 600 620

Wavelength (nm) 0

0.2 0.4 0.6 0.8

1 Experiment

MD

(a) Rolling angle 55˚

(b) Rolling angle 23˚

Figure 4.2 | Effect of rolling angle on absorption spectra. The left panels display the total absorption spectra calculated based on an MD trajectory of 20 ns. The right panels show the contributions from both walls. (a) A rolling angle of 55^±results in a total lack of oscillator strength in the high-energy region and too much oscillator strength in the two low-energy peaks polarized parallel to the tube’s axis. (b) A rolling angle of 23^±gives rise to too small an oscillator strength in the parallel polarized bands compared to the perpendicular (higher-energy) bands. In the spectra displayed here, the spectral broadening due to energy disorder and lifetime were not taken into account.

are strongly influenced by the rolling angle,

, of the lattice (Figure 4.1b). It was previously shown that the energy position of the optical band polarized parallel to the axis of tubular assemblies is not sensitive to variations in the rolling angle

.

By contrast, the energy of the optical band polarized perpendicular to the tube’s axis strongly depends on

. Moreover,

influences the relative oscillator strength of the parallel and perpendicular polarized bands. We thus performed a systematic study of structure and absorption properties in order to determine the rolling angle.

Importantly, we found that the rolling angle does not influence the stability of the

molecular aggregates and order in the molecular packing in the MD simulations

(see Figure 4.6). Consequently, the rolling angle can not be determined from the

MD simulations only, highlighting the importance of combining several theoretical

methods. Hence, relying on the absorption spectra calculated for the simulated

tubes, we found that a rolling angle around 30

ensures a good separation between

parallel and perpendicular transitions and gives rise to relative oscillator strengths in

4

agreement with those found in the experiment (Figure 4.1e,f). By contrast, Figure 4.2 shows examples of calculated spectra of proposed structures that were discarded as candidates. Structures with a rolling angle of 55

(Figure 4.2a) result in spectra that lack oscillator strength in the high-energy region (550-570 nm), where the experimental spectrum has considerable intensity.

A rolling angle of 23

(Figure 4.2b) results in too small an oscillator strength of the parallel polarized bands compared to the perpendicular ones.

Probing the effect of molecular disorder. Organic materials generally endure the presence of significant structural (interaction) and energy disorder, which can be detrimental or essential to their functional properties. In order to study how disorder impacts the (opto-)electronic properties, it is important to be able to estimate the amount of disorder. Our multiscale approach gives direct access to such information.

We obtain a qualitative picture of the structural disorder in the modeled aggregates by plotting the angle between pairs of chromophore transition dipole moments as a function of the distance between them. The obtained maps show how well-defined angle/distance distributions in the ideal (starting) structures (Figure 4.3a,b) are smeared out in the disordered structures upon MD simulation, while the overall peak structure is preserved (Figure 4.3c,d).

The quantitative impact of the structural disorder in the system on the absorption spectrum is obtained from the distribution of the sum of the intermolecular excitonic interactions per molecule,

, responsible for the collective excited states of the system (eq. 4.5). In Figure 4.3e, we show this distribution averaged over the snapshots of the MD trajectory in comparison to the distribution obtained from the ideal structure.

For the ideal structure, this distribution is very narrow with mean values of

cm

and

cm

for the IW and OW, respectively. These values are indicative for the red-shifts of the low-energy (parallel) absorption bands of both walls relative to the monomer transitions (J-aggregate), where for the IW this shift is larger than for the OW, as a result of the small difference in packing density (see below). The structural disorder acquired in the MD simulation, gives rise to broadening of the distribution and leads to a decrease of the mean value; in the spectra, this leads to broadening and blue-shifts of the spectral peaks relative to the spectra of the starting structures. The width of the distribution and the decrease of the mean value is slightly larger for the IW than the OW. This effect was seen for all the tubes and can be explained by the fact that the anionic sulfonate groups inside the IW are packed more closely than those outside the OW; the interactions between these groups in the IW thus lead to a larger stress, in turn leading to more structural disorder.

To probe energy disorder resulting from the local electrostatic environments

†This finding differs from the value of 50^±as best fit obtained for the phenomenological herringbone model in Ref.8; this apparent contradiction results from the fact that the orientations of the molecules within the unit cell of the optimal structure found here, differ from those in the phenomenological model (see Section 4.4.1of the Appendix).

4

(a)

(c)

(b)

(d)

(e)

(f)

Figure 4.3 | Quantitative presentation of the structural disorder for the optimal system. (a)-(d) Two- dimensional distributions of the angle between pairs of chromophore transition dipole moments (¡_{i j}) and the distance between them (r_{i j}) normalized to the number of pairs, shown separately for both walls. Results for the ideal (starting) structures (a)-(b) and for the structure obtained after 20 ns of MD simulation (c)-(d).

(e) Histogram of the sum of the intermolecular excitonic interactions per molecule^P_{i 6=j}J_{i j}(eq.4.5) for ideal and simulated structures, and table with the corresponding means and standard deviations, averaged over 20 ns.

around molecules in the IW and OW, we evaluate electrostatic and polarization (or induction) interaction energies between single molecules and their surroundings (Figure 4.4a,b) via microelectrostatic calculations. Molecules are treated classically and described by atomic point charges and atomic polarizabilities that differ for the molecular ground and excited states. The shift in the excitation energy induced by the molecular surroundings is calculated as the difference

between the interaction energies

and

of the charges and polarizabilities characterizing the molecular excited and ground states, respectively, with the charges and polarizabilities characterizing the surrounding molecules.

Figure 4.4c shows the distribution of the computed energy shifts for all the

molecules of a 10 nm thick fragment (taken from the middle) of the aggregate,

containing approximately 350 and 450 IW and OW molecules, respectively. As is seen,

molecules in both walls experience very similar electrostatic environments, reflected

by the fact that the shift distributions for both walls are practically identical, with mean

values around 350 cm

. The standard deviations of the Gaussian functions fitted to

these distributions is found to be 213 and 231 cm

for the IW and OW, respectively,

and is a measure of the energetic disorder. This disorder was used when calculating

the spectrum in Figure 4.1 by adding in the Frenkel exciton Hamiltonian for each C8S3

molecule a random shift relative to the gas-phase excitation energy taken from the

4

corresponding (IW or OW) Gaussian distribution.

Understanding the energy splitting between the two lowest-energy peaks. As- suming that the underlying lattice is the same for the IW and OW, it is not possible to explain the energy splitting between the lowest-energy optical bands of both walls, polarized parallel to the tube’s axis. This is because the energy position of the parallel component of the spectrum of a cylindrical aggregate is essentially independent of both the rolling vector

and the tube’s radius.

Moreover, the above calculations show that systematic differences in solvent shifts do not occur and also shifts induced by structural disorder can not explain the energy splitting.

Thus, we are left with the explanation that the energy splitting arises from a slightly different packing density in both walls. A higher density in the IW, as compared to the OW, would ensure stronger excitonic couplings

and therefore a larger red-shift of the lowest-energy absorption peak for the IW. To test this hypothesis, we construct starting aggregate structures assuming different densities of molecules in the two walls. It turns out that enlarging the in-plane unit cell lattice constants for the OW by 2.4

compared to the IW (see Table 4.1 in the Appendix) is enough to match the energy splitting in the experimental spectrum. Our MD simulations performed on multiple preassembled structures confirmed that such structures are stable and do not disintegrate. Only because in the above we ruled out relative shifts between both walls induced by surroundings and structural disorder, are we able to conclude on this small difference in density.

Our microelectrostatic calculations exclude the hypothesis of the energy splitting

being due to different excitation energies in the IW and OW as a result of different

electrostatic environments, used in the phenomenological modeling of Ref. 92. It is

instructive, however, to disentangle the different contributions to the distributions

of Figure 4.4c. In particular, removing the solvent (water and counterions) from the

microelectrostatic calculations results in a shift of both distributions toward more

stabilizing energies as compared to the full calculation (Figure 4.4c, dashed lines). A

relative shift between the IW and OW also appears, with the OW

less stabilizing

than the IW. Further investigations show that this shift between the two walls is

induced by a different electrostatic interaction with the polar heads of the molecules

(see Figure 4.8b). This is likely related to the inward vs. outward orientation of the

polar groups. This relative shift is compensated by the solvent. Interestingly, the effect

of the water is found to be the same for IW and OW molecules, implying that the

compensation of the shift between the walls observed without solvent comes from the

Na

counterions (see Figure 4.8a).

4

-1000 -500 0 500 1000 1500

∆ ϵ (cm⁻¹) 0

1 2 3 4 5

Probability Density Function

×10^-3

IW (no solvent), σ = 206 cm^-1 OW (no solvent), σ = 118 cm^-1 IW, σ = 213 cm^-1

OW, σ = 231 cm^-1

(a) (b) (c)

Figure 4.4 | Energetic disorder and energy shifts computed using microelectrostatic calculations. Based on the MD trajectory, microelectrostatic calculations have been performed on clusters of molecules positioned around a central molecule residing either in the outer (a, cyan) or inner (b, blue) walls. The solvent (water in red and white, and Na⁺in blue) is treated explicitly. The central molecule’s ground and excited state charge distributions interact differently with the environment. The difference between the ground and excited state electrostatic and polarization interaction energies, i.e., the solvent shift¢≤, is plotted in (c) for molecules residing in the inner (cyan) and outer (blue) walls. Notably, discarding electrostatic interactions with the solvent, a different¢≤in the two walls is induced by a different electrostatic interaction with the polar heads of the molecules as a result of the curvature (dashed lines). This effect cancels out, however, when also accounting for the interactions with the solvent (solid lines).

Conclusions

We have applied a multiscale approach to reveal structure—optical property relation- ships of a complex self-assembled structure of thousands of molecules, namely the double-walled tubular aggregates of C8S3 dye molecules. By combining (i) molecular dynamics simulations on preassembled structures to test the stability of the nanotubes and to retrieve realistic structural and corresponding interaction disorder, (ii) micro- electrostatic calculations to probe different electrostatic environments experienced by the molecules of the system and to retrieve solvent shifts and energy disorder, (iii) translating aggregate structures obtained from the molecular dynamics simulation into a Frenkel exciton Hamiltonian to model optical signatures, and (iv) using the comparison with experiment to guide modeling and validate the obtained structures, we found optimal structural parameters for the complex molecular packing of these systems. The absorption spectrum of the proposed optimal structure correctly predicts the optical features of the nanotube: namely the order of the optical bands, their relative oscillator strengths and polarization properties, and the optical bandwidths.

We explained the origin of the energy splitting between the inner and outer walls, which results from a slightly different molecular packing in both walls.

The proposed multiscale approach is a powerful tool to obtain understanding of the

structural and optical relationships that can be used to reveal design rules for creating

new materials. The obtained models with realistic structural parameters and disorder

can be further used for detailed studies of the exciton dynamics in these chlorosome

4

mimicking systems. The demonstrated methodology will prove useful for studying the structure and structure-property relationships for other complex supramolecular systems.

4.3. Methods

Molecular dynamics simulation of preassembled double-walled tubes. A set of preassembled double-walled structures was constructed by independently rolling a two-dimensional lattice, based on the adapted unit cell of a similar molecule, onto two cylindrical surfaces as described in Section 4.4.1 of the Appendix. MD simulations of preselected preassembled double-walled nanotubes were conducted using the Gromacs 5.1.4 simulation package

with the GAFF force field.

Preassembled systems were solvated in orthogonal boxes with periodic boundary conditions using the TIP3P water model.

Moreover, to neutralize the excess negative charge of C8S3 molecules, Na

ions were added. Initially, each system was minimized and equilibrated in NVT and NPT conditions followed by the production phase. The production phase was conducted in the NPT ensemble and lasted at least 20 ns or longer. It was found that it is important to initially place the Na

ions in a ratio close to 1/1 near the negatively charged C8S3 molecules, i.e., the number of Na

ions near each wall must be close to the number of C8S3 molecules in that wall. A charge imbalance tends to lead to considerable stress in the tube and may cause significant defects in the packing. Additional information for the MD parameters and the system preparation details are provided in Section 4.4.2 of the Appendix.

Microelectrostatic calculations. The energetic disorder and the relative energy shifts for the IW and OW were computed by means of a microelectrostatic scheme (see Ref. 62 for a recent review). Namely, the classical energy expression of the Direct Reaction Field (DRF) approach was used, as implemented in the DRF90 software.

In such a polarizable classical description, molecules were thus described by atomic point charges and atom-centered isotropic polarizabilities. The system can be thought of as being composed of two subsystems: a central molecule and its surroundings. The energy shift of a single C8S3 molecule was obtained by performing two calculations:

one where the C8S3 molecule is described by the charge distribution of its ground state,

and one where it is described by the charge distribution of its brightest excited state. In

both calculations, the C8S3 molecules were surrounded by a polarizable surroundings

comprising all the C8S3 and solvent (water, Na

) molecules within a radius of 3.5 nm

from the center of geometry of the C8S3 molecule (see Figure 4.4a,b). The surrounding

molecules were described by their ground state charge distributions. From these

two calculations, interaction energies between the central molecule represented by its

ground state (

) or excited state (

) charge distribution and its surroundings are

obtained. The energy shift,

, was then obtained as the difference (

). We refer to

4

Section 4.4.3 of the Appendix for further details on the atomic charges, polarizabilities, radius of the surroundings, and contributions to

.

Absorption spectra calculations. Electronic excited states of the system are described by a Frenkel exciton Hamiltonian

that takes into account dipole-dipole intermolecular excitation interactions.

The expression for the Hamiltonian is given in eq. (4.4). Exciton states are obtained by numerical diagonalization of the Hamiltonian. Homogeneous spectra are obtained from the ideal—cylindrically symmetric—structures assuming the same excitation energies for all molecules. Spectra of the disordered structures are obtained by averaging over the snapshots along the MD trajectory, where the structural parameters are directly translated into couplings in the Hamiltonian, and energy disorder is included by adding Gaussian random numbers to the gas-phase molecular excitation energies. The mean value and the width of the Gaussian distribution were estimated by microelectrostatic calculations as described above. For the final spectra, lifetime broadening was added in order to match with the experimental spectrum according to the procedure described in Ref. 8.

Further details are given in Section 4.4.4 of the Appendix.

Author contributions

I.P. and A.H.d.V. performed molecular dynamics simulations. A.S.B. performed optical

modeling. R.A. performed microelectrostatic calculations. P.C.T.S. provided key

insights for the molecular dynamics simulations. A.S.B., I.P., and R.A. analyzed the

data under the supervision of T.L.C.J., A.H.d.V., S.J.M, and J.K.

4

4.4. Appendix: Additional information

4.4.1. Obtaining the preassembled structures

Structural restraints

The available experimental data and previous modeling based on phenomenological studies allow us to set restraints for the construction of the initial set of tubes.

For the tube’s radius, we use values obtained from the cryo-TEM scan profile,

which are equal to 3.2

0.5 nm and 6.5

0.5 nm for the sulfonate groups of the inner and outer walls, respectively. For all the tubes, we use the optimal wall-to-wall thickness of 1.7 nm measured from the cores of the molecules, which is based on the systematic molecular dynamics (MD) study.

Moreover, we adhere to the herringbone arrangement proposed for the underlying molecular packing in this aggregate.

Using phenomenological modeling, a herringbone lattice with two molecules per unit cell results in a Davydov splitting of the optical bands, which explains well the complex optical spectrum of the C8S3 aggregate

and the optical changes for a modified version of the aggregate with an increased radius.

By contrast, a brickwork or staircase lattice with one molecule per unit cell cannot explain the multitude of peaks in the absorption spectra of these aggregates.

Thus, as a starting point we take a unit cell with the herringbone arrangement for a similar molecule—namely, TTBC—obtained crystallographically.

The unit cell has the following lattice parameters:

Å,

Å,

Å,

,

and

. Keeping the arrangement of the cores of the chromophores unchanged and replacing the side chains of TTBC for those of C8S3, we obtain a starting unit cell for our modeling (Figure 4.5a). Preliminary analysis of the arrangement of the C8S3 monomers showed that for optimal arrangement, an increase of the length of the unit cell (parameter

) was necessary. Consequently, we used

= 12.036 Å and one of the monomers was translated by 1.0 Å in the

-direction to accommodate the change in the unit cell.

Rolling of the herringbone lattice on the cylindrical surface

From the modified unit cell we construct two-dimensional herringbone lattices (Figure 4.5b). Full details are described in a separate publication,

but briefly, we used an in-house python script to take the atoms in the unit cell constructed as described in the previous section to eventually build a double-walled structure. To create a two-dimensional lattice, the unit cell is replicated along

and

directions

and

times, respectively. We wrap this lattice along the rolling vector

, by connecting its beginning and end points, onto a cylinder resulting in a tubular structure (Figure 4.1c).

The circumference of this tube is equal to the length of the rolling vector

, and is

related to the tube’s radius

via the relation

. The value of

, taken within

the experimental errors,

therefore, restricts the range for the values of the length

of the rolling vector. We construct a set of tubes with varying rolling angles (

), as

4

<latexit sha1_base64="IoELSitFJaTQ4WT4pr8f01q0csw=">AAAB7XicbVDLSgNBEJyNrxhfUY9eBoPgKexGQY9BLx4jmAckS+idzCZj5rHMzAphyT948aCIV//Hm3/jJNmDJhY0FFXddHdFCWfG+v63V1hb39jcKm6Xdnb39g/Kh0cto1JNaJMornQnAkM5k7RpmeW0k2gKIuK0HY1vZ377iWrDlHywk4SGAoaSxYyAdVKrNwQhoF+u+FV/DrxKgpxUUI5Gv/zVGyiSCiot4WBMN/ATG2agLSOcTku91NAEyBiGtOuoBEFNmM2vneIzpwxwrLQrafFc/T2RgTBmIiLXKcCOzLI3E//zuqmNr8OMySS1VJLFojjl2Co8ex0PmKbE8okjQDRzt2IyAg3EuoBKLoRg+eVV0qpVg4tq7f6yUr/J4yiiE3SKzlGArlAd3aEGaiKCHtEzekVvnvJevHfvY9Fa8PKZY/QH3ucPiDmPGQ==</latexit>

<latexit sha1_base64="EpwadKl79nuFFVBzWqBryCC0A38=">AAAB7HicbVBNS8NAEJ34WetX1aOXxSJ4KkkV9Fj04rGCaQttKJvtpl262YTdiVBCf4MXD4p49Qd589+4bXPQ1gcDj/dmmJkXplIYdN1vZ219Y3Nru7RT3t3bPzisHB23TJJpxn2WyER3Qmq4FIr7KFDyTqo5jUPJ2+H4bua3n7g2IlGPOEl5ENOhEpFgFK3k90KOtF+pujV3DrJKvIJUoUCzX/nqDRKWxVwhk9SYruemGORUo2CST8u9zPCUsjEd8q6lisbcBPn82Ck5t8qARIm2pZDM1d8TOY2NmcSh7YwpjsyyNxP/87oZRjdBLlSaIVdssSjKJMGEzD4nA6E5QzmxhDIt7K2EjaimDG0+ZRuCt/zyKmnVa95lrf5wVW3cFnGU4BTO4AI8uIYG3EMTfGAg4Ble4c1Rzovz7nwsWtecYuYE/sD5/AHFJI6o</latexit> <latexit sha1_base64="+wSBPeL8nxBdvzPXA2qswhGhfpg=">AAAB7XicbVBNS8NAEJ3Ur1q/qh69LBbBU0mqoMeiF48V7Ae0oUy2m3btZhN2N0IJ/Q9ePCji1f/jzX/jts1BWx8MPN6bYWZekAiujet+O4W19Y3NreJ2aWd3b/+gfHjU0nGqKGvSWMSqE6BmgkvWNNwI1kkUwygQrB2Mb2d++4kpzWP5YCYJ8yMcSh5yisZKrR6KZIT9csWtunOQVeLlpAI5Gv3yV28Q0zRi0lCBWnc9NzF+hspwKti01Es1S5COcci6lkqMmPaz+bVTcmaVAQljZUsaMld/T2QYaT2JAtsZoRnpZW8m/ud1UxNe+xmXSWqYpItFYSqIicnsdTLgilEjJpYgVdzeSugIFVJjAyrZELzll1dJq1b1Lqq1+8tK/SaPowgncArn4MEV1OEOGtAECo/wDK/w5sTOi/PufCxaC04+cwx/4Hz+AIzPjxw=</latexit>

c

b

a

cN

<latexit sha1_base64="JPLBV0RsTRj3PLGejJ64lbWd7/8=">AAAB63icbVBNS8NAEJ3Ur1q/qh69LBbBU0mqoMeiF09SwX5AG8pmu2mX7m7C7kYooX/BiwdFvPqHvPlv3KQ5aOuDgcd7M8zMC2LOtHHdb6e0tr6xuVXeruzs7u0fVA+POjpKFKFtEvFI9QKsKWeStg0znPZiRbEIOO0G09vM7z5RpVkkH80spr7AY8lCRrDJJHI/JMNqza27OdAq8QpSgwKtYfVrMIpIIqg0hGOt+54bGz/FyjDC6bwySDSNMZniMe1bKrGg2k/zW+fozCojFEbKljQoV39PpFhoPROB7RTYTPSyl4n/ef3EhNd+ymScGCrJYlGYcGQilD2ORkxRYvjMEkwUs7ciMsEKE2PjqdgQvOWXV0mnUfcu6o2Hy1rzpoijDCdwCufgwRU04Q5a0AYCE3iGV3hzhPPivDsfi9aSU8wcwx84nz/U544Z</latexit>

bN

_b

(a) (b)

C

<latexit sha1_base64="k/AN3iXODITolNbwP576Vc7DU5I=">AAAB7XicbVDLSgNBEOz1GeMr6tHLYBA8hd0o6DHoxWME84BkCbOT2WTM7Mwy0yuEkH/w4kERr/6PN//GyQPRxIKGoqqb7q4olcKi7395K6tr6xubua389s7u3n7h4LBudWYYrzEttWlG1HIpFK+hQMmbqeE0iSRvRIObid945MYKre5xmPIwoT0lYsEoOqnexj5H2ikU/ZI/BfkhwSIpwhzVTuGz3dUsS7hCJqm1rcBPMRxRg4JJPs63M8tTyga0x1uOKppwG46m147JqVO6JNbGlUIyVX9PjGhi7TCJXGdCsW8XvYn4n9fKML4KR0KlGXLFZoviTBLUZPI66QrDGcqhI5QZ4W4lrE8NZegCyrsQll5eJvVyKTgvle8uipXreRw5OIYTOIMALqECt1CFGjB4gCd4gVdPe8/em/c+a13x5jNH8Afexzems48t</latexit>

Figure 4.5 | Construction of the herringbone lattice. (a) Unit cell with two molecules. (b) Herringbone lattice.

this angle is not known from the experiment. The values of

and

acquire discrete values restricted by the dimensions of the unit cell and are defined by the the following equations:

R = 1 2º

q

(bN_b)²+ (cNc)²

(4.1)

cosµ =cNc

2ºR

(4.2)

From the obtained single-walled tubes, we construct a set of double-walled tubes (making sure the thickness is equal to

nm) with length

nm, resulting in aggregates containing about

C8S3 monomers. The rolling angle is varied between 20

and 55

.

Structural parameters of the inner and outer walls of the model tube used in Figure 4.1 are summarized in Table 4.1.

Table 4.1 | Structural parameters used for the model in Figure4.1.

a

, Å

, Å

, Å R, nm

,

IW 22.547 12.036 13.375 3.724 31.0

OW 22.547 12.324 13.696 5.442 32.7

4

4.4.2. MD Simulations on the preassembled structures

MD protocol

MD simulations of preselected preassembled double-walled nanotubes were performed using the Gromacs 5.1.4 simulation package

with the GAFF force field.

To neutralize the excess negative charge of C8S3 molecules, Na

ions were added. A random placement of ions in the simulation box (using, for example, the Gromacs tool gmx solvate) placed ions disproportionally outside the double-walled structure; the consequent charge imbalance tends to lead to considerable stress in the tube and may cause significant defects in the packing. Instead, Na

ions were distributed within cylinders, one for the ions neutralizing the IW, one for the ions neutralizing the OW, with a ratio of Na

/C8S3 close to 1/1 for each wall. The sizes of the cylinders were chosen to ensure that the Na

ions were near the C8S3 head groups. After counter ion placement, the systems were solvated (by using the Gromacs tool gmx solvate) in orthogonal boxes with periodic boundary conditions using the TIP3P water model.

After solvation, the systems were minimized with the Steepest Descent algorithm for 1000 steps and subsequently equilibrated in the NVT and NPT ensembles (with a time step of 1 fs) for 10 ps and 5 ns, respectively. During the equilibration, the aromatic core of the C8S3 monomers was restrained using harmonic potentials with a force constant of 1000 kJ mol

nm

. In this way, the side chains and the solvent were able to relax before the production phase. For the actual MD simulation, carried out in the NPT ensemble, the time step for integration of the equations of motion was set to 2 fs and snapshots were saved every 100 ps. To prevent the aggregates from crossing the periodic boundary conditions, the systems were translated and rotated to remove the center of mass motion of the nanotube every 100 steps. Production simulations lasted at least 20 ns or longer.

For all MD simulations, the temperature was set to 300 K and kept fixed by using the V-rescale algorithm

with a coupling constant of 0.1 ps. The Berendsen barostat

was used to maintain standard pressure in an isotropic pressure bath with a time constant of 1 ps and compressibility of 4.5 10

bar

. Electrostatic and van der Waals interactions were calculated by using the Verlet cut-off scheme with 1.4 nm cut-off.

The long range interactions were calculated with the reaction field method.

Stability of the tubes

Snapshots of the tubes with two different angles, 23

and 55

, after a 20 ns MD simulation are shown in Figure 4.6.

The stability of the tubes was measured by generating mass density plots of

the sulfur atoms in a direction perpendicular to the long axis of the tubes along

the trajectories. The sulfur density curves have a distinct pattern when the tubular

formation is maintained (Figure 4.7a). The pattern broadens or is even completely lost

when the tubes become disordered or they start to fall apart (Figure 4.7b).

4

(a) (b) (c)

Figure 4.6 | Final MD snapshots of the double-walled tubes after 20 ns. Tubes with rolling angle (a) µ= 23^±and (b)µ= 55^±are shown. (c) Close up of the arrangement of the C8S3 aromatic cores in the two different tubes for the 20 ns frame. Only the aromatic core of the C8S3 molecules (excluding hydrogens) is shown for clarity.

(b) (a)

Figure 4.7 | Mass density profile of the double-walled tubes as a function of simulation time. (a) A stable and ordered tube after 100 ns—the spectra given by this structure is the one reported in Figure4.1d.

(b) A tube with a high level of disorder after 15 ns; the simulation was terminated due to the tube’s disintegration after 15 ns. In both cases, the profiles were centralized and symmetrized along the X-axis.

4

4.4.3. Probing the energetic disorder

Microelectrostatic calculations

General details. The energetic disorder and the relative energy shifts for the IW and OW were computed by means of a microelectrostatic, or induced dipole, scheme.

We employed the classical energy expression of the Direct (or Discrete) Reaction Field (DRF) approach as implemented in the DRF90 software.

In such a polarizable classical description, molecules are described by atomic point charges and atom-centered isotropic polarizabilities (see below for details on the charges and polarizabilities used). Polarizabilities are described according to Thole’s method for interacting polarizabilities,

which avoids numerical instabilities by employing a distance- dependent damping function.

Definition of the energy shift in the excitation energy. The interaction of a C8S3 molecule with its environment generally will be different for its ground and excited states. This leads to a shift in its excitation energy away from the gas phase value (note that we are considering vertical excitations—i.e., no structural rearrangement takes place upon excitation—as those determine the absorption spectrum). We can think of the system as composed of two subsystems: a central molecule and its surrounding (Figure 4.4a,b). We represent the ground and excited states of the central molecule with a different set of charges (see Section 4.4.3). This means that what can change is 1) the electrostatic interaction energy, i.e., the (static) Coulomb interaction between the charges on the central molecules and the surrounding charges, and 2) the polarization interaction energy due to the mutual polarization of the central molecule and its surroundings. We compute the electrostatic (

) and polarization (or induction,

) interaction energies between a central molecule and its surrounding for both the ground state (

) and the (brightest) excited state (

) of the central molecule. The difference

E_es^(tot)° E^(tot)gs = (E^(pol)es + Ees^(el)) ° (Egs^(pol)+ E^(el)gs ) = ¢≤

(4.3) yields the shift (

) to be applied to the gas phase excitation energy.

Input parameters. Calculations are performed with DRF90—the program is avail-

able at http://www.marcelswart.eu/drf90/index.html—the classical implementation

of the DRF method.

The input requires the following parameters: i) the atomic

coordinates, which are obtained from the MD simulation; ii) the atomic charges for the

ground and excited states, which are obtained from quantum chemical calculations

as described in Section 4.4.3 below; iii) the atomic polarizabilities, which are taken

from the standard DRF set of (effective) atomic isotropic polarizabilities, that is a set

of atomic polarizabilities parametrized on the basis of a large set of experimental

molecular polarisabilities;

iv) the radius of the region surrounding the central

molecule, which was set to 30 Å for all the production calculations (see also below for

a discussion of the radius of the surrounding region).

4

In the calculations, the system is considered to be divided in two subsystems: the central molecule and its surroundings. A slab of a C8S3 nanotube of about 10 nm (from the middle) is considered. This comprises a total of approximately 400 and 500

“central” molecules in the inner and outer walls, respectively. Selection of the slab, selection of the surroundings of each central molecule, and preparation of DRF90 inputs are done via a Python script which makes extensive use of the MDAnalysis

library. The surrounding region is spherical, centered at the center of geometry of the central molecule (see Figure 4.4a,b). The effect of the radius of the surrounding region is investigated in Figure 4.9 and associated discussion. For each “central” molecule, two inputs are generated, with the central molecule bearing either the charges of the ground state or of the excited state. Molecules in the surrounding are described in both cases by the ground state charge distribution.

Contributions to the energy shift. Figure 4.4c shows a blue-shift of the two distribution upon solvation, but the question remains open whether the shift is caused by the water molecules or the Na

ions. To investigate this, we carried out additional calculations where the electrostatic interaction with the water was accounted for, but not with the Na

ions. Figure 4.8a shows the distribution of energy shifts for the same snapshot without any solvent, with water, and with the full solvent configuration, i.e., water and Na

ions. The results show that the effect of the water is the same for IW and OW molecules and that the compensation of the shift between the walls observed without solvent comes from the Na

counterions.

-1000 -500 0 500 1000 1500

∆ϵ(cm⁻¹) 0

1 2 3 4 5 6 7

Probability Density Function

×10^-3

IW (no solvent), σ = 206 cm^-1 OW (no solvent), σ = 118 cm^-1 IW (H₂O), σ = 248 cm^-1 OW (H₂O), σ = 238 cm^-1 IW, σ = 213 cm^-1 OW, σ = 231 cm^-1

(a)

-1000 -500 0 500 1000 1500

∆ϵ(cm⁻¹) 0

1 2 3 4 5 6 7

Probability Density Function

×10^-3

IW (core off), σ = 176 cm^-1 OW (core off), σ = 118 cm^-1 IW (SO₃^- off), σ = 120 cm^-1 OW (SO₃^- off), σ = 86 cm^-1 IW (no solv), σ = 206 cm^-1 OW (no solv), σ = 118 cm^-1

(b)

Figure 4.8 | Contributions to the energy shifts. (a) Energy shift distributions without any solvent (dashdotted line), with water (dashed), with both water and Na⁺ions (solid). (b) Distributions obtained excluding the solvent and when switching off the charges of the aromatic cores (dashdotted line), or of the polar head groups (dashed); the full calculation (excluding the solvent) is also shown (solid).

It remains to be seen what causes the initial shift which is present when excluding

the water and ions. To investigate this, the following calculations were performed: one

where the charges on the aromatic core of all but the central molecule were set to zero

4

(“core off”), and and one where the same was done for the polar head groups (“SO

off”). Figure 4.8b shows that when charges on the polar head groups of the molecules in the surrounding are set to zero, the relative shift of the IW and OW distributions is inverted. Thus, the stabilization of the IW excitation energies relative to the OW ones in the absence of the solvent is induced by an electrostatic interaction with the polar heads of the C8S3 molecules.

Convergence of microelectrostatic calculations. Figure 4.9a shows how the distri- bution and its Gaussian fit change upon varying the number of molecules considered in the calculations. It can be seen that already at 100 molecules the obtained Gaus- sian distribution is qualitatively converged if we take as reference the one with 350 molecules. In particular, we note how the standard deviation,

, is very robust and reasonably converged already at 100 molecules. However, the mean value of the distribution,

, is found to be very sensitive to little variations in the distribution, as can be seen from Figure 4.9a where practically overlapping fits give rise to changes in

¢≤

by 10-20 cm

. All the production calculations consider at least 350 molecules per distribution.

Figure 4.9b shows the effect of the radius of the surrounding region, varied from 20 to 30 Å. The distributions are very similar in all cases, and particularly between 25 and 30 Å. Again the

of the distributions converge quickly, while

is again found to be more sensitive. All production calculations have been run with a radius of 30 Å.

-1000 -500 0 500 1000 1500

∆ϵ(cm⁻¹) 0

0.5 1 1.5 2 2.5 3 3.5 4

Probability Density Function

×10^-3

IW, N=100, σ = 214 cm^-1 OW, N=100, σ = 233 cm^-1 IW, N=250, σ = 214 cm^-1 OW, N=250, σ = 233 cm^-1 IW, N=330, σ = 213 cm^-1 OW, N=330, σ = 234 cm^-1

(a)

-1000 -500 0 500 1000 1500

∆ϵ(cm⁻¹) 0

0.5 1 1.5 2 2.5 3 3.5 4

Probability Density Function

×10^-3

IW, 2.0 nm, σ = 217 cm^-1 OW, 2.0 nm, σ = 237 cm^-1 IW, 2.5 nm, σ = 213 cm^-1 OW, 2.5 nm, σ = 235 cm^-1 IW, 3.0 nm, σ = 213 cm^-1 OW, 3.0 nm, σ = 231 cm^-1

(b)

Figure 4.9 | Convergence of the Gaussian parameters extracted from the microelectrostatic calculations.

The convergence is shown as a function of (a) the number of molecules per wall for which the energy shift is computed and (b) the size of the surrounding region around the central molecule.

Electronic structure calculations

The only input parameters required by DRF90 not readily available are the atomic

charges of the ground and excited states. Charges which reproduce the electro-

static potential and the dipole moment of the molecule have been computed via the

CHELPG

scheme as implemented in Gaussian.

This is done for both the ground

4

state and the (vertical) excited state with the highest oscillator strength, as obtained via

B97XD

/6-31G(d,p)

density functional theory (DFT) and time-dependent DFT (TDDFT)

calculations, respectively. Further details on the TDDFT calculations and excited state analysis are reported in Table 4.2 and discussed below.

Excited state calculations. Table 4.2 collects the results for the TDDFT calculations in gas phase. The state-of-the-art

B97XD density functional, containing a long-range correction for the self-interaction error and a dispersion correction, is compared to results obtained with the B3LYP functional. Results obtained with

B97XD show that the brightest excited state, i.e., the excited state with the largest oscillator strength, is the lowest singlet excited state.

Table 4.2 | Results on the brightest excited state of C8S3 at the TDDFT level. The geometry was optimized at the B3LYP/6-31G(d,p) DFT level. The basis set employed for all the TDDFT calculations is 6-31G(d,p).µ is the transition dipole moment.

functional DFT Bright State Orbitals (nm) (cm Energy

) (au)

(D)

f. osc.

B3LYP S3 HOMO

2

LUMO 434.7 23,002 22.2741 12.0 1.56

!

B97XD S1 HOMO

LUMO 407.5 24,539 23.0113 12.2 1.71 HOMO

2

LUMO

Atomic charge distributions. CHELPG charge distributions reproducing the electrostatic potential and dipole moment have been obtained for the ground,

, and brightest excited state,

, at the

B97XD/6-31G(d,p) level. The difference between these two set of charges,

, is computed and plotted on the optimized molecular structure of C8S3 in Figure 4.10. As expected,

the excitation involves rearrangements of the electron density along the bridge connecting the two benzimidazole rings.

4.4.4. Absorption spectra calculation

The excitonic systems of inner and outer walls are assumed decoupled.

We use a Frenkel exciton Hamiltonian to describe the collective excited states of the excitonic system of either of the walls:

H =X^N

n=1

(≤₀+ ¢≤n)|nihn| + X^N

n6=m

J_nm|nihm|,

(4.4)

where

is a molecular gas-phase excitation energy;

is the fluctuation of the gas-phase excitation energy on molecule

obtained from a Gaussian distribution of width

, where the value of

is estimated from the microelectrostatic calculations;

|ni

is the state where only molecule

is excited and all other molecules are in the

ground state. We use

cm

.

is the intermolecular resonance (excitation

4

q (e)

Figure 4.10 | Difference in CHELPG charge distribution between the brightest excited state and the ground state. Ground and excited state charge distribtions have been obtained at the!B97XD/6-31G(d,p) (TD)DFT level via the CHELPG scheme. The charge shown on thei-th atom is obtained asqⁱ_es° qⁱg s, where qⁱ_esandqⁱ_{g s}indicate the charges on thei-th atom obtained by fitting the electrostatic potential around the molecule given by the density of the (brightest) excited state and of the ground state, respectively.

transfer) interaction between the molecules

and

.

To calculate

, we use the extended transition dipole model,

where the transition dipole moment of the molecule is considered as a dipole of length

with two point charges

and

such that the magnitude and orientation of the relevant molecular transition dipole moment is reproduced. In this approximation, the transfer interaction

is explicitly given by:

Jnm= Aµ² l²

∑ 1 r_nm⁺⁺° 1

r_nm^+°° 1 r_nm^°+ + 1

r_nm^°°

∏

,

(4.5)

with

r_nm^±±= |rnm± l(ˆen° ˆem) 2 |, r_nm^°+= |rnm° l(ˆen+ ˆem)

2 |, r_nm^+°= |rnm+ l(ˆen+ ˆem)

2 |,

(4.6)

where the transition dipole moment is related to the charge

and dipole length

(

), and

. We use the previously reported values

and

nm for the point charges and length that give rise to the magnitude of the transition dipole moment of 11.4

. The molecular position vector

and unit vector

that represents the direction of the transition dipole moment of molecule

,

, are obtained from the structure from each saved snapshot of the MD trajectory, where

en

is the unit vector between the carbon atoms of the imidazole rings of the C8S3

4

molecule (Figure 4.1a), and

is the coordinate of the carbon atom of one of the imidazole rings. The constant

cm

nm

Debye

comes from a conversion of units allowing one to express transition dipole moments, distances, and energies in Debye, nm, and cm

, respectively.

Linear absorption spectrum. After constructing the Hamiltonian of eq. 4.4, and numerically diagonalising it, we obtain the exciton states:

|qi =X

n 'qn|ni,

(4.7)

where

denotes the amplitude of the

th state on molecule

. The corresponding eigenvalue of the state is the energy

. This calculation is repeated for every snapshot of the MD simulation. For each snapshot, the molecular energy shifts are taken randomly for each molecule from a Gaussian distribution.

The general form of the linear absorption spectrum for an isotropic solution, obtained from the Fermi golden rule,

is given by

A(!) =DDX

q Oq±(! ° Eq)EE

,

(4.8)

with the oscillator strength of the state

given by

Oq=X

n,m'_qn'^§_qmh(µn· e)(µm· e)i.

(4.9) The double angular brackets

imply an average over snapshots from the MD trajectory, together with the single brackets

representing an average over the orientations of the cylinder relative to the electric polarization vector

of the incident linearly polarized electromagnetic wave. The averaged spectrum captures instantaneous fluctuations of the excitation transfer interactions

due to structural as well as energetic disorder and results in broadening of the modeled spectrum.

For Figure 4.1, the calculated transitions are additionally broadened with Lorentzian

lineshapes to facilitate the comparison with the experimental spectrum. For the lowest

energy Davydov split optical parallel and perpendicular polarized bands of the inner

wall, a homogeneous broadening with a full width at half maximum of 40 cm

and

1000 cm

, respectively, was used. For the outer wall, we use 110 cm

and 1000 cm

,

respectively for the lowest optical parallel and perpendicular polarized bands.