Communication: Molecular dynamics and ¹H NMR of n-hexane in liquid crystals

Conformational statistics of flexible alkanes partially oriented in liquid crystals (LCs) has long been a challenge to both theorists and NMR spectroscopists. It is fundamental to a host of physical phenomena ranging from the roles of flexible end chain(s) in the LC phase structures to functionality in biological systems. The ability to determine dipolar couplings from ¹H NMR of orientationally ordered chains is a powerful tool in these applications. However, the spectral complexity increases dramatically with the number of coupled spins, and until recently analysis was limited to rather small molecules. In recent work, we have demonstrated the use of evolutionary algorithms (EAs) for the analysis of complicated NMR spectra of solutes in LC phases.^1–5 In particular, we have used the covariance-matrix adaptation evolution strategy (CMA-ES)⁶ to fit spectra of n-butane⁵ and n-pentane.³ Analysis of the dipolar couplings obtained from these fits has yielded valuable information about the conformational averaging and the effect of the anisotropic environment on the conformer populations. We wish to extend our investigations to more complicated systems, in particular to n-hexane. However, the dramatic increase in complexity of the NMR spectrum poses problems. For the EA to succeed with such a complicated system, we need a reasonable set of starting dipolar coupling ranges. In the case of n-pentane, we successfully used models for orientational order and Molecular Dynamics (MD) simulations to give estimates of these ranges.⁷ However, in the case of n-pentane (of order 13 000 transitions with intensities greater than 0.3% of the most intense line) and n-butane (of order 2000 transitions >0.3%), the peaks in the NMR spectrum are predominantly separated, while in the more complex n-hexane case (of order 150 000 transitions >0.3%), there is a broad background containing many superimposed or overlapping peaks, as seen in the simulated spectra of Fig. 1 (left). In this paper, we discuss strategies to attack this challenging problem, including the use of models and MD simulations, in an attempt to solve the ¹H NMR spectrum of orientationally ordered n-hexane. Obtaining an accurate set of dipolar couplings will be key for a thorough study of the n-hexane conformational problem. Should the n-hexane spectrum in an oriented nematic (N) phase become solvable to yield precise dipolar couplings, the same strategy could open a new potential area of application for NMR spectroscopy, such as exploring the proton dipolar couplings of the alkane chain spacer in LC dimers.^8,9

In this paper, we report both the experimental proton NMR spectrum obtained for n-hexane as solute in the nematic LC Merck ZLI-1132 (1132) and the result of atomistic MD simulations. In the MD simulations, not only the interactions between molecules but also the solute and solvent internal motions such as vibrations and conformational changes are taken into account by appropriate molecular force fields (FFs). The MD simulation of solutes in a liquid-crystal solvent can be used to predict experimental observables, such as nuclear dipolar couplings, as has been demonstrated for a collection of small solutes¹⁰ and of biphenyl in 5CB.¹¹ Reasonable estimates of dipolar couplings have also been obtained for more complex and more flexible molecules, such as n-pentane in 5CB.⁷ Although we would have preferred to use 5CB as a solvent, the spectra obtained were not of the highest quality. Instead, therefore, we used 1132, a commercial mixture of p-butylcyclohexyl-p′-cyanophenyl (38.4%), p-pentylcyclohexyl-p′-cyanophenyl (34.2%), p-heptylcyclohexyl-p′-cyanophenyl (18.1%), and p-pentylcyclohexyl-p′-cyanobiphenyl (9.3%). It is not feasible at the moment to simulate solutes in such a multi-component mixture, particularly in view of the relatively small total number of molecules (of order a few thousand) that can be treated. On the other hand, the various components of the mixture have some similarity (a rigid core with terminal alkyl chain and cyano group) with p-n-alkyl-p′-cyanobiphenyls (nCB), including 5CB in particular. Thus, we assume that 5CB is a reasonable choice for a model solvent.

In this study, we use three different FFs for n-hexane as a solute in 5CB in order to explore the range of variability of the predicted dipolar couplings as a function of small variations of the solute FF parameters. Scheme 1 shows the molecule n-hexane with the proton labels.

Experimental ¹H-NMR spectra of n-hexane have been obtained in 1132, see Fig. 1. The spectral width is of order 20 kHz, with the lines in the wings being broader than the lines in the centre which have a linewidth of ca. 2 Hz. Although great pains were taken to equilibrate the sample, the broader lines are an experimental artifact due to remaining inhomogeneities in temperature and/or concentration across the sample tube.

In order for the CMA-ES to locate the global minimum (i.e., not get stuck in a false minimum), the following two criteria apply. First, we must choose search ranges for the NMR parameters (chemical shifts and dipolar couplings) that are not too large; ranges suggested in the earlier n-pentane work⁷ are the larger of ±50 Hz or ±15% of the trial value of a dipolar coupling. Tests with a simulated n-hexane spectrum confirm these limits and indeed mostly fail with ±100 Hz limits. Second, as also confirmed with the simulated spectrum, the parameter set that leads to the global minimum must be within the parameter ranges chosen. In order to make sensible guesses at sets of possible starting parameters, we utilize both MD simulations as well as predictions from what might be considered the best model specifically designed for estimating order parameters of hydrocarbon solutes, the so-called chord model.^12,13 This model describes alkanes as connected CH₂ building blocks. In its simple form, these blocks have axial symmetry and there is one adjustable parameter $w ̃ 0$⁠; the assumption of biaxiality leads to a second independent parameter $w ̃ 1$⁠. Here, we use the two-parameter version (named modified chord, or mod-chord for short) with the definitions of $w ̃ 0$ and $w ̃ 1$ as given in Refs. 3 and 12.

The idea is that at least one of the MD or mod-chord calculations might produce a “good” starting set of dipolar couplings for the ensuing fitting procedure with CMA-ES and that the targeted set of dipolar couplings searches for convergence between the calculated and experimental spectra without getting stuck in a false minimum (i.e., a failed search by CMA-ES). As already mentioned, computer simulations of n-hexane were carried out instead in the LC 5CB. The use of different LCs for experiment and simulation is at least partially justified because for the “magic solute” n-hexane, the anisotropic intermolecular interactions that are responsible for its orientational order are dominated by short-range size-and-shape effects, and longer-range interactions play only a minor role.^14–16 This equivalence has recently been demonstrated for a temperature-dependent study of n-pentane in both 5CB and 1132.¹⁷ 

The MD simulated sample consists of 1888 molecules of 5CB, 32 molecules of 1,3,5-trichlorobenzene as an orientational standard and 80 molecules of n-hexane. The 5CB solvent molecule is described at the United Atoms (UA) level of detail, using a FF that reproduces the experimental nematic-isotropic transition for n CB;¹⁸ while for 1,3,5-trichlorobenzene, a previously derived full atom parameterization is employed.¹⁰ The n-hexane conformations and their interactions with the solvent are simulated with three different potential energy functions: AMBER, CHARMM, and AMBER-QC. The first two are standard and of widespread use, namely, the AMBER-OPLS force field^19,20 and the CHARMM27 FF for alkanes.²¹ The third is named AMBER-QC because it still relies on AMBER-OPLS Lennard–Jones parameters, but with the atomic charges and the three torsions about hexane C–C bonds being parametrized with high level CCSD(T)/6-311G++(3df,3pd) calculations (performed with Gaussian09²² following the procedure described in Ref. 11). In principle, the AMBER-QC FF should provide a more accurate description of n-hexane conformations than do AMBER and CHARMM, where the torsional parameters and the atomic charges are not explicitly tailored for n-hexane. We ran simulations at constant number of particles N, pressure P and temperature T (NPT-MD) using the well-tested computer code NAMD²³ at the following thermodynamic conditions: P = 1 atm and T = 285 K. The typical MD simulation time was about 45 ns. We have analyzed the simulation results in various ways, and in particular, we have calculated atom-atom radial correlation functions between the various atoms of n-hexane and those of 5CB. These correlations allow us to investigate the preferred solute-solvent organization at short range. The radial distributions show a tendency for n-hexane to locate at the two ends of 5CB, these being local environments that are also present in the various components of 1132.

From simulation, we calculate the dipolar couplings between protons i and j using the relation

where r_ij is the internuclear distance, θ_ij describes the orientation of the vector connecting nuclei i and j with respect to the magnetic field B (in our case, coincident with the phase director), μ₀ = 4π10⁻⁷ T² J⁻¹ m³ is the magnetic permeability in vacuum, and γ_i = 2πg_iμ_N/h is the nuclear gyromagnetic ratio expressed in rads.

Hence, we use six possible starting sets of D_ij values (3 from MD with different FFs and 3 from mod-chord calculations with parameter $w ̃ 0 / w ̃ 1$ ratios of 0.753,³ 0.897,¹⁷ and 1.077²⁴) for CMA-ES to fit the experimental n-hexane spectrum (Fig. 1). Table I summarizes D_ij values obtained from the three different MD simulations and one mod-chord calculation. These sets of dipolar couplings were scaled to give simulated spectra of the same spectral width as the experimental spectrum (see Fig. 1, left). The different sets of couplings are quite different from each other with some smaller couplings having different signs. As mentioned above, because of the broader linewidth in the wings, the perceived experimental spectral width may be underestimated. Also listed in Table I are the D_ij values finally obtained from a successful CMA-ES fit (see below). The attempts to use MD simulations or mod-chord calculations were unsuccessful as (regardless of the scaling factor used) either one or more input dipolar coupling(s) was outside the search range unless it was set too wide.

At this point, we resort to earlier experimental data from the Pines group^25,26 of the D_ij for n-hexane in the N phase of p-pentyl-phenyl-2-chloro-4-(p-pentylbenzoyloxy)-benzoate (EK11650) as obtained by proton-proton double quantum 2D-NMR methods on partly deuterated n-hexane. The errors for the D_ij values obtained are of order 25 Hz, much greater than those expected from the analysis of a high-resolution spectrum. As pointed out above, n-hexane is a magic solute,^14–16 and these D_ij (in a different LC solvent) are a possible starting point for the spectral fitting. The first attempt with these couplings used search ranges of ±50 Hz or ±10% of the trial value of a dipolar couplings and a scaling factor of 0.60 to produce a spectrum that is consistent in width with the FF and mod-chord ones, but did not produce a fit. However, when we scaled these trial dipolar couplings up by 10% (original times 0.66, yielding a wider initial spectral width) and used the same search ranges, the EA was successful in yielding a final and precise set of n-hexane dipolar couplings. The observation of this fit was truly an “aha moment.”²⁷ Moreover, the present fit shows convincingly that even a huge amount of spectral overlap can be handled adequately by our methods. We note that scaling up of dipolar couplings worked only for the Pines couplings: in all other cases, the couplings were too far removed from the actual ones and the CMA-ES converged to a false minimum.

Table I lists the experimentally derived D_ij of n-hexane found by Pines et al.²⁶ in nematic EK11650, the original values scaled up by 0.66, and the precise values obtained from the CMA-ES fit (as well as the indirect spin-spin J_ij couplings and relative chemical shifts). In the CMA-ES, the initial values of J_ij were fixed to values for n-hexane at 300 K in chloroform;²⁸ in later calculations where the J_ij were fitted, their ranges were set to ±1 or 2 Hz for longer-range couplings up to about ±6 Hz from the literature value for two-bond couplings. Examination of the table shows that the Pines experimental dipolar couplings (scaled by 0.66) are much closer to the actual values and better than those predicted from MD or model calculations.

Now it is clear why the other predicted sets of D_ij (see Table I) were not suitable, as the criteria for CMA-ES to work properly were not met. Figure 1 (right) compares the background-removed simulated spectrum resulting from the CMA-ES fit with the background-removed experimental spectrum. Examination of the blow-up region at the bottom shows that the agreement between experimental and simulated spectra is indeed amazingly good. This agreement applies to all regions of the spectrum.

In conclusion, the precise set of D_ij values of n-hexane should prompt MD simulators to further improve the force field (note AMBER appears to predict the closest set of D_ij values) suitable for describing flexible alkanes. The work can now be extended to a full temperature study of n-hexane in liquid-crystal solvents which in principle will reveal information about its 27 conformations, 10 of which are symmetry-unrelated. It further points to other possibilities, such as unraveling the proton dipolar couplings in a partially deuterated LC dimer such as CBC7CB (CB_D8CD₂(CH₂)₅CD₂CB_D8) in the twist-bend nematic (N_tb) phase under deuterium decoupling,⁹ thus enhancing the possibility of investigating the important role that conformational freedom may play in that case.