Optimization of the GAFF force field to describe liquid crystal molecules: the path to a dramatic improvement in transition temperature predictions

Abstract

The physical properties and phase transitions of thermotropic liquid crystals are highly sensitive to small changes in chemical structure. However, these changes are challenging to model, as both the phase diagram and mesophase properties obtained from fully atomistic simulations are strongly dependent on the force field model employed, and the current generation of chemical force fields has not proved accurate enough to provide reliable predictions of transition temperatures for many liquid crystals. This paper presents a strategy for improving the nematic clearing point, T_NI, in atomistic simulations, by systematic optimization of the General Amber Force Field (GAFF) for key mesogenic fragments. We show that with careful optimization of the parameters describing a series of liquid crystal fragment molecules, it is possible to transfer these parameters to larger liquid crystal molecules and make accurate predictions for nematic mesophase formation. This new force field, GAFF-LCFF, is used to predict the nematic–isotropic clearing point to within 5 °C for the nematogen 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester, an improvement of 60 °C over the standard GAFF force field.

---

1 Introduction

Over the last few years, computer simulations techniques have proved immensely useful in the study of thermotropic liquid crystals.^1,2 Simulation studies have provided a series of insights into the structure and dynamics of liquid crystal phases over a wide range of time and length scales. These include director modelling, and coarse-grained and atomistic molecular models.^1,3 Although coarse-grained models are suitable for the study of large system sizes at reasonably low computational cost, fully atomistic simulations have the potential to link chemical details with the physical properties of a system.⁴ For example, the phase transition temperatures and the stability of a range of mesophases are particularly sensitive to small changes in chemical structure. Unfortunately, the phase diagram and mesophase properties derived from atomistic simulations are strongly dependent on the force field employed and its description of the molecular geometry and intermolecular interactions.⁵ Over many years, this has provided a major barrier to making accurate predictions of liquid crystal phase behaviour^6–9 and structure–property relationships.^10,11

In standard force fields, such as OPLS,¹² GROMOS,¹³ or AMBER,¹⁴ the parameters are derived from experimental and/or quantum mechanical (QM) values for a selected set of molecules. For the description of larger molecules, including liquid crystals, parameters describing atoms in small molecules are transferred to the larger entities, with the assumption that the chemical environment is similar.¹⁵ However, the use of standard force fields for the study of liquid crystal systems often produces only approximate results, indicating that the transferability hypothesis cannot always account for specific molecular details that may affect the stability of the mesophases and the macroscopic properties. For example, in initial test simulations, Tiberio et al. obtained a simulated nematic–isotropic transition temperature, T_NI ≈ 120 K above the experimental value using a combination of the AMBER and NERD (united atom alkane) force fields (prior to optimizing these potentials to obtain good transition temperatures) for the mesogen 5-alkyl-cyanobiphenyl (5CB) and homologues;^16,17 and recently Chami et al. simulated 8CB employing the all-atom General Amber Force Field (GAFF), producing a T_NI ≈ 61 K higher then the experimental value and a T_SN ≈ 33–53 K higher than the experimental transition temperature.¹⁸ Similarly, Kaprusevicius et al. obtained a T_NI ≈ 75 K higher than experimental for the same molecule with a force field based on OPLS-AA.¹⁹ In these cases, the high simulated T_NI values suggest that many modern force fields significantly over estimate the attraction between mesogenic molecules.

The objective for the current study is to optimize a general all-atom force field to make it suitable for the simulation of a range of liquid crystal molecules. We consider a series of fragment molecules (which when combined cover thousands of calamitic liquid crystal molecular structures in the literature) and produce refined parameters using GAFF as a starting point for the optimization work. Section 2 of the paper, provides a summary of the background to the development and testing of GAFF, Section 3 gives details of the computational techniques used. The results of the amended force field parametrization, GAFF-LCFF, are included in Sections 4.1 and 4.2, along with a discussion of these results. Section 4.3 presents the results from testing GAFF-LCFF parameters for a typical calamitic nematogen, 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester. Conclusions are drawn in Section 5.

2 Development and evaluation of a liquid crystal force field: initial choice of GAFF

One approach to improving the description of liquid crystal materials is to develop an original force field from QM calculations, as described by Cacelli et al.,^4,20–22 using the Fragmentation Reconstruction Method (FRM). FRM results for a number of liquid crystal systems have led to good agreement between experimental data and a number of calculated structural and thermodynamic properties. However, most standard force fields are based on effective two-body interactions, with three-body (and higher body) effects included in an average way. So even if the QM-derived force fields are sufficiently high quality to fully capture dispersion interactions, they are likely to suffer from a neglect of higher body interactions. While this is known to be a major problem for small molecules such as water,^23,24 the success of the Cacelli et al. approach suggests this may be less of a problem for larger molecules.²⁰

An alternative approach to obtaining a force field suitable for liquid crystal molecules is to focus on the refinement of standard force fields, via the amendment of some of the key parameters. This approach was demonstrated for a united atom force field by Tiberio et al.,¹⁶ who were able to tune the LJ parameters of the AMBER force field for the family of n-alkyl-cyanobiphenyls (n-CBs) so that the densities, phase transition temperatures, orientational order parameters and NMR residual dipolar couplings could be reproduced with good accuracy. In particular, the T_NI for 5CB was reproduced within ±4 K of the experimental value. Another study by Zhang et al., employing the TraPPE-UA force field with improved torsional potentials, involved reoptimizing a number of aromatic carbon LJ parameters of the biphenyl unit, with the aim of reproducing the experimental density of 5CB within 2%.²⁵ Results showed good agreement with experimental densities and the T_NI transition temperatures.

Among the standard force fields available for atomistic simulations, the General AMBER Force Field (GAFF) provides an attractive starting point for liquid crystals. GAFF is a force field developed with the objective of describing a wider range of molecules than those covered by the existing “AMBER” force fields, which were primarily developed for protein and nucleic acid systems.²⁶ Moreover, the Lennard-Jones 12:6 form for nonbonded parameters in GAFF is computationally cheap to simulate, and the force field itself can be easily employed within commonly used molecular dynamics programs such as AMBER,²⁷ GROMACS²⁸ and DL_POLY.²⁹

Although GAFF was designed as a general purpose force field with wide applicability,¹⁴ a number of recent attempts to reproduce accurate clearing temperatures for mesogens using GAFF have been unsuccessful.^18,19 We have shown also that simulations of a bent-core oxadiazole (ODBP) based mesogen and a linear phenyl ester mesogen (1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester as discussed below) result in T_NI values ≈110 K and ≈60 K higher than the experimental values respectively. Another study produced a T_NI temperature ≈200 K higher than experimental for a T-shaped benzothiazole mesogen.³⁰

This investigation attempts to improve GAFF for mesogenic molecules, via optimization and fitting of some key torsional angles (included a re-examination of existing conformational data for n-alkanes), refinement of Lennard-Jones parameters for mesogenic fragment molecules to improve the reproduction of experimental densities and heats of vaporization (Δ_vapH), and addition of extra parameters where simulations of fragment molecules show that GAFF parameters are not transferable.

We note in passing some recent systematic studies with GAFF, which have pointed to its potential if further optimized. Wang and Hou have tested the ability of GAFF in reproducing the bulk densities and Δ_vapH for 71 organic molecules. For densities, the average percent error is 4.43% but this can be improved greatly by focussing on molecules where GAFF performs poorly.²⁶ In another study Caleman et al.^31,32 have tested the ability of OPLS/AA and GAFF in reproducing some key properties of 146 small molecules. Their results showed that in general GAFF performs well but shows an overall slight underestimation of densities for the compounds studied, along with an over estimation of Δ_vapH for the majority of the compounds.

Finally, we note that previous attempts have been made to improve/optimize standard force fields for use with liquid crystals, notably with the OPLS-AA force field.³³ Here, it was thought that optimizing predictions of densities to within 2–3% would be sufficient for liquid crystals fragment molecules. In this study, we show that a criterion of better than 1% is ideally required for reasonable clearing temperature predictions.

3 Computational

3.1 Quantum chemical calculations

Calculations were performed using density functional theory (DFT), employing the B3LYP functional, or at the MP2 level. All calculations were carried out with the Gaussian09 suite of programs.³⁴ Initial structures were optimized at the B3LYP/6-31g(d,p), prior to higher level or larger basis set calculations. Dihedral scans were typically carried out at 6° intervals, and zero-point vibrational energies (ZPVE) have been neglected in all the calculations as these have been shown to contribute less than 1.0 kJ mol⁻¹ to the torsional potential profiles of conjugated systems.^35,36 The procedure used for the parametrization of the dihedral angles consisted of minimizing the squared difference, χ², between the molecular mechanics (MM) and the QM energieswhere E^QM(φ^j_i) and E^MM(φ^j_i) are the quantum mechanical and MM energies measured relative to the lowest energy conformation and N_pts represents the number of QM points for calculating the rotational profile of dihedral angle (φ^j_i).

The E^MM(φ^j_i) can be expressed as

E^MM(φ^j_i) = E_torsion(φ^j) + E^ff, (2)

where E^ff represents the other force field terms that contribute to the dihedral angle of interest in addition to the torsional terms, which are being fitted.

3.2 Atomistic simulations of fragment molecules

All calculations were performed using the GROMACS 4.5.5 package³⁷ using GAFF as the starting force field. The energy function employed in the MD simulations is given bywhere r_eq, θ_eq are respectively natural bond lengths and angles, K_r, K_θ and C_n are respectively bond, angle and torsional force constants, σ_ij and ε_ij are the usual Lennard-Jones parameters and q_i, q_j are partial electronic charges. Changes in E^MM arising from deviations in improper dihedral angles, ω, are represented by cosine functions using the force constants, k_d, the harmonic coefficients, n_d, and the phase angles ω_d. Throughout this work these improper dihedral angle parameters have been taken unchanged from the original GAFF force field. The standard Lorentz–Berthelot mixing rules, ε_ij = (ε_iε_j)^1/2 and σ_i = (σ_i + σ_j)/2, have been applied throughout this work. The Antechamber software from AmberTools 1.4 was used to generate GAFF topologies, with the point charges derived through the AM1-BCC method. The GAFF topologies and coordinate files were converted into the GROMACS format using the ACPYPE script.³⁸

Simulations of fragment molecules in the liquid phase were carried out using 1000 molecules at either 298 K or 293 K and a pressure of 1 bar (unless otherwise stated). A cutoff of 1.2 nm was used for short range nonbonded interactions, the Particle Mesh Ewald (PME) method was used for long-range electrostatics, and the simulations employed the usual corrections for the pressure and potential energies to compensate for the truncation of the vdW interactions. The Berendsen thermostat and barostat was used for initial simulation setups compressing, at 100 bar pressure, from low density random arrangements of molecules, followed by equilibration and production runs with a Nose–Hoover thermostat and Parrinello–Rahman barostat, once liquid state densities were reached. Bond lengths were kept fixed at their equilibrium values using the LINCS algorithm and a time step of 2 fs was employed. Typically, each production run was carried out for 20 ns.

The heat of vaporization was calculated using

Δ_vapH = (E_pot(g) + k_BT) − E_pot(l) (4)

where E_pot(g) represents the intramolecular energy in the (ideal) gas phase and E_pot(l) is the intermolecular energy in the liquid phase. The gas phase simulations, for 1 molecule, were performed using a stochastic dynamics (SD) integrator, which adds a friction and a noise term to Newton's equation of motion. Gas phase calculations were carried out over 200 ns.

We note that considerable work has already gone into the parametrization of GAFF, as documented in previous studies.^26,31 Hence, we concentrate here on key mesogenic fragments that are less than optimally represented within this force field. Results below are presented for the torsional potentials of phenylbenzoate, 2,5-diphenyl,1,3,4-oxadiazole and n-alkanes, and optimizations of Lennard-Jones parameters.

3.3 Atomistic simulations of mesogens

Simulations were also carried out for the phenylester-based mesogen, 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester. These simulations were carried out on 256 molecules, starting from an initial low density gas phase, compressing to a liquid at high temperature (using the procedure described above) and then cooling to a nematic phase. These simulations required long equilibration times and long production runs of ≈110 ns to see the formation of a nematic phase. The lengthy simulation times, coupled with cooling the system from disordered configurations, provides greater confidence in the results when observing the spontaneous onset of ordering in LC phases.^16,17,25

Molecular orientational order for the mesogens was monitored through the calculation of the orientational order parameter, S₂. This was calculated from the order of the long axis of the core, defined by a vector along the central phenyl ring, u⃑_i(t). The instantaneous average across molecular vectors in the simulations defines a director n⃑. n⃑ and S₂ are obtained by calculating the ordering tensor

where the sum runs over all N molecules. The largest eigenvalue of the Q tensor representswhere P₂ is the 2nd Legendre polynomial, and the associated eigenvector is the director n⃑(t). In practice, to minimize system size effects in locating the phase transition, we use −2× the middle eigenvalue of Q, which fluctuates about a value of zero in the isotropic phase but equals P₂(t) in the nematic phase. We define S₂ as the time average of P₂(t), over the final 80 ns of the production run simulations.

Structural information in the nematic and isotropic liquid phases was deduced by evaluating a set of pair distribution functions: g(r), calculates the most probable intermolecular distances between two particles irrespective of orientation

where r_ij is the vector between the centres of mass of two molecules i and j. The pairwise orientational distribution function, g₂(r), measures the average relative orientation of molecules separated by a distance r and has the form

g₂(r) = 〈P₂(cos θ_ij(r_ij))〉_ij, (8)

where θ_ij is the angle between vectors u⃑_i and u⃑_j for molecules i and j at separation r_ij. g₂(r) → 0 in the isotropic liquid and g₂(r) → S₂² in a nematic phase. Higher values of g₂(r) at short distances reflect strong local orientational ordering.^39,40 The anisotropic nature of liquid crystal phases requires additional distribution functions to differentiate between nematic and smectic phases. g_∥(r), and g_⊥(r) respectively measure the component of the g(r) parallel and perpendicular to the director.

4 Results and discussion

4.1 Optimization of torsional potentials

Phenylbenzoate. Many liquid crystal forming molecules contain a phenyl benzoate (PB) fragment as part of their structure. The position and conformation of the ester group is implicated in the development of spontaneous polarization in ferroelectric liquid crystal phases and spontaneous chiral segregation of bent-core liquid crystals as well as affecting the magnitude of the bend angle.^41–44 This behaviour is controlled by the two dihedrals, φ₁ and φ₂ (Fig. 1) associated with internal rotation around the C(=O)–C(phenyl) and C(phenyl)–O(ester) bonds. The torsion around the central C(=O)–O bond is generally considered to be rigid, with the associated dihedral assuming a fixed angle of 180°.⁴⁵

Fig. 1 Dihedrals φ₁ and φ₂ in phenylbenzoate and φ₁ in 2,5-diphenyl-1,3,4-oxadiazole.

Experimental and theoretical results (from this and previous work) for PB are summarized in Table 1. The torsional barrier for rotation around φ₁ from the theoretical calculations and gas-phase electron diffraction (GED) data show a reasonably high rotational barrier at 90° and a minimum energy dihedral angle, , close to 0°, indicating that π-conjugation is important in stabilizing a planar arrangement of the benzene ring and the C=O group. However, the GED barrier of 14.64 kJ mol⁻¹ suggests that the strength of π-conjugation may be less than that revealed by the DFT and MP2 calculations, and considerably smaller than that predicted by GAFF. It is generally recognized that molecules containing a double bond in conjunction with a benzene ring represent a special problem with respect to DFT calculations of the rotational barriers, resulting in over stabilization of the planar conformation and overestimated energy barriers.^45,47,48 However, the MP2/cc-pVTZ calculated barrier is similar to the most accurate DFT result, and so we have fitted the torsional profile to the MP2/cc-pVTZ result to give new GAFF-LCFF parameters (Table 2). It has been suggested that the inclusion of diffuse functions with some wave-function based methods, such as MP2, may provide a better description of delocalized electrons⁴⁹ and therefore it is likely that a larger aug-cc-pVTZ basis set may improve the results further.

Table 1 Lowest energy dihedral angles, and for the two key dihedral angles in phenyl benzoate (Fig. 1), and energy barriers for rotation

Method
a Results from ref. 45. b Results from ref. 46. c With bond length constraints in place. d Broad minimum energy well.
X-ray^a	−8.7		67.6
Exp. GED^a	0.0	14.64	64.0	5.02	0.13
GAFF	0.0	119^c	45.0	5.80	3.50
GAFF-LCFF (this work)	2.7	27^c	73.3^d	5.5	≈0.0
B3LYP/6-31G(d)^a	1.7	31.38	51.1	1.59	1.59
B3LYP/6-31G(d,p)			47.7	1.36	1.51
B3LYP/6-31+G(d,p)			63.7	2.83	0.26
B3LYP/6-311+G(d,p)	1.5	27.55	65.9	3.50	0.20
B3LYP/6-311+G(3df,3pd)	1.7	26.32	64.8	3.85	0.18
MP2/631+G(d)^b			71	12.52	0.57
MP2/631+G(d,p)^b			71	9.53	0.29
MP2/aug-cc-pVDZ			80.6	8.54	0.00
MP2/cc-pVTZ	1.3	26.41	66.4	5.99	0.18

---

Table 2 Optimized parameters for Ryckaert–Bellemans function (kJ mol⁻¹) obtained from fitting to ab initio data

Dihedral	Molecule	C0	C1	C2	C3	C4	C5
a Taken from ref. 50.
Dihedral φ1	Phenyl benzoate	7.335350	0.000000	−7.335350	0.000000	0.000000	0.000000
Dihedral φ2	Phenyl acetate	6.000000	0.000000	0.000000	0.000000	−6.000000	0.000000
Dihedral ϕ1	ODBP	4.000000	0.000000	−4.000000	0.000000	0.000000	0.000000
X–CH2–CH2–X^a	Heptane/butylbenzene	0.518587	−0.230192	0.896807	−1.491340	0.000000	0.000000

---

Both experimental values and theoretical results indicate that φ₂ is particularly flexible with a wide minimum energy region and a small main torsional barrier at 0° in addition to an insignificant barrier at 90° (see Fig. 2). The theoretical results are somewhat dependent on the level of theory and the basis sets employed in the calculations. Comparison with the GED data shows that the most accurate calculations for both the dihedral φ₂ angle and ΔE₀ barrier height are given by the B3LYP/6-311+G(3df,3pd) and MP2/cc-pVTZ calculations. The addition of diffuse functions in the DFT calculations increases the accuracy of the results. This has been attributed to the fact that electron lone pairs (e.g. the oxygen atom of dihedral φ₂) require the orbitals to occupy larger regions of space, which is better described with the addition of diffuse functions and has also been found to be particularly important when using DFT methods with the Pople basis sets.^45,51 The most accurate MP2 result in terms of the minimum energy structure and dihedral φ₂ torsional energy barriers is generated with the triple zeta, cc-pVTZ basis set. These are significantly closer to experiment than those seen previously in ref. 46, where the use of relatively small Pople basis sets show deviations from experimental values, in particular for the ΔE₀ torsional barrier. In contrast to the experimental values and the most accurate theoretical calculations, the GAFF minimum energy dihedral φ₂ angle is ≈45° as opposed to ≈65°, and there is a more significant torsional barrier at 90°. However, the main barrier at 0° is in reasonable agreement with the GED value. We have also calculated the torsional energy profile for rotation about the C(phenyl ring)–O(ester) bond in phenyl acetate at the MP2/cc-pVTZ level. As expected this is very similar to φ₂ with a minimum energy at 65° and a torsional barrier of just 1 kJ mol⁻¹ less than PB. New Ryckaert–Bellemans (RB) coefficients for GAFF-LCFF obtained by fitting to this angle are presented in Table 2.

Fig. 2 Experimental and calculated torsional energy profiles for dihedral φ₂ of PB.

In determining the structure of PB by GED, Tsuji et al.⁵² noted a relationship between mesogen core structure and the nematic–isotropic temperature. They compare PB with closely related mesogenic cores containing two phenyl rings but different linking units. The features that lowered the T_NI of PB relative to other closely related mesogens were: non-planarity of the core, a relatively large dihedral φ₂ angle of 64° for the minimum energy structure, a low torsional energy barrier for dihedral φ₂ and high flexibility of the phenyl ring attached to the ester oxygen atom, with the assumption that these structural features may be transferred into larger mesogens that contain the PB unit and may play a part in determining the T_NI. It is therefore important that GAFF should accurately represent these structural features of PB.

2,5-Diphenyl-1,3,4-oxadiazole (ODBP). The ODBP fragment (Fig. 1) represents the central core unit of oxadiazole based bent-core mesogens. These are a particularly interesting class of mesogens, displaying many unusual properties, including the (slightly controversial) possibility of forming a biaxial nematic phase.^8,53–57 There is no theoretical data available for the inter-ring (phenyl–heterocyclic ring) rotational energy barrier. However, there are numerous experimental structural studies of larger molecules containing the ODBP unit and these indicate a planar geometrical arrangement of the oxadiazole and phenyl rings is preferred.^58,59 For many bent-core mesogens, the flexibility of the central unit is important in determining the extent of local biaxial ordering,⁸ and has a major influence on the nematic–isotropic phase transition.

Our calculations show ODBP to prefer a planar geometry, with a torsional energy profile best fit by a single harmonic with barriers at 90° and 270°, reflecting the stabilization of planar conformation due to the π-conjugation between the phenyl and oxadiazole rings. The results (Table 3) show a barrier to rotation in the range 22.4–27.8 kJ mol⁻¹, and indicate a gradual decrease in the barrier height when the basis set is augmented with diffuse and multiple polarization functions. As mentioned previously, DFT calculations consistently over stabilize the planar conformation for π-conjugated systems. It would therefore seem reasonable to deduce that the barrier to rotation for the inter-ring dihedral is less than 23.47 kJ mol⁻¹, and therefore significantly lower than the predictions of the original GAFF force field. The MP2/cc-pVTZ calculation was used for fitting new RB coefficients (Table 2) to the calculated torsional profile.

Table 3 Rotational energy barriers calculated for dihedral ϕ₁ in the ODBP molecule

Method	ΔE90/kJ mol^−1
a Using bond length constraints.
GAFF	≈97.79^a
GAFF-LCFF (this work)	24.10^a
B3LYP/6-31G(d,p)	27.79
B3LYP/6-31+G(d,p)	24.17
B3LYP/6-311+G(d,p)	23.83
B3LYP/6-311+G(3df,3pd)	23.47
MP2/cc-pVTZ	22.39

---

n-Alkanes and n-butylbenzene. For flexible molecules such as n-alkanes, it is crucially important to correctly model the intramolecular interactions, as these contribute greatly to the flexibility of liquid crystal molecules, and therefore directly influence phase transition temperatures.

Many standard force fields, including GAFF, perform poorly in the reproduction of liquid properties and phase transitions of some n-alkanes.^50,60,61 There is some evidence that GAFF and OPLS-AA perform reasonably well for short-chain n-alkanes,^26,50,62,63 but significant deviations from experimental values are found for longer n-alkanes.^50,60 For example, the OPLS-AA force field results in liquid to gel phase transitions significantly higher than experimental values for longer n-alkanes.⁵⁰ Recently, a number of computer simulations employing the CHARMM27 and GROMOS force fields in the study of the structure and dynamics of lipid membranes have shown some disagreement with experimentally observed properties (e.g. area per lipid and NMR deuterium bond order parameters).⁶¹ These problems have largely been attributed to an under-prediction of the population of gauche states, leading to reduced flexibility, and in the case of lipid tails, an over estimation of the deuterium order parameters at the end of the chains.^60,61,64 Attempts to change the gauche and trans conformer ratio and rectify some of these problems have included refitting the torsional parameters of n-butane (as well as longer n-alkanes) to high quality ab initio data, or reducing the intramolecular vdW and electrostatic 1–4 scaling factors.^65,66 Noting that conformational energies are influenced by a combination of these factors, and to a much lesser extent many other terms in the force field, such as bending energies. This has led to continuous refinements of standard force fields for more accurate simulations of hydrocarbons and biomolecules (e.g. CHARMM27r, L-OPLS-AA GROMOS 43A2 and 45A3).^50,63,64,67 These findings prompted the current authors to re-examine the source data, both experimental and theoretical, used in the original parametrization of n-alkanes in standard force fields to look for some explanation for the discrepancies between the force field calculations and experimental values.

Experimental studies measuring the trans/gauche enthalpy difference (ΔH_tg) have predominately been performed in the gas phase using spectroscopic techniques. These have recently been summarized by Barna et al.⁶⁸ The results range from 2.08 to 4.58 kJ mol⁻¹. The lack of consistency in the results has largely been attributed to the complexity of the vibrational spectra of gaseous butane, and as early as 1991 doubts were raised by Murphy et al.⁶⁹ and then later by Herrebout (1995)⁷⁰ about previously calculated ΔH_g values, in particular those reporting larger values. Most recently, Balabin (2008) stated that the ratio of trans/gauche (t/g) conformer concentrations could only be predicted with an error margin of 40%. Despite these concerns, the two most accepted evaluations for ΔH_tg are Herrebout's value of 2.80 ± 0.40 kJ mol⁻¹ and Balabin's (2009) evaluation of 2.76 ± 0.09 kJ mol⁻¹, with the latter showing the least associated uncertainty.^70,71 These values have recently been revised by Barna et al. employing a more sophisticated statistical analysis of the original data and are presented in Table 4.

Table 4 Summary of the most reliable experimental and theoretical values for the trans/gauche energy and enthalpy difference for n-butane

ΔHg/kJ mol^−1	T/K	Method	Date^(ref.)
2.80 ± 0.09	133–196	Raman spectroscopy	2009^71
2.73 ± 0.52	223–297	Infrared spectroscopy	1995^70
2.71 ± 0.03	298	FPA/CCSD(T)	2012^68
2.83 ± 0.01	0	FPA/CCSD(T)	2012^68

---

ΔEg/kJ mol^−1
2.49 ± 0.01	0	FPA/CCSD(T)	2012^68

---

There are indications from experimental studies that the (t/g) liquid phase energy or enthalpy difference for n-butane is slightly less than in the gas phase, with the stability of the gauche conformer increasing by up to 0.42 kJ mol⁻¹.⁷¹ In terms of the rotational barriers of n-butane, it is not possible to measure these directly from spectroscopic data and methods are therefore based on estimates. These suggest that the t/g and cis barriers are comparable in energy, with one estimate giving values of 15.15 and 16.56 kJ mol⁻¹ for the t/g and cis barriers respectively.^72,73

The torsional energy about the dihedral C–C–C–C in n-butane involves 1–4 interactions between the methyl groups and it is therefore expected that the torsional energy about single bonds in longer n-alkanes may result in different values for the gauche energy and the rotational barriers.⁷³ There are indications from experiment that the gauche energy is slightly lower in longer n-alkanes compared to n-butane. For example, the conformational equilibration of n-pentane was studied by low-temperature gas-phase Raman spectroscopy by Balabin⁷¹ and resulted in a value of 2.59 kJ mol⁻¹ for the enthalpy difference between the trans–gauche and all trans states. This is slightly lower than his value of 2.76 kJ mol⁻¹ for n-butane. However, the author states that it is not clear whether the differences between the two alkanes are due to size differences or experimental uncertainty, and that further experimental values for various longer n-alkanes would be required to make general conclusions about the dependence of n-alkane size and conformer energies.

Similarly to n-butane, there is some evidence that the liquid phase energy or enthalpy difference between the gauche and trans conformers of longer n-alkanes is less than in the gas phase. For example, the results obtained from an infrared study gave an enthalpy difference of 2.08 ± 0.31 kJ mol⁻¹ for n-pentane, and an average value of 2.13 ± 0.21 kJ mol⁻¹ (ΔE_g) for liquid n-alkanes with n = 11–14 was obtained from a low frequency spectroscopy study.⁷⁴

The results of theoretical calculations of the t/g energy and enthalpy difference for butane are also somewhat inconsistent and appear to be dependent on the level of theory and specific ab initio technique used in the calculations.⁶⁸ The relative stability of the conformers of n-alkanes is thought to be strongly influenced by intramolecular dispersion interactions, with electrostatic effects playing a less important role.^75,76 However, the quantitative description of dispersion interactions still remains a great challenge for QM wave-function based methods, and in particular for density functional theories. The most accurate QM calculations are based on the coupled cluster (CC) single, doubles and perturbative triples {CCSD(T)} method. A number of sophisticated ab initio techniques have been developed to minimize the uncertainties in the calculations performed by improving the convergence of the electron correlation energy and addressing the problems arising from basis set incompleteness. These include focal point analysis (FPA) and Weizmann-n (Wn) methods, compound methods such as the Gaussian-2 (G2) method and complete basis set (CBS) methods.⁶⁸ Barna et al. claim that their study employing an improved ab initio method, with most of the energy contributions extrapolated to the CBS limit, currently provides the most reliable data for ΔE_g (see Table 4). It can be seen that their theoretical calculations for ΔH_g compare well with the experimental results. In terms of longer n-alkanes, there is some evidence from recent high level QM calculations that the t/g energy difference decreases with increasing chain length, at least up to n-octane.^63,73,77

Unlike experimental estimates, ab initio calculated energy differences for the rotational barriers of n-butane indicate that the cis barrier is significantly higher compared with the t/g barrier.^72,73 For example, Smith and Jaffe employed CCSD(T)/cc-pVTZ//MP2/6-311g(2df,p) level of theory, resulting in values of 13.85 ± 0.42 kJ mol⁻¹ for the t/g barrier and 22.93 ± 0.42 kJ mol⁻¹ for the cis barrier.⁷³ There is also some evidence that there is a small reduction in the t/g barrier compared to n-butane with increasing chain length.

The parameters for alkanes in standard force fields are generally obtained by fitting to different sets of QM and experimental data for n-butane. Table 5 shows a number of ΔE_g values for various force fields. With the exception of the CHARMM and GROMOS 43A2 force fields, all ΔE_g values are higher than the most recent QM and experimental enthalpy and energy values, considered to be the most reliable. For example, the MM3 parametrization of the conformational energetics of n-butane⁷⁸ was derived to be in accordance with the experimental results of Compton et al. (1980)⁷⁹ and Bartel et al. (1982)⁸⁰ that have since been superseded. The rather large values for AMBER99 and GAFF suggest that the ab initio conformational energies used for fitting alkane parameters were not of sufficiently high quality. The modified force fields CHARMM27r and GROMOS 43A2 show the best agreement with experiment and QM calculations.

Table 5 ΔE_g values for n-butane for a number of standard and modified force fields^63,78

Force field	ΔEg/kJ mol^−1
AMBER99	3.60
OPLS-AA	3.35
MM3	3.40
GAFF	4.50
CHARMM27	2.76
CHARMM27r	2.63
GROMOS 43A2	2.30

---

In light of these findings, and in particular the large ΔE_g value predicted by GAFF, it was considered that the n-alkane t/g conformer ratio needed to be amended, as a favouring of trans over gauche conformations in the alkane chains of LC molecules may be partly responsible for the high T_NI temperatures calculated with GAFF compared with the experimental values. Consideration was given to altering the 1–4 intramolecular scaling factors to achieve this. Reducing the electrostatic 1–4 interactions in particular, proved successful in tests on a number of n-alkanes, but raised considerable complexities when applied to the alkane fragments of larger LC molecules, as applying differential scaling factors was found to be a cumbersome approach. Instead, it was decided to reparametrize the torsional dihedrals of n-alkanes adopting the amended OPLS-AA torsional parameters specifically optimized for both short and long alkanes by Sui et al.⁵⁰ These were tested in GAFF-LCFF for n-heptane and butylbenzene, with the aim of reducing the t/g-energy difference. (Noting that butylbenzene is a common terminal structural component of many liquid crystal molecules.)

A comparison of the effective torsional profiles, obtained for the C₁–C₂–C₃–C₄ dihedral in n-heptane dihedral and n-butylbenzene, obtained by Boltzmann inversion from molecular dynamics torsional distribution functions, from the original GAFF and GAFF-LCFF are shown in Fig. 3. In both cases, the original GAFF force field yields a t/g-energy difference of ≈4.5 kJ mol⁻¹. The new parameters lead to a reduction in the t/g energy difference of ≈1.5 kJ mol⁻¹ to ≈3.0 kJ mol⁻¹ within the liquid phase, and consequently to a significant increase in chain flexibility. This in turn is likely to lead to reduced transition temperatures for many liquid crystal molecules.

Fig. 3 Effective torsional potentials obtained by Boltzmann inversion of dihedral angle distributions obtained from gas phase simulations of n-heptane and n-butylbenzene with GAFF and GAFF-LCFF.

Finally, we note that (specifically) for C₅ chains attached to a phenyl or cyclohexyl ring the effects of ΔE_tg are very subtle. A reduction in the ΔE_tg leads to a higher percentage of gauche conformations and hence more flexible chains. However, for terminal C₅ chains attached to the phenyl or cyclohexyl ring of a mesogen within a nematic phase, the most common conformations are ttt followed by tgt.^81–83 For alkyl (rather than alkoxy) chains an increase in tgt conformations enhances linearity of the mesogen, and hence enhances mesophase stability.^81,82 This balance of these two competing effects means that too high a value for ΔE_tg for C₅ chains has not been critical in past studies of common mesogens, such as CCH5^6,81 or 5CB.⁹ However, for chains longer than five carbons, the effects of increased chain flexibility are likely to be very significant.

4.2 Optimization of van der Waals parameters

Choice of fragments. The optimization of Lennard-Jones (LJ) parameters was initially guided by the work of Wang & Huo²⁶ and Caleman et al.³¹ (see Table S1 in the ESI† for a comparison of GAFF and experimental data for liquid crystal fragment molecules arising from their simulations). It was found that the predicted errors for ρ and Δ_vapH for a number of aromatic compounds could immediately be improved by tuning the aromatic carbon LJ parameters. Carboxylic acids containing highly polarizable groups with strong dipoles, such as C=O and OH groups, also produced very poor results with GAFF for ρ and Δ_vapH (Table S1, ESI†). To a lesser extent, the calculated properties for ester compounds displayed some significant deviations from experimental values. However, for aldehydes, ketones, alcohols and unbranched ethers, GAFF results show better agreement with experimental values. Given these considerations, the fragment esters, methylbenzoate, phenylacetate and methylformate were chosen for LJ parameter optimization as they regularly occur in calamitic mesogens. The fragments 1,3,4-oxadiazole and 2,5-diphenyl-1,3,4-oxadiazole, components of the central core unit of the ODBP bent-core mesogens, gave particularly poor density predictions with GAFF relative to the literature density predictions using the Advanced Chemistry Development (ACD/Labs) software.⁸⁴ (Here, we have used the ACD/Labs software predictions because very limited experimental data is available for these fragments, including no density data. Although there are uncertainties associated with the ACD calculations, examination of a number of calculations for small aromatic compounds for which there is available experimental data show good agreement with densities with a mean error of 0.38% – see ESI.†)

Caleman et al.³¹ have tested the GAFF force field for the heterocyclic compounds: pyrrole, 1,3-dioxalane, pyrimidine, morpholine and furan. Their calculated densities ranged from 3.7% to 9.8% greater than experimental values; and with the exception of pyrimidine, all calculated Δ_vapH were higher than experimental (+10.5 to +31.4% – see Table S1, ESI†). This suggests that the attraction between molecules is over-estimated for these compounds. The heterocyclic rings furan and pyrimidine share similar features with the 1,3,4-oxadiazole ring, for example π-electron density and heteroatoms with lone pairs. Optimization of LJ parameters for these compounds allowed the parameters to be reused for 1,3,4-oxadiazole and 2,5-diphenyl-1,3,4-oxadiazole.

GAFF performs well in predicting both ρ and Δ_vapH for small aliphatic molecules. For these compounds, we therefore used the original LJ parameters for GAFF in conjunction with updated torsional potentials. The updated RB-coefficients for n-heptane and n-butylbenzene, despite improving the t/g energy difference, make a negligible change to calculated ρ and Δ_vapH values for these compounds (Table 6).

Table 6 Density and heat of vaporization calculations for n-heptane and n-butylbenzene using GAFF and GAFF-LCFF with amended RB coefficients

Molecule	Property	T/K	Experiment	GAFF	New RB coefficients
n-Heptane	Density/g cm^−3	298	0.6788	0.6782 ± 0.0001	0.6783 ± 0.0001
	Heat of vap./kJ mol^−1	298	36.60	40.37 ± 0.03	40.12 ± 0.02
n-Butylbenzene	Density/g cm^−3	298	0.8559	0.8503 ± 0.0001	0.8505 ± 0.0001
	Heat of vap./kJ mol^−1	298	51.36	52.25 ± 0.02	51.95 ± 0.03

---

Aliphatic and aromatic esters compounds. Experimental ρ and Δ_vapH values for phenylacetate, methylbenzoate and methylformate are given in Table 7 along with those predicted by the original GAFF force field and GAFF-LCFF. The final parameters optimised for GAFF-LCFF are given in Table 8.

Table 7 Density and heat of vaporization calculations using the original and amended GAFF

Molecule	Property	T/K	Exp.	GAFF	% Diff.	New parameters	% Diff.
a Value taken from Advanced Chemistry Development (ACD/Labs) Software, V11.02, 1994–2014. b Value taken from Thermophysical Properties of Chemicals and Hydrocarbons, Carl L. Yaws, 2008, ch. 3 (pub. William Andrew). All other experimental values taken from ref. 31.
Phenylacetate	Density/g cm^−3	293	1.0739	1.1031	+2.7	1.0654 ± 0.0001	−0.8
	Heat of vap./kJ mol^−1	298	53.33	67.23 ± 0.02	+26.1	56.51 ± 0.02	+6.0
Methylbenzoate	Density/g cm^−3	298	1.0840	1.1105 ± 0.0001	+2.5	1.0796 ± 0.0002	−0.4
	Heat of vap./kJ mol^−1	298	54.28	65.41 ± 0.02	+20.5	54.98 ± 0.01	+1.3
Methylformate	Density/g cm^−3	298	0.9670	1.0467 ± 0.0001	+8.2	0.9880 ± 0.0001	+2.2
	Heat of vap./kJ mol^−1	298	30.59	38.89 ± 0.01	+27.1	32.14 ± 0.006	+5.1
Pyrimidine	Density/g cm^−3	298	1.0164	1.1022 ± 0.0003	+8.4	1.0246 ± 0.0002	+0.8
	Heat of vap./kJ mol^−1	298	49.81	48.63 ± 0.02	−2.4	46.87 ± 0.02	−5.9
Furan	Density/g cm^−3	298	0.9313	0.9495 ± 0.0003	+2.0	0.9379 ± 0.0005	+0.7
	Heat of vap./kJ mol^−1	298	27.46	29.27 ± 0.02	+6.6	27.67 ± 0.01	+0.8
1,3,4-Oxadiazole	Density/g cm^−3	293	1.1930^a	1.3093 ± 0.0002	+9.8	1.1959 ± 0.0002	+0.2
	Heat of vap./kJ mol^−1	298	37.1 ± 3.0^a	56.94 ± 0.04	+53.5	50.18 ± 0.01	+35.3
ODBP	Density/g cm^−3	293	1.1740^a	1.1950 ± 0.0003	+1.8	1.1695 ± 0.0002	−0.4
	Heat of vap./kJ mol^−1	563	49.964^b	80.26 ± 0.07	+60.6	57.97 ± 0.03	+16.0

---

Table 8 Amended vdW parameters and Ryckaert–Belleman parameters for selected dihedrals. Marked values * are the original GAFF parameters

Atom/description	σ/nm	ε ij /kJ mol^−1
O (carbonyl oxygen)	0.295992*	0.478608
C (aromatic carbon)	0.339967*	0.289824

---

Dihedral	C0	C1	C2	C3	C4	C5
C(Ar)–C(Ar)–C(sp^2)–O(sp^3) (in methyl benzoate)	7.335350	0.000000	−7.335350	0.000000	0.000000	0.000000
C(Ar)–C(Ar)–O(sp^3)–C(sp^2) (in phenyl acetate)	6.000000	0.000000	0.000000	0.000000	−6.000000	0.000000
C(sp^3)–C(sp^3)–C(sp^3)–C(sp^3)	0.518587	−0.230192	0.896807	−1.491340	0.000000	0.000000

---

Initial attempts to improve predicted properties for phenylacetate and methylbenzoate, focussed on the LJ parameter change for aromatic carbons reported by Wang and Hou.²⁶ These authors found that ρ and Δ_vapH were not sensitive to the radius parameter, σ_ij; but a well depth, ε_ij, that was slightly larger than the original, was required to reduce errors in the predicted values of both density and Δ_vapH. They also tested their new ε_ij for the aromatic compounds, phenol, m-cresol, aniline and fluorobenzene. Although the predicted properties for fluorobenzene were slightly better than those using the original GAFF parameters, for aniline and phenol only the Δ_vapH was improved, and for m-cresol both the predicted properties showed further deviations from the experimental data.

In our current study, using the Wang and Hou ε_ij parameter for the aromatic carbons of phenylacetate resulted in an increase in the errors obtained for both ρ and Δ_vapH. However, reducing the well depth of aromatic carbon by small increments resulted in an improvement in the predicted density and to a lesser extent Δ_vapH for phenylacetate, suggesting that LJ parameters for the aromatic carbons are not necessarily transferable to all aromatic compounds, i.e. there is some dependence of these effective pair potentials on chemical environment.

For both phenylacetate and methylbenzoate, the most accurate results were achieved by synchronously reducing the potential well depth of both the carbonyl oxygen atom as well as the aromatic carbons. New values of ε_ij = 0.289824 kJ mol⁻¹ for the aromatic carbons and ε_ij = 0.478608 kJ mol⁻¹ for the carbonyl oxygen atom produced the most accurate results for these bulk properties (Table 7). The new RB coefficients derived for the ester linkage of phenylacetate were also tested with the new LJ parameters. The new RB coefficients produce an insignificant change in the bulk properties, reinforcing the fact that ρ and Δ_vapH are almost exclusively associated with the non-bonded interactions.

The third fragment chosen for LJ parameter optimization, methylformate, showed significant deviations from experimental values for the bulk properties. An initial attempt to improve these predicted properties involved transferring the new ε value derived for the carbonyl oxygen atom of phenylacetate and methylbenzoate to this fragment. This resulted in a small reduction in the predicted errors with values of 1.0227 ± 0.002 g cm⁻³ for density and 35.38 kJ mol⁻¹ for Δ_vapH. However, reducing the well depth of the ester oxygen in addition to the carbonyl oxygen resulted in a significantly better agreement with the experimental data (Table 7).

The increased accuracy in the calculated bulk properties for these fragments does not necessarily imply any physical justification for the LJ parameter changes, but is most likely the result of a cancellation of errors due to the limitations of non-polarizable, atom-centred fixed charge force fields. The fragments studied in this work contain π-conjugated systems which are known to possess large quadruple moments, lone pairs as well as highly polarizable groups, such as C=O. Further improvements would probably require an improved description of intermolecular interactions using atomic multipoles and atomic polarizabilties.^85–87

Heterocyclic compounds. The original GAFF LJ parameters, give a significantly higher density for pyrimidine compared with experiment (Table 7), although the calculated Δ_vapH shows very good agreement with experiment. The nitrogens of pyrimidine in GAFF are assigned from earlier AMBER parameters derived for the basic nitrogens in adenine without further optimization. Additionally, the nitrogen well depth value of ε_ij = 0.71128 kJ mol⁻¹ in the original GAFF force field is significantly higher than those for analogous nitrogens (basic nitrogens with a lone pair) in other force fields. Literature values range from ε_ij = 0.4184 kJ mol⁻¹ (N in pyrimidine, Lopes et al.⁸⁸), ε_ij = 0.47369 kJ mol⁻¹ (basic N in ring structure, TraPPE⁸⁸), and ε_ij = 0.5857 kJ mol⁻¹ (amine nitrogen, Wang and Hou²⁶). In the current work, reduction of ε_ij by 0.1 kJ mol⁻¹ increments improved the density prediction error with ε_ij = 0.41128 kJ mol⁻¹ resulting in a density of 1.0598 ± 0.0002 g cm⁻³. However, this was still 4.7% higher than the experimental value and resulted in a significant deterioration in the calculated Δ_vapH (41.12 ± 0.02 kJ mol⁻¹ compared with an experimental value of 49.81 kJ mol⁻¹). Increasing the nitrogen σ_ij parameter by 0.01 nm along with this new ε_ij value improved the density prediction (1.0424 ± 0.0001 g cm⁻³) but had no effect on Δ_vapH which remained ≈−18% too low. This suggested that tuning the nitrogen LJ parameters alone was not sufficient to provide accurate predictions for both density and Δ_vapH.

GAFF assigns the same σ_ij and ε_ij parameters derived for aromatic carbons (ca) to all the carbon atoms of pyrimidine. However, the electronic nature of the C₂, C₄ and C₆ atoms, which are adjacent to the nitrogens, is different from the aromatic (ca) atoms of benzene, with a reduction in π-electron density and a tendency of the electrons to move towards the nitrogens. Retaining the original LJ parameters for nitrogen and increasing both the values of σ_ij and ε_ij for these carbons resulted in a very accurate density (1.0104 ± 0.0001 g cm⁻³) but a deterioration in the calculated Δ_vapH which was 14.5% higher than experiment. It was therefore decided to use the optimized ε_ij parameter for nitrogen described above and simultaneously tune the LJ parameters of the C₂, C₄ and C₆ atoms. It was found that a small increase in the values of σ_ij and ε_ij for these carbon atoms, along with the reduction in the nitrogen ε_ij parameter gave a significantly improved density and a reasonably good Δ_vapH. Although the latter was slightly worse than that with the original GAFF parameters, this was found to be the best compromise. The optimized LJ parameters for pyrimidine were, ε_NN = 0.41128 kJ mol⁻¹ and ε_CC = 0.42982 kJ mol⁻¹ for the carbons adjacent to the nitrogens, with ε_CC unchanged from the normal aromatic carbon for C₅.

In contrast to pyrimidine, the results for furan with the original GAFF parameters show a considerably smaller prediction error for ρ, but a larger prediction error for Δ_vapH. Furan, like benzene, is a π-electron rich aromatic compound and therefore it was decided to retain the original aromatic (ca) LJ parameters for the carbon atoms of the ring and focus on LJ parameter optimization for the oxygen atom. It was found that reducing the oxygen well depth (ε_ij = 0.61128 kJ mol⁻¹) brought both properties into good agreement with the experimental data.

In the absence of any experimental data on the bulk properties of the compounds 1,3,4-oxadiazole and 2,5-diphenyl,1,3,4-oxadiazole (ODBP), the predicted densities obtained with the ACD/labs software were used as references. These values suggest that the GAFF predicted densities are too high. Although there are uncertainties associated with the ACD calculations, examination of a number of calculations for compounds for which there is available experimental data show good agreement (see ESI† – Table S2). The 1,3,4-oxadiazole ring shares some features with furan (5-membered heterocyclic ring containing an oxygen atom). However, the electronic nature of the oxadiazxole ring is more closely related to that of pyrimidine, with relatively low π-electron density at carbon positions C2 and C5 and the presence of two basic nitrogen atoms that exert a withdrawal effect on the adjacent carbons. It was therefore decided to test the new LJ parameters derived for the nitrogens and carbons of pyrimidine and transfer these to the analogous atoms of the oxadiazole ring. This produced a density of 1.2028 g cm⁻³ which is close to the ACD calculation. In addition to these changes it was found that reducing the well depth of the oxygen atom of the oxadiazole ring, to that of furan described above, gave the best overall result for density when compared with the ACD result.

The original GAFF predicted density for the 2,5-diphenyl,1,3,4-oxadiazole fragment was slightly higher than the ACD result. Testing the new RB coefficients for the inter-ring dihedrals alone resulted in an insignificant increase to the GAFF predicted density. Adopting the new vdW parameters for the oxadiazole ring described above and retaining the new RB coefficients reduced the density to 1.1897 ± 0.0003 g cm⁻³. However, the best agreement with the ACD result was obtained through combining these changes with reducing the well depth of the carbon atoms of the phenyl rings to that derived for phenylacetate and methylbenzoate.

4.3 Testing the amended GAFF force field: simulation of 1,3-benzenedicarboxylic acid-1,3-bis(4-butylphenyl)ester mesogen

The mesogen 1,3-benzenedicarboxylic acid-1,3-bis(4-butylphenyl)ester (Fig. 4), containing aromatic and ester groups together with flexible chains, provides a particularly stringent test for GAFF. This mesogen shows a nematic phase with experimental phase transition temperatures of Cr–N 348 K, and N–I 452 K. One would expect the Cr–N transition to be very difficult to predict from simulation because of the ease of supercooling within standard simulations, together with sensitivity to system size. However, from previous simulation work, we would expect simulation to provide reasonable estimates for the nematic–isotropic phase transition temperature (T_NI).

Fig. 4 Structure of 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester.

A series of 80 ns runs at a range of temperatures, starting from a well-equilibrated isotropic system, were employed to obtain an approximate T_NI for the original GAFF and new GAFF-LCFF models. Close to the transition (as shown in Fig. 5 for four temperatures) large temporal fluctuations are seen in the nematic order parameter, S₂. Following Zhang et al.,²⁵ we use a value of S₂ > 0.4 to denote a stable nematic phase. For these lengths of simulation, we can reasonably achieve a prediction for T_NI of within ±5 K for individual model mesogens. (Closer than 5 K to the transition we expect whole simulation box fluctuations in orientational order to occur on a timescale longer than 80 ns. So a prediction of better than 5 K is not possible without a much larger simulation.)¹⁹

Fig. 5 Time evolution of the nematic order parameter, S₂, for 1,3-benzenedicarboxylic acid-1,3-bis(4-butylphenyl)ester at four different temperatures starting from a well equilibrated isotropic configuration at 550 K.

Fig. 6 plots the mean order parameter for the two models as the system is progressively cooled from the high temperature disordered liquid. The new GAFF-LCFF force field shows markedly different behaviour to the original GAFF. The combination of alkyl chains which are too stiff, ester groups which are not flexible enough and LJ interactions which are slightly too attractive lead to the original GAFF overestimating T_NI by ≈60 K. However, the re-optimized amended GAFF performs very well indeed exactly predicting the phase transition temperatures within the level of accuracy, ±5 K, possible for this size and length of simulation.

Fig. 6 Mean nematic order parameter S₂ as a function of temperature for the original and amended GAFF force fields. The vertical black line indicates the experimental transition temperature.

Radial distribution functions were used to check the identity of the low temperature phase. These are presented in Fig. 7 for selected temperatures. The standard radial distribution function g(r) exhibits liquid-like behaviour over the temperatures shown in Fig. 7, with a characteristic peak at short range, ≈5 to 7 Å, followed by convergence to a value of one at longer range. At the lower temperatures of 400 K and 360 K, the main peak is split into two subsidiary peaks. Additionally, the magnitude of the first peak increases with decreasing temperature, suggesting stronger local correlations between neighbouring molecules. Examination of the orientational correlation function confirms that the phase is isotropic at 460 K with g₂(r) decaying to zero at long range above this temperature. However, at 450 K, g₂(r) converges to a value of ≈0.26 at large r distances, consistent with an S₂² = 0.51² = 0.2601 and commensurate with an orientationally ordered phase at this temperature. The lack of structure in g₂(r) is consistent with a nematic phase.

Fig. 7 Radial distribution functions, top: g(r); middle: g₂(r); and bottom: g_⊥(r); calculated as a function of distance between the centre of mass of 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester for a series of temperatures.

To rule out the possibility of any translational ordering of the systems at the temperatures expected to be nematic, the contribution to the radial distribution function parallel, g_∥(r) and perpendicular g_⊥(r) to the director were also examined. g_∥(r) (not shown) gives a value of one at all temperatures examined, indicating no layering; and with the exceptions of a small peak at short range, g_⊥(r) displays minimal structure at 450 K, 400 K and 360 K (Fig. 7) confirming the nematic nature of the phase at these temperatures.

In optimizing GAFF, we noted the particular sensitivity of T_NI to both the density of the system and to molecular flexibility. This is particularly noticeable for elongated calamitic molecules. In these cases, although the systems are thermotropic with both anisotropic attractive and anisotropic repulsive interactions, the molecular shape is sufficiently rod-like for repulsive interactions to exert a dominant influence on the phase transition. It is interesting to compare such calamitic systems with hard colloidal rods. In the latter the density of rods controls the balance between the competing effects of rotational and translational entropy. Hence longer rods form nematics at lower densities.⁸⁹ For thermotropic calamitics, if the molecule is made more flexible and therefore less rod-like, we see an immediate reduction in the phase stability and hence the transition temperature. Likewise if LJ parameters are made slightly less attractive, leading to a small decrease in density, the transition temperature is also reduced. This suggests that the strength of anisotropic attractive interactions alone, as used in Maier and Saupe theory,⁹⁰ are unlikely to be sufficient to explain the changes in nematic stability induced by subtle changes in intermolecular interactions.

Previously, the relative roles of attractions and repulsions in determining the isotropic to nematic phase transition have been assessed by considering the quantity

which measures the relative dependence of the orientational order parameter on density and temperature.⁹⁰ Considering hard-core interactions only, in an athermal system, γ = ∞, whereas with angle-dependent attractions only, γ = 1. Measurements on a real mesogen for example, para-azoxyanisole (PAA) show γ to be ≈4 which is in accordance with estimates provided by combined models. Additionally, theoretical calculations show that γ is very sensitive to the packing fraction, validating the dominant role of hard-core interactions at high density.⁹⁰ This seems entirely consistent with the results of the current study.

The influence of molecular shape on the location of T_NI was also examined in the current study. An indication of the molecular dimensions, and hence overall shape, can be obtained from the averaged principle moments of inertia 〈I₁〉, 〈I₂〉 and 〈I₃〉. These values enable the average length 2a, width 2b and breadth 2c of a mesogen to be calculated (Table 9), using and cyclic permutations for b and c.⁸¹ Both force fields show a small increase in molecular length with decreasing temperatures. However, the main difference occurs between the two force fields, with new GAFF-LCFF resulting in a decrease in length of ≈0.23 Å (in the isotropic phase at 550 K) as well as a small increase in molecular width and breadth compared with GAFF. This small change in molecular length with force field is significant. Tiberio et al.⁵ in their investigations of the linear T6 mesogen, found a decrease in the average molecular length of less than 0.2 Å in going from the isotropic to the nematic phase. It is probable that the amended torsional potentials introduced into the phenylester-LC have increased its flexibility, enabling the molecular structure to sample a broader range of configurations. This is likely to reduce the length : breadth ratio and contribute to a lower T_NI temperature, as expected from results of DPD simulations of semi-rigid mesogens⁹¹ and semi-flexible chains of hard spheres.⁸⁹

Table 9 The average length, 2a, width, 2b, and breadth, 2c, of the phenylester-LC molecule at the simulated temperatures of 480 K, 500 K and 550 K for original and the new GAFF-LCFF force fields

	〈2a〉/Å	〈2b〉/Å	〈2c〉/Å	a/(b + c)
Original GAFF
550 K (isotropic)	29.63	5.86	3.63	6.24
500 K (nematic)	29.85	5.70	3.54	6.46
480 K (nematic)	29.93	5.62	3.51	6.55
New GAFF-LCFF
550 K (isotropic)	29.50	5.85	3.68	6.19
500 K (isotropic)	29.62	5.75	3.65	6.30
480 K (isotropic)	29.68	5.70	3.63	6.36

---

To test the relative effects of the changes in LJ parameters and torsions, we repeated the transition temperature predictions using GAFF-LCFF torsions together with the original GAFF LJ parameters. The new value of T_NI was found to be 15 K lower than found with GAFF, indicating that the torsional parameters have a significant effect, even though the main improvement in T_NI arises from optimized LJ parameters. For some liquid crystal systems we expect the accuracy of torsional parameters to be even more important than for 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester. In recent work on quinquephenyl, Olivier et al.⁹² have shown that the polydispersity in aspect ratio, which arises in quinquephenyl from the inter-ring torsional angles (together with a small amount of bond bending), is important in stabilizing the nematic phase relative to smectic phases. It will be important to test further the transferability of the improved torsions of GAFF-LCFF on a range of other liquid crystals.

Finally, we note that in this work, it has proved necessary in some cases to use different Lennard-Jones parameters for situations where the same type of atom appears in different molecular environments (specifically aromatic carbons). This implicitly recognises the limitations of effective two-body potentials, where the influence of higher-body effects are averaged into the two body potentials on fitting. As a consequence transferability is reduced when the surrounding environment changes. In our small molecule testing of GAFF-LCFF this wider range of Lennard-Jones parameters has clearly improved the transferability of the force field in comparison to the original GAFF parameter set. Again, it will be important to test the transferability of GAFF-LCFF further using a wide range of different mesogens.

5 Conclusions

In summary, we have amended the GAFF force field by:

• careful tuning of a selected number of LJ parameters of component fragments of standard calamitic mesogens with the aim of reproducing the experimental properties, density and Δ_vapH, in particular, to obtain a density deviation of less than 1% from experimental values;

• re-parametrization of a number of torsional potentials for fragment molecules using high-level quantum chemical calculations, with the aim of improving the description of the overall ‘shape’ and flexibility of the mesogen.

MD simulations employing the new amended, GAFF-LCFF, force field provide a very good estimate of the experimental T_NI for the 1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)ester, reducing the original GAFF prediction of the T_NI temperature by ≈60 K.

GAFF-LCFF is being tested further in our laboratory on a number of other systems with promising results. For C5-Ph-ODBP-Ph-OC12 and C4-Ph-ODBP-Ph-C7 (members of the bis-(phenyl)oxadiazole family of bent core mesogens), GAFF-LCFF predicts transition temperatures within 10 K of experiment. This prediction is within typical system size dependency errors for 256 molecules, and represents >100 K improvement on the over-estimated T_NI values obtained through use of the standard GAFF parameter set. It will be interesting to test the general applicability of GAFF-LCFF on a wider range of liquid crystal systems in future work.

Acknowledgements

The authors would like to thank Durham University for providing computer time on its high performance computer system, Hamilton; and are grateful for the use of the S.H.E.D. facility in Middleton-in-Teesdale.

References

1. M. R. Wilson, Int. Rev. Phys. Chem., 2005, 24, 421–455.
2. M. R. Wilson, Chem. Soc. Rev., 2007, 36, 1881–1888.
3. C. M. Care and D. J. Cleaver, Rep. Prog. Phys., 2005, 68, 2665–2700.
4. I. Cacelli and G. Prampolini, J. Chem. Theory Comput., 2007, 3, 1803–1817.
5. A. Pizzirusso, M. Savini, L. Muccioli and C. Zannoni, J. Mater. Chem., 2011, 21, 125–133.
6. M. R. Wilson and M. P. Allen, Mol. Cryst. Liq. Cryst., 1991, 198, 465–477.
7. M. J. Cook and M. R. Wilson, Mol. Cryst. Liq. Cryst., 2001, 363, 181–193.
8. J. Peláez and M. R. Wilson, Phys. Rev. Lett., 2006, 97, 267801.
9. J. Peláez and M. R. Wilson, Phys. Chem. Chem. Phys., 2007, 9, 2968–2975.
10. D. L. Cheung, S. J. Clark and M. R. Wilson, J. Chem. Phys., 2004, 121, 9131–9139.
11. D. L. Cheung, S. J. Clark and M. R. Wilson, Chem. Phys. Lett., 2002, 356, 140–146.
12. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
13. C. Oostenbrink, A. Villa, A. E. Mark and W. F. Van Gunsteren, J. Comput. Chem., 2004, 25, 1656–1676.
14. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
15. I. Cacelli, A. Cimoli, P. R. Livotto and G. Prampolini, J. Comput. Chem., 2012, 33, 1055–1067.
16. G. Tiberio, L. Muccioli, R. Berardi and C. Zannoni, ChemPhysChem, 2009, 10, 125–136.
17. M. F. Palermo, A. Pizzirusso, L. Muccioli and C. Zannoni, J. Chem. Phys., 2013, 138, 204901.
18. F. Chami, M. R. Wilson and V. S. Oganesyan, Soft Matter, 2012, 8, 6823–6833.
19. E. Kuprusevicius, R. Edge, H. Gopee, A. N. Cammidge, E. J. L. McInnes, M. R. Wilson and V. S. Oganesyan, Chem. – Eur. J., 2010, 16, 11558–11562.
20. I. Cacelli, C. F. Lami and G. Prampolini, J. Comput. Chem., 2009, 30, 366–378.
21. I. Cacelli, G. Prampolini and A. Tani, J. Phys. Chem. B, 2005, 109, 3531–3538.
22. I. Cacelli, A. Cimoli, L. De Gaetani, G. Prampolini and A. Tani, J. Chem. Theory Comput., 2009, 5, 1865–1876.
23. E. M. Mas, R. Bukowski and K. Szalewicz, J. Chem. Phys., 2003, 118, 4386–4403.
24. E. M. Mas, R. Bukowski and K. Szalewicz, J. Chem. Phys., 2003, 118, 4404–4413.
25. J. Zhang, J. Su and H. Guo, J. Phys. Chem. B, 2011, 115, 2214–2227.
26. J. Wang and T. Hou, J. Chem. Theory Comput., 2011, 7, 2151–2165.
27. D. A. Case, T. E. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang and R. J. Woods, J. Comput. Chem., 2005, 26, 1668–1688.
28. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
29. I. Todorov and W. Smith, The DL POLY 4 user manual, 2011.
30. M. R. Wilson, unpublished work.
31. C. Caleman, P. J. van Maaren, M. Hong, J. S. Hub, L. T. Costa and D. van der Spoel, J. Chem. Theory Comput., 2012, 8, 61–74.
32. R. Doe, http://virtualchemistry.org, 2014.
33. D. L. Cheung, S. J. Clark and M. R. Wilson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 051709.
34. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ã. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian09 Revision D.01, Gaussian Inc., Wallingford, CT, 2009.
35. J. C. Sancho-Garcia, J. L. Bredas and J. Cornil, Chem. Phys. Lett., 2003, 377, 63–68.
36. J. C. Sancho-Garcia and J. Cornil, J. Chem. Phys., 2004, 121, 3096–3101.
37. D. van der Spoel, E. Lindahl, B. Hess, R. van Buuren, E. Apol, P. Meulenhoof, D. Tieleman, A. Sijbers, K. Feenstra, R. van Drunen and H. Berendsen, Gromacs User Manual, version 4.0.
38. A. W. S. da Silva and W. F. Vranken, BMC Res. Notes, 2012, 5, 367.
39. M. R. Wilson, J. Chem. Phys., 1997, 107, 8654–8663.
40. R. Pecheanu and N. M. Cann, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041704.
41. R. Korlacki, M. Steiner, A. J. Meixner, J. K. Vij, M. Hird and J. W. Goodby, J. Chem. Phys., 2007, 126, 224904.
42. H. Kurosu, M. Kawasaki, M. Hirose, M. Yamada, S. M. Kang, M. Sone, H. Takezoe and J. Watanabe, J. Phys. Chem. A, 2004, 108, 4674–4678.
43. T. Imase, S. Kawauchi and J. Watanabe, J. Mol. Struct., 2001, 560, 275–281.
44. R. Y. Dong and A. Marini, J. Phys. Chem. B, 2009, 113, 14062–14072.
45. R. Wrzalik, K. Merkel and A. Kocot, J. Mol. Model., 2003, 9, 248–258.
46. G. Cinacchi and G. Prampolini, J. Phys. Chem. A, 2005, 109, 6290–6293.
47. R. J. Meier, Comput. Mater. Sci., 2003, 27, 219–223.
48. R. J. Meier and E. Koglin, Chem. Phys. Lett., 2002, 353, 239–243.
49. J. C. Sancho-Garcia and A. J. Perez-Jimenez, J. Chem. Phys., 2003, 119, 5121–5127.
50. S. W. Siu, K. Pluhackova and R. A. Böckmann, J. Chem. Theory Comput., 2012, 8, 1459–1470.
51. A. D. Boese, J. M. L. Martin and N. C. Handy, J. Chem. Phys., 2003, 119, 3005–3014.
52. T. Tsuji, H. Takeuchi, T. Egawa and S. Konaka, J. Am. Chem. Soc., 2001, 123, 6381–6387.
53. V. Gortz, C. Southern, N. W. Roberts, H. F. Gleeson and J. W. Goodby, Soft Matter, 2009, 5, 463–471.
54. O. Francescangeli, V. Stanic, S. I. Torgova, A. Strigazzi, N. Scaramuzza, C. Ferrero, I. P. Dolbnya, T. M. Weiss, R. Berardi, L. Muccioli, S. Orlandi and C. Zannoni, Adv. Funct. Mater., 2009, 19, 2592–2600.
55. O. Francescangeli and E. T. Samulski, Soft Matter, 2010, 6, 2413–2420.
56. O. Francescangeli, F. Vita, F. Fauth and E. T. Samulski, Phys. Rev. Lett., 2011, 107, 207801.
57. L. A. Madsen, T. J. Dingemans, M. Nakata and E. T. Samulski, Phys. Rev. Lett., 2004, 92, 145505.
58. A. Carella, A. Castaldo, R. Centore, A. Fort, A. Sirigu and A. Tuzi, J. Chem. Soc., Perkin Trans. 2, 2002, 1791–1795.
59. C. H. Mao, Q. M. Wang, R. Q. Huang, L. Chen, J. Shang and H. B. Song, Acta Crystallogr., Sect. E: Struct. Rep. Online, 2004, 60, O1823–O1825.
60. C. J. Dickson, L. Rosso, R. M. Betz, R. C. Walker and I. R. Gould, Soft Matter, 2012, 8, 9617–9627.
61. A. P. Lyubartsev and A. L. Rabinovich, Soft Matter, 2011, 7, 25–39.
62. W. L. Jorgensen, D. S. Maxwell and J. TiradoRives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
63. J. Klauda, B. Brooks, A. MacKerell, R. Venable and R. Pastor, J. Phys. Chem. B, 2005, 109, 5300–5311.
64. L. D. Schuler and W. F. van Gunsteren, Mol. Simul., 2000, 25, 301–319.
65. X. Ye, S. Cui, V. F. de Almeida and B. Khomami, J. Mol. Model., 2013, 19, 1251–1258.
66. C.-J. Hogberg, A. M. Nikitin and A. P. Lyubartsev, J. Comput. Chem., 2008, 29, 2359–2369.
67. L. D. Schuler, X. Daura and W. F. Van Gunsteren, J. Comput. Chem., 2001, 22, 1205–1218.
68. D. Barna, B. Nagy, J. Csontos, A. G. Csaszar and G. Tasi, J. Chem. Theory Comput., 2012, 8, 479–486.
69. W. Murphy, J. Fernandezsanchez and K. Raghavachari, J. Phys. Chem., 1991, 95, 1124–1139.
70. W. Herrebout, B. Vanderveken, A. Wang and J. Durig, J. Phys. Chem., 1995, 99, 578–585.
71. R. M. Balabin, J. Phys. Chem. A, 2009, 113, 1012–1019.
72. M. J. Hafezi and F. Sharif, THEOCHEM, 2007, 814, 43–49.
73. G. Smith and R. Jaffe, J. Phys. Chem., 1996, 100, 18718–18724.
74. J. Scherer and R. Snyder, J. Chem. Phys., 1980, 72, 5798–5808.
75. J. M. L. Martin, J. Phys. Chem. A, 2013, 117, 3118–3132.
76. J. Miao, S. Hua and S. Li, Chem. Phys. Lett., 2012, 541, 7–11.
77. J. Klauda, R. Pastor and B. Brooks, J. Phys. Chem. B, 2005, 109, 15684–15686.
78. N. L. Allinger, Y. H. Yuh and J. H. Lii, J. Am. Chem. Soc., 1989, 111, 8551–8566.
79. D. Compton, S. Montero and W. Murphy, J. Phys. Chem., 1980, 84, 3587–3591.
80. R. Heenan and L. Bartell, J. Chem. Phys., 1983, 78, 1270–1274.
81. M. R. Wilson and M. P. Allen, Liq. Cryst., 1992, 12, 157–176.
82. M. Wilson, Liq. Cryst., 1996, 21, 437–447.
83. C. W. Cross and B. Fung, J. Chem. Phys., 1994, 101, 6839.
84. Advanced Chemistry Development (ACD/Labs) Software V11.02, 1994–2014.
85. C. Kramer, P. Gedeck and M. Meuwly, J. Chem. Theory Comput., 2013, 9, 1499–1511.
86. A. Holt, J. Bostrom, G. Karlstrom and R. Lindh, J. Comput. Chem., 2010, 31, 1583–1591.
87. P. Cieplak, F.-Y. Dupradeau, Y. Duan and J. Wang, J. Phys.: Condens. Matter, 2009, 21, 333102.
88. N. Rai and J. I. Siepmann, J. Phys. Chem. B, 2007, 111, 10790–10799.
89. M. R. Wilson, Mol. Phys., 1994, 81, 675–690.
90. W. M. Gelbart, J. Phys. Chem., 1982, 86, 4298–4307.
91. Z. Zhang and H. Guo, J. Chem. Phys., 2010, 133, 144911.
92. Y. Olivier, L. Muccioli and C. Zannoni, ChemPhysChem, 2014, 15, 1345–1355.

---

Footnote

† Electronic supplementary information (ESI) available: Fig. S1. The time dependence of the orientational order parameter at different temperatures for the GAFF force field using the original Lennard-Jones potentials and torsions amended in this study. Table S1. List of the experimental densities and heats of vaporization used to assess GAFF for liquid crystal fragment molecules. Table S2. Comparison of ACD Labs calculated densities and experimental values. See DOI: 10.1039/c5cp03702f

---