Locally ordered junctions govern diffusion in triglycerides: insights from molecular dynamics simulations

Abstract

Triglycerides (TGs) exhibit diverse flow behaviours that are essential for applications ranging from food and pharmaceuticals to renewable fuels and advanced lubricants. Despite sharing nearly identical atomic compositions and similar bulk thermodynamic properties, different TGs display markedly different viscosities, and the microscopic origins of these distinct behaviours remain unclear. Herein, we performed all-atom molecular dynamics simulations to investigate three representative TGs with different chain lengths and degrees of unsaturation: trioctanoin (8 : 0), triolein (18 : 1), and trilinolenin (18 : 3). Our results qualitatively reproduced the experimental viscosity ordering, with 18 : 1 being the most viscous and 8 : 0 showing a slightly higher viscosity than 18 : 3. Neither the cohesive energy density, melting point correlation, nor molecular size could account for these differences. Likewise, single-molecule conformational statistics and mesoscale cluster morphologies revealed no decisive distinctions that explain the observed viscosity variations. Instead, structural analyses revealed that the local parallel alignment of C-chain segments produced junctions that formed extended network structures. These networks were abundant in 18 : 1, intermediate in 8 : 0, and absent in 18 : 3, directly correlating with the observed viscosities. Overall, fine-scale packing appears to govern the dynamics of TGs.

---

1 Introduction

Triglycerides (TGs), consisting of three fatty acids esterified to a glycerol backbone, exhibit remarkable structural diversity in both the carbon-chain length and degree of unsaturation. This structural variability governs their physical and chemical behaviour, enabling a wide range of applications, including energy storage and conversion (e.g., vegetable oils and animal fats as renewable fuels^1–3), advanced lubricants in tribological systems,^4–7 and food and pharmaceutical crystallisation,^8,9 as well as biological functions such as tear-film stability,¹⁰ adipocyte energy storage,¹¹ and cutaneous protection via sebum.¹²

TGs are commonly referred to as edible oils and have been extensively used in human diets for centuries. According to the U.S. Department of Agriculture Economic Research Service, global vegetable oil production and consumption are projected to exceed 220 million tons in the 2025/26 fiscal year.¹³ In addition to conventional edible oils, functional products with health benefits—such as medium-chain TG oil—are now commercially available and widely marketed. The flow behaviour of these oils is determined by their viscosity, which directly affects food processing in homes, industries, and restaurants,^14,15 and may influence aspects of nutritional utilisation.

As a general feature of oils, several thermodynamic properties, such as density and cohesive energy, are markedly similar among TGs because they share the same atomic composition. In contrast, their dynamic properties significantly vary. Empirical trends observed across fatty acid and edible oil viscosity data suggest that viscosity increases with chain length and decreases with the number of double bonds.^14,16 These observations are partially consistent with the behaviours of the TGs examined in the present study. For instance, trilinolenin (18 : 3) exhibits a viscosity of around 24.64 mPa s, which is lower than that of triolein (18 : 1, 84.18 mPa s), in line with the expected fluidising effect of multiple double bonds. However, despite its substantially longer carbon chains, trilinolenin also shows a lower viscosity than the medium-chain trioctanoin (8 : 0, 26.65 mPa s), as confirmed by our rheological measurements (Table 1). Such a trend cannot be rationalised by conventional empirical correlations. This discrepancy indicates that viscosity cannot be explained solely by chain length or unsaturation; rather, a molecular-level framework based on collective organisation and intermolecular energetics is required to capture the distinct dynamics of TGs.

Table 1 D, η, and structural/thermodynamic properties of the TGs. Values in parentheses for η are experimental

	8 : 0	18 : 1	18 : 3
D^∞ (10^−12 m^2 s^−1)	5.2	0.68	3.6
η (mPa s)	100 (26.65)	345 (84.18)	92 (24.64)
Rg (nm)	0.52	0.88	0.84
Rh (nm)	0.42	3.2	0.61
E^dc (J cm^−3)	−315	−312	−317
ρ (g cm^−3)	0.955	0.917	0.940

---

Molecular dynamics (MD) simulations have been employed to elucidate various TG phenomena, such as the crystallisation of tripalmitin and triolein at physiological temperatures,¹⁷ coarse-grained modelling and crystallisation of systems of tridecanoin and mixed stearin,^18,19 and interfacial assembly of tri-cis-6-hexadecenoin at air–water interfaces,²⁰ among others.^{2,3,7,9,21–27} These studies have clarified the packing, phase behaviour, and interfacial ordering in TGs. However, a direct connection between molecular-scale structure and viscosity has not been established, and thus the microscopic origin of their distinct dynamic behaviour remains elusive.

In this work, we employ all-atom MD simulations to investigate three homologous TG molecules differing in the chain length and degree of unsaturation, namely trioctanoin (8 : 0), triolein (18 : 1), and trilinolenin (18 : 3) (Fig. 1). By systematically comparing their diffusion and viscosity, we demonstrate that bulk thermodynamic properties are insufficient to explain the observed trends. By examining the alkyl chain conformational distributions, local alignment, and collective network organisation, we observe how fine-scale structural features govern the viscosity of TGs. This molecular-level framework provides insight for the design of oils and lubricants with tailored flow properties.

Fig. 1 Chemical structures of the TGs denoted as 8 : 0, 18 : 1, and 18 : 3.

2 Methods

2.1 Diffusion and shear viscosity coefficients

The self-diffusion coefficient, D, was obtained from the Einstein relation:where 〈⋯〉 denotes an ensemble average, r_i(t) is the centre-of-mass position of molecule i at time t, and N is the number of molecules. Experimentally, D was evaluated from the slope of the linear regime of the mean-squared displacement (MSD) plotted against time.²⁸ It is well known that the D value depends on the simulation cell length L. Therefore, we calculated D^∞, which is consistent with experimental values, using the correction formula: D^∞ = D + 2.83729k_BT/(6πηL), where k_B, T, and η are the Boltzmann constant, temperature, and viscosity, respectively.²⁹

The shear viscosity, η, was computed using the Green–Kubo relation:²⁸

where P_αβ(t) is an off-diagonal component (α ≠ β ∈ {x, y, z}) of the pressure tensor, V is the system volume, k_B is the Boltzmann constant, and T is the temperature. Zhang et al. proposed the time decomposition method to remove long-time tails in the correlation function.³⁰ Accordingly, for the viscosity calculations, 300 ns trajectories were analysed using 31 overlapping segments; each segment was 50 ns in length and initiated every 10 ns. MSDs and η were calculated using the gmx msd and gmx energy tools in GROMACS 2022 and 2024, respectively.

2.2 Cohesive energy density

The cohesive energy, E_c, of a liquid is the energy required to separate it into independent gaseous molecules. Following previous work,^31,32 we defined

E_c = E_system − E_{isolated, total} = E_inter, (3)

where E_system is the potential energy of the simulated system, E_{isolated, total} is the sum of the potential energies of the N isolated molecules, and E_inter is the intermolecular potential energy. For isolated molecules, single-molecule calculations were performed in vacuum using the same simulation conditions employed in the condensed-phase simulations. Because cohesive energy scales with system size, we report the cohesive energy density (CED) as:where V_cell is the simulation cell volume.

2.3 Definition of TG aggregates

TG aggregates were identified using Voronoi polyhedral analysis.^33,34 A Voronoi region around a given atom is defined as the set of all points in space closer to that atom than to any other. The analysis was performed on the glycerol backbone atoms, and molecules whose glycerol atoms shared a Voronoi face were assigned to the same aggregate. Periodic boundary conditions were accounted for by replicating the simulation cell when constructing neighbour lists for the Voronoi analysis. Voronoi polyhedral analysis was carried out using the Qhull package.³⁵

2.4 Distribution of the C-chain alignment

We calculated the joint distribution function f (r, cos θ), which gives the probability of finding two C-chain units separated by a distance r between their central carbon atoms, with a relative orientation cos θ. Here, a “C-chain unit” refers to a contiguous segment of seven carbons in 8 : 0 or eight carbons in 18 : 1 and 18 : 3 (see Fig. 1). Accordingly, each fatty acid chain of 8 : 0 and 18 : 3 contains one C-chain unit, whereas each fatty acid chain of 18 : 1 contains two C-chain units. A pronounced feature at r ≈ 0.5 nm and cos θ ≈ 1 corresponds to parallel alignment of the C-chain segments.

2.5 Bond correlation function

The bond correlation function between C-chain units, B(t), was defined analogously to the hydrogen-bond correlation function:²⁶

B(t) = 〈b(t)b(0)〉, (5)

where b(t) is bonded to another C-chain unit at time t and equals 0 otherwise. A bond was considered to exist if all the following criteria were satisfied:

• The distance between the centre atoms of the two C-chain units was less than 0.6 nm (contact condition).

• The end-to-end length of the C-chain was larger than 0.8 nm for 8 : 0 or 1.0 nm for 18 : 1 and 18 : 3 (straight condition).

• The relative orientation satisfied cos (θ) > 0.95 (orientation condition).

The angle brackets denote an ensemble average.

2.6 Molecular dynamics calculations

All-atom MD simulations were performed using GROMACS 2022.³⁶ The initial configurations of 4096 molecules were generated by random packing using Packmol,^37,38 followed by energy minimisation. Periodic boundary conditions were applied in all three directions. Long-range electrostatics were treated using the smooth particle–mesh Ewald method,³⁹ with a grid spacing of 0.10 nm and an Ewald alpha of 2.6 nm⁻¹. The real space cut-off for Lennard–Jones interactions and short-range electrostatics was 1.2 nm. These parameters ensure energy conservation to within five significant digits, a commonly accepted criterion in MD simulations. Long-range dispersion corrections were applied to account for energy and pressure contributions. The temperature was maintained at 298.15 K using the velocity-rescaling thermostat⁴⁰ with a coupling constant of 1.0 ps, and the pressure was controlled at 1.0 bar using the Parrinello–Rahman barostat⁴¹ with a coupling constant of 1.0 ps. The integration time step was 1 fs. Bonds involving hydrogen atoms were constrained using the LINCS algorithm.⁴² The OPLS-AA force field, which is well-established for organic molecules, was employed.⁴³ Each system was simulated for 500 ns in the NPT ensemble. The first 200 ns were treated as equilibration, and the subsequent 300 ns were used for the analysis of dynamic properties.

2.7 Experimental viscosity measurements

Viscosities of the TG samples were measured using a stress-controlled rheometer (MCR 102, Anton Paar GmbH, Graz, Austria) equipped with a 25 mm diameter cone–plate geometry; the cone–plate gap was set to 0.103 mm. During each measurement, the shear rate was varied from 1 to 30 s⁻¹ over a period of 50 s to obtain flow curves.

Data were analysed using RheoCompass software (Anton Paar GmbH, Graz, Austria). After confirming that the samples exhibited Newtonian behaviour, the viscosity was determined as the average of the values obtained from 6 to 50 s during the measurement. Temperature control was achieved using a Peltier system with a circulating fluid bath, and all measurements were conducted at 20 °C.

3 Results and discussion

3.1 System equilibration

The time evolution of the densities for 8 : 0, 18 : 1, and 18 : 3 is shown in Fig. 2. The density of 18 : 1 continued to gradually increase until approximately 200 ns and then equilibrated, whereas the other systems stabilised earlier. Accordingly, only data collected after 300 ns were used in the subsequent analyses; only the first 300 ns of post-equilibration data are shown in the following figures, except for Fig. 3 and 9. Profiles of the potential energy, temperature, and pressure as functions of time are shown in the SI and confirm that all three systems reached equilibrium under the present conditions. The stabilised values are listed in Table 1, along with the other calculated data.

Fig. 2 Densities of 8 : 0 (red), 18 : 1 (blue), and 18 : 3 (green) as a function of time.

Fig. 3 Viscosities of 8 : 0 (red), 18 : 1 (blue), and 18 : 3 (green) as a function of time.

3.2 Viscosity and diffusion coefficients

Fig. 3 shows the viscosities calculated using eqn (2). The stabilised values of 8 : 0, 18 : 1, and 18 : 3 were approximately 100, 345, and 92 mPa s, respectively (Table 1; experimentally determined values are also shown in parentheses). The calculated viscosities are roughly four times larger than the experimental values, but the relative ordering among the three TGs is consistent with the experiment. This discrepancy is comparable to that reported in [ref. 27] investigating twelve TGs, where the calculated viscosities were several times higher than the experimental values. Among the TGs, 18 : 1 exhibited the highest viscosity, and the value of 8 : 0 was slightly larger than that of 18 : 3, suggesting that 18 : 1 molecules have the greatest resistance to motion.

Fig. 4 shows the MSD and the corresponding logarithmic plot for 8 : 0, 18 : 1, and 18 : 3. The MSDs were fitted in their linear regimes to obtain the diffusion coefficients, D^∞ (Table 1). As with the viscosities, the diffusion coefficients of 8 : 0 and 18 : 3 were nearly identical, whereas that of 18 : 1 was the lowest, further indicating that 18 : 1 molecules are the most hindered in their motion. Hydrodynamic radii R_h were calculated using the Stokes–Einstein relation and compared with the radii of gyration R_g, whose distribution is shown in Fig. 5. The geometric R_g of 8 : 0 was the smallest among the three TGs, whereas the R_g of 18 : 1 was nearly identical to that of 18 : 3, despite the presence of three double bonds in 18 : 3. Notably, 18 : 1 had the largest R_h, but 8 : 0 and 18 : 3 exhibited nearly the same value (differing by only 0.19 nm). This mismatch between R_g and R_h signifies that molecular size alone cannot account for the observed diffusion behaviour, necessitating energy and structural analyses.

Fig. 4 (a) MSD of 8 : 0 (red), 18 : 1 (blue), and 18 : 3 (green), and (b) the corresponding logarithmic plots as a function of time. The dashed lines represent linear fits obtained between 75 and 150 ns with a slope of 1.

Fig. 5 Probability distributions of the radius of gyration R_g for 8 : 0 (red), 18 : 1 (blue), and 18 : 3 (green).

3.3 Cohesive energy density and melting point expectations

We analysed the CED using eqn (4). As shown in Table 1, the CED values of 8 : 0, 18 : 1, and 18 : 3 differ by less than 2%. Similarly, the equilibrium densities are close in value. Using the Hansen solubility parameters listed in the Database of Safer Solvents⁴⁴ for 8 : 0, 18 : 1, and 18 : 3, the solubility parameters (SP) of each triglyceride TG were estimated to be 17.1, 16.8, and 17.1, respectively. Applying these SP values to the Hildebrand equation⁴⁵ to calculate the CED, the resulting values for 8 : 0, 18 : 1, and 18 : 3 were −292, −282, and −292 J cm⁻³, respectively. These values show good agreement with MD simulations, indicating that the three TGs have nearly identical CED values. Thus, bulk cohesion and packing density cannot explain the observed differences in diffusion.

Another widely used empirical descriptor is the melting temperature (T_m), which correlates with viscosity via Vogel–Fulcher–Tammann-type relations. A similar temperature dependence, scaled by T_m, is found in Andrade's model, an empirical relation widely used to describe viscosity as a function of temperature and to rationalise structure–property trends. In both cases, the viscosity is governed by how far the system is from the melting point. Here, we considered the melting points of the three TGs: 8 : 0 (T_m ≈ 9–10 °C),⁴⁶ 18 : 1 (T_m ≈ −5.5 °C),⁴⁷ and 18 : 3 (T_m ≈ −23–24 °C).⁴⁸ Based on these melting points, one would expect the viscosity to decrease in the order 8 : 0 > 18 : 1 > 18 : 3 because lower melting points generally correspond to more fluid systems. However, our experiments and simulations reveal that 18 : 1 exhibits a markedly higher viscosity than both 8 : 0 and 18 : 3, which have similar values, contradicting this empirical expectation. This discrepancy indicates that neither CED nor melting point arguments can account for the distinct dynamic behaviour.

3.4 Hierarchical structures of the TGs

3.4.1 Single-molecule conformations. The molecular conformation of the TGs plays a critical role in determining their properties. We classified TG conformers into three planar states based on the relative orientation of their three fatty acid chains. Specifically, when all three inter-chain angles were less than 90°, the conformation was assigned to State 1; when all three angles exceeded 90°, the conformation was assigned to State 2; and all other conformations were classified as State 3. Across all three species, State 2 was rare (<5%), whereas States 1 and 3 accounted for approximately 30% and 65% of the population, respectively, with no significant differences observed among the TGs. Notably, State 1 adopts a tripodal conformation, whereas State 3 exhibits a trident-like conformation, suggesting that TG molecules predominantly adopt asymmetric geometries. Supporting movie data for a single TG molecule over a 100-ns further simulation reveals that its conformation repeatedly changes over time, indicating that individual TGs do not remain in a single state. Taken together, these results demonstrate that the overall conformational distribution is essentially independent of TG type.

The three TGs studied here display similar overall molecular conformations (Fig. 6). Thus, the pronounced differences in viscosity and hydrodynamic radii among 8 : 0, 18 : 1, and 18 : 3 cannot be attributed to molecular shape effects.

Fig. 6 Probabilities of the conformational states for 8 : 0 (a), 18 : 1 (b), and 18 : 3 (c). Red, blue, and green denote States 1, 2, and 3, respectively. Representative molecular snapshots corresponding to each state are also shown.

3.4.2 Collective aggregation structure. We conducted a cluster analysis based on Voronoi tessellation to characterise the collective aggregation of the TGs. Fig. 7(a) shows a representative snapshot in which molecules belonging to the same cluster are coloured identically, with only the glycerol moieties shown for clarity. In the 8 : 0 system, a single large cluster was observed, whereas multiple smaller clusters (i.e., discrete aggregates) formed in the 18 : 1 and 18 : 3 systems (Fig. 7(b and c)). The size distributions of the aggregates exhibited an approximately exponential decay. Furthermore, the 18 : 3 system formed a greater number of smaller aggregates than the 18 : 1 system, likely because the kinks introduced by the three double bonds hinder tight packing of the hydrocarbon chains. These structural differences suggest that chain length and unsaturation influence mesoscale morphology.

Fig. 7 (a) Snapshots of the 8 : 0 (left), 18 : 1 (middle), and 18 : 3 (right) aggregates (fatty acid chains are omitted for clarity). (b) Number of 18 : 1 (blue) and 18 : 3 (green) aggregates. (c) Size distribution of the 18 : 1 (blue) and 18 : 3 (green) aggregates.

However, the cluster morphology of 18 : 1 resembled that of 18 : 3, despite their viscosities differing by a factor of four. Conversely, although the cluster morphology of 8 : 0 was markedly different from that of 18 : 3, their viscosities were comparable. Therefore, cluster structures alone cannot account for the viscosity differences among the TGs, and more local analyses are required.

3.4.3 Network structures via C-chain alignment. We analysed the local arrangements between C-chain units to identify potential network junctions. Fig. 8 shows the distributions of the distances between the centre atoms of the C-chain units and of the angles between pairs of C-chain units. For 8 : 0 and 18 : 1, pronounced peaks appear at r = 0.5 nm and cos θ = 1, corresponding to two C-chain units aligned in parallel. In contrast, only a weak peak was observed for 18 : 3. These locally aligned C-chain units act as junction points, which, when connected through the glycerol backbones, generate extended network structures in 8 : 0 and 18 : 1. Considering the number of C-chain segments, 8 : 0 contains three C-chain units per molecule, whereas 18 : 1 contains six. Note that 18 : 3 also has three C-chain units, but the large kinks introduced by the double bonds adjacent to these segments hinder their parallel packing, thereby preventing effective junction formation.

Fig. 8 Distributions of the C-chain alignment for 8 : 0 (a), 18 : 1 (b), and 18 : 3 (c). The right panels show representative snapshots at r = 0.5 nm and cos θ = 1, corresponding to different contact pairs.

Our findings suggest that network formation is responsible for the markedly reduced diffusion—higher viscosity—observed in 18 : 1. Although the compact geometry of 8 : 0 would typically imply a much higher diffusion coefficient, its mobility is comparable to that of 18 : 3 because of network formation. Quantitatively, the peak intensities in Fig. 8 indicate that the density of aligned junctions is highest in 18 : 1, moderate in 8 : 0, and negligible in 18 : 3, correlating with the observed order of viscosity.

Fig. 9 presents the bond correlation function B(t), which characterises the lifetime of the junction structures. For 8 : 0 and 18 : 3, B(t) reached zero at roughly 100 ns, indicating that the local alignments between C-chain units are relatively short-lived compared with those in 18 : 1. Conversely, for 18 : 1, a significant fraction of correlations persisted at 100 ns, demonstrating that C-chain–C-chain alignments are more long-lived. These results support the view that the anomalously high viscosity of 18 : 1 arises not only from the formation of network junctions but also from their enhanced temporal stability.

Fig. 9 Bond correlation function of 8 : 0 (red), 18 : 1 (blue), and 18 : 3 (green) as a function of time.

Furthermore, we calculated the relaxation times of B(t) for 8 : 0, 18 : 1, and 18 : 3 by fitting the data to a double exponential function: , where T₁ and T₂ represent the relaxation times associated with fast vibrational bond fluctuations and long-lived C-chain alignment, respectively. T₁ was found to be 2.5 ps for all three TGs, whereas the T₂ values for 8 : 0, 18 : 1, and 18 : 3 were determined to be 3, 14, and 5 ns, respectively. From the perspective of the star polymer relaxation analogy,⁴⁹ each TG molecule contains three fatty acid chains, and for the entire molecule to relax and diffuse, all three chains must relax first (Fig. 10). Thus, the presence of aligned C-chains delays molecular relaxation and leads to increased viscosity.

Fig. 10 Schematic illustration of TG diffusion mechanics based on the star polymer relaxation analogy.

In summary, the viscosities of TGs can be rationalised by the emergence—or absence—of network structures mediated by the parallel alignment of C-chain units. The ease of forming such networks is determined by the balance between the number of available C-chain units and the disruption due to kinks introduced by double bonds.

4 Conclusion

We employed all-atom MD simulations to elucidate the diffusion and viscosity of three homologous TGs—trioctanoin (8 : 0), triolein (18 : 1), and trilinolenin (18 : 3). The simulations quantitatively reproduced the ordering; 18 : 1 was the most viscous, whereas 8 : 0 and 18 : 3 exhibited nearly identical viscosities.

Bulk thermodynamic descriptors, including CED, density, and melting point correlations of fatty acids, could not explain the viscosity trends among the three TGs. Similarly, single-molecule conformations were essentially comparable across the three systems and thus cannot account for the observed differences. Mesoscale cluster morphologies also failed to clarify the trend; although 18 : 1 and 18 : 3 exhibited similar aggregation patterns, their viscosities differed by a factor of four, and 8 : 0 displayed a distinct morphology but a comparable viscosity to 18 : 3. The decisive factor was the presence or absence of locally aligned C-chain units that acted as network junctions. Such junctions were abundant in 18 : 1, intermediate in 8 : 0, and negligible in 18 : 3. Consequently, the viscosity of 18 : 1 was markedly larger, and that of 8 : 0—which was expected to be much lower owing to the small molecular size—was elevated by network formation, yielding a value comparable to or slightly higher than that of 18 : 3.

These findings demonstrate that fine-scale structural motifs, rather than bulk thermodynamics or molecular size alone, govern the dynamic properties of TGs. The molecular-level framework developed herein provides new guidelines for tailoring the flow behaviour of oils and lubricants through the design of chain length, unsaturation, and alignment propensity.

Author contributions

Yuki Kitamura: formal analysis, investigation, visualisation, writing – original draft. Yusuke Yasuda: conceptualisation, writing – review & editing. Junya Metoki: investigation, data curation, validation, writing – review & editing. Shogo Tsujino: validation, writing – review & editing. Kazuma Yoshimura: validation, writing – review & editing. Takahide Watanabe: validation, writing – review & editing. Kazushi Fujimoto: conceptualisation, supervision, project administration, funding acquisition, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data supporting the findings of this study are available within the article and its supporting information (SI). See DOI: https://doi.org/10.1039/d5ra07973j.

Acknowledgements

This work was supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) through the Program for Promoting Research on the Supercomputer Fugaku (“Computational Research on Materials with Better Functions and Durability Toward Sustainable Development,” Grant No. JPMXP1020230325). We used computational resources (Fugaku supercomputer) provided by the RIKEN Center for Computational Science (Project IDs: hp230205, hp240224, hp250227). Calculations were also performed using supercomputers at the Research Center for Computational Science in Okazaki, Japan (Projects: 24-IMS-C502, 24-IMS-C602, 25-IMS-C118, 25-IMS-C502). We also thank Robert Ireland, PhD, from Edanz (https://jp.edanz.com/ac) for editing a draft of this manuscript.

References

1. K. D. Maher and D. C. Bressler, Bioresour. Technol., 2007, 98, 2351–2368.
2. Y. Zhang, X. Wang, Q. Li, R. Yang and C. Li, Energy Fuels, 2015, 29, 5056–5068.
3. Z. Zhang, K. Yan and J. Zhang, J. Mol. Model., 2014, 20, 1–9.
4. N. H. Jayadas, K. Prabhakaran Nair and G. Ajithkumar, Tribol. Int., 2007, 40, 350–354.
5. A. Adhvaryu and S. Z. Erhan, Ind. Crops Prod., 2002, 15, 247–254.
6. M. P. Schneider, J. Sci. Food Agric., 2006, 86, 1769–1780.
7. W.-D. Hsu and A. Violi, J. Phys. Chem. B, 2009, 113, 887–893.
8. P. Lonchampt and R. W. Hartel, Eur. J. Lipid Sci. Technol., 2004, 106, 241–274.
9. S. S. Rane and B. D. Anderson, Adv. Drug Delivery Rev., 2008, 60, 638–656.
10. A. H. Rantamäki, T. Seppänen-Laakso, M. Oresic, M. Jauhiainen and J. M. Holopainen, PLoS One, 2011, 6, e19553.
11. S. Martin and R. G. Parton, Nat. Rev. Mol. Cell Biol., 2006, 7, 373–378.
12. A. Pappas, Derm.-Endocrinol., 2009, 1, 72–76.
13. E. R. S. U.S. Department of Agriculture, Oil Crops Outlook: May 2025, U.S. Department of Agriculture, Economic Research Service, 2025.
14. S. N. Sahasrabudhe, V. Rodriguez-Martinez, M. O'Meara and B. E. Farkas, Int. J. Food Prop., 2017, 20, 1965–1981.
15. M. Flores, V. Avendaño, J. Bravo, C. Valdés, O. Forero-Doria, V. Quitral, Y. Vilcanqui and J. Ortiz-Viedma, Int. J. Food Sci., 2021, 2021, 7105170.
16. H. Noureddini, B. C. Teoh and L. Davis Clements, J. Am. Oil Chem. Soc., 1992, 69, 1189–1191.
17. A. Hall, J. Repakova and I. Vattulainen, J. Phys. Chem. B, 2008, 112, 13772–13782.
18. A. Brasiello, S. Crescitelli and G. Milano, Phys. Chem. Chem. Phys., 2011, 13, 16618–16628.
19. A. Pizzirusso, A. Brasiello, A. De Nicola, A. G. Marangoni and G. Milano, J. Phys. D: Appl. Phys., 2015, 48, 494004.
20. A. S. Tascini, M. G. Noro, R. Chen, J. M. Seddon and F. Bresme, Phys. Chem. Chem. Phys., 2018, 20, 1848–1860.
21. P. Campomanes, J. Prabhu, V. Zoni and S. Vanni, Biophys. Rep., 2021, 1, 100034.
22. J. Telenius, A. Koivuniemi, P. Kulovesi, J. M. Holopainen and I. Vattulainen, Langmuir, 2012, 28, 17092–17100.
23. A. Bacle, R. Gautier, C. L. Jackson, P. F. J. Fuchs and S. Vanni, Biophys. J., 2017, 112, 1417–1430.
24. G. Henneré, P. Prognon, F. Brion and I. Nicolis, Chem. Phys. Lipids, 2009, 157, 86–93.
25. O. H. S. Ollila, A. Lamberg, M. Lehtivaara, A. Koivuniemi and I. Vattulainen, Biophys. J., 2012, 103, 1236–1244.
26. I. Chandrasekhar and W. F. van Gunsteren, Eur. Biophys. J., 2002, 31, 89–101.
27. A. K. Sum, M. J. Biddy, J. J. de Pablo and M. J. Tupy, J. Phys. Chem. B, 2003, 107, 14443–14451.
28. M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, publisher = Oxford University Press, 2023.
29. I.-C. Yeh and G. Hummer, J. Phys. Chem. B, 2004, 108, 15873–15879.
30. Y. Zhang, A. Otani and E. J. Maginn, J. Chem. Theory Comput., 2015, 11, 3537–3546.
31. M. Wu, G. Xu, Y. Luan, Y. Zhu, T. Ma and W. Zhang, Constr. Build. Mater., 2022, 333, 127403.
32. J. K. Maranas, M. Mondello, G. S. Grest, S. K. Kumar, P. G. Debenedetti and W. W. Graessley, Macromolecules, 1998, 31, 6991–6997.
33. S. Kawada, M. Komori, K. Fujimoto, N. Yoshii and S. Okazaki, Chem. Phys. Lett., 2016, 646, 36–40.
34. Y. Ito, Y. Yasuda, Y. konishi, T. Shiba, K. Dokko and K. Fujimoto, ChemRxiv, 2025,  DOI:10.26434/chemrxiv-2025-1gkf0.
35. C. B. Barber, D. P. Dobkin and H. Huhdanpaa, ACM Trans. Math. Softw., 1996, 22, 469–483.
36. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
37. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
38. J. M. Martínez and L. Martínez, J. Comput. Chem., 2003, 24, 819–825.
39. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
40. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126.
41. M. Parrinello and A. Rahman, Phys. Rev. Lett., 1980, 45, 1196–1199.
42. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
43. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
44. Database of Safer Solvents, https://doss.turi.org/, accessed October 16, 2025.
45. W. Zeng, Y. Du, Y. Xue and H. L. Frisch, in Physical Properties of Polymers Handbook, ed. J. E. Mark, Springer New York, New York, NY, 2007, pp. 289–303,  DOI:10.1007/978-0-387-69002-5_16.
46. ChemicalBook, https://www.chemicalbook.com/ChemicalProductProperty_EN_CB8170985.htm, accessed October 6, 2025.
47. ChemicalBook, https://www.chemicalbook.com/ChemicalProductProperty_EN_CB2715065.htm, accessed October 6, 2025.
48. ChemicalBook, https://www.chemicalbook.com/ChemicalProductProperty_US_CB5178728.aspx, accessed October 6, 2025.
49. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, 2003.

---