NIST Technical Note 1805

Mass diffusion of organic fluids: a 

molecular dynamics perspective 

Alex Smolyanitsky

Andrei F. Kazakov

Thomas J. Bruno

Marcia L. Huber
NIST Technical Note 1805 

Mass diffusion of organic fluids: a 

molecular dynamics perspective 

Alex Smolyanitsky

Andrei F. Kazakov

 Thomas J. Bruno

Marcia L. Huber

Applied Chemicals and Materials Division 

Materials Measurement Laboratory 

May 2013

U.S. Department of Commerce 

Rebecca Blank, Acting Secretary

National Institute of Standards and Technology 

Patrick D. Gallagher, Under Secretary of Commerce for Standards and Technology and Director 
Certain commercial entities, equipment, or materials may be identified in this

document in order to describe an experimental procedure or concept adequately.

Such identification is not intended to imply recommendation or endorsement by the

National Institute of Standards and Technology, nor is it intended to imply that the

entities, materials, or equipment are necessarily the best available for the purpose. 

National Institute of Standards and Technology Technical Note 1805 

Natl. Inst. Stand. Technol. Tech. Note 1805, 16 pages (May 2013) 

CODEN: NTNOEF
iii

Mass diffusion of organic fluids: a molecular dynamics perspective 

Abstract

The Optimized Potential for Liquid Simulation All-Atom (OPLS-AA) forcefield, a wellestablished interaction potential, within the MD framework was used to determine self- and 

mutual diffusivity of several fluids, in the high-pressure and near- and supercritical regimes. The 

test fluids were benzene, heptane and carbon dioxide. The Green-Kubo method of calculating 

mass diffusivity yielded considerably overestimated values of self-diffusivity; better results were 

obtained with use of the Stokes-Einstein method. In addition, a number of relatively significant 

problems with the default OPLS-AA parameterization were identified. There exist ranges of 

pressure and temperatures, in which the results for the selected pure fluids show good agreement 

with experiment, and thus reparameterization of OPLS-AA may not be required. However, the 

results for mixtures all show significant deviation from experiment. Based on the observed 

results, the current predictive power of OPLS-AA for mass diffusivity of the fluids in question is 

not satisfactory. We are therefore convinced that a systematic reparameterization of OPLS-AA 

both in terms of the intermolecular (van der Waals, electrostatic) and intramolecular (bonding) 

interactions is needed to increase its accuracy. In addition, OPLS-AA may not be generally 

appropriate for the gaseous state and therefore near-critical systems may benefit from using a

different forcefield.
1

1 Introduction

Knowledge of the fluid transport properties is crucial for our understanding of the way a 

pure fluid or a mixture responds to changes in pressure, or temperature. Of great 

importance in the design and optimization of mass transfer processes is the diffusion 

coefficient. A recent application of particular interest is the concept of clean diesel 

combustion that incorporates the use of supercritical carbon dioxide to reduce emissions 

and improve engine efficiency. Accurate values of the diffusion coefficients of 

hydrocarbons in supercritical carbon dioxide are essential to the understanding and 

optimal design of this process.

Determining fluid diffusion properties in the supercritical regime is challenging from the 

experimental standpoint [1-4], which often reduces the practicality of direct 

measurement. Current analytical calculations based on kinetic theory can provide an 

adequate description of the self-diffusivity of pure fluids, but yield significantly less 

accurate results for mixtures [5]. As a result, molecular dynamics (MD) simulations have 

been utilized to determine the mutual mass diffusivity of mixtures (as well as pure fluids) 

with relative success [5-7].

The aim of this investigation is to determine the feasibility of using one of the most 

established interaction potentials, the Optimized Potential for Liquid Simulation AllAtom (OPLS-AA) force-field [8] within the MD framework, to determine self- and 

mutual diffusivity of near-critical and supercritical carbon dioxide, as well as organic 

fuels, at high pressure. For test fluids we chose benzene, heptane, and carbon dioxide due 
to their relevance as representative fluids for clean diesel combustion as well as the 

availability of some experimental data to validate our calculations. We report the 

calculated mass density and diffusion data obtained with OPLS-AA, as a function of 

pressure and temperature. We discuss the current performance of OPLS-AA in terms of 

accuracy of calculations, as well as possible modifications necessary to obtain better 

agreement with known experimental data, and therefore increase the general predictive 

power of the MD-based approach.

2 Description of the simulated system

The fluid systems simulated here represented carbon dioxide, benzene, and heptane; each 

pure fluid system consisted of a total of 4000 molecules, unless stated otherwise. All 

mixtures studied here contained 4000 molecules of solvent and 40 molecules of solute, in 

order to mimic conditions close to infinite dilution. The simulations were performed with 

use of the GROMACS 4.5.6 package [9] configured with double precision. The pressures 

and temperatures were appropriately maintained by baro- and thermostatics, respectively

[9]. Each system was pre-relaxed with temperature and pressure control in a two-step 

procedure for a total of 1 ns, followed by a production simulation. The production 

simulations of pure fluids were performed at constant volume and temperature and run 

for a total of 5 ns; simulations of mixtures were run for 15 ns. The simulated time 

intervals were selected to ensure overall numerical convergence and thus statistical 

significance of the calculated properties. All statistical averages were calculated in the 

production simulations during the second half of the simulated time interval. 
The interatomic interactions for heptane and benzene were described by the 

corresponding default OPLS-AA models [8]. Carbon dioxide was represented by a rigid 

linear system with two virtual atomic sites of appropriate mass, two oxygen atoms, and 

one carbon atom. All atoms were fully represented in terms of the van der Waals and 

electrostatic interactions [9]. The electrostatic interactions were calculated from the 

particle-particle sum and a long-range particle-mesh correction obtained from the 

appropriate Ewald sum [9]. The van der Waals and short-range electrostatic cut-off radius 

used here was 1 nm. Within the OPLS-AA formulation, the non-bonded short-range van 

der Waals interactions were represented by the Lennard-Jones pairwise potential:

φ(r) = 4ε ��

σ

r

�

12

− �

σ

r

�

6

�, (1)

where ε and σ are the depth and position of the pairwise energy well. The Lennard-Jones 

parameters along with the partial electric charges q are provided in Table 1, and the interspecies mixing rules were implemented according to OPLS-AA, as follows: σAB =

√σAσB. εAB = √εAεB.

ε (K) σ (Å) Q (|e|)

Carbon dioxide 

(CO2)

C-C: 28.13

O-O: 80.51

C-C: 2.76

O-O: 3.03

C: +0.652

H: -0.326

Benzene (C6H6) C-C: 35.24

H-H: 15.10

C-C: 3.55

H-H: 2.42

C: +0.115

H: -0.115

Heptane (C7H16) C-C: 33.23

H-H: 15.10

C-C: 3.50

H-H: 2.50

C (terminus): -0.18

C (backbone): -0.12

H: +0.06

Table 1. Lennard-Jones parameters and partial electric charges used in the simulations. 

The carbon dioxide parameters are from literature [10, 11], while the parameters for 

benzene and heptane are the OPLS-AA’s defaults for the corresponding compound.
The mass density was calculated directly from the ratio of the total mass of the system to

its volume. The self-diffusion coefficients for the pure fluids were calculated with use of

the Stokes-Einstein relation:

DA = 1

6τ 〈∑ ∆r2〉|τ, (2)

where ∆r is the molecular displacement and the averaging is calculated for all molecules 

of species A, integrated over the period τ. For mixtures, the corresponding expression is 

more complex and is provided elsewhere [6]. However, for the case of infinite dilution, 

the expression for mutual mass diffusivity conveniently reduces to the solute’s selfdiffusivity. Note that for numerical comparison in some cases, we also calculated the 

diffusion coefficients with use of the Green-Kubo method [5].

3 Results

Carbon dioxide

Shown in Fig. 1 is the carbon dioxide mass density data obtained for various 

temperatures, compared with the NIST Refprop [12] data. Good agreement is obtained 

for the sub-critical regime (< 50 bar) at T = 298 K, while a significant discrepancy (25 %) 

is observed for the supercritical point at ≈ 140 bar. For the temperatures of 323 K and 348 

K, the agreement between Refprop data and the MD results also decreases with 

increasing pressure. The observed discrepancies strongly suggest that the compressibility 

of carbon dioxide with the current van der Waals and electrostatic interaction parameters 

is tailored within default OPLS-AA toward lower pressures simulated here. 
Figure 1. Carbon dioxide mass density results obtained with OPLS-AA force-field for the 

temperatures of 298 K (a), 323 K (b), and 348 K (c), compared with NIST Refprop data. 

As we see further, the effect of incorrect compressibility on the carbon dioxide’s selfdiffusivity is not significant. Shown in Fig. 2 is the mass self-diffusivity obtained for the 

same set of temperatures, compared with experimental data [4].

Figure 2. Carbon dioxide mass self-diffusivity obtained with OPLS-AA force-field for 

the temperatures of 298 K (a), 323 K (b), and 348 K (c), compared with experimental 

data [4]. 

As we can see from Fig. 2, the self-diffusivity results obtained with the Stokes-Einstein 

expression are in good agreement with the experimental data for all temperatures, 

indicating that an incorrect compressibility at high pressures does not significantly affect 

the simulated mass diffusion for the carbon dioxide molecule. Moreover, the agreement 

with experiment is considerably better, compared with other interatomic potentials [13]. 
For all temperatures in Fig. 2, the Green-Kubo method yields considerably overestimated 

values of self-diffusivity. Such discrepancies, resulting from general numerical 

convergence issues, are not surprising and have been reported previously [14]. Because 

of this discrepancy, the diffusivities further in the text were only calculated with use of 

Stokes-Einstein method.

Benzene

Shown in Fig. 3 is the benzene mass density data obtained for the temperatures of 298 K 

and 328 K, compared with the NIST Refprop [12] data. We see that for T = 298 K the 

simulated density agrees with Refprop within 5 %. For T = 328 K, the difference between 

MD and Refprop at low pressures is approximately 6 %, improving somewhat toward 

higher pressures. Such behavior may indicate a slight offset in the van der Waals, as well 

as electrostatic parameters of OPLS-AA; also, the internal rigidity of each benzene ring 

may be overestimated, which affects intermolecular packing.

Figure 3. Benzene mass density results obtained with OPLS-AA force-field for the 

temperatures of 298 K (a), 328 K (b), compared with NIST Refprop data. 
In Fig. 4 we report the benzene self-diffusion results, obtained for the same set of 

temperatures, compared with experimental data [1]. For T = 298 K, the average 

discrepancy with experiment is 20 %, while for T = 328 K, the discrepancy steadily 

decreases from 25 % to 3 % at the lowest and the highest pressure, respectively.

Figure 4. Benzene mass self-diffusivity obtained with OPLS-AA force-field for the 

temperatures of 298 K (a) and 328 K (b), compared with experimental data [1]. 

This result indicates that for the model of benzene, both the compressibility and mass 

diffusivity are affected by a relatively severe lack of proper parameterization in OPLSAA. It also appears that the discrepancy between calculated and experimentally obtained 

self-diffusivity decreases strongly with increasing pressure.

Heptane

Shown in Fig. 5 are the results for mass density obtained for pure heptane at temperatures 

of 293 K and 353 K. The results show considerable discrepancy from Refprop and 

experimental data at lower pressures, with 10 % to 5 % deviation at P = 1 bar for 353 K 

and 293 K, respectively. Toward higher pressures, the agreement improves and is on 
average 2 % to 3 %. These results suggest that the default OPLS-AA parameterization for 

heptane (at the temperatures considered) is tailored toward pressures significantly above 

atmospheric.

Figure 5. Heptane mass density results obtained with OPLS-AA force-field for the 

temperatures of 298 K (a) and 353 K (b), compared with NIST Refprop data and 

experimental values [15]. 

Shown in Fig. 6 are the self-diffusivity results for heptane, compared with experimental 

data [16]. Similar to the results for benzene in Fig. 4 (b), for both temperatures 

considered, the discrepancy with experiment decreases from 25 % at the atmospheric 

pressure to 2 % at the highest pressure of ≈ 2450 bar. It appears that for the highly 

flexible heptane molecule, the default OPLS-AA representation provides an excess of 

bending and rotational degrees of freedom, resulting in overestimated self-diffusivity. It 

is therefore not surprising that at high pressures these shortcomings of the OPLS-AA 

default parameterization become negligible and the self-diffusivity approaches 

experimental values. A comprehensive reparameterization of OPLS-AA energies 

associated with molecule bending and rotation may solve the observed issue. At the same 
time, if one is interested in the self-diffusivity of heptane for pressures above 1000 bar, 

the current OPLS-AA performance is already within 10 % of the experimental data.

Figure 6. Heptane mass self-diffusivity obtained with OPLS-AA force-field for the 

temperatures of 293 K (a) and 353 K (b), compared with experimental data [16]. 

Binary mixtures

Shown in Fig. 7 are the mutual diffusivity data of benzene and heptane infinitely diluted 

in supercritical carbon dioxide, compared with available experimental data. The relatively 

narrow range of pressures selected for simulations was due to the fact that the density of 

carbon dioxide, which dominates the environment of the solute in this case, agreed with 

experiment within 15 % only for pressure values below 120 bar (Fig. 1). As we can see 

from Fig. 7, despite the carefully selected range of simulated pressures, the mutual 

diffusivity for both benzene and heptanes is considerably overestimated.
Figure 7. Mutual diffusivity obtained with OPLS-AA force-field for benzene (a) and 

heptane (b) infinitely diluted in supercritical carbon dioxide. The experimental data for 

benzene and heptane are from [17] and [18], respectively.

We believe the discrepancy observed here is due to a combination of the underestimated 

density of the solvent, together with incorrectly described internal degrees of freedom of 

the solute molecules, which results in greatly overestimated values of mutual diffusivity.

While the same phenomenon does not strongly affect the self-diffusivity of carbon 

dioxide (although its self-diffusivity was also somewhat overestimated; see Fig. 2), the 

effect for more complex (in terms of the number of degrees of freedom) solute molecule, 

such as benzene and heptane is significant, resulting in overestimation of mutual 

diffusivity by at least 100 %, compared to experiment. For isotherms at higher pressures, 

we expect the discrepancy to grow further. 

4 Conclusions, recommendations, and future work

In this study, we performed a variety of MD simulations with use of OPLS-AA forcefield

and reported on the accuracy of density and mass diffusivity for pure near-critical and 
supercritical fluid carbon dioxide, high-pressure benzene, and heptane, as well as (carbon 

dioxide – benzene) and (carbon dioxide – heptane) mixtures. 

A number of significant problems with the default OPLS-AA parameterization were

identified. For carbon dioxide, the mass density is strongly underestimated at higher 

pressures, indicating an incorrect high-pressure compressibility. At the same time, the 

self-diffusivity of carbon dioxide is calculated within 5 % to 10 % of the established 

experimental data, which indicates that the issues with compressibility here do not 

significantly affect self-diffusion. For benzene and heptane, it was observed that the mass 

densities are consistently underestimated at pressures close to atmospheric (reaching a 

discrepancy with experiment of at least ≈ 6 %), while the self-diffusion data is 

overestimated at lower pressures, compared to experiment. We believe that the observed 

discrepancies arise from an incorrect default OPLS-AA parameterization; thus, a 

systematic reparameterization could be of great value. For mixtures, we observe a 

definite combined effect, in which the mutual diffusivities are overestimated by as high 

as 100 %. We believe that is due to an incorrect compressibility of carbon dioxide in 

OPLS-AA, combined with the lack of accuracy of OPLS-AA in describing the internal 

degrees of freedom for benzene and heptane. Finally, because OPLS-AA is designed for 

the simulation of liquids, it may not be generally appropriate for simulating systems in 

the selected states, and thus other molecular dynamics forcefields could be used.

Given these observations, great care must be exercised when simulating self- and mutual 

diffusivity of organic fuels, using the default OPLS-AA parameterization. There exist 

ranges of pressure and temperatures, in which the results for pure fluids show good 
agreement with experiment, and thus reparameterization of OPLS-AA may not be 

required (as seen, for example, for heptane in Fig. 6). However, all results for mixtures 

show significant deviation from experiment. Based on the observed results, the current 

predictive power of OPLS-AA for mass diffusivity of organic fluids in selected states is 

not satisfactory. We therefore recommend that a systematic reparameterization of OPLSAA, both in terms of the intermolecular (van der Waals, electrostatic) and intramolecular 

(bonding) interactions, should be undertaken in order to increase the accuracy of OPLSAA representation of fluid transport properties calculated within the molecular dynamics 

framework.

References

1. McCool, M.A., A.F. Collings, and L.A. Woolf, Pressure and temperature 

dependence of the self-diffusion of benzene. J. Chem. Soc., Faraday Trans., 1972. 

1(68): p. 1489-1497.

2. Polzin, B. and A. Weiss, Transport Properties of Liquids. VIII. Molar Volume and 

Selfdiffusion of Organic Liquids at Pressures up to 200 MPa. Berichte der 

Bunsengesellschaft für physikalische Chemie, 1990. 94(7): p. 746-758.

3. Moore, J.W. and R.M. Wellek, Diffusion Coefficients of n-Heptane and n-Decane 

in n-Alkanes and n-Alcohols at Several Temperatures. J. Chem. Eng. Data, 1974. 

19(2): p. 136-140.

4. Etesse, P., J.A. Zega, and R. Kobayashi, High pressure nuclear magnetic 

resonance measurement of spin-lattice relaxation and self-diffusion in carbon 

dioxide. Journal of Chemical Physics, 1992. 97(3): p. 2022-2029.

5. Chae, K. and A. Violi, Mutual diffusion coefficients of heptane isomers in 

nitrogen: A molecular dynamics study. Journal of Chemical Physics, 2011. 134.

6. Zhang, L., Y.-C. Liu, and Q. Wang, Molecular dynamics simulation of self- and 

mutual diffusion coefficients for confined mixtures. Journal of Chemical Physics, 

2005. 123.

7. Zabaloy, M.S., V.R. Vasquez, and E.A. Macedo, Description of self-diffusion 

coefficients of gases, liquids and fluids at high pressure based on molecular 

simulation data. Fluid Phase Equilibria, 2006. 242: p. 43-56.

8. Jorgensen, W.L., D.S. Maxwell, and J. Tirado-Rives, Development and Testing of 

the OPLS All-Atom Force Field on Conformational Energetics and Properties of 

Organic Liquids. J. Am. Chem. Soc., 1996. 118(45): p. 11225–11236.

9. Apol, E., et al., GROMACS 4.5.6 manual, 2010.
10. Chatzis, G. and J. Samios, Binary mixtures of supercritical carbon dioxide with 

methanol. A molecular dynamics simulation study. Chemical Physics Letters, 

2003. 374(1): p. 187-193.

11. Yoo, J.-H., et al., Diffusion coefficients of supercritical carbon dioxide and its 

mixtures using molecular dynamic simulations. Korean Journal of Chemical 

Engineering, 2012. 29(7): p. 935-940.

12. Lemmon, E.W., M.L. Huber, and M.O. McLinden, REFPROP, NIST Reference 

Fluid Thermodynamic and Transport Properties Database 23, Version 9.0, 2010, 

National Institute of Standards and Technology: Gaithersburg, MD.

13. Coelho, L.A.F., J.V. Oliveira, and F.W. Tavares, Dense fluid self-diffusion 

coefficient calculations using perturbation theory and molecular dynamics.

Brazilian Journal of Chemical Engineering, 1999. 16(3).

14. Wei-Zhong, L., C. Cong, and Y. Jian, Molecular Dynamics Simulation of SelfDiffusion Coefficient and Its Relation with Temperature Using Simple LennardJones Potential. Heat Transfer - Asian Research, 2008. 37(2): p. 86-93.

15. Sagdeev, D.I., et al., Experimental Study of the Density and Viscosity of nHeptane at Temperatures from 298 K to 470 K and Pressure up to 245 MPa.

International Journal of Thermophysics, 2013. 34(1): p. 1-33.

16. Arkhipov, V.P., Experimental study of the diffusion coefficient in liquid nparaffins at high pressures. Izvestiia vysshikh uchebnykh zavedenii. Neft' i gaz, 

1982: p. 34-44.

17. Funazukuri, T. and N. Nishimoto, Tracer diffusion coefficients of benzene in 

dense CO2 at 313.2 K and 8.5-30MPa. Fluid Phase Equilibria, 1996. 125(1-2): p. 

235-243.

18. Umezawa, S. and A. Nagashima, Measurement of the Diffusion Coefficients of 

Acetone, Benzene, and Alkane in Supercritical CO2 by the Taylor Dispersion 

Method. The Journal of Supercritical Fluids, 1992. 5(4): p. 242-250.