Genetic algorithm-accelerated computational discovery of liquid crystal polymers with enhanced optical properties

Abstract

Liquid crystal polymers with exceptional optical properties are highly promising for next-generation virtual, augmented, and mixed reality (VR/AR/MR) technologies, serving as high-performance, compact, lightweight, and cost-effective optical components. However, the growing demands for optical transparency and high refractive index in advanced optical devices present a challenge for material discovery. In this study, we develop a novel approach that integrates first-principles calculations with genetic algorithms to accelerate the discovery of liquid crystal polymers with low visible absorption and high refractive index. By iterating within a predefined space of molecular building blocks, our approach rapidly identifies reactive mesogens that meet target specifications. Additionally, it provides valuable insights into the relationships between molecular structure and properties. This strategy not only accelerates material screening but also uncovers key molecular design principles, offering a systematic and scalable alternative to traditional trial-and-error methods.

---

1 Introduction

Liquid crystal polymers (LCPs) have played an important role in the development of advanced optical components. By leveraging the optical anisotropy of the liquid crystal molecules and their flexible alignment capabilities, liquid crystal polymer enable many thin-film based planar optical elements that are high-efficiency, compact, lightweight, and cost-effective, with examples including photochromic films, optical fibers, and polarization holographic gratings and lenses.^1–3 With the rapid advancement of display technologies and virtual, augmented, and mixed reality (VR/AR/MR) systems, liquid crystal polymer-based spatial light modulators and waveguides have attracted significant interest due to their critical functionality in these devices.^4–6

A key advantage of liquid crystal polymers in optical applications is their combination of low optical absorption in the visible range and high refractive indices, which improves light efficiency while minimizing component weight and volume.⁷ The fundamental building blocks of liquid crystal polymers, known as reactive mesogens (RMs), consist of a rigid mesogenic core, alkyl chain spacers, and polymerizable terminal groups, as depicted in Fig. 1(a). The mesogenic core retains liquid crystal properties due to its rod-like conjugated structure, which results in molecular anisotropy. Strong intermolecular interactions between mesogenic cores promote self-assembly and alignment in the liquid crystal phase. The length of the alkyl spacer can be adjusted to modify the liquid crystal phase transition temperature. Additionally, reactive mesogens act as monomers through polymerizable end groups, allowing polymerization by thermal or UV curing to form stable optical films.⁸ The unique structures of reactive mesogens enable tunability of both optical absorption and refractive index, providing a molecular design strategy for enhancing optical properties.

Fig. 1 (a) Structures of reactive mesogen. (b) Alignment of reactive mesogens in the nematic phase, followed by polymerization into a liquid crystal polymer network.

Over the past few decades, design and synthesis of reactive mesogens with desired optical properties remain a significant challenge due to both the theoretical complexity of molecular design and the inefficiency of material discovery.^9,10 Mesogenic cores with longer conjugation generally exhibit increased polarizability, enhancing the refractive index. However, this also tends to increase light absorption at longer wavelengths, as the delocalized electrons are more easily excited by lower energy.^11,12 The experimental process often relies on empirical methods and modifications of existing molecular structures, making it time-consuming and resource-intensive.^13–15 Although computational methods can accelerate molecular design compared to experimental approaches, they also face challenges. Unlike crystal materials with periodic structures, simulating the optical properties of liquid crystal polymers traditionally requires modeling multi-molecular systems, making it similarly time-consuming and unsuitable for large-scale screening. Variations in intramolecular and intermolecular interactions in liquid crystal polymers can significantly influence their optical properties. Therefore, accurate modeling largely relies on capturing the molecular ordering of the liquid crystal polymer system.¹⁶

To address these challenges, it is crucial to develop a pipeline that enables the rapid and accurate simulation of the optical properties of liquid crystal polymers, specifically focusing on absorption and refractive index. This pipeline would facilitate extensive screening of candidate structures, providing insights into structure–property relationships and opening new avenues in the molecular design space. Recent advancements in optical material exploration, particularly through machine learning and generative models, have provided valuable insights for liquid crystal polymer development.^17–19 However, a comprehensive liquid crystal material database with detailed optical data is still lacking, underscoring the need to explore the liquid crystal polymer space using advanced search and generative methods.

In this study, we develop a first-principles-based computational framework for large-scale screening of liquid crystal polymers with optical transparency and refractive index. By approximating the polymer network configuration using dimer conformations of reactive mesogens in nematic phase, we significantly reduce the computational cost of multi-molecular system simulations while retaining key molecular interactions. Further validation through molecular dynamics (MD) simulations and experimental data confirm the accuracy of this approach, supporting the dimer-based modeling strategy. We then integrate a genetic algorithm to screen, optimize, and generate novel liquid crystal polymer candidates. The resulting molecules demonstrate low absorbance in the visible wavelengths and exceptionally high refractive indices, offering valuable guidance for molecular design and material discovery. Furthermore, the computational pipeline developed in this research is adaptable to other organic molecule discovery platforms, particularly for polymer materials with specific optical properties. This adaptability holds benefits for a wide range of applications in optics, robotics, and other fields, providing a versatile tool for advancing material design and discovery.

2 Computational pipeline methodology and validation

To facilitate high-throughput computational molecular screening, it is essential to develop a streamlined model that incorporates necessary approximations, supported by robust scientific evidence and experimental validation. In this section, we first introduce liquid crystal polymer networks. The configuration of these networks is simplified as several dimer conformations, which are then used for absorption and refractive index calculations using time-dependent density functional theory (TD-DFT) and density functional theory (DFT), respectively.

2.1 Formation of liquid crystal polymer network

The formation of a liquid crystal polymer network, depicted in Fig. 1(b), begins by heating reactive mesogens to the liquid crystal phase to achieve alignment. This is followed by a curing step to form a cross-linked network, effectively locking the molecular alignment of reactive mesogens and preserving their conformations in the liquid crystal polymer. Consequently, the optical properties of the liquid crystal polymer are largely determined by the packing of the reactive mesogens.¹ Although minor differences in the order parameter and density between reactive mesogens and corresponding liquid crystal polymers arise due to spacer shrinkage before and after cross-linking process, these differences have a smaller impact compared to the influence of molecular arrangement on optical properties.^8,20 In this work, we use the conformational states of reactive mesogens in the nematic phase as an approximation for the polymer network when simulating the optical properties of liquid crystal polymers, as discussed in the following sections.

2.1.1 Dimer simplifications. The arrangement and packing of reactive mesogens are crucial in determining the optical properties of liquid crystal polymers. To improve computational efficiency while preserving accuracy, we further simplify the full assembly of reactive mesogens to a dimer system with various conformations. Dimers represent the minimal aggregate size that can capture nearest-neighbor π–π interactions and reproduce the characteristic spectral shifts, while avoiding the rapidly increasing computational cost and convergence difficulties of larger aggregates such as trimers or tetramers.^21–23 These dimer conformations represent different degrees of π–π stacking, and the dimer system energies adhere to the Boltzmann distribution. This approach is based on the Lebwohl–Lasher lattice model.²⁴ This simple spin lattice model, rooted in classical physics, is particularly suitable for high-throughput screening due to its computational efficiency. The model utilizes the Monte Carlo method to simulate phase behavior by stochastically altering the orientations of individual molecules, with the Boltzmann distribution governing the process. Studies have shown that it effectively captures the first-order phase transition between ordered and disordered states.²⁴

Initially, 50 reactive mesogen conformers were generated from the SMILES of input molecular structure via RDKit,^25,26 followed by structural optimization using semi-empirical GFN2-xTB method²⁷ with Crest.^28–30 Molclus isostat tool³¹ was used to merge identical conformers by applying a geometry threshold of 0.25 Å and an energy threshold of 0.5 kcal mol⁻¹. Next, 200 dimers were generated from the selected conformers, with inter-dimer distances controlled to within 2 van der Waals radii, and random relative spatial orientations were assigned using Molclus genmer tool.³¹ Dimers were further structurally optimized using the GFN2-xTB method in Crest, followed by the calculation of dimer probabilities at 400 K based on the Boltzmann distribution. Subsequently, the eight lowest energy dimers were identified among those with energy differences of less than 3 kcal mol⁻¹, serving as simplified computational units. This approach efficiently selects low energy conformers and captures the impact of intermolecular interactions on molecular geometry. Finally, quantum chemical calculations were performed to obtain the UV-vis spectra and refractive index of the selected dimers, and the averaged properties were used to approximate the absorption and refractive index of the polymer film. The dimer generation workflow is illustrated in Fig. 3, with the key steps highlighted in blue.

2.2 UV-vis spectra computation

2.2.1 Spectra calculation and experimental validation. TD-DFT has proven to be a reliable method for predicting the optical properties of organic molecules.^32–34 By applying traditional DFT to time-dependent systems, TD-DFT allows for accurate modeling of excited states, which are crucial for understanding electronic transitions in organic materials.³³ To evaluate the absorption of liquid crystal polymers using the dimer pipeline, UV-vis spectra for five commercial reactive mesogens, HCM-008 (RM257), HCM-009 (RM82), HCM-020 (RM23), HCM-021 (RM006/RM105), and HCM-083, were calculated with TD-DFT at the PBE0/6-31G(d,p) level. Concurrently, experimental measurements were conducted using UV-vis spectroscopy. The results show that the computed values deviate by an average of only 7 nm from the experimental data, indicating a high degree of consistency with experimental observations. All calculations were performed using the Gaussian 16.c software package,³⁵ and the simulated spectra are presented alongside the experimental data in Fig. S3.

As shown in the schematic drawing in Fig. 2(a), the conjugated mesogenic cores are distributed in a liquid crystal polymer network connected by non-conjugated spacers. This arrangement naturally gives rise to dimer conformations with varying degrees and geometries of π–π stacking. For clarity of discussion, we group these conformations into three representative types. Type A corresponds to configurations with negligible π–π overlap, typically resembling edge-to-face or highly displaced slip-stack arrangements. Type B represents dimers with partial orbital overlap, indicative of weak π–π interactions. Type C corresponds to nearly cofacial stacking with significant overlap between conjugated segments, resulting in the strongest π–π interactions. We emphasize that this classification is introduced only as a schematic framework to illustrate the diversity of possible packing geometries, rather than a strict energetic definition applied in our simulations.

Fig. 2 (a) Schematic illustration of a liquid crystal polymer network featuring three distinct dimer conformations, each representing varying degrees of π–π stacking interactions. (b) Structure of HCM-008 and its HOMO–LUMO orbitals. (c) UV absorption spectra of HCM-008 from experimental monomer solution, polymer thin film and computational results of different dimer types.

As observed in both computational and experimental data, different degrees of π–π stacking are critical in influencing optical responses, often leading to pronounced red-shifted absorption spectra or notable spectral broadening. To illustrate this effect in detail, we take HCM-008 as an example, as shown in Fig. 2(b). In HCM-008, the calculated orbital contributions indicate that the dominant allowed transition occurs at 250 nm, corresponding to a HOMO–LUMO transition identified as a π → π* excitation within the molecule. Fig. 2(c) illustrates that Type A, characterized by negligible π–π overlap, exhibits spectra closely resembling those of isolated monomers. By contrast, Type C, with nearly cofacial stacking, shows the most significant red shift in its absorption spectrum. The combined contributions from these diverse dimer configurations lead to the broadening of absorption peaks in polymer films. These findings highlight the importance of dimer conformations in determining the optical properties of liquid crystal polymers and support our approach of approximating UV-vis spectra through ensemble averaging over multiple dimer conformations of the mesogenic core.

Building on the averaged spectra strategy, we developed an automated computational pipeline (Fig. 3) to perform TD-DFT calculations on the top-ranked conformer candidates identified earlier, generating an averaged spectrum. This framework enables efficient, high-throughput molecular absorption calculations, facilitating subsequent genetic algorithm-based optimization. Beyond streamlining the computational process, this approach provides deeper insights into how mesogenic core interactions in the nematic phase influence the optical properties of liquid crystal polymers.

Fig. 3 Computational pipeline for dimer generation and optical property calculations.

2.2.2 Validation via molecular dynamics. To further validate the accuracy of our computational pipeline based on dimer simplification, we conducted additional MD simulations, which are known for their ability to capture complex molecular interactions and phase transitions.^36–38 The phase behavior of five reactive mesogens, HCM-008 (RM257), HCM-009 (RM82), HCM-020 (RM23), HCM-021 (RM006/RM105), and HCM-083, was captured using MD simulations using Schrödinger's Materials Science Suite with the OPLS4 force field (details in SI, Section 1).^39–41 In this process, all molecules successfully reproduced the nematic liquid crystal phase, validating the robustness of the protocol Fig. S3.

We extracted dimer conformations by randomly selecting 200 dimers from the initialized liquid crystal phase (details in SI, Section 1), ensuring an intermolecular distance of less than 4 Å based on the nearest-neighbor criterion. Single-point energy calculations were then performed on these dimers using GFN2-xTB in Crest. The energy distribution of the dimers was found to approximate a Gaussian distribution, providing additional evidence of the distinct intermolecular interactions characteristic of the liquid crystal phase (Fig. S2).

The average UV-vis spectra were compared with experimental data (Table S2). The results showed a strong correlation between the computed and experimental values, with an average deviation of only 6.2 nm. Similarly, TD-DFT calculations on the eight dimers generated via the designed dimer pipeline also yielded average UV-vis spectra that closely matched the experimental results, with the same average error of 9.4 nm.

Compared with the maximum absorption wavelength obtained from the MD-extracted dimer conformations, our designed dimer pipeline exhibits a red-shifted maximum absorption wavelength. This difference arises because the pipeline selects dimers based on their lowest energy conformations, resulting in stronger π–π stacking and lower energy than MD-generated dimers, which follow a Gaussian distribution. Although this discrepancy arises from the pipeline's inability to fully capture the MD-based Gaussian distribution, the error is acceptable given the substantial acceleration of spectral calculations achieved through dimer simplification. These dimers will also serve as a basis for further quantum chemical calculations of refractive indices, providing valuable support for subsequent material screening efforts.

2.3 Refractive index computation

The macroscopic refractive indices of liquid crystals are typically derived from the Lorentz–Lorenz equation. Due to the anisotropic nature of extraordinary (n_e) and ordinary (n_o) refractive indices, Vuks⁴² introduced a semi-empirical equation that incorporates the isotropic local field, relating microscopic molecular polarizability to the macroscopic refractive indices in anisotropic media.^43,44where N is the number of molecules per unit volume, α_e and α_o are the respective molecular polarizabilities, and 〈n²〉 is defined as:

The ratio of normalized polarizability can be expressed as:

Therefore, by calculating the molecular anisotropic polarizability and number density, the refractive index can be accurately predicted. The following sections provide a detailed explanation of how these results are calculated through first-principles methods, complemented by a calibration against experimental data.

2.3.1 Number density. The number density is typically determined from the density of MD simulation trajectory. However, this approach is both computationally expensive and time-consuming, posing challenges for high-throughput material screening. To address this, current methods that substitute molecular volume for MD-calculated density have demonstrated strong correlations with experimental results for amorphous polymers and liquid systems, offering a more cost-effective alternative to traditional MD simulations.⁴⁵

To validate that the molecular volume can be applied to the density calculations of liquid crystals, we investigated the temperature-dependent density behavior of the five reactive mesogens. The results showed that while density decreased with temperature, the slope of density change remained nearly constant (Fig. S5). This consistency arises from the highly uniform ellipsoidal shapes of reactive mesogen molecules, which, despite variations in molecular weight, exhibit similar packing coefficients in the nematic phase. This finding supports the feasibility of using molecular volume to accelerate number density calculations.

We employed Bader's definition of molecular volume in the condensed phase, where molecular volume is defined using the electron density isosurface as the van der Waals surface, accounting for volume deformation due to electronic effects.^46,47 The molecular volume was computed using the marching tetrahedra method on grid data, as implemented in Multiwfn software.^48,49 A comparison of the densities from the molecular volumes of monomers and dimers for the five reactive mesogens revealed a less than 2% difference in the densities obtained from wavefunctions (Table S3). Thus, using the molecular volume of monomers significantly enhances computational efficiency, while maintaining desired accuracy.

Subsequently, we compared the nematic phase densities at 298 K of 77 literature-reported liquid crystal molecules obtained via MD simulations with densities derived from molecular volumes calculated using DFT wavefunctions (Fig. S6). The results showed a strong correlation (R² = 0.92) between the two methods. This indicates that using the molecular volume of reactive mesogens can greatly simplify and accelerate the number density calculations.

2.3.2 Polarizability computation. Polarizability (α) of liquid crystal molecules can be derived from wave function-based methods and DFT.^50,51 By using the linear response theory, polarizability is computed by analyzing the system's reaction to an external perturbation, such as a constant electric field.⁵² The first-order response of the zero-field energy E(0) is the molecular dipole moment, reflecting the electron density distribution in polar molecules. The second-order response is the static polarizability tensor (α), which quantifies the deformation of the electron cloud under the external field. Due to the rod-like structure and π-conjugated backbone in liquid crystal molecules, the polarization along the molecular long axis is typically greater than that along the short axis. This structural anisotropy results in a directionally dependent dipole response to external electric fields. Therefore, the extraordinary polarizability α_e and ordinary polarizability α_o are defined as:

α_e = α_xx, (4)

where α_xx, α_yy, α_zz represent the polarizability components in the x, y, and z directions, respectively.

The polarizability is a frequency-dependent property. However, only the static polarizability was calculated to derive the refractive index under zero-field conditions. Given that the target liquid crystal molecules are transparent with negligible absorption in the visible region and lack low energy excited states, the frequency-dependent polarizability remains nearly constant or slightly monotonically decreases. It ultimately converges to the static polarizability value.⁵³ Consequently, the relative refractive indices of different molecular structures remain unaffected by frequency, thus not impacting the subsequent material screening process. Moreover, the computational cost of determining static polarizability is significantly lower than that of frequency-dependent polarizability, enabling the use of static values to approximate the experimentally measured refractive indices at specific wavelengths.

Accurately calculating the polarizability is a computationally intensive task that imposes strict requirements on the dispersion functions of the basis set and the selection of DFT functionals.^54–56 To reduce computational expense and improve efficiency, the ZPol basis set has been rigorously validated as a reliable choice to achieve moderate accuracy in the evaluation of molecular linear electrical properties.^57–59 It has demonstrated competitiveness with larger and more general basis sets while maintaining reliability in polarizability estimates for organic molecules.^60,61 To further enhance computational efficiency, we calculated the polarizability of dimers for the five reactive mesogens and their constituent monomers separately, with results summarized in Table S3. The findings show that the averaged polarizability of the two monomers deviates by less than 2% from that of the corresponding dimer, demonstrating that monomers can serve as a computationally efficient alternative for polarizability calculations. Furthermore, conformational differences between dimers and monomers contribute to an error of less than 2%. Therefore, we select the lowest energy dimer conformation, compute the polarizability of its individual reactive mesogens separately, and use the averaged value.

We used the M06-HF functional⁶² along with the ZPol basis set (containing both polarization and diffuse functions) to calculate the polarizabilities of 77 liquid crystal small molecules reported in the literature, and compared the results with experimental values, as shown in Fig. S7. All calculations were performed using the Gaussian 16.c software package.³⁵ The n_avg calculated by α_avg showed excellent agreement with the experimental results, while the calculated n_e was slightly higher and the n_o was slightly lower than the experimental values. This discrepancy arises because DFT calculations for individual molecules do not account for the influence of order parameter on n_e and n_o in liquid crystal polymer films. Consequently, in the subsequent genetic algorithm-based structure exploration, we used n_avg as the criterion to identify candidates with high refractive index.

In the high-throughput refractive index calculation pipeline, DFT calculations are performed on monomers derived from the lowest energy dimer configuration. Molecular volume is obtained from wavefunctions, and density is simulated using a calibrated packing coefficient. By integrating these parameters with polarizability calculations, the polymer film's refractive index is approximated. This approach optimally balances computational efficiency and accuracy, making it particularly suitable for high-throughput screening and the systematic exploration of liquid crystal materials.

3 Genetic algorithm design

A major challenge in designing novel materials for optical applications is identifying optimal molecular structures within the vast chemical space. To overcome this challenge, we use an inverse design strategy based on a genetic algorithm for global multi-objective optimization – a widely recognized approach for accelerating the discovery of organic materials with desired optical and electronic properties.^63–65 Compared to machine learning, genetic algorithms are computationally efficient, as they do not require large training datasets, and they offer a significant speed advantage over traditional experimental workflows. During multiple iterations of the evolutionary algorithm, the population of monomer structures undergoes computational crossover and mutation operations (mimicking biological processes), such as the exchange or substitution of molecular fragments. Ultimately, candidates that pass the initial screening are subjected to further property computations with first-principles methods in our computational pipeline. We note https://github.com/hutchisonlab/hutchison-polymer-GA as the basis for the implementation and adaptation of the genetic algorithm for this work.

The initial batch of generated reactive mesogen structures begins with the iterative process of the genetic algorithm. Using an established computational pipeline, UV-vis spectra and theoretical refractive index values are obtained for each molecule. A fitness function ranks the molecules, selecting those that best balance transparency and refractive index for the next iteration. New molecular structures are generated through crossover and mutation, and the process continues until the generated molecules meet the optical requirements. The iteration then terminates, yielding the desired reactive mesogen structures (Fig. 4(a)).

Fig. 4 (a) Design of reactive mesogen structures and the five steps of genetic algorithm for evolutionary iterations: initialization, selection, crossover, mutation, and termination. (b) Distribution of maximum absorption wavelength and n_e of molecules obtained from Fitness functions 1 and 2.

The highlight of our pipeline is the molecular generation from building blocks, which is unique to liquid crystals, which will be discussed in more detail in the following sections.

3.1 Molecular building blocks

The molecular building blocks include the mesogen core, an acrylate-linked alkyl chain, and a terminal substituent, as illustrated in Fig. 4(a). Each moiety types was identified and extracted from a liquid crystal molecule literature.

The database comprises four polymerizable alkyl chain building blocks incorporating different heteroatoms, 15 substituent building blocks positioned at the opposite end of the mesogen, and 52 mesogen core building blocks, including 16 linear bridges and 36 aromatic rings (Fig. S12). The length of the mesogen core is constrained to 3–7 units, as shorter cores exhibit insufficient anisotropy between their long and short axes, hindering liquid crystal phase formation. However, excessively long cores exhibit extended conjugation, which significantly reduces optical transparency. The unique design of these building blocks for reactive mesogens precisely captures key molecular structural features, enhancing the efficiency of the genetic algorithm by accelerating convergence toward target properties. The liquid crystalline phase behavior of the generated molecules is further validated through MD simulations, as discussed in subsequent sections.

3.2 Fitness functions

According to the sum-over-states (SOS) equations, achieving both high optical transparency and high polarizability in a single monomer is inherently limited by electronic delocalization. Molecules with higher transparency per unit volume typically exhibit lower polarizability, as the delocalized electrons are more easily excited by shorter wavelengths. Given these inherent trade-offs, where low absorption and high refractive index often conflict, optimizing these characteristics simultaneously becomes a complex challenge. Therefore, the fitness function plays a critical role in ranking, selecting, and filtering target structures, ensuring that the selected molecules are aligned with the desired target properties.

To optimize both optical transparency and refractive index, we designed two specialized fitness functions to guide the iterative selection process. The fitness functions consist of three components: the transparency score, the average refractive index value score, and the anisotropic refractive index difference score. By adjusting the weights of each component, we refine the ranking of candidates at each iteration, ensuring that only the most promising candidates progress to the next stage.

The transparency score is determined by the maximum absorption wavelength and the ratio of the absorbance of a molecule at its maximum absorption wavelength to its absorbance at 460 nm. Molecules with shorter maximum absorption wavelengths and lower absorbance at 460 nm receive higher transparency scores. The average refractive index score prioritizes molecules with higher average refractive indices, promoting the selection of high refractive index candidates. The anisotropic refractive index difference is included to avoid deviations from a rod-shaped molecular shape. Randomly generated reactive mesogen core building blocks may contain multiple bridge linkers, which might reduce molecular rigidity. This can hinder nematic phase packing or even prevent the molecule from being considered a liquid crystal monomer. Including the anisotropic refractive index difference in the fitness function helps favor rod-like molecules, making it more likely that they will go to the next iteration.

Fitness function 1 gives more weight to the transparency score. It starts with transparent molecules, especially those with a maximum absorption wavelength shorter than 460 nm, selecting candidates with higher refractive indices in subsequent iterations. Fitness function 2 focuses on molecules with refractive indices greater than 2.0, filtering for transparency in later stages. The specific formulas for both fitness functions can be found in SI Section 3. Although the two functions have different weightings, both prioritize molecules that achieve high refractive index and transparency, making it more likely for them to progress to the next iteration. As shown in Fig. 4(b), molecules optimized with Fitness function 1 tend to cluster in regions with low absorption and low n_e, while those optimized with Fitness function 2 are found in areas with higher n_e and closer to the visible absorption range. This highlights the challenge of balancing transparency and refractive index, and demonstrates how our dual fitness functions effectively guide molecules toward the desired optical properties. With each iteration, the molecules move closer to our target area.

4 Results and discussion

4.1 Genetic algorithm convergence

For both fitness functions, we performed five independent runs of the genetic algorithm, starting from a random initial seed of 32 molecules, each run lasting 30 epochs. To evaluate the convergence behavior of the genetic algorithm, we analyzed the evolution of fitness scores, refractive indices, and transparency scores over epochs for each run. We report the individual and combined behaviors of the runs shown in Fig. 5 and 6.

Fitness function 1 showed minimal variability across the five runs, and all runs demonstrated a clear convergence toward maximizing the fitness score. The maximum, minimum, and average scores of the population consistently increased, indicating that the genetic algorithm successfully identified molecules with both high transparency and refractive index (Fig. 5(a)). Transparency scores within the population also showed significant convergence toward the maximum value of 10, with high consistency across runs (Fig. 5(d)). Moreover, both n_e and n_avg showed an increase across epochs in all five runs, with the average values improving from 1.80 to 1.95 for n_e and from 1.51 to 1.58 for n_avg (Fig. 5(b) and (c)). In contrast, n_o, representing the refractive index along the short axis of the reactive mesogen, remained unchanged during the molecular design process (Fig. 5(e)). Given the competing relationship between transparency and refractive index, Fitness function 1 prioritizes maximizing molecular transparency. As a result, the improvement in transparency scores during the iterations far outpaced the increases in n_e and n_avg, which aligns with the expected outcome for the target molecules.

Fig. 5 Score and property distributions across all five independent genetic algorithm runs: (a and f) show fitness score, (b and g) show extraordinary refractive index n_e, (c and h) show average refractive index n_avg, (d and i) show transparency score, (e and j) show ordinary refractive index n_o for Fitness functions 1 and 2, respectively.

Fitness function 2 showed greater variability in the five independent runs compared to Fitness function 1. Both n_e and n_avg demonstrated increases during the initial ten epochs before convergence, with n_e increasing from 1.95 to 2.03 and n_avg increasing from 1.62 to 1.65 (Fig. 5(g) and (h)). Compared to Fitness function 1, Fitness function 2 achieved higher values for both n_e and n_avg. Obvious variability was observed in the fitness function score, with the mean value remaining stable at approximately 205 throughout iterations, while the maximum score in one run reached as high as 260 (Fig. 6(e)). Similar variability was observed for the transparency score, where the average value remained nearly constant, but the maximum value improved significantly between epochs. These distinct variations in transparency scores can be attributed to the constraint imposed by Fitness function 2, requiring n_avg to exceed 1.6. This restriction effectively narrows the search space, making advancements in transparency highly reliant on the stochastic discovery of top-performing candidates. As a result, the variability between different runs became more pronounced, with the success of the process relied on the inherent randomness of crossover and mutation in the genetic algorithm. Moreover, the design of the molecule along the short axis remained unchanged, so the n_o value showed little fluctuation (Fig. 5(j)). Despite the considerable variability in the optimization outcomes across different runs with Fitness function 2, the average refractive index between 5 runs consistently reached values above 2.0 for all runs. Furthermore, the refractive index demonstrated an upward trend with fluctuations in all three runs, aligning with our expectations.

Fig. 6 Metrics across each of the five independent genetic algorithm runs: (a and e) show fitness score, (b and f) show extraordinary refractive index n_e, (c and g) show average refractive index n_avg, (d and h) show transparency score for Fitness functions 1 and 2, respectively.

4.2 Building block frequencies

We analyzed the high-frequency building blocks in reactive mesogens optimized by the two fitness functions. In Fitness function 1, A17 (benzoxazole) and A20 (quinoline) consistently exhibited high occurrence in all five runs, whereas other aromatic building blocks showed high frequencies in only one run (Fig. 7(a)). In Fitness function 2, A12, A16, and A25 were dominant in at least one run, with frequencies exceeding 50% and remaining prevalent in others (Fig. 7(e)). Despite the differing selection criteria of the two fitness functions, both favored aromatic rings with multiple heteroatoms and fused ring systems. The presence of heteroatoms modifies the electron distribution, influencing dipole interactions and electronic transitions, which can enhance polarizability. Similarly, fused rings with extended π-conjugation facilitate electron delocalization, further improving polarizability. However, excessive ring fusion or the introduction of strongly electron-withdrawing heteroatoms can shift absorption into the visible range, reducing optical transparency. These high-frequency building blocks provide a balance between polarizability and transparency.

Fig. 7 Building block frequencies of each of five genetic algorithm runs with Fitness functions 1 and 2, including (a and e) aromatic rings, (b and f) linear bridges, (c and g) end groups, (d and h) tails.

For the bridge linkers in the mesogen core, C≡C triple bonds were highly frequent and consistent in all five runs in both fitness functions, especially in Fitness function 1, which prioritizes transparency. Their rigid structure stabilizes molecular geometry, creating a more structured backbone that enhances the packing efficiency of reactive mesogens. This improved packing promotes intermolecular interactions, contributing to greater order in the nematic phase. Additionally, the high bond energy and localization of π-electrons in triple bonds widen the HOMO–LUMO gap, shifting absorption into the UV region and avoiding visible absorption. As a result, C≡C triple bonds improve both transparency and structural rigidity, making them a preferred choice over other linear linkers. In contrast, Fitness function 2 favored azo (–N=N–) linkages and halogen-containing bridge linkers. Halogen substitution modifies electron density distribution and can enhance polarizability via induced dipole effects. However, azo bonds are not strictly linear and can undergo cis–trans isomerization under light exposure, potentially disrupting the nematic phase stability. Consequently, they were less effective in Fitness function 1, which prioritizes optical transparency. This highlights C≡C triple bonds as an optimal choice, offering a balance between high polarizability, rigidity, and excellent optical transparency.

The –N=C=S end-group, prevalent across all optimization runs, significantly improves both polarizability and transparency. The electron-withdrawing sulfur and nitrogen atoms, by reducing electron density, promote electron delocalization within the conjugated system, thereby enhancing the polarizability. Furthermore, conjugation between the –C=N and –C=S groups further facilitates electron delocalization, modifying the molecular electronic structure. This change in electronic structure shifts the absorption to the UV region, as the modified energy levels make electronic transitions more likely to occur in the UV range. Consequently, the –N=C=S end-group not only enhances polarizability but also maintains transparency, making it an ideal choice for the design of liquid crystal molecules.

For flexible chains with different heteroatoms, the results across all five runs were consistent in both Fitness functions 1 and 2. Molecules substituted with oxygen and nitrogen contributed more to polarizability, while sulfur played a more significant role in enhancing transparency. Sulfur, with its larger size and lower electronegativity, participating in conjugation through electron delocalization. This electron donation helps distribute electron density more evenly across the molecule, avoiding electronic transitions in the visible range. As a result, sulfur-containing flexible chains are more transparent, as their electronic transitions tend to shift into the UV range. In contrast, carbon atoms in the flexible side chain had the weakest impact on both polarizability and transparency across both fitness functions, likely due to their lower electron-donating ability and smaller size, which make them less effective at stabilizing the structure or shifting absorption toward shorter wavelengths.

These outcomes may be attributed to the inherent stochasticity of genetic algorithms, which require a greater number of iterations and repeated experiments to avoid local minima. The current limited iterations and repetitions may have overlooked top candidates. However, within the limited number of runs and epochs, we have still identified several consistent candidates, which have been highly valuable for the molecular design.

4.3 Generated material validations

To verify the existence of the nematic phase in the generated reactive mesogen structures, we selected candidates that simultaneously satisfied two criteria: n_e > 2.0 and absorption intensity below 0.1% at 460 nm relative to the maximum peak. After removing duplicates and candidates with identical main-chain building blocks, 10 unique candidates with distinct main-chain cores remained for MD validation, ensuring representative coverage (Fig. S10). The results showed that six candidates exhibited a nematic phase within the temperature range of 200 K to 700 K. Candidate 1 did not undergo a nematic phase transition and maintained an order parameter below 0.2. This result can be attributed to the elongated aspect ratio of the mesogen core, which is entirely composed of C≡C triple bond bridge linkers (Fig. S10). Compared to cores formed by aromatic rings, linear bridges produce a much smaller radius along the n_o direction, hindering the formation of ordered structures during the cooling process. In contrast, Candidates 8–10 exhibited phase transition temperatures that exceeded 700 K. These candidates consistently incorporated the fused ring building block and the mesogen core composed of four building blocks, imparting greater rigidity and enhanced π–π stacking interactions. Consequently, the increased molecular rigidity and improved stacking efficiency substantially elevated the transition temperatures of these candidates.

Among the six molecules exhibiting liquid crystalline phases, Candidate 2 was selected for polarizability density analysis, as its transition temperature is the lowest and falls within an experimentally accessible range. The polarizability density of Candidate 2 is primarily concentrated near the triple bonds and the thiocyanate group at one end of the molecule, as well as around the sulfur atom at the opposite end. This localization arises from the high electron density and polarizability of the triple bonds and the strong dipole moment of the thiocyanate group, which enhances the molecule's responsiveness to external electric fields. Similarly, the sulfur atom contributes significantly to the polarizability due to its large atomic radius and high polarizability, highlighting its role in influencing the optical properties.

Subsequently, we employed the TD-DFT method to calculate the excitation energies and oscillator strengths of the first 20 excited states of Candidate 2. The maximum absorption peak was observed at 365 nm (Fig. 8). A strong absorption band is attributed to the doubly degenerate HOMO−1 → LUMO and HOMO → LUMO+1, which correspond to the π → π* transition, while all other excitations contribute negligibly to the absorption spectrum (Fig. S11). As no significant absorption band is observed in the visible region, the polymer formed by Candidate 2 is expected to be colorless.

Fig. 8 (a) Order parameter as a function of temperature for Candidate 2 during the MD cooling process. (b) Polarizability density distribution of Candidate 2, with the green and blue isosurfaces representing positive and negative values, respectively. (c) Average and individual spectra of eight dimers of Candidate 2 selected through the dimer generation pipeline.

5 Conclusion

In this study, we developed a material discovery pipeline for liquid crystal molecules that integrates the semi-empirical GFN2-xTB method with first-principles (TD-)DFT calculations. This integration enables rapid iterations within a genetic algorithm framework, allowing for efficient exploration of the chemical search space. Our approach successfully identifies potential reactive mesogen structures characterized by high transparency and high refractive indices.

The dimer generation pipeline incorporates the influence of liquid crystal molecular packing. We have validated the accuracy and efficiency of this pipeline by comparing experimental and calculated values for five commercial reactive mesogens. These comparisons demonstrate the pipeline's capability to reliably predict material properties, thereby supporting its use in the discovery of new liquid crystal materials.

The genetic algorithm demonstrated its effectiveness as a methodology for optimizing material absorption and refractive index by efficiently exploring and targeting specific regions within the chemical search space. Despite the inherent randomness of genetic algorithms, it consistently identified recurring high-frequency building blocks and specific structural arrangements across independent runs. This consistency suggests that increasing the number of iterations and independent runs will further enhance the robustness and reliability of the results, solidifying the effectiveness of this design strategy.

Our computational analysis revealed that the 10 screened candidates exhibit less than 0.1% absorption at 460 nm compared to the maximum absorption peak absorbance, with n_e values exceeding 2.0. The highest candidate shows an n_e value exceeding 2.1, and all candidates were further validated with MD simulations to confirm the existence of the nematic phase. In cases with lower transparency constraints or additional side-chain modifications, we anticipate even higher n_e values for our candidates.

Future research could focus on improving the computational efficiency of the pipeline and advancing the prediction of nematic phase transitions. These improvements would enable deeper exploration of the chemical space for materials with enhanced optical properties. In addition, our physics-driven workflow, validated at the quantum chemical level, can serve as a complementary source of high-quality datasets for training and benchmarking machine learning/deep learning models. Such integration of physics-based validation with data-driven approaches holds promise for expanding the potential applications of liquid crystal-based optical materials and accelerating the development of high-performance optical devices.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: additional figures, tables, and experimental details supporting this study. The liquid crystal refractive index dataset referenced in the main text is accessible at attach file. See DOI: https://doi.org/10.1039/d5ra04477d.

References

1. V. Shibaev, A. Bobrovsky and N. Boiko, Prog. Polym. Sci., 2003, 28, 729–836.
2. Y. Li and D. Luo, Opt. Mater. Express, 2016, 6, 691–696.
3. P. Chen, B.-Y. Wei, W. Hu and Y.-Q. Lu, Adv. Mater., 2020, 32, 1903665.
4. J. Kim, Y. Li, M. N. Miskiewicz, C. Oh, M. W. Kudenov and M. J. Escuti, Optica, 2015, 2, 958–964.
5. Y.-H. Lee, K. Yin and S.-T. Wu, Opt. Express, 2017, 25, 27008–27014.
6. J. Xiong and S.-T. Wu, eLight, 2021, 1, 3.
7. X. Zhang, J. Wang, Z. Han, Y. Huang, M. Wang, X. Feng, L. Lu and B. Silverstein, SID Symposium Digest of Technical Papers, 2024, pp. 1163–1166.
8. D. Liu and D. J. Broer, Langmuir, 2014, 30, 13499–13509.
9. Q. Wei, R. Pötzsch, X. Liu, H. Komber, A. Kiriy, B. Voit, P.-A. Will, S. Lenk and S. Reineke, Adv. Funct. Mater., 2016, 26, 2545–2553.
10. J.-g. Liu and M. Ueda, J. Mater. Chem., 2009, 19, 8907–8919.
11. K. Sasagane, F. Aiga and R. Itoh, J. Chem. Phys., 1993, 99, 3738–3778.
12. J. Li, C.-H. Wen, S. Gauza, R. Lu and S.-T. Wu, J. Disp. Technol., 2005, 1, 51.
13. C. Sekine, K. Iwakura, M. Minai and K. Fujisawa, Liq. Cryst., 2001, 28, 1505–1512.
14. R. He, E. Oh, Y. Ye, P. Wen, K.-U. Jeong, S. H. Lee, X.-D. Li and M.-H. Lee, Polymer, 2019, 176, 51–59.
15. J. D. Allen, K. Adlem and M. Heeney, J. Mater. Chem. C, 2021, 9, 17419–17426.
16. M. R. Wilson, Chem. Soc. Rev., 2007, 36, 1881–1888.
17. H. S. Stein, D. Guevarra, P. F. Newhouse, E. Soedarmadji and J. M. Gregoire, Chem. Sci., 2019, 10, 47–55.
18. L. Simine, T. C. Allen and P. J. Rossky, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 13945–13948.
19. G. Chen, Z. Shen, A. Iyer, U. F. Ghumman, S. Tang, J. Bi, W. Chen and Y. Li, Polymers, 2020, 12, 163.
20. E. J. Cooper, M. Reynolds, T. Raistrick, S. R. Berrow, E. I. Jull, V. Reshetnyak, D. Mistry and H. F. Gleeson, Macromolecules, 2024, 57, 2030–2038.
21. B. Tirri, G. Mazzone, A. Ottochian, J. Gomar, U. Raucci, C. Adamo and I. Ciofini, J. Comput. Chem., 2021, 42, 1054–1063.
22. M. Breza, Comput. Theor. Chem., 2018, 1143, 1–8.
23. C. Allolio, T. Stangl, T. Eder, D. Schmitz, J. Vogelsang, S. Hoger, D. Horinek and J. M. Lupton, J. Phys. Chem. B, 2018, 122, 6431–6441.
24. P. A. Lebwohl and G. Lasher, Phys. Rev. A: At., Mol., Opt. Phys., 1972, 6, 426.
25. D. Weininger, J. Chem. Inf. Comput. Sci., 1988, 28, 31–36.
26. G. Landrum, RDKit: Open-Source Cheminformatics, http://www.rdkit.org.
27. C. Bannwarth, S. Ehlert and S. Grimme, J. Chem. Theory Comput., 2019, 15, 1652–1671.
28. P. Pracht, F. Bohle and S. Grimme, Phys. Chem. Chem. Phys., 2020, 22, 7169–7192.
29. P. Pracht, S. Grimme, C. Bannwarth, F. Bohle, S. Ehlert, G. Feldmann, J. Gorges, M. Müller, T. Neudecker and C. Plett, et al., J. Chem. Phys., 2024, 160, 114110.
30. S. Grimme, J. Chem. Theory Comput., 2019, 15, 2847–2862.
31. T. Lu, Molclus program, Version 1.12, http://www.keinsci.com/research/molclus.html.
32. M. E. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem., 2012, 63, 287–323.
33. C. Adamo and D. Jacquemin, Chem. Soc. Rev., 2013, 42, 845–856.
34. R. Chen, Z. An, W. Wang, X. Chen and P. Chen, Liq. Cryst., 2017, 44, 1515–1524.
35. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford, CT, 2016.
36. M. P. Allen, Mol. Phys., 2019, 117, 2391–2417.
37. M. R. Wilson, Int. Rev. Phys. Chem., 2005, 24, 421–455.
38. C. McBride and M. R. Wilson, Mol. Phys., 1999, 97, 511–522.
39. C. Lu, C. Wu, D. Ghoreishi, W. Chen, L. Wang, W. Damm, G. A. Ross, M. K. Dahlgren, E. Russell and C. D. Von Bargen, et al., J. Chem. Theory Comput., 2021, 17, 4291–4300.
40. Schrödinger, LLC, Schrödinger Release 2024-3: Materials Science Suite, 2024.
41. Schrödinger, LLC, Schrödinger Release 2024-3: Force Fields, 2024.
42. M. Vuks, Opt. Spectrosc., 1966, 20, 361.
43. M. Zakerhamidi and H. Rahimzadeh, J. Mol. Liq., 2012, 172, 41–45.
44. Z. Alipanah, M. S. Zakerhamidi and A. Ranjkesh, Sci. Rep., 2022, 12, 12676.
45. M. A. F. Afzal, C. Cheng and J. Hachmann, J. Chem. Phys., 2018, 148, 241712.
46. R. F. Bader, M. T. Carroll, J. R. Cheeseman and C. Chang, J. Am. Chem. Soc., 1987, 109, 7968–7979.
47. Z. Liu, T. Lu and Q. Chen, Carbon, 2021, 171, 514–523.
48. T. Lu and F. Chen, J. Mol. Graphics Modell., 2012, 38, 314–323.
49. T. Lu, J. Chem. Phys., 2024, 161, 082503.
50. T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen and K. Ruud, Chem. Rev., 2012, 112, 543–631.
51. D. Hait and M. Head-Gordon, Phys. Chem. Chem. Phys., 2018, 20, 19800–19810.
52. M. Stott and E. Zaremba, Phys. Rev. A: At., Mol., Opt. Phys., 1980, 21, 12.
53. B. Jansik, D. Jonsson, P. Sałek and H. Ågren, J. Chem. Phys., 2004, 121, 7595–7600.
54. A. Baranowska-Łączkowska, B. Fernández, A. Rizzo and F. Pawłowski, Mol. Phys., 2015, 113, 1786–1793.
55. A. L. Hickey and C. N. Rowley, J. Phys. Chem. A, 2014, 118, 3678–3687.
56. M. A. F. Afzal and J. Hachmann, Phys. Chem. Chem. Phys., 2019, 21, 4452–4460.
57. A. Baranowska, M. Siedlecka and A. J. Sadlej, Theor. Chem. Acc., 2007, 118, 959–972.
58. Z. Benkova, A. J. Sadlej, R. E. Oakes and S. E. Bell, J. Comput. Chem., 2005, 26, 145–153.
59. Z. Benkova, A. J. Sadlej, R. E. Oakes and S. E. Bell, Theor. Chem. Acc., 2005, 113, 238–247.
60. S. Lee, S. S. Park and F. Hagelberg, Chem. Phys., 2012, 394, 40–45.
61. A. Baranowska-Łączkowska, K. Z. Łączkowski and B. Fernandez, J. Chem. Phys., 2012, 136, 024302.
62. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2006, 110, 13126–13130.
63. C. Suh, C. Fare, J. A. Warren and E. O. Pyzer-Knapp, Annu. Rev. Mater. Res., 2020, 50, 1–25.
64. D. C. Hiener and G. R. Hutchison, J. Phys. Chem. A, 2022, 126, 2750–2760.
65. O. D. Abarbanel and G. R. Hutchison, Phys. Chem. Chem. Phys., 2023, 25, 11278–11285.

---

Footnotes

† Work done during Meta internship.

‡ Work done while at Meta.

§ Joint last author.

---