Microstructural Model of Indacenodithiophene-co-benzothiadiazole Polymer: π-Crossing Interactions and Their Potential Impact on Charge Transport

Morphological and electronic properties of indacenodithiophene-co-benzothiadiazole (IDTBT) copolymer with varying molecular weights are calculated through combined molecular dynamics (MD) and quantum chemical (QC) methods. Our study focuses on the polymer chain arrangements, interchain connectivity pathways, and interplay between morphological and electronic structure properties of IDTBT. Our models, which are verified against GIWAXS measurements, show a considerable number of BT-BT π–π interactions with a (preferential) perpendicular local orientation of polymer chains due to the steric hindrance of bulky side chains around IDT. Although our models predict a noncrystalline structure for IDTBT, the BT-BT (interchain) crossing points show a considerable degree of short-range order in spatial arrangement which most likely result in a mesh-like structure for the polymer and provide efficient pathways for interchain charge transport.


S1. Force Field Parametrisation and Model Details
Force field parameters of the backbone and sidechains of the polymer were obtained separately.OPLS force field was used for bonded and non-bonded parameters of the side chains.The backbone parameters were derived from optimised structure of IDT-BT monomer with two methyl groups attached to each sidechain connection point of IDT (in total four methyl groups representing the place of four side chains that will be attached later) and the monomer was capped with IDT and BT fragments (see Figure S1, the capping fragments are marked with blue and red ellipses).Note that all DFT calculations were performed by Gaussian 16 by using B3lyp hybrid functional and cc-pVTZ basis set for all steps of force field calculations.The monomer optimisation process is summarised in "monomer optimisation".The backbone non-bonded and bonded parameterisation methods are detailed in "non-bonded parameters" and "bonded parameters", respectively.The sidechain attachment procedure and the backbone-sidechain connection point parametrisation are explained in section "sidechain attachment procedure"."Polymer models" summarises the resulting IDT-BT polymers.

Monomer Optimisation
For the model monomer shown in Figure S1, all possible Cis-Trans isomers have been made (see Figure S2) and their total energies of optimised structures were calculated and the lowest energy structure (7) was used as a starting point for the other component of the force field generation.

Non-bonded Parameters Lennard Jones parameters
The Lennard jones parameters of atoms in the polymer backbone were directly taken from their analogues in OPLS force field.

Atomic Charge of Backbone (Conjugated Part)
The atomic charges were calculated by using CHELPG scheme (developed by Breneman and Wiberg, J. Comp.Chem.1990, 11, 361) on the optimised structure of the monomer shown in Figure S1.The resulting point charges from CHELPG computation for one repeat unit (RU, i.e., the monomer structure after removing the capping (BT and IDT) fragments, see Figure S3) are shown in table S1.The sum of point charges for the whole monomer was zero and, after removing the capping BT and IDT fragments, the total charge of one RU was calculated -0.004 e.The extra charge of the RU was redistributed on all atoms (equally) so that the total charge of each RU in the polymers equals to zero.Table S1.Atomic charges as calculated by CHELPG method and as used for the repeat unit force field (RU-FF).Note that the total charge of each repeat unit is set to zero by redistribution of the total excess charge calculated by CHELPG.

Bonds and Angles
Harmonic potentials in the forms of equations SE1 and SE2 were used for implementing bond and angle interactions, respectively, in the force field.
where Vb and Va are bond and angle potentials, k b and K θ are the force constants representing the stiffness of the bond and angle, and bij and θijk are the bond and angle equilibrium values taken from DFT-optimised monomer for each combination of bonded atoms.It is important that the distances and angles are consistent with DFT calculations because the MD simulations will be used to generate models for DFT calculations.Note that employing force field like GAFF that ignore the pi-bonding may introduce errors.
The stiffness of all bonds (k ij b ) and angles (k ijk θ ) were set to 320,000.0 kJ/mol/nm 2 and 500.0kJ/mol/rad 2 .

Intra Fragment Torsional Potentials
The torsional potentials for internal dihedral angles of IDT and BT fragments in the RU were implemented in form of a Ryckaert-Bellemans function (equation SE3) for all dihedral angles around sp2 hybridised heavy atoms.These fragments are flat and rigid due to their (partial) double bond nature.A representative and well-studied example for a flat and sp2 hybridised molecule is the benzene ring; thus, the OPLS constants of Ryckaert-Bellemans function for carbon atoms of benzene molecule was used for these dihedral potentials.It is worth noting that benzene is an ideal molecule if one wants a single torsional potential for conjugated fragments as each bond is 50% double and 50% single.The Ryckaert-Bellemans reads as:

SE3
where Vrb is the torsional potential of the dihedral angle between the planes of ijk and jkl atoms.Note that ψ=ϕ-180 and C n are the six constants of the function.The Cn [kJ/mol] values of the C-C-C-C dihedral of the benzene ring from OPLS force field are shown in Table 2.For the dihedral angles in which the sp3 hybridised carbons (the connecting carbon to the sidechain) take part, equation S4 was used (proper dihedral type 1 in GROMACS).Note that all equilibrium dihedral angle values (  in equation S4) were directly taken from the optimised monomer.
k ϕ of 10 kJ/mol was used for all dihedrals around the sp3 carbons connecting the backbone to the sidechains.Note that ϕ= 0 is corresponding to the cis configuration (i.e., ijk and jkl on the same side).

Inter Fragment Torsional Potential (Connecting Point between IDT and BT)
Torsional potential (V DFT ) of dihedral angle between IDT and BT (IDT-BT) in the monomer structure (marked with red circles in Figure S4 a) was calculated by DFT (via B3LYP/cc-pvtz) through a dihedral scan with 10° spacing (in total 37 scans from -180° to 180°).Figure S4 b shows the calculated torsional potential as a function of dihedral angle.It should be noted that V DFT shows the total potential energy of the monomer structure at each of 37 dihedral angle points.Therefore, after implementing this torsional potential correctly as a force field parameter, the total potential energy of the monomer at each dihedral angle as calculated by the force field should match the V DFT (IDT-BT).To this end, we used the parametrisation scheme as explained below.
The DFT-optimised structures at each IDT-BT (i.e., -180, -170, …, 170, 180) were obtained and used as the input coordinate file.Then, for each structure, an energy minimisation based on steepest descent algorithm using the generated force field parameters (excluding the torsional potential of IDT-BT) with a stiff dihedral restraint (30,000 kJ/mol/rad), which ensures IDT-BT remains reasonably constant (< 1° fluctuations) during minimisation, was performed.The total energy (V FF ), excluding the energies of the restrained dihedral, for each structure (in total 37 values) were calculated after energy minimisation.Figure S5 a shows V FF values for all 37 structures.Accordingly, the torsional potential for the force field will be V CORR = V DFT -V FF .Figure S5 b shows the correct torsional potential, as it is given to the force field in the format of a table potential for MD simulations.The dihedral table was provided as a three-column table.The first column shows angles (-180 to and including 180, 1-degree spacing is recommended), the second column is the potential value (a cubic spline was fitted on 10-degree spline and the potential for every degree was calculated), and the third column represents the negative value of the first derivates of the potential (i.e., force) as obtained from the cubic spline fit.As a final check, we calculated the total energy of the monomer, but this time the force filed torsional potential (V CORR ) was also included.Figures S6 shows the comparison between the total energy (excluding the dihedral restraint energy) from force field (V tot ) and the DFT calculated torsional potential (V DFT ) for IDT-BT.As shown, the total potential energy of the monomer after imposing the implemented torsional potential (V CORR ) is matching the DFT calculated values (V DFT ).
Figure S6.Comparison between DFT calculated (V DFT ) and force field calculated (V tot ) total potential energy at each dihedral angle of ϕIDT-BT.

Sidechain Attachment Procedure
The optimised RU structure (coordinates) and force field parameters have been obtained so far.The next step is to attach side chains to the repeat unit and add the force field parameters accordingly.As earlier mentioned, OPLS force field parameters were used for sidechains.However, the interactions between the backbone and the sidechain at the connection points should be treated correctly.
First, the sidechains are assumed to have a total zero net charge.Each methylene (-CH2-) and the end methyl (CH3) group for each sidechain are parametrised by the united atom parameters of OPLS.Therefore, the total charge of each united atom (and accordingly, the charge of each sidechain) is zero.However, to attach a side chain to each of four methyl groups in the RU, one hydrogen atom of the methyl group should be removed, and its charge should be redistributed on the remaining two hydrogen atoms, see Figure S7.In this way, the total charge of the repeat unit with sidechain (RU-SC) will remain zero.The bonded potentials defined around the connection points also need a reasonable treatment.Figure S8 shows the RU-SC structure.The atom labeling around one of the two connection points are shown.The sp2 hybridised carbon of IDT are shown in green numbers, the sp3 carbon and hydrogen atoms are represented by blue numbers, and the sidechain united atoms are marked with red numbers.
After removing one hydrogen of each methyl group, all the bonded parameters in which this hydrogen was involved were removed from the force field.Then the new bonded potentials were added to the force field as explained here.Based on the atom labeling shown in Figure S8, the bond potential for 16-24 bond was taken from OPLS and added to the force field.Also, all angle potentials for any newly formed angle in which any sidechain united atom exists (e.g., 51-24-16, 24-16-26, 24-16-25, 24-16-10, etc.) were taken from OPLS and added to the force field.In a same way, the potential of all the newly formed dihedrals in which any sidechain united atom exists (e.g., 51-24-16-26, 24-16-10-17, etc.) excluding the ones that sp2 hybridised carbons of the IDT also exist (i.e., 24-16-10-11, and 24-16-10-4) were taken from OPLS and added to the force field.Excluding the two dihedral potentials for each sidechain in which one sidechain united atom (e.g., 24) and one sp2 hybridised carbon atom (e.g., 4 and 11) take part is due to avoiding any distortion in the flatness of the conjugated fragment (IDT) by imposing additional dihedral potentials on it.

Polymer models
For all RU-SC structures as shown in Figure S9, a polymer with DP (degree of polymerisation) =5 was made, i.e., IDT-BT16-5, IDT-BT8-5, IDT-BT4-5, and IDT-BT1-5.For IDT-BT16, polymers with DP=10 and 20 were also made, i.e., IDT-BT16-10 and IDT-BT16-20.The force field of polymers are simply made by repeating the force field of each RU-SC.Note that all polymers were capped with one united atom carbon (with zero charge) at both ends and improper dihedral potentials were used to keep it coplanar with the backbone of the polymer.

S2. Simulation and Analyses Details Simulation Procedure and Parameters
Fully stretched IDTBT chains were randomly inserted into the simulation box and after energy minimisation a rapid contraction of the box were done (under NPT and with P = 1000 bar) to quickly pack the box to the correct density.Then several annealing cycles, as shown in the main manuscript, were performed to relax the structures.Relaxation simulations were performed under NPT condition with time step of 3 fs by using GROMACS.A 1.0 nm cutoff for Lennard Johns and electrostatic interactions was used and all nonbonded interactions for 1-2 and 1-3 bonded pairs were excluded and a scaling factor of 0.5 was used for 1-4 bonded pairs.V-rescale thermostat and C-rescale barostat were used for packing steps and Nose-Hoover thermostat and Parrinello-Rahman barostat were employed for equilibration runs.Verlet cut-off scheme was employed for non-bonded interactions and Particle-mesh Ewald was used for long-range electrostatic interactions.Examples of coordinate, topology, and run files can be found in https://github.com/HMakkiMD/IDTBT.

T g evaluation and annealing procedure setup
IDTBT shows a glass transition temperature Tg of around 510 K in our MD simulations (see Figure S10); therefore, equilibration of initial configurations at service temperature (around 300 K) is not feasible.Besides, there is a high chance of being kinetically trapped in unfavorable configurations if one performs equilibrations well above Tg followed by a (direct) rapid cooling to 300 K, similar to many polymer glasses 1 .Therefore, we used a "sub-Tg relaxation" protocol, as was previously used for similar polymers in refs 2,3 .Thus, we use an intermediate annealing at 500 K (just below Tg) between the equilibration well above Tg (at 900 K) and at 300 K.A schematic of this procedure is depicted in Figure 1b (in the main manuscript).

π-π interaction analysis
Similar to our previously implemented method 2 , we defined the plane spanned by the heavy atoms of each IDT, and BT fragments (through a least squares fitting) and recognized π-π interaction between them that satisfy the following criteria: 1) the angle between two planes' normal vectors is smaller than 10° (parallel fragments), 2) π-π distance (marked as Dπ-π in Figure S11a) is smaller than 0.5 Å, and 3) the horizontal distance between the centre of geometry of each parallel pair (marked as HCOG in Figure S11a) is smaller than 5.0 Å. Figure S11b shows the distribution of Dπ-π as calculated for 150 π-π coupled BT fragments obtained from 5 independent simulations snapshots (each extracted from a different annealing cycle).
It should be note that possible π-π interactions between the central six-carbon ring of IDT (not the whole fragment) was also explored to check possible crossings between IDT fragments and the number of π-π interactions between them were as rare as the IDT-IDT π-π interactions (<0.005 of total IDT fragments in the simulation box).The code which identifies and calculates π-π interactions and examples of required coordinate and index files for this calculation are provided in https://github.com/HMakkiMD/IDTBT.

Sidechain effect
To confirm the origin of BT-BT crossings, we generated IDTBT chain models with different sidechain lengths, i.e., IDTBT8-5, IDTBT4-5, and IDTBT1-5 (see monomer structures in  We calculated the torsional angle distribution and the chain orientation (against z axis) for equilibrated models of IDTBT16-5 and IDTBT1-5, see Figures S13 and S14, respectively.As shown, the change in the relative orientation of chains (as a result of shortening the sidechains) does not influence the distribution of torsional angles, however, the packing of chains lead to an obvious deviation from random orientation-a common observation for semicrystalline polymers.It is worth noting that according to the reports in the literature 4 , it is not possible to synthesize IDTBT polymer with sidechains shorter than 4 carbon length (IDTBT4) due to the solubility problem of the monomers during polymerisation, and our analysis on the hypothetical polymers (e.g., IDTBT1) has been only done to unravel the explicit role of sidechains on IDTBT morphology.

S3. X-ray Scattering Pattern Calculation
We use Debye scattering equation (SE5) 5 to calculate the X-ray scattered intensity of polymer models.

SE5
where Iq is the scattering intensity, f is atomic form factor, q is the magnitude of the scattering vector, and rij is the distance between atoms i and j.We employed a distance-histogram approximation [4] to avoid calculation of expensive sine function for each atomic pair.Therefore, we split computations into two steps: (i) a histogram of distances for each pair of atoms (excluding hydrogen atoms) is calculated and (ii) the Debye formula treats distances in each histogram bin at once so that the Debye formula (assuming a mono atomic system) can be written as SE6 (derivation can be found in https://debyer.readthedocs.io/en/latest/).
where N is the number of atoms (excluding hydrogen) in the simulation box, nk and rk are the are the number of pairs and the distance corresponding to the k-th bin.
The code by which the scattering patterns were calculated is provided here https://github.com/HMakkiMD/IDTBT.

Sampling Procedure for Electronic Structure Calculations
Chain conformations for electronic structure calculation were extracted periodically throughout molecular dynamics simulation as snapshots, which were deemed to be sufficiently uncorrelated from one another at a simulation time difference of one annealing cycle (25 ns total in the NPT ensemble with 10 ns at 900 K, 10 ns at 500 K, 2 ns at 300 K, and 3 ns for temperature ramping).5 snapshots for each system were considered, with equivalent sampling for chains of varying lengths determined according to the following equation, where TS is the total sampling, NC is the number of chains in the simulation box, L is the chain length (in number of monomers), and NS is the number of snapshots.If we set TS to 2500, and NS is 5, then the values of NC for IDTBT16-5, IDTBT16-10, and IDTBT16-20 are 100, 50, and 25, respectively.Some approximations were necessary for efficient computation of the electronic structure.Firstly, the -C16H33 side-chains modelled in MD simulation were truncated to -CH3 groups in order to significantly reduce computational expense with negligible effect on the electronic structure due to the dominant contribution from the backbone.A second approximation was the treatment of the electronic structure of individual chains separately, as the overall influence of coupling between chains is expected to be minimal.To consider the electrostatic environment of each chain, all monomers containing a charged atom within 20 Å of a charged atom on the central chain were included as point charges in the calculation, where the charges on each atom were the same as those used for the force field and therefore balanced.The B3LYP/3-21G level of theory was used as implemented in Gaussian16 software for these calculations.

Density of States
The density of states (DOS) was computed for all one-electron states according to the abovementioned sampling procedure.We define the DOS here similarly as the per chain per monomer DOS, where Ei (m) is the energy of molecular orbital m for chain i, M is the number of monomers in the chain, and g denotes a Gaussian approximation for the Dirac delta function.The bulk DOS, ρ b (E) (shown as DOS(E) in the manuscript) is then an average over the per chain per monomer DOS for all chains, where NC is the total number of chains considered.A Gaussian broadening parameter of 0.025 eV was used for computation of the DOS.

Transfer Integral Calculations
Transfer integrals were computed as a measure of the coupling between 150 BT crossing points, which we deem to be an interchain BT pair in a stacking configuration.We define a stacked pair as one where both the vertical and horizontal distance between the centres of geometry of each fragment is smaller than 5 Å and angle between normal vectors of planes fitted to the heavy atoms of each is smaller than 10°.BT fragments were isolated by cutting at the carbons connected at either side (inclusive), and made whole by replacing these atoms with hydrogen.
Transfer integral calculations were performed using in-house Python code according to the method presented in ref 6 , using the definition, where φ i and φ j are the unperturbed HOMO orbitals of the monomers and F ̂ is the Fock operator of the molecular dimer.Electronic structure calculations were carried out at the B3LYP/3-21G* level of theory as implemented in Gaussian16.

Figure
Figure S1.IDT-BT monomer structure.The blue (BT) and red (IDT) highlighted fragments were used to cap the repeat unit.

Figure S2 .
Figure S2.All possible Cis-Trans isomers of IDT-BT.Structure 7 shows the minimum energy value among all.

Figure 3 .
Figure 3. Atom labels as used for force field parameterisation of one IDT-BT repeat unit (RU).

Figure
Figure S4.(a) IDT-BT monomer capped with BT (left-end) and IDT (fight-end) molecules; the red circle shows the place of rotation for dihedral scan (ϕIDT-BT).(b) DFT calculated torsional potential by B3LYP/cc-pvtz.

Figure
Figure S5.(a) Total energy (excluding the restraint energy imposed to keep ϕ IDT-BT constant) based on the force field (excluding the ϕIDT-BT torsional potential) for energy minimised structures.(b) Force field torsional potential (V CORR = V DFT -V FF ).

Figure S7 .
Figure S7.One hydrogen of each methyl group attached to IDT is removed (red circle) and its charge is redistributed on the two remaining hydrogens (blue circles).

Figure S8 .
Figure S8.Atom labeling of the RU-SC (repeat unit with side chain).Green, blue, and red colored numbers represent sp2 carbons, non-sp2 hybridised atoms, and sidechain united atoms, respectively.In this way, we constructed RU-SC structures with different side chain lengths.RU-SC16, RU-SC8, RU-SC4, and RU-SC1 are shown in figureS9.For instance, RU-SC8 is the repeat unit of IDT-BT polymer with four side chains with a length of 8 carbon atoms.

Figure S10 .
Figure S10.Volume-temperature curves for the three polymer models.The intersection of fitted lines on the data well-above (700-800 K) and well-below (200-300 K) of the transition is used to estimate Tg.The cooling rate for all simulations is 5 K/ns.The standard deviation of volume at each temperature step is shown with vertical error bars.

Figure
Figure S11.(a) Illustration of Dπ-π and HCOG for one example of BT-BT in π-π interactions.(b) distribution of Dπ-π as calculated for 150 π-π coupled BT fragments in 5 independent simulations snapshots (each extracted from a different annealing cycle).

Figure
Figure S12 a), and equilibrated these polymer models in a similar way as IDTBT16-5.Then, the number of BT-BT π-π interactions and θij for all models were calculated.As shown in S12 b and c, the number of π-π interactions increases and the distribution of θij tends towards 0 ° as the sidechain length is shortened.This reconfirms that the reason behind relative perpendicular orientation of BT-BT crossings is the steric hinderance of the bulky sidechains.

Figure
Figure S12.(a) IDTBT monomer models with varying sidechain lengths.(b) Number of BT-BT π-π interactions for different polymer models with varying sidechain lengths.(c) Distribution of θij for different polymer models with varying sidechain lengths.

Figure S14 .
Figure S14.The distribution of angles between chain directors (a vector connecting the two end of chains) with the z axis for IDTBT16-5 (left) and IDTBT1-5 (right) models.The normalised sin θ function is shown evaluate the randomness of distributions.
Figure5cof the main paper shows a bulk partial density of states, PDOS(E), for IDTBT, which is a projection of the contributions of IDT and BT fragments to the total DOS.The per chain per monomer PDOS for a particular fragment is calculated as follows,

Figure S15 .
Figure S15.Electronic bandgap of the three polymer models calculated by our QC/MD method.

Table 2 .
The OPLS force field constants of Ryckaert-Bellemans function for C-C-C-C dihedral angle of benzene.