Fully Atomistic Molecular Dynamics Simulation of a TIPS-Pentacene:Polystyrene Mixed Film Obtained via the Solution Process

Organic thin-film transistors using small-molecule semiconductor materials such as 6,13-bis(triisopropylsilylethynyl)pentacene (TIPS-P) have been recently studied for the production of flexible and printed electronic devices. Blending a semiconductor with an insulating polymer, such as polystyrene, is known to improve the device performance; however, its molecular-level structure remains unknown. In this study, we performed molecular dynamics (MD) simulations on a mixed system of TIPS-P and atactic polystyrene (aPS) with fully atomistic models to understand the structure of the mixed thin film at the molecular level and the influence on the device properties. To reproduce the deposition from the solution, we gradually reduced the number of toluene molecules in the simulation. The dynamic characteristics of the system, mean squared displacement, diffusion coefficient, density profile, and P2 order parameter were analyzed. Some of the simulated systems reached the equilibrium state. In these systems, the simulated structures suggested the presence of more TIPS-P molecules on the surface than inside the bulk, even at the low molecular weight of aPS, where phase separation was not observed experimentally. The results of the fully atomistic MD simulations are also a basis for the coarse-grained model to increase the speed of the MD simulation.


Introduction
Flexible and printed organic electronic devices have been extensively studied over the last three decades [1] because of their unique mechanical features, such as their ability to bend, stick, and be disposed of [2][3][4], which are different from those of conventional electronic devices based on rigid inorganic semiconductors. These unique features make flexible electronic devices suitable for healthcare applications such as biosensing [5,6] and electronic organs [7][8][9] (bioelectronics). The two main trends in the research field of flexible electronic devices are (i) reducing the thickness of films integrating or constituting the electronic components [10] and (ii) developing new materials suitable for producing flexible electronic devices [11][12][13]. One of the molecules used as organic semiconductors in flexible electronic devices is pentacene, which exhibits a high carrier mobility [14,15] (~1.2 cm 2 /Vs) required for applications in electronic devices. However, the material requires a vacuum deposition process, which notably increases the cost [16,17]. An alternative to pentacene is 6,13-bis(triisopropylsilylethynyl)pentacene (TIPS-P), which has been developed as a highly soluble semiconductor for flexible electronic devices [18,19]. In particular, TIPS-P is used as the active layer of organic thin-film transistors as a small-molecule semiconductor owing to its solubility and conductivity [17,20]. Owing to its high solubility, TIPS-P can be deposited via solution processes at room temperature at a lower cost than pentacene and is suitable for large-area electronic products [16,18,21]. TIPS-P also has a high carrier mobility, exceeding 1 cm 2 /Vs [17]. Ohe et al. experimentally reported that blending TIPS-P with an insulating polymer, poly-α-methylstyrene (PαMS), can improve the performance and homogeneity of organic thin-film transistors [17]. In their experiments, spin-coating a solution of TIPS-P and PαMS in toluene formed a tri-layer structure of a TIPS-P layer/mixed layer of TIPS-P and PαMS/TIPS-P layer via spontaneous phase separation if the weight-average molecular weight (M w ) of PαMS was larger than 28,200 [17]. The blend of TIPS-P and PαMS with a low M w of 2200 resulted in lower device performance without phase separation. However, the reason why the insulating polymer affects the device performance, as well as its interfacial structures, has not been revealed to date. The technique of blending small-molecule semiconductors with insulating polymers has been extended to different combinations of materials [22].
In this study, we performed fully atomistic molecular dynamics (MD) simulations of a TIPS-P:atactic polystyrene mixed film to understand its structure and interactions at the molecular level. The thin-film formation from the solution was mimicked by gradually reducing the toluene content in the system. Using MD simulations with atomistic models makes it possible to reproduce the molecular state at the surface in more detail, which is otherwise difficult to accomplish by using experimental approaches such as atomic force microscopy and X-ray diffraction. To date, trial-and-error experiments based on device features have been conducted for the selection of polymer materials. Elucidating the mechanism from the results of MD simulations will guide material selection. The results of the fully atomistic MD simulation are also a basis for the coarse-grained (CG) model to increase the speed of the MD simulation.

Atomistic Models
In this study, atomistic models of atactic polystyrene (aPS, C 8n H 8n+2 ) [23], TIPS-P (C 44 H 54 Si 2 ) [24], and toluene (C 7 H 8 ) [25] were used. Figure 1 shows the molecular formulas of aPS, TIPS-P, and toluene. For all the considered models, OPLS-AA force field parameters were used. A short oligomer of aPS (M w = 1056, C 80 H 82 ), composed of ten repeating units, was selected because of the limitation of the calculation cost and to have access to the characteristic relaxation time of aPS occurring at different timescales, which can be achieved by full-atomistic simulations. The aPS model, which is able to correctly reproduce the structural and dynamic properties of the polymer bulk, was taken from the literature [23,26]. The TIPS-P model, proposed by Steiner et al., was used to study the crystalline structure of TIPS-P [24]. Nanomaterials 2023, 13, x 2 of 10 than pentacene and is suitable for large-area electronic products [16,18,21]. TIPS-P also has a high carrier mobility, exceeding 1 cm 2 /Vs [17]. Ohe et al. experimentally reported that blending TIPS-P with an insulating polymer, poly-α-methylstyrene (PαMS), can improve the performance and homogeneity of organic thin-film transistors [17]. In their experiments, spin-coating a solution of TIPS-P and PαMS in toluene formed a tri-layer structure of a TIPS-P layer/mixed layer of TIPS-P and PαMS/TIPS-P layer via spontaneous phase separation if the weight-average molecular weight (Mw) of PαMS was larger than 28,200 [17]. The blend of TIPS-P and PαMS with a low Mw of 2200 resulted in lower device performance without phase separation. However, the reason why the insulating polymer affects the device performance, as well as its interfacial structures, has not been revealed to date. The technique of blending small-molecule semiconductors with insulating polymers has been extended to different combinations of materials [22]. In this study, we performed fully atomistic molecular dynamics (MD) simulations of a TIPS-P:atactic polystyrene mixed film to understand its structure and interactions at the molecular level. The thin-film formation from the solution was mimicked by gradually reducing the toluene content in the system. Using MD simulations with atomistic models makes it possible to reproduce the molecular state at the surface in more detail, which is otherwise difficult to accomplish by using experimental approaches such as atomic force microscopy and X-ray diffraction. To date, trial-and-error experiments based on device features have been conducted for the selection of polymer materials. Elucidating the mechanism from the results of MD simulations will guide material selection. The results of the fully atomistic MD simulation are also a basis for the coarse-grained (CG) model to increase the speed of the MD simulation.

Atomistic Models
In this study, atomistic models of atactic polystyrene (aPS, C8nH8n+2) [23], TIPS-P (C44H54Si2) [24], and toluene (C7H8) [25] were used. Figure 1 shows the molecular formulas of aPS, TIPS-P, and toluene. For all the considered models, OPLS-AA force field parameters were used. A short oligomer of aPS (Mw = 1,056, C80H82), composed of ten repeating units, was selected because of the limitation of the calculation cost and to have access to the characteristic relaxation time of aPS occurring at different timescales, which can be achieved by full-atomistic simulations. The aPS model, which is able to correctly reproduce the structural and dynamic properties of the polymer bulk, was taken from the literature [23,26]. The TIPS-P model, proposed by Steiner et al., was used to study the crystalline structure of TIPS-P [24].

Calculation for the Pristine TIPS-P Crystal
To validate the force field for TIPS-P, an MD simulation of pristine TIPS-P crystals without PS and toluene was performed. The 10 × 8 × 4 supercell of the experimental structure [27] was used as the initial configuration. Energy minimization was performed using 5000 steps of the steepest descent algorithm. A time step of 2 fs for this system and a cut-off

Calculation for the Pristine TIPS-P Crystal
To validate the force field for TIPS-P, an MD simulation of pristine TIPS-P crystals without PS and toluene was performed. The 10 × 8 × 4 supercell of the experimental structure [27] was used as the initial configuration. Energy minimization was performed using 5000 steps of the steepest descent algorithm. A time step of 2 fs for this system and a cut-off distance for the Lennard-Jones non-bonded interaction of 1.2 nm were set. Periodic boundary conditions were considered in all directions (x, y, and z). For the production run, the temperature was kept constant using a velocity-rescaling algorithm [28] with a characteristic relaxation time τ T = 0.02 ps and a fixed volume. Simulated annealing was performed from 0 to 200 K every 100 ps in 50 K increments, with a simulation time of 1 ns. Chemical structures of all molecules used in the MD simulations are reported in Figure 1.

Calculation for the Mixed Solution System
The free software GROMACS [29] (Ver. 2016.4) was used for all simulations, and the temperature was kept constant using a velocity-rescaling algorithm [28] with a characteristic relaxation time τ T = 0.02 ps. The pressure was kept constant using the Berendsen weakcoupling scheme, and the box was scaled isotropically with a characteristic relaxation time τ P = 0.3 ps. Periodic boundary conditions were used in all directions (x, y, and z). A time step of 2 fs was set for all systems. The cut-off distance for the Lennard-Jones non-bonded interactions was set to 1.2 nm. The electrostatic interactions were computed using the Ewald particle mesh (mesh spacing in Fourier space: 0.12 nm) [30,31]. All the bonds involving hydrogen atoms were constrained using the LINCS algorithm [32]. For all the initial coordinate sets, the energy of the systems was minimized by performing 5000 steps of the steepest descent algorithm. On each initial set of coordinates, an NPT simulation of the bulk system was performed until the total mass density reached equilibrium. As a further validation model test, we performed an additional run in the NPT ensemble by using a different barostat algorithm, the Parrinello-Rahman algorithm (with same target temperature and pressure) [33]. For the additional run, we used the last wellequilibrated configuration obtained from the set-up with a Berendsen barostat. In Figure 2, the comparison of the time evolution of the mass density of the system calculated for the two barostat algorithms is reported. As can be seen from the trends in Figure 2, we can assume that the key interactions reproduced by the Berendsen set-up are in reasonable agreement with the Parrinello-Rahman set-up (which is known to be accurate in the reproduction of the canonical ensemble) [34,35]. For both algorithms, the mass density converges to the same equilibrium value. distance for the Lennard-Jones non-bonded interaction of 1.2 nm were set. Periodic boundary conditions were considered in all directions (x, y, and z). For the production run, the temperature was kept constant using a velocity-rescaling algorithm [28] with a characteristic relaxation time τT = 0.02 ps and a fixed volume. Simulated annealing was performed from 0 to 200 K every 100 ps in 50 K increments, with a simulation time of 1 ns. Chemical structures of all molecules used in the MD simulations are reported in Figure 1.

Calculation for the Mixed Solution System
The free software GROMACS [29] (Ver. 2016.4) was used for all simulations, and the temperature was kept constant using a velocity-rescaling algorithm [28] with a characteristic relaxation time τT = 0.02 ps. The pressure was kept constant using the Berendsen weak-coupling scheme, and the box was scaled isotropically with a characteristic relaxation time τP = 0.3 ps. Periodic boundary conditions were used in all directions (x, y, and z). A time step of 2 fs was set for all systems. The cut-off distance for the Lennard-Jones non-bonded interactions was set to 1.2 nm. The electrostatic interactions were computed using the Ewald particle mesh (mesh spacing in Fourier space: 0.12 nm) [30,31]. All the bonds involving hydrogen atoms were constrained using the LINCS algorithm [32]. For all the initial coordinate sets, the energy of the systems was minimized by performing 5000 steps of the steepest descent algorithm. On each initial set of coordinates, an NPT simulation of the bulk system was performed until the total mass density reached equilibrium. As a further validation model test, we performed an additional run in the NPT ensemble by using a different barostat algorithm, the Parrinello-Rahman algorithm (with same target temperature and pressure) [33]. For the additional run, we used the last wellequilibrated configuration obtained from the set-up with a Berendsen barostat. In Figure  2, the comparison of the time evolution of the mass density of the system calculated for the two barostat algorithms is reported. As can be seen from the trends in Figure 2, we can assume that the key interactions reproduced by the Berendsen set-up are in reasonable agreement with the Parrinello-Rahman set-up (which is known to be accurate in the reproduction of the canonical ensemble) [34,35]. For both algorithms, the mass density converges to the same equilibrium value. After the NPT run, we performed an NVT simulation by taking the last configuration of the NPT simulation and extending the box side length up to 20 nm to create the interface under vacuum. The toluene content was gradually reduced from 50 to 30, 20, 10, and 5 w/w%. A constant number of TIPS-P molecules and aPS chains was used (Table 1). First, we prepared a 50 w/w% toluene system and randomly set all molecules in the simulation After the NPT run, we performed an NVT simulation by taking the last configuration of the NPT simulation and extending the box side length up to 20 nm to create the interface under vacuum. The toluene content was gradually reduced from 50 to 30, 20, 10, and 5 w/w%. A constant number of TIPS-P molecules and aPS chains was used (Table 1). First, we prepared a 50 w/w% toluene system and randomly set all molecules in the simulation box without a vacuum phase. Minimization and NPT simulations were performed for the bulk, followed by an NVT simulation, in which the interface of the bulk with vacuum was included. Considering the experimental conditions [17], the simulations were performed at room temperature (300 K) and the annealing temperature (333 K). A trajectory of 1 µs was accumulated for all systems. Starting from a set of coordinates (taken at equilibrium) of the system at a higher toluene content, some toluene molecules were removed from the system according to the new (lower) concentration. The same protocol, minimization, and NPT run were then applied, and a long NVT production run (1 µs with the bulk/vacuum interface) was performed. This procedure was performed until a system with a toluene content of 5 w/w% was simulated. The compositions and total times of the simulated systems are listed in Table 1.  box without a vacuum phase. Minimization and NPT simulations were performed for the bulk, followed by an NVT simulation, in which the interface of the bulk with vacuum was included. Considering the experimental conditions [17], the simulations were performed at room temperature (300 K) and the annealing temperature (333 K). A trajectory of 1 µs was accumulated for all systems. Starting from a set of coordinates (taken at equilibrium) of the system at a higher toluene content, some toluene molecules were removed from the system according to the new (lower) concentration. The same protocol, minimization, and NPT run were then applied, and a long NVT production run (1 µs with the bulk/vacuum interface) was performed. This procedure was performed until a system with a toluene content of 5 w/w% was simulated. The compositions and total times of the simulated systems are listed in Table 1.

Pristine TIPS-P Crystal for Validation of the Force Field
Different planar views of a representative equilibrium configuration of the simulated TIPS-P crystal (Figure 3a [27]. These results indicate the validity of the adopted force field to simulate the crystalline structure of TIPS-P.  The mean squared displacement (MSD) [37][38][39] can provide an approximate estimation of the capability of sampling the entire phase space of the system for each component.  Figure 4 shows the time-dependence of the MSDs in the direction normal to the bulk/vacuum interface. At toluene content of 50 w/w%, the MSDs of TIPS-P and aPS saturated at 15 and 10 nm 2 , respectively. The thickness T of the ternary system was approximately 10 nm at all times. If each molecule has enough time to diffuse within the thickness, the MSD is expected to be where s and t are the initial and final positions of a molecule, and 1/T gives the uniform probability density of the position between 0 and T. This formula and T = 10 nm give MSD = 16.7 nm 2 . Therefore, the saturation of MSDs at 15 and 10 nm 2 indicates that the TIPS-P and aPS had enough time to diffuse within the thickness in the systems with toluene content of 50 w/w%. However, the MSDs in the systems with toluene content of 30 w/w% or less did not saturate within the simulated time, which indicates that systems with lower toluene content require simulation times longer than 1000 ns. The diffusion coefficient of each component was also calculated using the least-squares method of linear fitting to the time-MSD curve. The diffusion constants for all toluene contents are listed in Table 2.
ray diffraction patterns of the TIPS-P crystal.

Mean Squared Displacement
The mean squared displacement (MSD) [37][38][39] can provide an approximate estimation of the capability of sampling the entire phase space of the system for each component. Figure 4 shows the time-dependence of the MSDs in the direction normal to the bulk/vacuum interface. At toluene content of 50 w/w%, the MSDs of TIPS-P and aPS saturated at 15 and 10 nm 2 , respectively. The thickness T of the ternary system was approximately 10 nm at all times. If each molecule has enough time to diffuse within the thickness, the MSD is expected to be where s and t are the initial and final positions of a molecule, and 1/T gives the uniform probability density of the position between 0 and T. This formula and T = 10 nm give MSD = 16.7 nm 2 . Therefore, the saturation of MSDs at 15 and 10 nm 2 indicates that the TIPS-P and aPS had enough time to diffuse within the thickness in the systems with toluene content of 50 w/w%. However, the MSDs in the systems with toluene content of 30 w/w% or less did not saturate within the simulated time, which indicates that systems with lower toluene content require simulation times longer than 1000 ns. The diffusion coefficient of each component was also calculated using the least-squares method of linear fitting to the time-MSD curve. The diffusion constants for all toluene contents are listed in Table 2.  (2.5 ± 0.5)×10 −6 (3 ± 2)×10 −6 (7.2 ± 0.5)×10 −8

Polymer Configuration
To evaluate the effectiveness of sampling the polymer configurations, we computed the end-to-end relaxation time for the aPS chains. The autocorrelation function ( ) =

Polymer Configuration
To evaluate the effectiveness of sampling the polymer configurations, we computed the end-to-end relaxation time for the aPS chains.
The autocorrelation function C(t) = (R(τ)R(τ+t)) (R 2 ) was used to compute the relation time ( Figure 5), where R is the length of the end-to-end distance vector, and the brackets indicate the time average. Then, the relaxation time (τ PS ) was computed by integrating the stretched exponential function fitted to the autocorrelation functions, as described in Equation (1).
The τ PS values are listed in Table 3. As the toluene content decreased from 50 to 5 w/w%, the relaxation time of aPS increased from~0.6 to 200 ns. In addition, shorter relax-Nanomaterials 2023, 13, 312 6 of 10 ation times were observed at higher temperatures. The total simulation time was approximately 4-5 times the relaxation time τ PS at 5 w/w% toluene. From these data, we can conclude that the sampling of chain conformations is sufficient for all the toluene contents.  Figure 5. Autocorrelation functions of the end-to-end distance for the systems of 50 w/w% toluene at 300 K and 333 K (a,b) and 5 w/w% toluene at 300 K and 333 K (c,d).

Figure 5.
Autocorrelation functions of the end-to-end distance for the systems of 50 w/w% toluene at 300 K and 333 K (a,b) and 5 w/w% toluene at 300 K and 333 K (c,d).

Density Profile
From the results of MSD in Figure 4, the system with 50 w/w% toluene was considered under the equilibrium after 500 ns. The density profile of each component was calculated in 50 ns increments from 0 to 1000 ns. Figure 6 shows the density profiles, computed along the normal direction of the bulk/vacuum interface, of the system with 50 w/w% toluene at 50, 500, and 1000 ns. Figure 7 shows the configurations of the same system. The density profile indicates that the concentration of TIPS-P is higher within 1 nm of the surface. However, experimental studies have reported no phase separation when the M w of PαMS is 2000 or lower [17]. This is because the experimental method does not have a high spatial resolution of 1 nm. Our simulation reveals, for the first time, that the surface contains more TIPS-P molecules than the film, even at a low M w of 1056. Such phenomena can be considered a precursor of phase separation into a tri-layer structure at the high molecular weight of aPS.
puted along the normal direction of the bulk/vacuum interface, of the system with 50 w/w% toluene at 50, 500, and 1000 ns. Figure 7 shows the configurations of the same system. The density profile indicates that the concentration of TIPS-P is higher within 1 nm of the surface. However, experimental studies have reported no phase separation when the Mw of PαMS is 2000 or lower [17]. This is because the experimental method does not have a high spatial resolution of 1 nm. Our simulation reveals, for the first time, that the surface contains more TIPS-P molecules than the film, even at a low Mw of 1056. Such phenomena can be considered a precursor of phase separation into a tri-layer structure at the high molecular weight of aPS. Figure 6. Density profiles of aPS (black curve) and TIPS-P (green curve) together with the P2 order parameters of TIPS-P for the systems with 50 w/w% toluene at 300 K.

P2 Order Parameter
To investigate the orientation in different regions in the bulk and at the bulk/vacuum interface, we calculated the order parameter P2 [40]. The order parameter P2 is defined as , where is the angle between the molecular vector of TIPS-P and the direction normal to the bulk/vacuum interface. P 2 = −0.5 if the molecular vector is perpendicular to the interface, and P 2 = 1 if the molecular vector is parallel to the interface. If the molecules are randomly oriented, P 2 = 0. The TIPS-P molecular vector was defined using the positions of the two silicon atoms.
The order parameter P2 was calculated at each slice (thickness: 2 nm) parallel to the surface for a system with 50 w/w% toluene at 300 K at 50, 500, and 1000 ns (see Figure 6). As can be seen from the profiles, the P2 values were almost equal to 0 under all these Figure 6. Density profiles of aPS (black curve) and TIPS-P (green curve) together with the P 2 order parameters of TIPS-P for the systems with 50 w/w% toluene at 300 K.
puted along the normal direction of the bulk/vacuum interface, of the system with 50 w/w% toluene at 50, 500, and 1000 ns. Figure 7 shows the configurations of the same system. The density profile indicates that the concentration of TIPS-P is higher within 1 nm of the surface. However, experimental studies have reported no phase separation when the Mw of PαMS is 2000 or lower [17]. This is because the experimental method does not have a high spatial resolution of 1 nm. Our simulation reveals, for the first time, that the surface contains more TIPS-P molecules than the film, even at a low Mw of 1056. Such phenomena can be considered a precursor of phase separation into a tri-layer structure at the high molecular weight of aPS.

P2 Order Parameter
To investigate the orientation in different regions in the bulk and at the bulk/vacuum interface, we calculated the order parameter P2 [40]. The order parameter P2 is defined as , where is the angle between the molecular vector of TIPS-P and the direction normal to the bulk/vacuum interface. P 2 = −0.5 if the molecular vector is perpendicular to the interface, and P 2 = 1 if the molecular vector is parallel to the interface. If the molecules are randomly oriented, P 2 = 0. The TIPS-P molecular vector was defined using the positions of the two silicon atoms.
The order parameter P2 was calculated at each slice (thickness: 2 nm) parallel to the surface for a system with 50 w/w% toluene at 300 K at 50, 500, and 1000 ns (see Figure 6). As can be seen from the profiles, the P2 values were almost equal to 0 under all these

P 2 Order Parameter
To investigate the orientation in different regions in the bulk and at the bulk/vacuum interface, we calculated the order parameter P 2 [40]. The order parameter P 2 is defined as P 2 = 3 cos 2 θ −1 2 , where θ is the angle between the molecular vector of TIPS-P and the direction normal to the bulk/vacuum interface. P 2 = −0.5 if the molecular vector is perpendicular to the interface, and P 2 = 1 if the molecular vector is parallel to the interface. If the molecules are randomly oriented, P 2 = 0. The TIPS-P molecular vector was defined using the positions of the two silicon atoms.
The order parameter P 2 was calculated at each slice (thickness: 2 nm) parallel to the surface for a system with 50 w/w% toluene at 300 K at 50, 500, and 1000 ns (see Figure 6). As can be seen from the profiles, the P 2 values were almost equal to 0 under all these conditions, which indicates the random orientation of the investigated molecules. This is consistent with the experimental result for PαMS with M w ≤ 2000 [17].

Surface Profile
Some indication of the roughness of the surface can be obtained from the profile decay. The roughness of the surface can be measured using the root mean squared roughness R q [41], where R q is defined as To this end, the system was divided into 10 × 10 bins in the direction parallel to the bulk/vacuum interface (x, y). For each bin, the largest height h(x, y) was calculated by considering the atom coordinates of both TIPS-P and aPS. The term h in Equation (2) is the average height of the surface, computed by considering all 10 × 10 bins. The R q values computed for all systems are compared in Figure 8.
Some indication of the roughness of the surface can be obtained from the profile decay. The roughness of the surface can be measured using the root mean squared roughness (R q ) [41], where R q is defined as To this end, the system was divided into 10 × 10 bins in the direction parallel to the bulk/vacuum interface (x, y). For each bin, the largest height h(x,y) was calculated by considering the atom coordinates of both TIPS-P and aPS. The term ℎ ̅ in Equation (2) is the average height of the surface, computed by considering all 10 × 10 bins. The R q values computed for all systems are compared in Figure 8. The main result emerging from the roughness calculation confirms that the TIPS-P and aPS molecules were arranged in a manner that minimizes R q . R q was the highest with 50 w/w% toluene (0.31 nm). R q decreased as the toluene content decreased, and finally, R q was 0.17 nm with 5 w/w% toluene. The range of agreed with the experimental range (0.2-0.5 nm) [42]. The increase in the amount of toluene perturbed the roughness of both molecules. The temperature did not seem to affect the roughness significantly.

Conclusions
We performed a fully atomistic MD simulation for a blend system of TIPS-P and aPS in toluene. The evaporation of toluene was mimicked by a gradual reduction in the number of toluene molecules. The MSDs and diffusion coefficients for all systems indicated that equilibrium states were reached under two conditions: 50 w/w% at 300 K and 50 w/w% at 333 K. The density profile revealed that the concentration of TIPS-P was higher at the interface with vacuum, whereas no clear phase separation was observed because of the low molecular weight of aPS. The high concentration of TIPS-P at the surface is advantageous to forming conductive paths in organic thin-film transistors. The random orientation of the TIPS-P molecules was observed, which is consistent with the experimental results. The analysis of the surface roughness indicated that the TIPS-P molecules were arranged in a manner that produced a smooth surface, which is compatible with the high ordering of the molecules. Further studies with a higher Mw would clarify the film structure with improved mobility in the experiments. However, this requires a CG model The main result emerging from the roughness calculation confirms that the TIPS-P and aPS molecules were arranged in a manner that minimizes R q . R q was the highest with 50 w/w% toluene (0.31 nm). R q decreased as the toluene content decreased, and finally, R q was 0.17 nm with 5 w/w% toluene. The range of R q agreed with the experimental range (0.2-0.5 nm) [42]. The increase in the amount of toluene perturbed the roughness of both molecules. The temperature did not seem to affect the roughness significantly.

Conclusions
We performed a fully atomistic MD simulation for a blend system of TIPS-P and aPS in toluene. The evaporation of toluene was mimicked by a gradual reduction in the number of toluene molecules. The MSDs and diffusion coefficients for all systems indicated that equilibrium states were reached under two conditions: 50 w/w% at 300 K and 50 w/w% at 333 K. The density profile revealed that the concentration of TIPS-P was higher at the interface with vacuum, whereas no clear phase separation was observed because of the low molecular weight of aPS. The high concentration of TIPS-P at the surface is advantageous to forming conductive paths in organic thin-film transistors. The random orientation of the TIPS-P molecules was observed, which is consistent with the experimental results. The analysis of the surface roughness indicated that the TIPS-P molecules were arranged in a manner that produced a smooth surface, which is compatible with the high ordering of the molecules. Further studies with a higher M w would clarify the film structure with improved mobility in the experiments. However, this requires a CG model instead of a fully atomistic model to reduce the calculation cost and to gain access to the intrinsically faster dynamics of the CG models. The present results can be used to develop CG models and to approach the problem with a multiscale strategy, reintroducing atomic details via backmapping methodology where needed.