EPR spectroscopy and molecular dynamics modelling: a combined approach to study liquid crystals

ABSTRACT

This review outlines the recent theoretical and computational developments for the prediction of motional electron paramagnetic resonance spectra with introduced spin probes from molecular dynamics simulations. The methodology is illustrated with applications to thermotropic and lyotropic liquid crystals at different phases and aggregate states.

Graphical Abstract

KEYWORDS:

• Electron paramagnetic resonance (EPR)
• molecular dynamics (MD) simulations
• MD-EPR simulation methodology
• nitroxide spin probes
• thermotropic liquid crystals
• lyotropic liquid crystals
• discotic liquid crystals
• phase transition
• aggregate states

Introduction

Electron paramagnetic resonance (EPR) spectroscopy with paramagnetic spin probes offers several advantages in the studies of soft matter systems such as liquid crystals (LCs). First, because of high sensitivity of EPR, only a small amount of introduced spin probes is required (typically 1–2 mol %), making it generally a non-invasive technique. Second, EPR signals vary strongly with the orientation of the spin probe relative to the external magnetic field, making it an ideal technique to study the motions and orientations of the host LC mesogens on different temporal and spatial scales. Application of continuous-wave (CW) EPR with spin probes to LCs has been pioneered in the group of Professor Geoffrey Luckhurst from the University of Southampton with the prominent contributions by Professor Claudio Zannoni who was the group member at the time. Their early papers established theoretical foundations for the analysis of the angular variation of EPR line shapes and spin relaxation of spin probes in anisotropic environments with applications to radical, biradical and triradical spin probes in LCs [1–3]. Later at the University of Bologna, Claudio continued to unravel the potential of EPR spectroscopy in LC research alongside other spectroscopic techniques. Claudio’s theoretical interests and mathematical vigour in the use of 3D rotation group theory and the associated algebra of rotational operators for the analysis of EPR made these and his proceeding papers an instrumental and inspiring tool for later developments in the field by other researchers. Over the years, Claudio Zannoni’s group has produced eminent contributions to the study of various LC phenomena by EPR. Examples of applications include the study of poly-liquid crystalline states [4], determination of the helix pitch and handedness of induced cholesteric mesophases, including the ones induced by a chiral probe [5,6], analysis of the effects of dispersed hydrophobic aerosol particles and solutes such as azobenzene on the order and dynamics of nematic LCs [7,8] and determination of director configuration in the twist-bend nematic phase [9] as well as applications to lyotropic LCs [10,11].

Figure 1(a) shows two nitroxide spin probes commonly used in LC studies, namely 3-DOXYL-5-cholestane (CSL) and 5-DOXYL stearic acid (5DSA) alongside two nematic 4-cyano-4-n-pentylbiphenyl (nCB) LC mesogens. In the probes, the unpaired electron is largely localised within the π-orbital of the N-O bond. The states of the system are described by the following spin Hamiltonian (SH) [12]:

(1) $\hat{H}=\beta \widehat{\mathit{S}g}\vec{\mathit{B}}+\widehat{\mathit{S}A}\hat{\mathit{I}}$(1)

Figure 1. (Colour online) (a) Structures of 5CB and 8CB LC molecules and CLS and 5DSA nitroxide spin probes employed in the studies of nematic LCs. Magnetic axes of the nitroxide head group of the probes are shown. (b) Schematic representation of the three EPR hyperfine coupling resonance lines in the absence of motions with the orientations of the nitroxide moiety corresponding to different resonance positions. (c) and (d) Illustrate the sensitivity of X-band EPR line shapes to the motions and orientational constraints of the nitroxide spin probe for isotropic rotation with different correlation times and anisotropic rotations along different magnetic axes, respectively.

where g and A are the magnetic tensors, with g defining the interaction of the spin (S = 1/2) with the external magnetic field (B) and A the hyperfine coupling to the single ¹⁴N nuclear spin I = 1, respectively (β is the Bohr’s magneton). Both tensors are anisotropic, leading to a strong dependence of the EPR resonances on the direction of magnetic field relative to the principle magnetic axes of the probe. Figure 1(b) shows schematically the distribution of the resonance positions for different orientations of the nitroxide plane relative to the external magnetic field. As a result, the EPR line shapes are highly sensitive to both the dynamics and orientation of the probe (see examples with simple rotational motions in Figure 1(c,d)). Note that in EPR measurements, the first derivative of microwave absorption is reordered.

An important feature of EPR is that it is a ‘fast’ spectroscopy able to resolve the details of molecular re-orientational dynamics over time scales between 10⁻¹¹ and 10⁻⁷ s [13]. CW EPR at X-band frequency (9.5 GHz) is particularly sensitive to the nanosecond molecular motions that are characteristic to a wide range of LC systems. The lines are sensitive to both local dynamics and the order of the probe, and their shapes can range from completely motionally averaged three sharp lines in the case of fast isotropic motion to broad asymmetric features in the case of the ‘frozen’ limit (Figure 1(c))). As such, when spin probes are doped inside LCs, the following pieces of information can be extracted from the analysis of EPR line shapes: (i) rotational diffusion of the probe, (ii) motional restrictions imposed on the probe by the surrounding host mesogens and (iii) global order imposed by the LC state.

Although motional EPR spectra of LCs are potentially highly informative, they require careful analysis in order to extract reliably the information about molecular motions and local order as well as director distribution. A typical strategy is the fitting of EPR line shapes using a model of a rigid rod undergoing a stochastic rotation in the asymmetric restricted environment combined with a stochastic Liouville equation (SLE)-based EPR simulation approach pioneered by Freed and co-workers [14–17]. Although being a highly successful approach, the requirement of using multiple adjustable parameters in many cases prevents it from reaching a unique fitting and unambiguous interpretations of the results. With the recent advances in computer power and the increasing accuracy of molecular dynamics (MD) simulations, it becomes certainly desirable to be able to predict EPR spectra directly and completely from MD simulations of actual structures of LC aggregates. This review outlines the most recent theoretical developments and application of the our EPR simulation methodology to different types of LCs.

MD-EPR simulation methodology

Trajectory-based approach using Liouville von Neumann equation

The trajectory-based approach employed for EPR simulations employs the Liouville von Neumann equation in the semi-classical approximation, often called the Langevin form of the SLE [12,18]. In this approach, the SLE is transformed into a system of coupled stochastic differential equations with explicit time dependence of the Liouvillian:

(2) $\frac{d\rho \left(t\right)}{dt}=-i\widehat{\stackrel{ˆ}{L}}\left(t\right)\rho \left(t\right)$(2)

Here ρ is a spin density matrix whose evolution in time is described by Equation (2) and $\widehat{\stackrel{ˆ}{L}}$ is a superoperator of the interaction Hamiltonian (1), expressed in the units of ℏ, defined by its action on the density matrix ρ as: $\widehat{\stackrel{ˆ}{L}}\rho =\left[\hat{H},\rho \right]$. In the case of the nitroxide spin label of the electron spin coupled to a nuclear spin of ¹⁴N, the SLE in the Langevin form is reduced to a system of nine coupled differential equations. The form (2) requires an explicit computational model for the stochastic process to define the time propagation of the spin density matrix. This is achieved by the use of stochastic dynamical trajectories (DT) required for the explicit time variation of $\widehat{\stackrel{ˆ}{L}}\left(t\right).$ A DT can be generated by different means and level of modelling. For instance, using an appropriate Langevin dynamical equation [19,20] applied to a rigid rod undergoing rotational tumbling in an ordering potential, a Brownian dynamics (BD) trajectories can be generated. From direct integration of the Langevin equation, the stochastic trajectory of the orientations of the magnetic axes of the probe defined by three Euler angles $\mathrm{\Omega }=\left(\alpha ,\beta ,\gamma \right)$ is obtained and employed in Equation (2). Alternatively, DTs can be generated from MD simulations of actual molecular structures at either fully atomistic (FA) or coarse-grained (CG) levels. Once propagation of the spin density matrix is completed, EPR line shapes are obtained by applying a Fourier-Laplace transform to the time-dependent averaged transverse magnetisation into the frequency of field domain according to the following equation:

(3) $M_{+}\left(t\right)\propto Tr\left({\hat{S}}_{+}\rho \left(t\right)\right)$(3)

In Equation (3), ${\hat{S}}_{+}=1/\sqrt{2}\left({\hat{S}}_x+i{\hat{S}}_y\right)$ and Tr is the trace of the matrix product. Equation (2) for the propagation of spin density matrix has the following formal solution:

(4) $\rho \left(t\right)=〈exp\left(-i{\displaystyle {\int }_0^t\widehat{\widehat{L}}\left(\tau \right)d\tau }\right)〉\rho \left(0\right)$(4)

where angle brackets indicate the ensemble average and $\rho \left(0\right)$ is the density matrix at the thermal equilibrium. For an appropriately short discrete time step $\mathrm{\Delta }t$, the propagation is achieved numerically and recursively using the following approximation:

(5) $\rho \left(t+\mathrm{\Delta }t\right)\approx \mathrm{e}\mathrm{x}\mathrm{p}\left(-i\widehat{\stackrel{ˆ}{L}}\left(t\right)\mathrm{\Delta }t\right)\rho \left(t\right)$(5)

In Hilbert space, Equation (5) is equivalent to the following unitary transformation of the density matrix occurring during $\mathrm{\Delta }t$ time interval:

(6) $\rho \left(t+\mathrm{\Delta }t\right)\approx \mathrm{e}\mathrm{x}\mathrm{p}\left(-i\hat{H}\left(t\right)\mathrm{\Delta }t\right)\rho \left(t\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(i\hat{H}\left(t\right)\mathrm{\Delta }t\right)$(6)

According to the properties of the Laplace-Fourier transform for an accurate spectral simulation, the following two conditions, known as Nyquist criteria, should be satisfied: $T\ge 1/\delta \omega$ and $\delta t\le 1/\mathrm{\Delta }\omega$, where $\delta \omega$ and $\delta t$ are resolutions in the frequency and time domains, respectively. $\mathrm{\Delta }\omega$ is the width of the entire spectrum in the frequency domain and T is the total propagation (sampling) time. To achieve a reasonable resolution of, say, $\delta B$ ~0.1 Gauss (0.28 MHz between frequency points) after Fourier transformation, a trajectory should be T >3 μs [21]. On the other hand, the second condition of the Fourier transformation imposes the limit on the discrete propagation time step $\mathrm{\Delta }t$ used in Equation (5), namely, $\mathrm{\Delta }t\le \delta t\le 1/\mathrm{\Delta }\omega$. The need for the application of relatively short incremental times and numerical calculation of matrix exponential in Equations (5–6) represents a significant computational challenge for a numerical propagation of the density matrix when the propagation is carried out for the entire sampling time T and when a large number of trajectories are required for statistical averaging of Equation (4). Generally, for relatively small intervals, the matrix exponential can be estimated by either the use of Chebyshev polynomial expansion [22] or decomposing stochastic Liouvillian into a finite sum of matrices followed by the application of a symmetrised version of Trotter’s formula [23]. Such methods, however, are prompt to accumulated numerical errors, particularly for long propagations. An alternative to numerical propagation of the full Liouvillian is to employ simplified forms of the SH (1) where certain terms (mainly non-secular) are neglected allowing analytical solutions. Several reports have utilised approximations to the SH (1), taking into account that the hyperfine coupling of the electron with the nuclear spin is much smaller compared to its interaction with the magnetic field [19,24,25]. The electron spin is, therefore, mainly quantitised along the z axis, whereas the x- and y-contributions of $\hat{S}$ can be neglected (so-called intermediate field [IF] treatment). The Hamiltonian is then given by:

(7) $\hat{H}=\beta {\hat{S}}_z{g}_{zz}B+{\hat{S}}_z{A}_{zz}{\hat{I}}_z+{\hat{S}}_z{A}_{zx}{\hat{I}}_x+{\hat{S}}_z{A}_{zy}{\hat{I}}_y$(7)

For each orientation $\mathrm{\Omega }\left(\mathrm{t}\right)$, the appropriate elements of the g- and A-tensors are determined from:

(8) $\mathbf{g}=R\left(\mathrm{\Omega }\left(\mathrm{t}\right)\right){\mathbf{g}}_{\mathit{d}}{R}^T\left(\mathrm{\Omega }\left(\mathrm{t}\right)\right)\mathbf{A}=R\left(\mathrm{\Omega }\left(\mathrm{t}\right)\right){\mathbf{A}}_d{R}^T\left(\mathrm{\Omega }\left(\mathrm{t}\right)\right)$(8)

where g_d and A_d are the tensors diagonalised in the frame of the nitroxide and $R\left(\mathrm{\Omega }\left(\mathrm{t}\right)\right)$ is the Euler matrix for the rotation from the nitroxide to laboratory frame [26]. According to Steinhoff and Hubbell under this approximation, the Liouvillian in Equation (2) can be effectively presented in a diagonal form, leading to three uncoupled hyperfine coupling lines [19]:

(9) ${\omega }^m\left(t\right)=\left(g_{zz}\left(\mathrm{\Omega }\left(t\right)\right)\beta B+A_{zz}\left(\mathrm{\Omega }\left(t\right)\right)m\right)/\hslash$(9)

As shown in [19] in the slow motional limit, Equation (9) takes into account for the effects of pseudo-secular terms by introducing additional correlation function into the expression for the eigenvalues in ad hoc manner. The transverse magnetisation is then just a sum of three independent scalar contributions:

(10) $⟨M_{+}^m\left(t\right)⟩=⟨exp\left(-\left[i\bar{\omega }t-i{\int }_0^t{\omega }^m\left(\tau \right)d\tau \right]\right)⟩M_{+}^m\left(0\right)⟩$(10)

In Equation (10), the term i\bar \omega t comes from rotation invariant contributions to the SH (7) [21,27]. This approach speeds up the propagation of the density matrix employed in the EPR simulation programmes. Another approach is to use explicit analytical expressions for both the eigenvalues and eigenvectors of the SH (7) [24]. In such a case, analytical propagation is achieved by using:

(11) $exp\left(-i\hat{H}\left(t\right)\mathrm{\Delta }t\right)=\hat{V}\left(t\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(i\hat{\mathrm{\Lambda }}\left(t\right)\mathrm{\Delta }t\right){\hat{V}}^{+}\left(t\right)$(11)

where $\hat{V}$ and $\hat{\mathrm{\Lambda }}$ are the matrix of the eigenvector columns and diagonal matrix of the eigenvalues, respectively, of SH (7) for orientations $\mathrm{\Omega }\left(t\right)$. Although the IF approximation does not include terms ${\hat{S}}_{x/y}{\hat{I}}_i$ that contribute to T₁ relaxation effects, overall the simulation results confirm that in most cases it is an adequate approach for the simulation of CW EPR spectra of nitroxide spin probes at X-band and higher frequencies. In MD-EPR simulation methodology, the EPR spectral line shapes of nitroxide spin probes are determined entirely by the variation with time of the angles that define the orientation of the applied magnetic field to the principal axis of the nitroxide group. Therefore, the orientational history of the magnetic axes in the fixed frame of the MD simulation box is calculated and processed. The z axis of the nitroxide ring (coincident with the direction of the p_z-orbital of N) is calculated from the cross-product of the unit vectors of two N-C bonds of the ring (see Figure 1(a)). The x axis is calculated as a projection vector of the N-O bond on the nitroxide plane (defined by the C-N-C atoms) and the y axis is taken as a cross-product of the z and x vectors [24,27].

Simulation of EPR spectra completely from truncated MD trajectories

In many cases, the rotational dynamics associated with nitroxide spin labels and probes falls within the 0.1–100 ns interval. Also, EPR line shapes of a widely used X-band are not sensitive to a slower dynamics. This timescale, however, is much shorter compared to the sampling time of T >3 μs required to the propagation of the density matrix in EPR spectral simulations. Oganesyan has shown [27,28] that, in fact, the propagation of the density matrix along the entire sampling time T is not strictly necessary, and an accurate simulation can be achieved from a truncated MD trajectory until the point when the autocorrelation function of the re-orientational motion of a spin probe has completely relaxed. The evolution of the spin density matrix for the remaining working length can then be predicted from the truncated trajectory [27,28]. The reported simulation method employs the following: (i) the propagation of the spin density matrix in the Liouville space for this initial time interval ($\tilde{T}\sim 10\tau$, where τ is the effective correlation time for the re-orientational motion of the spin probe) is carried out numerically according to Equations (5–6); (ii) after the autocorrelation function of the spin probe motions has completely relaxed, the rest of the magnetisation trajectory ($t>\tilde{T}$) for the sampling time is calculated using well-defined parameters (see for details [12,27,28]). For example, when the calculation of transverse magnetisation based on Equation (10) is justifiable, calculation of only two parameters for each of the three hyperfine coupling lines, namely ${\bar{\omega }}^m=\frac{1}{\tilde{T}}{\int }_0^{\tilde{T}}⟨{\bar{\omega }}^m\left(\tau \right)⟩d\tau$ and ${\lambda }^m={\int }_0^{\tilde{T}}⟨\mathrm{\Delta }{\omega }^m\left(0\right)\mathrm{\Delta }{\omega }^m\left(\tau \right)⟩d\tau$, is required. Here, $⟨\mathrm{\Delta }{\omega }^m\left(0\right)\mathrm{\Delta }{\omega }^m\left(\tau \right)⟩$ are three autocorrelation functions of $\mathrm{\Delta }{\omega }^m\left(\tau \right)={\omega }^m\left(\tau \right)-{\omega }^m$. The de-phasing of magnetisations is caused by the modulation of $\mathrm{\Delta }{\omega }^m\left(\tau \right)$ due to the re-orientational dynamics of the spin probe. The rest of the evolution of the transverse magnetisation curve can be calculated with high accuracy using the following Equation (12) [27]:

(12) $⟨M_{+}^m\left(t\right)⟩=exp\left(-\left[i\left(\bar{\omega }+{\bar{\omega }}^m\right)+{\lambda }^m\right]\left(t-\tilde{T}\right)\right)⟨M_{+}^m\left(\tilde{T}\right)⟩$(12)

The autocorrelation functions $\mathrm{\Delta }{\omega }^m\left(0\right)\mathrm{\Delta }{\omega }^m\left(\tau \right)$ can be readily used to provide an estimate of the initial propagation time. This methodology is applicable to DTs of any kind (e.g. generated from either BD or MD calculations). The simulation approach does not require numerical propagation of the spin density matrix for the significant part of the sampling time T and thus relies on relatively short single DTs. This offers two major advantages. First, this substantially reduces the amount of computational time for the simulation of EPR spectra. Second, the method relies on the EPR simulation from relatively short DT trajectories over timescales that are within the reach of current MD computational methods [27,28]. This approach has paved the way for the first-time prediction of motional EPR spectra directly and completely from MD trajectories of actual molecular structures. It has been successfully applied to proteins [29,30], DNA [31], lipid bilayers [32] and also LCs [33–35]. A computer code developed by the author [27,28], which includes both options, namely the propagation of the trajectories along the entire sampling time and the use of truncated trajectories combined with autocorrelation function method, has been employed for the simulation of EPR spectra of LC systems reported in this review.

Simulation of EPR spectra indirectly from MD trajectories using SLE in the Fokker-Planck (F-P) form

An alternative approach for the calculation of EPR spectra from the outputs of MD simulations is to solve the SLE in the F-P form using parameters for diffusional tensor and ordering potential extracted from MD outputs [12,36–38]. The application of the SLE in the F-P form has been originally developed by Kubo in the early 1960s [39] and adapted for the simulation EPR spectra by Freed and co-workers using fitting routines with adjustable parameters [14,16,17,40]. SLE in the F-P employs the stochastic generalised diffusional operator $\hat{\mathrm{\Gamma }}$ describing stochastic motion of the probe:

(13) $\frac{\partial \rho }{\partial t}=-i\left[\widehat{\widehat{L}}\left(\Omega \right)-i\widehat{\Gamma }\right]\rho \left(\Omega ,t\right)$(13)

Equation (13) describes the time evolution of the density matrix, $\rho \left(\mathrm{\Omega },t\right)$, in the spin-space basis [14,39–41]. The Quantum Mechanical (QM) part of Equation (13) is represented by the Liouville superoperator $\widehat{\stackrel{ˆ}{L}}\left(\mathrm{\Omega }\right)$ of the system which does not depend explicitly on time in the F-P formulation. The space states of the system are defined by stochastic rotational coordinates, $\mathrm{\Omega }$, that evolve in space under the action of $\hat{\mathrm{\Gamma }}$. The stochastic diffusive operator can be presented in the following form [26,42]:

(14) $\hat{\mathrm{\Gamma }}=-{\hat{\mathbf{J}}}^{+}\mathbf{D}{P}_{eq}\left(\mathrm{\Omega }\right)\hat{\mathbf{J}}{P}_{eq}^{-1}\left(\mathrm{\Omega }\right)$(14)

which is expressed in terms of a diffusion tensor D and $P_{eq}\left(\mathrm{\Omega }\right)$, the orientational equilibrium distribution. $\hat{\mathbf{J}}$ is the angular momentum operator with the components defined in the molecular frame. For instance, in the case of a single rotational diffusion, $\hat{\mathrm{\Gamma }}$ reduces to the following well-known form $\hat{\mathrm{\Gamma }}=D_{xx}{\hat{J}}_x^2+D_{yy}{\hat{J}}_y^2+D_{zz}{\hat{J}}_z^2$, where $D_{ii}$ are the principle values of the rotational diffusion tensor describing this motion and ${\hat{J}}_i^{\hspace{0.17em}}$ are the components of $\hat{\mathbf{J}}$. $P_{eq}\left(\mathrm{\Omega }\right)$in Equation (14) is predetermined by the probe’s ordering potentials [42]. In many cases of LCs, the restriction of the motion can be well approximated by the axially symmetric potential in the Maier–Saupe form, namely:

(15) $V\left(\theta \right)=-\epsilon \left(3\cos {(\theta )}^2-1\right)/2$(15)

with $\theta$ representing the azimuthal angle of the potential. The strength of the potential $\epsilon$ is directly related to the order parameter S of the spin probe [43]. EPR spectrum is then obtained by Fourier-Laplace transform of the SLE F-P equation:

(16) $I\left(\omega -{\omega }_0\right)=1/\pi Re〈〈u\text{|}\left(i(\omega -{\omega }_0\right)+i\widehat{\widehat{L}}\left(\Omega \right)+\widehat{\Gamma }{)}^{-1}\text{|}u{P}_{eq}〉〉$(16)

where $\omega$ and ${\omega }_0$ are the sweep and Larmour frequencies, respectively, $|u⟩$ defines the QM spin transition states in the Liouville space [15,17,42] and angular brackets represent the ensemble average over the classical degrees of freedom which in our case are the Euler angles $\mathrm{\Omega }$. As have been recently shown by Prior and Oganesyan [38], the parameters describing rotational diffusion motion of the spin probe in LCs can be readily extracted from MD trajectories using the Model-Free (MF) approach of Lipari and Szabo [44]. A detailed analysis has been reported to different aggregate states of lyotropic LCs sodium dodecyl sulphate (SDS) and dodecyltrimethylammonium chloride (DTAC) [38]. In addition, it has been demonstrated that the simulation of motional EPR line shapes can be considerably simplified if the fast internal motion of the probe falls within fast motional regime with correlation times below 1 ns. In such cases, the effects of fast local motions can be accounted for prior to solving the SLE by partial averaging of the magnetic parameters A and g in the local frame. Then, the solution of the SLE equation requires only the motional parameters describing the global dynamics of the molecular system.

Calculation of magnetic tensors g and A

In the last decade, accurate Density Functional Theory (DFT) methods for the calculation of the principle values and directions of both A and g tensors have been developed. Among them, DFT calculations with the N07D [45,46] basis set developed by Barone and co-workers, as implemented in the Gaussian 09 software package, are particularly popular choice[47]. Structures are normally first optimised in the gas phase followed by further optimisation in water using implicit continuum solvent schemes. Calculation of magnetic parameters from first principles makes the entire MD-EPR simulation protocol fully predictive.

Experimental measurements

All experimental EPR spectra reported in this review were measured using an X-band Bruker EMX spectrometer equipped with the digital temperature control system (ER4131VT) for high-temperature measurements using a heated flow of nitrogen gas. For each temperature, samples were equilibrated for several minutes before taking the measurement. The following conditions were used: microwave frequency of 9.55 GHz; microwave power of 2 mW; modulation frequency of 100 kHz; and modulation amplitude of 1.0 G.

Application to LC systems

Thermotropic nematic LCs

Our MD-EPR simulation methodology has been applied for the first time to 5CB LC system doped with CSL probe using truncated MD trajectories with autocorrelation function-based method (Equations 5–6 and 12) [33]. A hybrid MD simulation model was used where the spin probe was represented at FA level and the host LCs molecules were treated at a CG level. Simulations of EPR spectra were performed directly and entirely from generated MD trajectories and were found in excellent agreement with the experimental spectra, showing strong sensitivity to both the local order parameter of the probe and director distribution.

The first prediction of EPR spectra for a nematic LC with a doped spin probe from MD simulations at the FA level has been reported in [35]. MD simulations were performed on a system of 5CB consisting of c.a. 140 randomly arranged (head to tail) molecules in a cubic box of 4.5 nm dimensions with cubic periodic boundary conditions.

The EPR predictions and analysis focused particularly on the LC states along the nematic–isotropic (N-I) phase transition region and tested against the measured EPR spectra.

MD simulations were carried out using GROMACS suite [48] with Optimized Potentials for Liquid Simulations - All atom (OPLS-AA) force field. Force field parameters for the spin probe have been developed using a combination of the OPLSS-AA force field and density functional theory calculations at the B3LYP/6-31G** level using Gaussian 03 software [47]. For each temperature, the system, initially composed of 5CB molecules, was equilibrated for 50 ns followed by the doping of the probe at the central position of the box. The system was further equilibrated for another 30 ns followed by production runs of ~100 ns for EPR simulations. EPR spectra were simulated from MD using a direct propagation method combined with autocorrelation function approach. The simulation snapshot of the 5CB with doped CSL is shown in the right top corner of Figure 2. As mentioned above, at X-band, the sensitivity of EPR line shape is predetermined mainly by the anisotropy of hyperfine coupling of ¹⁴N with strong dependence of the resonances on the orientation of magnetic field relative to the principle axes of the hyperfine coupling tensor. In EPR experiment, the magnetic field is set along the orientation of the director of 5CB. Figure 2(a)) shows the comparison between predicted and experimental EPR spectra for a pure nematic (N) phase. It is characterised by three closely positioned narrow lines. In the N phase, the magnetic field on average is lying in the plane of nitroxide ring. For such orientations of the field, the distance between the outer resonance field positions is approaching 2A_xx/yy. In addition, the axial director distribution around the direction of magnetic field in this state was only ~15° consistent with the observed narrow shapes of the lines. Upon increasing the temperature, the system gradually becomes less ordered with the distances between the hyperfine lines increasing and the shapes broadening (see comparison between predicted and experimental spectra in Figure 2(b)). Drastic changes of the shape of EPR spectrum were observed at the phase transition (Figure 2(e)) indicating increasing contributions to the resonance field positions from off-plane orientations of the magnetic field. EPR measurements have shown that this characteristic shape maintains only for a very narrow temperature interval (∆T = 1 K). An associated MD trajectory is shown in the bottom left corner of Figure 2 as a red line. One can see that at the start of the MD run, the system is partially ordered with the order parameter S_LC ~0.6. At approximate time, t = 15 ns, the system clears to the isotropic state with an order parameter, S_LC ~0.15, and at t = 70 ns, a nematic phase is re-grown. Thus, the system is meta-stable and is characterised by a slow dynamical exchange between two distinct phases I and II highlighted in Figure 2(e)). This observation was in full agreement with the bi-modal distribution of instantaneous values of the order parameter P₂ calculated for a range of nCB (n = 5–8) LCs as reported by Zannoni and co-workers [49]. This dynamical exchange process of approximately two events per 100 ns is, however, too slow on the EPR timescale, resulting in the total EPR spectra being a superposition of the line shapes associated with the partially ordered (I) and disordered (II) meta-states. An excellent agreement with the experimental spectrum at the critical point confirms a relatively slow exchange rate between the meta-states. The fluctuations between partially ordered and disordered states can be rationalised based on the use of Isothermal–isobaric (NPT) ensemble employed in MD simulations. NPT represent the state variables of the Gibbs’s free energy whose variation is determined by both the heat exchange ($\mathrm{\Delta }H$) and the entropy change ($\mathrm{\Delta }S$) according to:

(17) $\mathrm{\Delta }G=\mathrm{\Delta }H-T\mathrm{\Delta }S$(17)

Figure 2. (Colour online) Subplots a)–e) show comparison between simulated (solid lines) and measured (dash-dotted lines) EPR spectra at different temperatures. In a)-d) T-T_NI values are indicated in each subplot. Top right corner: Atomistic snapshot of 5CB with the doped CSL spin probe (yellow) positioned in the middle of the simulation box. Bottom left corner: time evolution of the LC order parameter at simulation temperatures T = 380K (blue line) and 385 K (red line). Experimental and simulated spectra are reproduced from [35] with permission from WILEY.

Under the equilibrium condition $\mathrm{\Delta }G=0$, the absorption of the heat from the surrounding LC domains leads to the increase of the entropy, while the loss of heat leads to the decrease of the entropy. Decreases and increases in the entropy are associated with the changes in the order parameter clearly seen in Figure 2. For each LC domain, the frequency of such an exchange would depend on its size and the temperature [50].

At simulation temperature 387 K, both the observed and calculated from MD trajectories EPR spectra have a specific shape which is explained as being a result of a partially restricted local motion of the probe (Figure 2(c)). The local alignment of the probe is strongly averaged out by slow diffusional dynamics of the surrounding LC molecules. Also, at this temperature, the order parameter of the LC is substantially reduced and the director distribution function becomes almost isotropic. When the system reaches completely isotropic state at T ~ 395 K (Figure 2(d)), the fast dynamics of the probe averages out the hyperfine coupling almost completely, and the distance between the lines reaches its maximal values.

FA MD simulations and EPR spectroscopic studies have been performed also on 8CB LC [34]. In this work, AMBER MD simulation suite [51] was used in combination with the General AMBER Force Field (GAFF), and partial charges from CM4 scheme led to the improvement of the predicted value of T_N-I transition temperature as well as observation of both nematic and smectic-A LCs and reproduction of experimental layer spacing in Sm-A [34]. Overall, this study confirmed that the behaviour of 8CB in terms of MD and order in the phase transition region was fully consistent with the one reported earlier for 5CB. Spectra predicted for CSL spin probe doped in 8CB are compared with experimental ones in top panels of Figure 3. As expected, MD simulations reveal the behaviour of CSL in 8CB being similar to the one observed in 5CB across the phase transition region including the dynamical fluctuations at the critical point. This is confirmed by the similarity of the EPR line shapes observed between 5CB and 8CB. It is beneficial to apply structurally different spin probes to the same LC systems with anticipation that their EPR spectra would be different but provide highly complementary information about molecular motions and order. For that reason, 8CB was also doped with 5DSA spin probes and experimental spectra measured at the same temperatures. The predicted EPR spectra are compared to experimental ones in Figure 3 bottom panels and are significantly different from the ones corresponding to CSL probe. In the case of 5DSA probe in the nematic state of 8CB, the averaged orientation of the magnetic z-axis lies along the director (see Figure 1(a)), leading to the shift of resonances towards the outer peak positions of the left and right hyperfine coupling lines (B ~ 3366 Gauss and B ~ 3415 Gauss, respectively). This results in a significant broadening of the spectrum (line a) in the bottom panels of Figure 3). At the same time, the intensity of the inner peaks, which are dominant in the case of CLS probe, is significantly reduced. The slow tilting motions of the z-axis of 5DSA lead to a wide spread of resonances across the hyperfine coupling lines. The resulting EPR line shape is, therefore, a sensitive measure of such motions of the z-axis of 5DSA probe imposed by the alignment of the surrounding host 8CB molecules in the nematic phase. Upon the temperature increase, one would expect the correlation times and the amplitude of this motion to become faster and wider, respectively, causing the EPR spectrum to change significantly. Indeed, at T = 370 K, the line shapes of each of the hyperfine coupling lines become much narrower (line b) in the bottom panels of Figure 3), in agreement with the experiment. The observation of two meta-states I and II at the critical point of the phase transition was much less pronounced in the case of 5DSA (line c)) compared to CSL. Indeed, the line shapes corresponding to the partially ordered (I) and partially disordered (II) states are expected to be close to the ones just above and below line c) (lines b) and d), respectively). Contrary to the case of CLS spin probe, where these two line shapes are drastically different, in the case of 5DSA, the difference between them is less pronounced, leading to a much uniform EPR spectrum observed at T = 375K (lines c)). Heating of the system up to 400 K (lines e)) results in a completely isotropic state of the 8CB system with the fast dynamics of 5DSA averaging out almost completely the hyperfine coupling. As such, the resulting spectrum is very close in the shape to the one observed for CSL probe at the same temperature.

Figure 3. (Colour online) Comparison between experimental and predicted EPR spectra at different values of T-T_NI. a) T-T_NI = −15 K, b) T-T_NI = −5 K, c) T-T_NI = 0 K, d) T-T_NI = 5 K and e) T-T_NI = 25 K for CSL (left panels) and 5DSA (right panels) spin probes in 8CB. Spectra for CSL are reproduced from [34] with permission from the Royal Society of Chemistry. Spectra for 5DSA were calculated from MD simulations of F. Chami and V. S. Oganesyan (unpublished). Contributions from meta-stable states I and II are indicated.

The autocorrelation function of each magnetic axis $\vec{v}\left(t\right)$ is calculated from an MD trajectory according to the following expression:

(18) $C_v\left(t\right)=\underset{T\to \mathrm{\infty }}{lim}\frac{1}{T}{\int }_0^T⟨P_2\left(\vec{v}\left(\tau \right)\cdot \vec{v}\left(t+\tau \right)\right)⟩d\tau$(18)

where $P_{2\left(x\right)}$ is the second-order Legendre polynomial and the bracket denotes the average taken over the orientation distribution, time and the number of available molecules. A ‘sliding time window’ approach is conventionally used for time averaging [27]. The effective correlation times are then calculated by time integration of the associated correlation function in accordance with:

(19) ${\tau }_v=\frac{{\int }_0^{\mathrm{\infty }}[C_v\left(t\right)-C_v\left(\mathrm{\infty }\right)]dt}{\left[C_v\left(0\right)-C_v\left(\mathrm{\infty }\right)\right]}$(19)

With the increasing temperature, both the order parameter and the correlation time of the z-magnetic axis of 5DSA are progressively reduced. In particular, the order parameter $S_z$ is reduced from 0.11 to 0.03, while the correlation time τ_z is reduced from 1.28 to 0.23 ns at T = 360 K and T = 400 K, respectively.

Figure 4 compares the variation with temperature of the calculated order parameters and effective correlation times among 8CB, CSL and 5DSA. Parameters for the main molecular axis of 8CB, y magnetic axis of CSL probe and z magnetic axis of 5DSA probe are shown by black, red and blue circles, respectively. Also, for illustration, the autocorrelation functions of 5DSA spin probe calculated from single MD trajectories at different temperatures are presented in Figure 5. Strong correlation among all three order parameters is observed with S_LC being the largest since it is related to the core vector of the 8CB molecules, while S_y and S_z are associated with the flexible nitroxide head group of the CLS and 5DSA spin probes, respectively. The effective rotational correlation times of 8CB and both probes are distinctly different at the lowest simulation temperature 320 K but become comparable upon approaching the phase transition temperature and above it.

Figure 4. (Colour online) Panels (a) and (b) show the variation with temperature of the calculated order parameters and effective correlation times, respectively, of 8CB, CSL and 5DSA. Parameters for the main axis of 8CB, y magnetic axis of CSL probe and z magnetic axis of 5DSA probe are shown by black, red and blue circles, respectively.

Figure 5. (Colour online) Variation with temperature of autocorrelation functions of 5DSA spin probe in 8CB calculated from single MD trajectories with calculated effective correlation times indicated.

Application to discotic LCs using BD trajectories

Columnar discotic LCs (DLCs) currently attract considerable interest for their potential technological applications as one-dimensional conductors and in sensors, field-effect transistors and photovoltaic solar cells [52,53]. The stacked poly-aromatic cores provide an efficient structure for charge transport along the columns and considered as one-dimensional pathways for an efficient charge and energy transfer. Thus, understanding the dynamics and microscopic phase behaviour of columnar discotics at the molecular scale is crucial for designing devices with desired functionalities. Hexakis(n-hexyloxy)triphenylene (HAT6) is the most commonly known columnar LC based on hexasubstituted triphenylenes. A purpose-designed novel spin probe has been synthesised to study HAT6 [54] (see Figure 6). The design of the probe was based on two important conditions: (i) the probe should resemble the host matrix molecules as closely as possible to favour intercalation and cause minimal disruption and (ii) the nitroxide fragment must be rigid with respect to the discotic core. HAT6 can adopt three phases, namely the isotropic (I), hexagonal columnar (Col) and crystalline (Cr), with temperatures ~373 and ~335 K corresponding to I–Col and Col–Cr phase transitions, respectively.

Figure 6. (Colour online) Top: Structures of HAT6 and discotic rigid core nitroxide spin probe. Bottom: Schematic diagrams representing isotropic, columnar and crystalline domain distributions of HAT6 molecules. Spin probes are shown in red.

Experimental EPR spectra of HAT6, doped with the discotic spin probe between 420 and 360 K, are presented in Figure 7. EPR spectra demonstrate high sensitivity of EPR line shapes to the states of HAT6 with prominent changes observed. The outer right peak position (inset) is particularly indicative of all four states of HAT6, namely I, Mixed, Col and Cr. For temperatures 378 K < T < 420 K of the isotropic phase, EPR spectra have characteristic line shapes corresponding to slow anisotropic rotational diffusion. Upon approaching the critical point, the spectra undergo dramatic change showing contribution from the Col phase in the sample. In the presence of magnetic field, the columns of this LC phase are oriented perpendicular to the field as schematically shown in Figure 6.

Figure 7. (Colour online) Panel (a) shows experimental VT EPR spectra of HAT6 doped with discotic rigid core probe recorded every 3 K between 420 and 360 K on cooling. Inset shows a magnified image of the high-field edge of EPR spectra. Vertical lines indicate the features in EPR spectra corresponding either to I or to Col states of HAT6. ‘Mix’ in the inset indicates the high-field feature associated with the mixed state. Panel (b) shows EPR spectra between 360 and 320 K on cooling. Inset indicates a magnified image of the high-field edge of EPR spectra. Vertical lines indicate the features in EPR spectra corresponding either to Col or to Cr phases of HAT6. Panels (c) and (d) correspond to the simulated EPR spectra for selected temperatures representing most characteristic line shapes. I, Col and Cr shapes are shown by green, red and purple lines. The I–Col phase transition critical point is shown by blue line. Figures are reproduced from [54] with some modification and permission from WILEY.

Spectra measured upon further cooling of the sample through the Col–Cr phase transition are shown in Figure 7(b). The Col–Cr transition is characterised by an abrupt change of the EPR line shape which occurs within a narrow temperature interval (<3 K).

Relatively slow re-orientational dynamics of HAT6 molecules would require generation of long MD trajectories, making the application of all-atom MD simulations at the time too computationally expensive and therefore impractical. Instead, EPR simulations were performed using BD trajectories generated from the solution of Langevin equation for the rotational diffusion of the spin probe in the presence of an ordering potential described by Equation (15). The results presented in bottom panels of Figure 7 demonstrated good agreement with experiment for all phases. The adjusted values of rotational diffusion coefficients $D_{||}^{Col}$ and $D_{\mathrm{\perp }}^{Col}$ for the Col phase were in good agreement with those reported previously [55]. Also, EPR simulations have provided the estimation of the director distribution of the columnar axes from their average plane of alignment which was calculated using a normal distribution with a bandwidth of ~42°. Importantly, the sensitivity of variable temperature EPR spectra to the order and rotational diffusion of the probe molecule in different states has confirmed the columnar domain distribution of HAT6 molecules in the presence of a magnetic field. In the Col phase, the magnetic field on average is lying in the plane of the nitroxide ring (see schematic diagram in Figure 6). As a result, the line shape of a pure Col phase (at 360 K) is characterised by a substantial increase of the resonance intensity corresponding to the orientation of the magnetic field in the xy-plane (inner edges of hyperfine coupling lines at 3377 Gauss and 3392 Gauss, respectively). If the columns would be oriented along the applied magnetic field, the field would be oriented perpendicular to the nitroxide plane, resulting in significantly increased resonance intensity at the outer edge positions of the spectra (3361 Gauss and 3422 Gauss, respectively). This was not observed experimentally. An abrupt change in the line shape upon going from Col to powder Cr is also easily understood from the view point of the changes in the averaged orientational distribution of the spin probe. Cr phase is associated with the even distribution, resulting in the decreased number of resonances in the xy-plane and increased ones in the z orientations (outer edges of the lines) with the noticeable broadening of the total spectrum (Figure 7, right panels). Overall, the Variable Temperature (VT) EPR spectra showed high sensitivity to the HAT6 phase composition, molecular rotational dynamics and columnar order as well as the director distribution. Simulations of the EPR line shapes using a BD simulation model provided accurate estimate of the motional and order parameters at different temperatures along both I–Col and Col–Cr phase transitions.

Application to lyotropic LCs

Surfactant–water systems are highly important in both nature’s (cell membranes) and industrial (e.g. detergents, lubricants, colloidal agents and drug delivery systems) processes and applications and have been the subject of intense study for many years [56,57]. SDS and DTAC are classical examples of anionic and cationic surfactants, respectively.

Recently, Prior and Oganesyan have reported the first simulation of CW EPR spectra of lyotropic LCs at different aggregate states doped with 5DSA spin probe from the results of all-atom MD modelling [38]. Large-scale MD simulations using the GAFF with AMBER suite have been performed on pre-micellar, micellar, rod and lamellar aggregates of SDS and DTAC LCs. The resulting MD trajectories were used to predict and interpret the EPR spectra [38]. A comparative analysis between the results for SDS and DTAC has been undertaken. The purpose of this work was twofold. First, simulation of EPR line shapes from the results of FA MD modelling and comparing them to experiment provided an ultimate test bed for the force fields currently employed to model lyotropic systems. Second, and of equal importance, such a combined MD-EPR methodology allowed the authors to test the validity of the application of the MF approach. MF approach aims to parameterise the time autocorrelation functions of molecular rotations without specific assumptions of the models employed [44]. The approach first introduced by Lipari and Szabo has been used widely in the interpretation of Nuclear Magnetic Resonance (NMR) relaxation data [58–60] and is popular due to the ease of use and clarity of parameters estimated and interpreted. Its formulation is based on the assumption that local and global motions are statistically independent from each other. It was shown in [38] that both the rotational diffusion coefficients and the order parameters associated with the restricted motion of the spin probe can be extracted from the autocorrelation functions and subsequently employed in the SLE-FP EPR simulation procedure using Equations (13–16) in combination with the available computer software [43,61]. In the case when two motional contributions, namely the local restricted dynamics of the probe and the global dynamics of the environment of the LC aggregation, are well separated on the motional timescale, a useful simplification can be employed. Provided that the local motions of the probe are within the so-called motional narrowing limit (τ ≤ 1 ns), the magnetic tensors associated with the nitroxide moiety can be partially averaged prior to the application of the SLE-FP procedure. Then, only the global diffusion of an appropriate symmetry and restriction is taken into consideration by the SLE-FP simulation through the use of an appropriate rotational operator. Importantly, the advantage of the MF-based MD-EPR simulation approach is that its application is relatively straightforward and fast. This is because the simulation parameters estimated by the fitting of autocorrelation curves generated from MD trajectories using the MF approach usually require much shorter trajectory lengths compared to the ones used in the direct propagation methods that require long sampling times [38].

The results of the application of direct propagation technique and indirect simulation using the MF approach were compared to each other and experimental spectra. By performing the MD analysis on the host LC molecules, the relationship between the motions of the host and the probe was also demonstrated [38].

The results showed perfect agreement between predicted and experimental EPR spectra for both micro-aggregate and micellar phases of SDS and DTAC LCs. A good agreement has also been achieved for hexagonal phases of SDS and DTAC and for the lamellar phase of SDS. The disagreement between the simulated and experimental EPR spectra for the DTAC lamellar phase was attributed to the overestimation by GAFF of the packing and order of the carbon acyl chains of the host lyotropics. Contrary to SDS, in DTAC lamellae, the nitroxide head group of the probe penetrates deep inside the environment of the carbon chains, resulting in its increased motional restraint and consequent broadening of the EPR spectrum. This was in agreement with recent modelling studies on non-ionic chromonic LCs that have confirmed that GAFF also overestimates the interaction between the hydrophilic ethyleneoxy chains of the molecules, resulting in their aggregation into compact disorganised clusters instead of ordered stacks [62].

In the current review for illustrative purposes, we will focus on the application of the MD-EPR simulation approaches to the dynamics of SDS micelles at different temperatures. Figure 8 shows the equilibrated structures of 5DSA (purple) doped in SDS micelle composed of 60 surfactant molecules [63]. It also shows snapshots of SDS micelle at different times highlighting the effect of surface diffusion using a selection of molecules shown in red. Snapshots indicate that relatively fast molecular surface diffusion occurs, leading to the loss of coordinated motions among them on the scale of a few nanoseconds. The autocorrelation plots of the magnetic and molecular axes of the probe and host molecules, respectively, are presented in the bottom panels of Figure 8 and were fitted using the following MF expressions:

(20) $\begin{array}{c}{C}_z(t)=e^{-6D_{G}t}\cdot \left(S_L^2+(1-S_L^2)\cdot e^{-t/{\tau }_{\mathrm{\perp }}^L}\right)\\ C_{xy}(t)=e^{-6D_{G}t}\cdot \left(\frac{1}{4}\left(S_L^2+(1-S_L^2)\cdot e^{-t/{\tau }_{\mathrm{\perp }}^L}\right)+\frac{3}{4}\left(e^{-\frac{2}{3}t/{\tau }_{||}^L}\left(S_L^2+(1-S_L^2\right)\cdot e^{-\frac{1}{3}t/{\tau }_{\mathrm{\perp }}^L}\right)\right)\end{array}$(20)

Figure 8. (Colour online) Equilibrated structure of N60 SDS and snapshots of SDS micelles showing the effect of surface diffusion using a selection of molecules represented in red. In all images, water molecules are excluded for clarity; (a) and (b) represent rotational autocorrelation functions of 5DSA spin probe and SDS, respectively. Functions corresponding to magnetic axes x, y and z of the probe are shown by green, blue and red lines, respectively. Functions corresponding to SDS molecular axes X/Y and Z are shown by blue and red lines, respectively. Autocorrelation functions are reproduced from [38] with permission from WILEY.

where the superposition of the axial restrained local motion with two-dimensional isotropic surface diffusion is represented by ${\tau }_{||}^L$ and ${\tau }_{\mathrm{\perp }}^L$ correlation times and, $D_G$, isotropic rotational diffusion coefficient, respectively. $D_G\mathrm{i}\mathrm{s}$ related to the associated correlation time according to $D_G=1/\left(6{\tau }_G\right)$. As one can observe in each case, the decay curves are consistent with the overall axial dynamics which has distinct bi-exponential character, representing the superposition of fast restricted internal and unrestricted global motions. For the SDS micelle, the local order parameter and global diffusion coefficient ($S_L$= 0.36 and $D_G$ = 3.8 × 10⁷ s⁻¹ at 300 K) were found in good agreement with that of Westlund and co-workers who fitted the EPR spectra with the two-dynamic model ($S_L$= 0.43 and $D_G$ = 5.0 × 10⁷ s⁻¹ at 298 K) [64].

The magnetic tensors A and g were partially averaged out in response to the local axial order, represented by the order parameter$S_L$, and in accordance with the equations derived by Gaffney and McConnell for their averaged values in the axially symmetric local environment [65]. The EPR spectra predicted from MD data using the MF approach are shown in Figure 9 in red and are in perfect agreement with the experimental ones shown in blue. They were also in full agreement with the results of EPR simulations using direct propagation technique for the entire sampling time (300 ns), although evidently much shorter trajectory lengths were required for the simulation using the MF approach (~30 ns).

Figure 9. (Colour online) Comparison between experimental (blue) and predicted (red) EPR spectra of SDS micelle doped with 5DSA at three selected temperatures. Spectra are reproduced from [38] with permission from WILEY.

The global tumbling of the micelles is on the timescale of nanoseconds. As mentioned above, EPR line shapes at X-band are particularly sensitive to the motions occurring on the 0.1–100 ns timescale, making the application of our MD-EPR simulation methodology highly feasible for testing water force field models. Thus, additional simulation of SDS micelle using the same force fields for surfactants but combined with TIP3P water model, instead of SPC/E, was undertaken. The result for T = 320 K is presented in Figure 10 and is compared to both the prediction using SPC/E water model and experiment. The predicted from MD data EPR spectrum using TIP3P has small but noticeable narrowing of the high field line in the EPR spectrum compared to both the experimental spectrum and the one simulated using SPC/E. The observed differences are attributed to the fact that contrary to SPC/E, which is known for an adequate representation of water diffusion [66], TIP3P tends to underestimate the viscosity of water molecules [66] and consequently overestimate the rotational diffusion of micellar aggregates, resulting in the narrowing of the hyperfine coupling line in the EPR spectrum. Indeed, the extracted global diffusion coefficient in the cases of SPC/E and TIP3P were $D_G$ = 6.7 × 10⁷ s⁻¹ and $D_G$ = 9.9 × 10⁷ s⁻¹, respectively. This behaviour has been confirmed in most recent MD-EPR simulation study of double- and single-strand spin-labelled DNA where more pronounced difference in EPR spectra was observed between SPC/E and TIP3P models because of relatively slow rotational diffusion of DNA fragments in water [31].

Figure 10. (Colour online) Comparison between experimental EPR spectrum of SDS micelle doped with 5DSA at 320 K (blue lines) and those predicted from MD (red lines) using two different force field models for water.

Conclusions

A novel MD-EPR simulation methodology bridges the gap between an advanced spectroscopic technique and MD modelling, increasing substantially the informational content generated from the analysis of EPR measurements on LCs. The main advantage of our approach is that motional contributions are explicitly accounted in MD trajectories of actual structures without the need for evoking highly parametrised models with multiple adjustable parameters. Such an approach not only greatly simplifies the interpretation and analysis of experimental results but also provides a rigorous test bed for the state-of-the-art molecular force fields. Recent applications to both thermotropic and lyotropic LCs show that this strategy has a power of providing unambiguous conclusions about molecular structure and dynamics in different states and across phase transition regions of LCs. In particular, an accurate estimate of the molecular rotational correlation times, local order and global director distribution can be achieved. With the fast improvements in the available computing power and parallelisation methods (including novel cluster architectures with Graphics processing units/Message Passing Interface (GPU/MPI)), the major focus in the near future is likely to be on achieving longer trajectory lengths over realistic times and the possibility of all-atom MD-EPR studies on larger and more complex LC systems. This will entail the possibility of applying direct propagation methods for the entire sampling time as well as using larger bulk system for introduction of multiple spin probes for improved statistical calculation and the increased accuracy of the predicted spectra as a result. Further advances are expected in the development and use of hybrid MD models where the spin probe is modelled at all-atom level, but the bulk is simulated at CG level and thus allowing significantly longer trajectories over reasonable time.

Acknowledgements

I would like to acknowledge the financial support of the LC work in my lab by EPSRC (Grants EP/H020411/1, EP/L001322/1 and EP/P007554/1). I would like to thank all the past and present members of my group involved one way or another in the studies of LCs and, in particular, Dr Fatima Chami and Dr Chris Prior for their contributions and for carrying out additional MD simulations. Finally, I would like to express special gratitude to my long-standing collaborators Prof Mark Wilson from Durham University and Prof Andrew Cammidge from the University of East Anglia and for their contributions to the work reported in this review.

Disclosure statement

No potential conflict of interest was reported by the author.

Additional information

Funding

This work was supported by the Engineering and Physical Sciences Research Council [EP/H020411/1, EP/L001322/1, EP/P007554/1].