Theoretical Evaluation of the Influence of Molecular Packing Mode on the Intramolecular Reorganization Energy of Oligothiophene Molecules

Accurate determination of the relationships among packing mode, molecular structure and charge transfer mobility for oligothiophene analogues has been significantly impeded, due to the lack of crystal structure information. In the current study, molecular dynamics (MD) were used to investigate the packing mode of non-, methyl- and ethyl-substituted poly(3-alkylthiophenes) (P3ATs). Obvious conformational changes were observed when comparing the packed and isolated oligothiophene molecules, indicating the important influence of packing mode on the geometric structures of these materials. Considering the crucial role played by reorganization energy (RE) in the charge transfer process, both quantum mechanics (QM) and quantum mechanics/molecular mechanics (QM/MM) were performed to examine the impact of different conformations on energy. Our simulations revealed that the geometric structures have distinct effects on the RE. Our data suggest that MD could give a reliable packing mode of oligothiophene analogues, and that QM/MM is indispensable for precisely estimating RE.


Introduction
Thiophene-based materials have attracted much attention in organic photovoltaic cells, due to the large π-conjugation and excellent chemical and physical properties [1]; regioregular poly(3-hexylthiophene) (P3HT) is even used as a model polymer for research in organic solar cells. [2,3]. Owing to the excellent charge transfer properties and structural simplicity, poly(3-alkylthiophene) (P3AT) based oligothiophene analogues have been considered ideal targets to investigate mechanisms of charge transfer processes [4][5][6][7][8]. Although these mechanisms are still not fully understood, Marcus theory has been used extensively to study the factors that influence the charge motilities [9], and two key parameters, i.e., the transfer integral (intermolecular electronic coupling), and the reorganization energy λ (RE, the local electron-phonon coupling) have been shown to determine self-exchange charge transfer rates. The reorganization energy of self-exchange charge transfer reaction in a hole-hopping material is defined as the sum of the geometrical relaxation energies of a molecule upon going from the neutral state to the charge state geometry, and the neighboring molecules upon going from the inverse Software Inc., San Diego, CA, USA) [25]. Optimization of the structures was achieved with the Gaussian09 package (Gaussian, Inc., Wallingford, CT, USA) [26] at the HF/6-31G* theoretical level. Electrostatic potentials (ESP) were then generated with Merz-Singh-Kollman van der Waals parameters [27]. Fitting of the charges to the ESP was conducted with the restrained ESP (RESP) program (University of California, San Francisco, CA, USA) [28] of the AMBER package [29]. GAFF [29] force field parameters and RESP partial charges were assigned by using the ANTECHAMBER module (University of California, San Francisco, CA, USA). The starting aggregation models were established with a systematic density of 2 g/cm 3 on the amorphous cell module of the Material Studio suite (Accelrys Software Inc., San Diego, CA, USA), based on the force field parameters in a consistent valence force field (CVFF) [30].
To make the simulation more consistent with real cases, each system consisted of amorphous structures containing 20 oligothiophene molecules ( Figure 1). The aggregation of 20 methyl-substituted 5T P3ATs was taken as an example. The MD simulations were carried out with the sander module provided in the dynamic software package Amber14 (University of California, San Francisco, CA, USA) [29]. With periodic boundary conditions, the 18 systems were range-limited in periodic water boxes of 36 [31] was adopted to calculate long-range electrostatic interactions. Energy optimization of each system was made within 1500 steps in total. In the former 500 steps, the steepest descent method was used and after the 500 steps, a conjugate gradient method was utilized. The systems were then heated to 298 K within 100 ps, so as to carry out the MD simulation under constant temperature and pressure in 1 fs for 100 ns. In this process, temperature and pressure were controlled by a weak-coupling algorithm and an isotropic position scaling method [32], respectively. The data processing and analysis was made using cpptraj [33] module in Amber14 (University of California, San Francisco, CA, USA).
Polymers 2018, 10, 30 3 of 11 the Gaussian09 package (Gaussian, Inc., Wallingford, CT, USA) [26] at the HF/6-31G* theoretical level. Electrostatic potentials (ESP) were then generated with Merz-Singh-Kollman van der Waals parameters [27]. Fitting of the charges to the ESP was conducted with the restrained ESP (RESP) program (University of California, San Francisco, CA, USA) [28] of the AMBER package [29]. GAFF [29] force field parameters and RESP partial charges were assigned by using the ANTECHAMBER module (University of California, San Francisco, CA, USA). The starting aggregation models were established with a systematic density of 2 g/cm 3 on the amorphous cell module of the Material Studio suite (Accelrys Software Inc., San Diego, CA, USA), based on the force field parameters in a consistent valence force field (CVFF) [30].
To make the simulation more consistent with real cases, each system consisted of amorphous structures containing 20 oligothiophene molecules ( Figure 1). The aggregation of 20 methyl-substituted 5T P3ATs was taken as an example. The MD simulations were carried out with the sander module provided in the dynamic software package Amber14 (University of California, San Francisco, CA, USA) [29]. With periodic boundary conditions, the 18 systems were range-limited in periodic water boxes of 36 [31] was adopted to calculate long-range electrostatic interactions. Energy optimization of each system was made within 1500 steps in total. In the former 500 steps, the steepest descent method was used and after the 500 steps, a conjugate gradient method was utilized. The systems were then heated to 298 K within 100 ps, so as to carry out the MD simulation under constant temperature and pressure in 1 fs for 100 ns. In this process, temperature and pressure were controlled by a weak-coupling algorithm and an isotropic position scaling method [32], respectively. The data processing and analysis was made using cpptraj [33] module in Amber14 (University of California, San Francisco, CA, USA).

QM and QM/MM Calculations of Reorganization Energy
RE λ was computed as half the energy difference between the vertical excitation and florescence energies as described before [16,34]. For the pure quantum mechanics calculations of λ, two separate optimizations of the structure of the isolated molecular cation in the ground state (S0) and in the first singlet excited state (S1) were firstly performed. Single-point calculations of the S0 → S1 excitation energy of the cation at the ground state geometry ΔE(S0), and at the S1 excited state geometry ΔE(S1) were then performed, i.e., λ = [ΔE(S0) − ΔE(S1)]/2. (1) All the geometry optimization computations were carried out using the Gaussion09 program at the B3lyp/6-31G and B3lyp/6-311G** levels, respectively.
The packing effect on the RE was estimated using the QM/MM method by the ONIOM module [35] in the Gaussian 09 program (Gaussian, Inc., Wallingford, CT, USA). In the QM/MM approach,

QM and QM/MM Calculations of Reorganization Energy
RE λ was computed as half the energy difference between the vertical excitation and florescence energies as described before [16,34]. For the pure quantum mechanics calculations of λ, two separate optimizations of the structure of the isolated molecular cation in the ground state (S 0 ) and in the first singlet excited state (S 1 ) were firstly performed. Single-point calculations of the S 0 → S 1 excitation energy of the cation at the ground state geometry ∆E(S 0 ), and at the S 1 excited state geometry ∆E(S 1 ) were then performed, i.e., All the geometry optimization computations were carried out using the Gaussion09 program at the B3lyp/6-31G and B3lyp/6-311G** levels, respectively.
The packing effect on the RE was estimated using the QM/MM method by the ONIOM module [35] in the Gaussian 09 program (Gaussian, Inc., Wallingford, CT, USA). In the QM/MM approach, the QM was only applied to the center molecule, while the MM computations of the surrounding molecules were conducted with the general Amber force field (GAFF) [29], as depicted in Figure 2. During the calculations, all the ambient oligothiophene molecules were fixed, and the center molecule was free to move. The same methodology and theory levels for the QM section were used as described in the pure QM calculations for precise comparisons. the QM was only applied to the center molecule, while the MM computations of the surrounding molecules were conducted with the general Amber force field (GAFF) [29], as depicted in Figure 2.
During the calculations, all the ambient oligothiophene molecules were fixed, and the center molecule was free to move. The same methodology and theory levels for the QM section were used as described in the pure QM calculations for precise comparisons.

Packing Mode of Different Length and Different Substituted P3AT
The packing mode of non-substituted P3AT by our MD simulations has been reported elsewhere [22], only the results of methyl-and ethyl-substituted P3ATs are discussed here. Figure 3 shows the variations of the root mean square deviations (RMSDs) of all the systems during the dynamic processes. As seen, the curves for each of the 12 systems stabilized after a certain time-span of fluctuations, indicating the molecular packing processes from the starting amorphous mode to a certain arrangement pattern.

Packing Mode of Different Length and Different Substituted P3AT
The packing mode of non-substituted P3AT by our MD simulations has been reported elsewhere [22], only the results of methyl-and ethyl-substituted P3ATs are discussed here. Figure 3 shows the variations of the root mean square deviations (RMSDs) of all the systems during the dynamic processes. As seen, the curves for each of the 12 systems stabilized after a certain time-span of fluctuations, indicating the molecular packing processes from the starting amorphous mode to a certain arrangement pattern. the QM was only applied to the center molecule, while the MM computations of the surrounding molecules were conducted with the general Amber force field (GAFF) [29], as depicted in Figure 2.
During the calculations, all the ambient oligothiophene molecules were fixed, and the center molecule was free to move. The same methodology and theory levels for the QM section were used as described in the pure QM calculations for precise comparisons.

Packing Mode of Different Length and Different Substituted P3AT
The packing mode of non-substituted P3AT by our MD simulations has been reported elsewhere [22], only the results of methyl-and ethyl-substituted P3ATs are discussed here. Figure 3 shows the variations of the root mean square deviations (RMSDs) of all the systems during the dynamic processes. As seen, the curves for each of the 12 systems stabilized after a certain time-span of fluctuations, indicating the molecular packing processes from the starting amorphous mode to a certain arrangement pattern. In contrast to data generated from our previous theoretical models [24] and the crystal structures [17][18][19]36] of the herringbone packing pattern of the non-substituted P3AT chains, the methyland ethyl-substituted oligothiophene molecules were arranged in parallel during the dynamic equilibrium, which ultimately led to a stacking plate packing mode, as depicted in Figure 4. Compared with the non-substituted structure [19], both the methyl-and ethyl-substituted 3T P3AT were less closely-packed, which was caused by side-chain replacement. However, with an increase in the main poly-thiophene chain length, all of the oligothiophene molecules appear to have a close-packed arrangement, although side-chain substitution is unfavorable for close-packing, the negative influence from which could be overcome by the increase in the main chain length. From the figure, it is also clear that the methyl-and ethyl-substituted poly-thiophene ring main chain mainly shows a flattened conformation, which is consistent with our simulation results of the non-substituted P3AT molecules [24].  In contrast to data generated from our previous theoretical models [24] and the crystal structures [17][18][19]36] of the herringbone packing pattern of the non-substituted P3AT chains, the methyl-and ethyl-substituted oligothiophene molecules were arranged in parallel during the dynamic equilibrium, which ultimately led to a stacking plate packing mode, as depicted in Figure 4. Compared with the non-substituted structure [19], both the methyl-and ethyl-substituted 3T P3AT were less closely-packed, which was caused by side-chain replacement. However, with an increase in the main poly-thiophene chain length, all of the oligothiophene molecules appear to have a close-packed arrangement, although side-chain substitution is unfavorable for close-packing, the negative influence from which could be overcome by the increase in the main chain length. From the figure, it is also clear that the methyl-and ethyl-substituted poly-thiophene ring main chain mainly shows a flattened conformation, which is consistent with our simulation results of the non-substituted P3AT molecules [24].  The radial distribution functions (RDFs) of selected carbon atoms in the center of the main-chain of the non-, methyl-, and ethyl-substituted 3T, 5T, and 8T P3ATs were computed to quantitatively describe the atomic packing and thiophene ring arrangement. In Figure 5, the non-substituted P3ATs interchain C-C RDFs exhibit strong peaks that indicate well-ordered nanostructure of these chains. It is clear that the peak intensity gets weaker and weaker in the methyl-and ethyl-substituted chains in all the 3T, 5T, and 8T cases, which means less compact packing modes, and is consistent with the qualitative results in Figure 4. The radial distribution functions (RDFs) of selected carbon atoms in the center of the main-chain of the non-, methyl-, and ethyl-substituted 3T, 5T, and 8T P3ATs were computed to quantitatively describe the atomic packing and thiophene ring arrangement. In Figure 5, the non-substituted P3ATs interchain C-C RDFs exhibit strong peaks that indicate well-ordered nanostructure of these chains. It is clear that the peak intensity gets weaker and weaker in the methyl-and ethyl-substituted chains in all the 3T, 5T, and 8T cases, which means less compact packing modes, and is consistent with the qualitative results in Figure 4.

Different Conformations from QM, MD, and QM/MM Simulations
To precisely estimate the RE, the structures from either X-ray experiments or MD simulations had to first be optimized. Single isolated molecules were optimized by QM without any restrictions, and well packed molecules were managed via QM/MM, as described in the Section 2.2. In the case of QM optimizations, all the poly-thiophene main chains appear to have various degrees of curvature. The conformation of ethyl-substituted P3ATs is shown as an example in Figure 6. The torsion angles between the neighboring thiophene rings and the degrees of curvature of all the systems are given in Table 1 to illustrate the conformational deviation from planarity. As discussed, in both the crystallography data and the MD simulation results, the conformation of the oligothiophene molecules in the well packed groups was found to be planar, as shown in Figure 4 and listed in Table  1. From both the figure and the table, it is obvious that the curvature of the isolated candidates from the QM method is much larger than that of the packed structure from the QM/MM calculations. Taking the non-, methyl-and ethyl-substituted 8T analogues as an example, the curvature from the QM method was 0.02687, 0.07546, and 0.10528, respectively, whilst that from the QM/MM method was 0.00008, 0.00854, and 0.019833, respectively. The smaller curvature relates to less curved conformation.

Different Conformations from QM, MD, and QM/MM Simulations
To precisely estimate the RE, the structures from either X-ray experiments or MD simulations had to first be optimized. Single isolated molecules were optimized by QM without any restrictions, and well packed molecules were managed via QM/MM, as described in the Section 2.2. In the case of QM optimizations, all the poly-thiophene main chains appear to have various degrees of curvature. The conformation of ethyl-substituted P3ATs is shown as an example in Figure 6. The torsion angles between the neighboring thiophene rings and the degrees of curvature of all the systems are given in Table 1 to illustrate the conformational deviation from planarity. As discussed, in both the crystallography data and the MD simulation results, the conformation of the oligothiophene molecules in the well packed groups was found to be planar, as shown in Figure 4 and listed in Table 1. From both the figure and the table, it is obvious that the curvature of the isolated candidates from the QM method is much larger than that of the packed structure from the QM/MM calculations. Taking the non-, methyl-and ethyl-substituted 8T analogues as an example, the curvature from the QM method was 0.02687, 0.07546, and 0.10528, respectively, whilst that from the QM/MM method was 0.00008, 0.00854, and 0.019833, respectively. The smaller curvature relates to less curved conformation. From the torsion angle values in Table 1, it can be seen that with increasing side chain length, the distortion between the neighboring rings also increases. For example, in the case of 8-thiophene-ring chains, QM and QM/MM torsion angles for the non-, methyl-, ethyl-substituted molecules were 149.82°, 122.32°, 114.23°, and 175.74°, 176.07°, 174.78°, respectively. The QM/MM results reveal much less effect of the side chain on the planarity of the poly-thiophene derivatives. The straight and planar conformations could be further confirmed by the consistency with the previous report about the non-substituted P3AT and P3HT main chain conformation [17][18][19]37,38]. The phenomenon also proves the influence of the packing on the conformation of the oligothiophene structures.   From the torsion angle values in Table 1, it can be seen that with increasing side chain length, the distortion between the neighboring rings also increases. For example, in the case of 8-thiophene-ring chains, QM and QM/MM torsion angles for the non-, methyl-, ethyl-substituted molecules were 149.82 • , 122.32 • , 114.23 • , and 175.74 • , 176.07 • , 174.78 • , respectively. The QM/MM results reveal much less effect of the side chain on the planarity of the poly-thiophene derivatives. The straight and planar conformations could be further confirmed by the consistency with the previous report about the non-substituted P3AT and P3HT main chain conformation [17][18][19]37,38]. The phenomenon also proves the influence of the packing on the conformation of the oligothiophene structures.

Molecular Packing Influences on the Reorganization Energy
The RE of the 6 non-substituted oligothiophene molecules was firstly calculated by both the QM and QM/MM methods to explore the effects originating from the curved and flat conformations.
The energy values are listed in Table 2. It can be seen from the first two columns of Table 2 that the molecular packing has a great influence on the calculated intramolecular RE. When the molecular packing was considered, the RE was about 15-20% lower than that of the isolated candidates. The RE of the methyl-and ethyl-substituted candidates was then calculated when the packing mode was included, to investigate the influence of conformational variations on the λ values. As can be observed in the table, with increasing side chain length, the λ values were also shown to increase. The increment between the methyl-and non-substitutions, and the ethyl-and non-substitutions is 1.0-4.5 kcal/mol and 1.4-5.3 kcal/mol, respectively. This is also consistent with the data in Figure 5, in which less close packing modes of the bigger side chain replacements always leads to lower RDFs. At the same time, it is evident that with the elongation of the poly-thiophene main chain, the RE decreases. However, it is observed that the presence of substitutions increases the RE by a greater amount for short systems (3T) compared to longer systems (8T). By taking a closer look at the data in Figure 5, the averaged RDF value of the non-substituted 3T, 5T, and 8T chains is 0.041, 0.030, and 0.020, respectively; that of the methyl-and ethyl-substituted 3T, 5T, and 8T molecules is 0.027, 0.024, 0.013 and 0.023, 0.014, 0.009, respectively. The reduction of the averaged RDFs from non-to ethyl-substituted 3T, 5T, and 8T molecules is 0.018, 0.014, and 0.011, respectively.

Discussion
Although pure poly-thiophene rings normally adopt a herringbone packing pattern, the simulation results clearly demonstrate that all the substituted P3ATs aggregate in a stacking plate form, indicating the important influence of the substitution on the packing mode. Furthermore, the conformation of the chains results in very different geometric structures when the packing mode is taken into consideration. Compared to the bowed conformations from the QM calculations, the QM/MM optimizations show a straight and planar conformation. From Table 1, it is obvious that, for the QM optimized isolated molecules, the bigger the substitution group, the more twisted the torsion angle between the neighboring thiophene rings, i.e., the bigger the substitution, the bigger influence on the planarity of the oligothiophene molecules. Considering the RE calculations as shown in the first two columns in Table 2, our data strongly suggests that when calculating systems with bigger substitutions, more emphasis should be given to choosing the most appropriate theoretical methodology, to obtain conformations that are as close as possible to the experimental packing mode. According to the Marcus theory, smaller RE leads to larger hopping rates, which means the charge mobility could be underestimated if the molecular packing is overlooked.
When the environment of the isolated molecules was considered, all the optimized structures are near planar, and the data in Table 1 clearly demonstrate that, irrespective of the planarity or the torsion angles, with the elongation of the poly-thiophene main chain, the RE decreases. These changes reflect the fact that a larger π-conjugation system is always more stable, and benefits the charge transport.

Conclusions
In summary, the packing modes of the non-, methyl-, and ethyl-substituted P3ATs were firstly simulated by MDs, and the RE of the P3ATs was calculated in the cases of single isolated and packed molecules. Our simulations show significant influence of the packing on the conformation of all the structures. When the packing pattern was considered, all the molecules show a flat conformation. Compared with the curved and twisted geometries, the flat conformation gives a smaller RE, and would eventually affect the charge mobility estimation. A longer main chain with more monomer units expands the π-electron delocalization, and decreases RE as well. Our simulations also reveal that the bigger substitution results in larger RE. Further investigations, therefore, need be done to find the optimal combination of a suitable size of substitution on a suitable length of main chain for an optimal RE value.
Our work proves that MD simulation is a reliable tool to locate proper packing modes of the oligothiophene molecules, and that QM/MM calculations are necessary and precise for the consideration of molecular packing, especially when big substitution groups are attached on the main chains.