Automatic Determination of $n$-Cyanobiphenyl Elastic Constants from Molecular Simulation

TL;DR

This paper presents an automated molecular simulation method to accurately determine the elastic constants of n-cyanobiphenyl liquid crystals, aligning well with experimental data and enabling future molecular engineering.

Contribution

It applies and extends existing techniques to the nCB series, providing an automated approach for calculating elastic properties of liquid crystals from molecular models.

Findings

Strong agreement with experimental elastic constants for nCBs.

Automated process enables efficient elasticity calculations for new molecules.

Path established for computational design of liquid crystal materials.

Abstract

New applications of liquid crystalline materials have increased the need for precise engineering of elastic properties. Recently, Sidky et al. presented methods by which the elastic coefficients of molecular models with atomistic detail can be accurately calculated, demonstrating the result for the ubiquitous mesogen 5CB. In this work, these techniques are applied to the homologous series of nCB materials, focusing on the standard bend, twist, and splay deformations, using an entirely automated process. Our results show strong agreement with published experimental measurements for the nCBs and present a path forward to computational molecular engineering of liquid crystal elasticity for novel molecules and mixtures.

Equations

$$f=\frac{1}{2}{k}_{11}{\left(\nabla \cdot \hat{\mathbf{n}}\right)}^2$$

$$+\frac{1}{2}\left(k_{22}+k_{24}\right)\left[\mathrm{Tr}{\left(\nabla \hat{\mathbf{n}}\right)}^2-{\left(\nabla \cdot \hat{\mathbf{n}}\right)}^2\right].$$

Full text

Introduction

Liquid crystals (LC) are a class of fluids which exhibit long-range orientational ordering and a resistance to spontaneous deformation. Kleman and Lavrentovich (2003) While traditionally used in display technologies Gray and Kelly (1999), advances in both our understanding of elasticity and material fabrication has opened up many new and interesting applications where the unique properties of LCs can be exploited. Chemoresponsive LCs have been designed to respond to targeted chemical species, which can be quantified through polarized optical microscopy. Szilvási et al. (2017) LC elastomeric materials can function as elements in soft robots and sensors. Babakhanova et al. (2018) Topological defects and inclusion formation in LCs have also been used as templates for molecular self-assembly. Wang et al. (2015) Many of these novel applications rely on the interplay between external, boundary and elastic restoring forces, which is dictated by both the absolute magnitudes and ratios of the elastic constants.

Advances in atomistic simulations have made it possible to gain an unprecedented level of detail into the molecular underpinnings of liquid crystalline behavior which can be rigorously verified against experiment. Palermo et al. (2013); Ramezani-Dakhel et al. (2017) By understanding the role of molecular features such as conformational flexibility and charge distribution, new LCs can be engineered in silico to have highly specific thermophysical properties necessary for biosensing, templating, and other emerging technologies. Until very recently Sidky et al. (2018a), it was not possible to economically predict the elastic properties of molecular LCs and thus made it extremely difficult to engineer novel LCs for these applications.

The work of Sidky, et al. Sidky et al. (2018a) was a crucial first step in this process, as there it was demonstrated how molecular simulations may be used as a tool to predict the bulk elastic moduli of the molecular liquid crystal models, focusing on the ubiquitous 4-pentyl-4’-cyanobiphenyl, more commonly known as 5CB. There, it was also shown that the saddle-splay elastic constant $k_{24}$, difficult to obtain in experiments, could be directly measured, and this revealed the underlying positive definite nature of the mode in contradiction to recent experimental suggestions based on indirect observation. Allender et al. (1991); Sparavigna et al. (1994); Polak et al. (1994); Pairam et al. (2013); Davidson et al. (2015) The stability of 5CB and its homologues, in addition to their achromic appearance and strong positive dielectric anisotropy, has made them widely used in research Wang et al. (2016); Eimura et al. (2016) and industrial applications. Gray and Kelly (1999); Kleman and Lavrentovich (2003) Importantly, their elastic ˘properties have been well-characterized through experiment. Madhusudana and Pratibha (1982); Hakemi et al. (1983); Chen et al. (1989); Chatopadhayay and Roy (1994) Molecular models have also been parameterized to reproduce basic thermodynamic properties such as density, orientational order, and phase transitions. Tiberio et al. (2009); Boyd and Wilson (2015); Cacelli et al. (2007, 2009)

The focus of this work is on applying the methodology proposed in Ref. Sidky et al., 2018a to the $n$-cyanobiphenyls ($n$CBs) 6CB, 7CB, and 8CB, which are illustrated in Figure 1. Although there are minimal structural differences between this series of molecules, the extended alkyl chain has a strong effect on ordering, which gives rise to the “odd-even” effect as the homologous series is ascended. Madhusudana and Pratibha (1982); Hakemi et al. (1983) This phenomenon, where the molecular shape anisotropy and corresponding nematic transition temperature increases only when $n$ goes from even to odd, has been observed in both experiments and simulations.Tiberio et al. (2009); Gray and Harrison (1971) While the number $n$ is arbitrary, for $n<5$ a nematic phase is not observed, while for $n>7$, the LCs admits additional smectic behavior. Here, we seek to predict the splay, twist, and bend elastic moduli of $n$CB nematogens through an automated version of the prior algorithms Sidky et al. (2018a) in order to test the limits and robustness of our methods. This represents a step towards the the ultimate goal of computationally-guided LC engineering, and extends the characterization of the elastic properties of existing LC forcefields. We restrict ourselves to the three bulk elastic constants because they are directly measurable in experiment and their behavior is generally well-understood.

Methods

The computational methods used here follow Ref. Sidky et al., 2018a closely. For completeness, we will provide a brief summary and highlight the differences. The linear-order elastic free energy density of a uniaxial nematic, absent of any molecular polarity or chirality, may be represented in the Frank-Oseen form Stewart (2004)

[TABLE]

This expression contains the most commonly used elastic terms: splay ($k_{11}$), twist ($k_{22}$) and bend ($k_{33}$). The additional term, referred to as ‘saddle-splay’, depends on $k_{24}$ and penalizes bidirectional deformations. While the prior work demonstrated this can be measured, the process requires many more molecules to set up the appropriate anchoring and boundary conditions for the problem, it is thus significantly slower. Notably, no measurements of $k_{24}$ exist for $n$CB molecules other than 5CB. We thus focus here on the standard modes of deformation for which data can be readily obtained.

The $n$CB molecules are represented using the united atom forcefield parameterized by Tiberio et al. Tiberio et al. (2009) After preparing an initial configuration for each molecule type, the process of obtaining elastic constant estimates is entirely automated. Each cyanobiphenyl system contains 400 molecules, and is initially prepared by running 100 nanosecond NPT simulations at 1 atm in Gromacs 5.1.3 Abraham et al. (2015). The temperatures examined are between $T_{\mathrm{NI}}$ and $T_{\mathrm{NI}}-$20\text{,}\mathrm{K}$$. The transition temperature is independently determined by performing an initial simulation sweep centered on the values reported in Ref. Tiberio et al., 2009, and determining when the equilibrium value of the nematic order parameter $S$ exceeds 0.1; the measured $T_{\rm NI}$ correspond to those reported in the previous paper. The equilibration time is considerably less than the 400 ns used for 5CB in our previous work, though a reasonable length is necessary to allow for elastic constant prediction at this scale. A Langevin thermostat and Parrinello–Rahman barostat with $\tau_{\mathrm{p}}=5$ picoseconds are used with a time step of 2 femtoseconds for all simulations.

In accordance with Ref. Sidky et al., 2018a, the average volume at the each temperature is obtained in the initial round of simulations, and configurations with this volume are automatically drawn from the completed $NPT$ trajectories and used to initiate a second round of $NVT$ simulations with edge restrictions applied using a harmonic restraint having modulus $k=10^{5}$ kJ/mol which serves to align the molecules along the $\hat{\mathbf{z}}$ axis. Four uncorrelated instances are generated at each simulated temperature to function as independent walkers during the elastic measurement, and enhance the convergence of the free energy measurement. The simulations at this step were carried out for $400\text{\,}\mathrm{ns}$, considerably less than the $1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ timescale of the previous study Sidky et al. (2018a). We arrived at this number through experimentation, but feel that it is entirely appropriate for most calamitic-type nematogens. Importantly, this step may be trivially parallelized to increase the speed of convergence to the underlying free energy landscape and improve the accuracy of the calculations.

The final step involves estimating the elastic constants from a biased simulation. With edge restrictions in place aligning the molecules along the $\hat{\mathbf{z}}$ axis, a deformation $\xi$ is applied in the central region according to the elastic mode desired. We take $\xi$ to be $\partial n_{x}/\partial x$ for splay, $\partial n_{y}/\partial x$ for twist, and $\partial n_{z}/\partial z$ for bend. Parabolic free energy profiles obtained via basis function sampling Whitmer et al. (2014) (BFS) as implemented in SSAGES 0.6 Sidky et al. (2018b) are used to extract the final elastic constant as $f=\frac{1}{2}\gamma k_{ii}\xi^{2}$ where $\gamma$ is a geometric factor accounting for the finite restriction and deformation regions. To assist readers in reproducing our results, we have posted all scripts and Gromacs and SSAGES runfiles in a free online repository located at https://github.com/hsidky/atomistic_elastics.

Results and Discussion

While data for $n$CB molecules is more sparse than that for 5CB, a few measurements are available in the literature. We compare all of our calculations to experimental data from Refs. Madhusudana and Pratibha, 1982,and Chatopadhayay and Roy, 1994; other measurements such as those reported in Ref. Hakemi et al., 1983 were not used due to substantial deviations from other reported measurements. Figure 2 shows the elastic constant predictions for 6CB. The estimated values are in good agreement with experiment across the reported temperature range. Values for $k_{33}$ and $k_{11}$ fall within the typical range of variability of experimental measurements reported for 6CB. Chen et al. (1989); Madhusudana and Pratibha (1982) Twist elastic constants show a systematic positive drift away from experimental data as temperature decreases, but remain significantly smaller than the other moduli, as expected, which results in a negative deviation of $k_{ii}/k_{22}$. For this measurement, there is minimal noise observed in our data for $k_{ii}$ with decreasing $T-T_{\rm NI}$, and the trendline is reasonably smooth and monotonic.

The elastic constants for 7CB are shown in Figure 3. They are in excellent agreement with experiment across a broader temperature range than 6CB, with only a minor negative deviation of $k_{11}$ at lower temperature. Compared to 6CB, the $k_{22}$ is exhibits a similar temperature dependence, with significantly higher $k_{33}$ and $k_{11}$ as temperature decreases. The deviations between experimental and simulation behavior at low temperatures could be due to a form of entropy–enthalpy tradeoffChaimovich and Shell (2010) within the model of Ref. Tiberio et al., 2009. These deviations are neverthless small compared to the magnitude of the elastic constants $k_{ii}$. We note that although quantitative agreement appears to be better than our measurements for 6CB, the measured elastic constants are also noisier, with data variability exceeding the uncertainty bounds associated with each individual measurement. This can be understood by clarifying precisely how the uncertainty estimates were generated, the inherent limitations, and how they an be addressed.

As previously mentioned, each free energy estimate and corresponding elastic constant is obtained using BFS which projects a running estimate of the partition function onto a truncated basis set. For our purposes, four independent walkers all contribute to the same free energy estimate to accelerate convergence. Here the uncertainty is estimated from higher order terms truncated in the projection, which uses the second order Legendre polynomial, $P_{2}$. In other words, the uncertainty is a measure of the “goodness of fit” of the discrete biased data to $P_{2}$. BFS, like many other free energy sampling methods, is a convergent algorithm, meaning that later statistics added into the free energy estimate contribute exponentially less than early estimates. For slow-relaxing collective variables such as the elastic deformation used in this work, noisy early estimates can overemphasize early, noisy data that becomes more difficult to correct over time. It is possible that this issue can be alleviated utilizing new, neural network-based methods Sidky and Whitmer (2018) which obtain an initial estimate of the surface very quickly from a small number of counts, or by using “forgetting” methods Whitmer et al. (2014) which discard a specified fraction of the previous bias history to limit the influence of early state visits, though these options were not explored in this study.

Ideally, one can also run independent free energy simulations with different starting configurations and seeds, and to average the resulting statistics into a single measurement exhibiting upper and lower bounds. At this stage, running a sufficient number of copies to be confident about the average and error estimate makes these simulations prohibitively expensive—the “wall” time incurred in the simulations as described in the Methods section will increase from $\approx 1$ month of computing time to $\approx 1$ year or more. An additional potential source of error is that our elastic constant measurements are carried out in the NVT ensemble, which introduces additional, albeit less apparently significant, error, due to our inability to relax the box dimensions in response to applied orientational stresses. This is a limitation of the deformation collective variable which does not admit a simple virial expression and thus cannot be run in NPT simulations. This could be addressed using augmented, joint density-of-states methods where the box dimensions are additionally sampled via a coupled Monte-Carlo algorithm. Shell et al. (2002) Finally, the cumulative simulation time for systems in this work were a factor of 5 less than those for 5CB Sidky et al. (2018a), which appears to be a significant factor in overcoming the initial noisy estimate mentioned. Average wall time was still substantial, considering the amount of data reported in this work; an cumulative average of 45 days was needed to obtain an elastic constant estimate for each $k_{ii}$ at each temperature on 24-core Intel Haswell nodes.111Computers used in this study are 24-core systems comprised of two 12-core Haswell E5-2680v3 processors equipped with a total of 32GB DRAM. Computational time should be significantly less on updated hardware. Nonetheless, some of the measurement uncertainty is captured in the reported data, and quantitative accuracy and trends are not compromised.

Figure 4 shows the elastic moduli for 8CB. It is well known that 8CB exhibits a smectic-A–nematic transition at approximately $7\text{\,}\mathrm{K}$ below $T_{NI}$. which results in a divergence of $k_{22}$ and $k_{33}$. Madhusudana and Pratibha (1982) The formation of layers along the nematic axis permits only deformations which are curl-free. de Gennes (1972) This gives rise to a power-law scaling relation with a critical exponent of $\nu\approx 0.67$. Hakemi (1989) Additionally, the model used in this study has been shown to accurately capture the both transition temperatures and smectic layering behavior Palermo et al. (2013), although with strong finite size dependence. Our results across the nematic phase, shown in the Figure 4 inset, are an excellent match to the experimental results. However, as the experimental $k_{22}$ begins to diverge through the smectic-nematic transition, we see that our simulated results do not exhibit the same behavior. Plotted across the entire range of experimental measurements, it is clear that our measurements are not indicative of the presence of a smectic phase within our small samples.

This apparent discrepancy is due to a number of effects. The first is a strong finite size dependence which has already been reported previously. Palermo et al. (2013) Compared to larger systems, those containing fewer than $N=1000$ mesogens, which is our case, exhibit a considerable discrepancy in smectic ordering. Palermo et al. (2013) Similar behaviors were previously observed in combined phase and elastic property simulations of Gay –Berne models.Joshi et al. (2014) The second is that we force the alignment of the 400 mesogens along the $\hat{\mathbf{z}}$ axis which may not be the preferred orientation of layers in a fixed-size periodic box. Thirdly, deformations to the mesogens are applied without considering the presence of smectic layering. Finally, simulations are carried out in the NVT ensemble which does not allow collective orientational relaxation. All of these factors contribute to the suppression and disruption of smectic layering which would manifest in the divergent elasticities absent from our predictions. To properly accommodate and characterize natural smectic layering a semi-isotropic barostat would be necessary, allowing the box dimension ratios to adopt those corresponding to layer separation. The difficulties associated with this within the context of our elastic measurements have been discussed above. In principle, however, the methods we present may be extended to study this intriguing behavior in systems on the order of $N=3000$ mesogens Palermo et al. (2013).

Conclusion

In this work we have estimated the elastic constants of the cyanobiphenyls 6CB, 7CB, and 8CB from biased molecular simulation using an entirely automated procedure. We show that the simulated elastic moduli are generally in good agreement with reported experimental values. Difficulties arise in generating meaningful uncertainties and dealing with the onset of the smectic phase for 8CB. Both improving the prediction accuracy and applying this methodology on a large scale are only possible if the cost of elasticity measurements is substantially reduced. Currently, up to 4 $\mu s$ of total simulation runtime is required to generate a estimate of a single elastic mode at a single temperature.

Developing newer or better optimized sampling algorithms may significantly reduce this time and open the door to in silico liquid crystal engineering. This would necessitate some means of attenuating early-time noise present in elastic measurements, and perhaps a simplification of the order-parameter calculation which is computationally expensive and involves third-order tensor algebra. As presented, the methodology developed is capable of producing accurate elastic constants and merits further study to elucidate the precise limit of efficiency while retaining measurement fidelity. Further, a crucial next step involves moving beyond cyanobiphenyls and examine atomistic models of the myriad liquid crystal phases, from molecules of historical interest, such as PAA and MBBAHaller (1972); Jeu et al. (1976) to more exotic bent-core molecules and models exhibiting chiralDierking (2014) or twist-bend phases.Chen et al. (2013) Improving the breadth, speed, and accuracy of these techniques will enable their use in computational screening of new mesogenic compounds.

Code Availability

All scripts and information necessary to run the examples contained in this article are posted at https://github.com/hsidky/atomistic_elastics.

Acknowledgement

HS acknowledges support from the National Science Foundation Graduate Research Fellowship Program (NSF-GRFP). JKW acknowledges the support of MICCoM, the Midwest Center for Computational Materials, as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, for the development of algorithms and codes used within this work. HS and JKW additionally acknowledge computational resources at the Notre Dame Center for Research Computing (CRC).