Brown’s characteristic curves deﬁne lines on the thermodynamic surface where special thermodynamic conditions hold. These curves are an important tool for the development of thermodynamic models of ﬂuids. Yet, practically no experimental data for Brown’s characteristic curves is available. In this work, a rigorous and generalized method for determining Brown’s characteristic curves based on molecular simulation was developed. As multiple thermodynamic equivalent deﬁnitions apply for the char- acteristic curves, diﬀerent simulation routes were compared. Based on this system- atic approach, the most favorable route for determining each characteristic curve was identiﬁed. The computational procedure developed in this work combines molecu- lar simulation, molecular-based equation of state, and the evaluation of the second virial coeﬃcient. The new method was tested on a simple model system (the clas- sical Lennard-Jones ﬂuid) and diﬀerent types of real substances (toluene, methane, ethane, propane, and ethanol). It is thereby shown that the method is robust and yields accurate results. Moreover, a computer code implementation of the method is presented. 1 Introduction Brown’s characteristic curves deﬁne lines on the thermodynamic equilibrium surface. Each of the characteristic curves has – for a given thermophysical property – the same behavior as the ideal gas.1The four characteristic curves are the Zeno, Amagat, Boyle, and Charles curve. These curves are located within a large pressure and temperature range. For a molecular ﬂuid, Brown’s characteristic curves are known to exhibit certain thermodynamic features.1 Therefore, Brown’s characteristic curves have become an important tool for the assessment of the extrapolation behavior of new equations of state (EOS).2–4This test has also been included in the IUPAC guidelines for publishing new EOS.5Nevertheless, due to the extreme conditions, there is practically no experimental data available on the characteristic curves. Molecular simulation based on force ﬁelds provides an attractive alternative to classical laboratory experiments, which is explored in this work. The four characteristic curves are also known by other names. The Zeno curve is at times simply called ideal curve . The Amagat curve is also called Joule inversion curve , and the Charles curve is often called Joule-Thomson inversion curve . In this work, the original names introduced by Brown1are used. Brown’s characteristic curves are usually discussed in the double logarithmic pressure- temperature projection as schematically shown in Fig. 1. In this representation, all four characteristic curves exhibit a concave dome shape. Thereby, each characteristic curve ex- hibits a single pressure maximum. Moreover, Brown postulated three speciﬁc points where the characteristic curves contact each other:1the Zeno and Charles curve intersect each other in the pressure maximum of the Zeno curve, the Zeno and Boyle curve converge into each other in their zero-density limit →0at the Boyle temperature (zero-crossing point of second virial coeﬃcient), and the Zeno and Amagat curve intersect each other at low temperatures. The Zeno and Amagat curve are usually truncated by the solid-liquid equilibrium at low temperatures. The Amagat, Boyle, and Charles curve must not have an intersection point. Instead, the Amagat curve surrounds the Charles curve and the Charles curve surrounds 2 the Boyle curve in the p-Tprojection, cf. Fig. 1. The Boyle and Charles curve originate on the vapor pressure curve. More speciﬁcally, the Charles curve originates on the binodal and the Boyle curve on the liquid-side spinodal. Moreover, all four characteristic curves circulate around the critical point in the p-Tprojection – even though the Boyle curve passes the critical point in close proximity, cf. Fig. 1. Furthermore, the temperature of the zero-density limit →0of the characteristic curves is directly linked to speciﬁc points of the second virial coeﬃcient, cf. Fig. 1. Brown deduced based on rational thermodynamics that one 0th-order and three 1st-order characteristic curves exist for a simple molecular ﬂuid. For state points on the 0th-order curve, the compressibility factor Zis unity, i.e. Z=pv TR=1; (1) where Rindicates the universal gas constant and the speciﬁc volume is v=1/slash.left. The 0th-order characteristic curve is called Zeno curve. For state points on the three 1st-order characteristic curves, the temperature, pressure, and density derivative of the compressibility factor, respectively, correspond to the ideal gas values, i.e. are zero.1For the Amagat curve, the temperature derivative of Zat isochoric conditions is zero /parenleft.alt3@Z @T/parenright.alt3 v=0: (2) For the Boyle curve, the pressure derivative of Zat isothermal conditions is zero /parenleft.alt3@Z @p/parenright.alt3 T=0: (3) For the Charles curve, the speciﬁc volume derivative of Zat isobaric conditions is zero /parenleft.alt3@Z @v/parenright.alt3 p=0: (4) 3 Several thermodynamic equivalent conditions formulated in other variables can be derived from the conditions (1)-(4).3,6Table 1 gives an overview of the most relevant conditions based on the molar enthalpy h=H/slash.leftN, the molar internal energy u=U/slash.leftN, and molar Helmholtz energy a=A/slash.leftN. For the Helmholtz energy derivatives, the following short-hand notation is used ~akl=(1/slash.leftT)kl@k+l~a @(1/slash.leftT)k@l; (5) where ~a=a/slash.left(kBT)indicates the reduced Helmholtz energy with kBbeing the Boltzmann constant. For the evaluation of the characteristic curves from an algebraic thermodynamic model, e.g. an EOS, the diﬀerent conditions (cf. Table 1) are practically equivalent. This does however not hold for molecular simulations, which are computationally and numerically more demanding. For evaluating the diﬀerent thermodynamic conditions (cf. Table 1) from a molecular simulation, diﬀerent algorithms and simulation methods have to be applied for the sampling that go hand in hand with diﬀerent statistical and systematic uncertainties as well as required computer time. There is presently no information available on the advantages and disadvantages of the diﬀerent thermodynamic sampling routes for determining characteristic curve state points, cf. Table 1. The Charles curve, a.k.a. Joule-Thomson inversion curve, has been at times studied using molecular simulations.7–13Molecular simulation data for all four characteristic curves has been reported only rarely and only for simple models.6,14In most cases, no statistical uncertainties were reported. In this work, we developed a generalized approach for determin- ing Brown’s characteristic curves of a given force ﬁeld of a molecular ﬂuid. We demonstrate the general applicability of the computational procedure by determining Brown’s charac- teristic curves for diﬀerent types of substances (model ﬂuid, alkane, alcohol, and aromatic component). The new approach combines molecular-based equation of state calculations for determining initial guess values for the characteristic curves and three sequential molecular simulation steps. In the ﬁrst simulation step, the second virial coeﬃcient is determined, 4 which provides the zero-density limit state points of the characteristic curves. In the second simulation step, simulations are systematically carried out in the vicinity of the initial guess values, which provides the characteristic curve state points of the given force ﬁeld. In the third simulation step, an additional simulation is carried out at each determined character- istic curve state point. For the second simulation step, diﬀerent thermodynamic conditions (cf. Table 1) were systematically compared in this work. Thereby, the most favorable route was identiﬁed for each characteristic curve. The rigorous approach proposed in this work has several advantages: (a) it can be applied to both simple model ﬂuids and complex molecular ﬂuids; (b) statistical uncertainties of the characteristic curve state points are estimated; (c) the thermodynamic consistency of the results is assessed by the method itself; (d) characteristic curve state points are also obtained for the vapor-liquid and solid-liquid metastable regions; (e) high-accurate results for the characteristic curves are appended in the zero-density limit (obtained from the virial route). For applying the computational procedure, we provide a source code implementation of the entire simulation workﬂow that can be straightforwardly used in combination with a given simulation engine. In this work, the simulation engine ms215was used in combination with the MolMod database.16 The method was applied to six substances to demonstrate its general applicability: the classical Lennard-Jones (LJ) model ﬂuid14(used for validation with results from the lit- erature); toluene, methane (using the Lennard-Jones ﬂuid results), ethane, propane, and ethanol as real substance ﬂuids; and the Lennard-Jones truncated and shifted (LJTS) ﬂuid, which is an important model ﬂuid in computational physics and physical chemistry.17–23 This work is structured as follows: ﬁrst, the method and tools are introduced. Then, results for the ﬂuids are presented and discussed (in parts in the Supporting Information). All primary simulation data is reported in the electronic Supporting Information. 5 Development of Computational Procedure In the following, the method for determining Brown’s characteristic curves using molecular simulation is presented. This section is structured as follows: ﬁrst, an overview of the computational procedure is given. Then, the central elements are discussed in more detail following the sequential order of the algorithm, i.e. the determination of the initial guess values, the setting-up of the simulation points, details on the molecular simulations, and the evaluation of the characteristic curve state points and their statistical uncertainties as well as the assessment of their thermodynamic consistency. Moreover, the Matlab implementation of the algorithm is presented. Overview Fig. 2 gives an overview of the method. It can be applied to an arbitrary molecular ﬂuid described by a force ﬁeld. At the core of the method, molecular simulations are carried out in three steps that are mediated by molecular-based equation of state calculations. The algorithm is structured in four parts. First, initial guess values for all three characteristic curves are calculated using a molecular- based EOS in combination with second virial coeﬃcient data (purple part in Fig. 2). The latter is directly obtained from the force ﬁeld. Thus, all four characteristic curves as well as the vapor-liquid equilibrium are computed from the EOS. Using the virial coeﬃcient data of the force ﬁeld, the zero-density limit of the EOS characteristic curve data is corrected such that the limit of the initial guess values converge to the correct limit. The second virial coeﬃcient is computed from Mayer’s f-function by averaging over numerous random orientations at diﬀerent radii based on a Monte Carlo approach.24Alternatively, for simple (e.g. spherical) force ﬁelds, the integration of Mayer’s f-function can be carried out directly by integration of the interaction potential. Both routes yield high accurate virial coeﬃcient data and are part of the implementation of the algorithm (see below). 6 Second, based on the virial coeﬃcient-corrected EOS initial guess values, homogeneous state molecular simulations are carried out (green part in Fig. 2). Therefore, a simulation state point set is generated as a grid in the vicinity of the initial guess values for each char- acteristic curve. For the evaluation of the diﬀerent thermodynamic conditions (cf. Table 1), diﬀerent state point grids are required. A state point grid consists of multiple sets of simu- lation points; each set of simulation points is used to determine one state point on a given characteristic curve. Depending on the thermodynamic condition, the state points within a set of simulation points are deﬁned at a constant thermodynamic variable X, e.g. the temperature (see details below). Depending on the considered ﬂuid, the used force ﬁeld, and the reliability of the initial guess values, the width of the simulation state point grid should be adapted. For each state point on the grid characterized by a second thermody- namic variable Y, a molecular simulation is carried out to sample the values required for the evaluation of the applied condition of a thermodynamic property W, e.g. the compressibility factor (cf. Table 1). For determining a characteristic curve state point, each simulation set is evaluated independently. The respective thermodynamic condition is applied by ﬁtting a second-order polynomial function to the results from each simulation set. The polynomial function is evaluated such that the thermodynamic variable Y, e.g. the density, specifying the characteristic curve state point is determined. As a simulation set is speciﬁed with one thermodynamic variable Xbeing constant (and known), the characteristic curve state point is known at this point by two independent variables X;Y (two out of pvT). Third, the statistical uncertainties of the characteristic curve state points are estimated from the statistical uncertainties of the simulation results of each individual state point of a given set using error propagation (red part in Fig. 2). Finally (blue part in Fig. 2), an additional simulation is carried out for each characteris- tic curve state point to obtain further thermodynamic properties and in particular the state variable Z(third out of pvT). The state points obtained for the four characteristic curves from the molecular simulation sets are moreover extended by the zero-density limit charac- 7 teristic curve state points obtained from the virial route. Based on all compiled results, the thermodynamic consistency of the results is assessed. In this work, diﬀerent thermodynamic conditions (cf. Table 1) and diﬀerent types of sim- ulation state point grids (cf. Fig. 3) were used and compared regarding their computational eﬃciency, robustness, and accuracy. Here, the simulations were carried out in the NVT or NpT ensemble (but also other ensembles could in general be used). The method outlined here is in general compatible with any molecular simulation engine. Also, both MD and MC sampling can in general be applied. In this work, the simulation code ms2 was used15,24,25for testing the method. The code ms2 can be used beneﬁcially in the computational procedure proposed in this work as it automatically samples a large number of thermodynamic properties, including the Helmholtz energy derivatives with respect to density and inverse temperature using the Lustig formalism.26–28Thereby, a large number of thermodynamic characteristic curve conditions can be tested and consistency tests be applied. The computational procedure proposed in this work was tested on diﬀerent types of substances. For determining the most favorable sampling method, the Lennard-Jones ﬂuid and toluene were used. Then, the most favorable sampling route was moreover tested on ethane, propane, ethanol, and the LJTS ﬂuid. For all substances, rigid united atom force ﬁelds were used, i.e. no internal degrees of freedom were considered and hydrogen atoms in CH, CH 2, and CH 3groups were implicitly modeled as fused interaction sites. For toluene, the force ﬁeld proposed by Guevara-Carrion et al.29was used, which consists of seven Lennard- Jones interactions sites – each modeling a CH xgroup – and six point quadrupole moments that model the -orbital of the aromatic component. For ethane, the force ﬁeld of Vrabec et al30was used, which consists of two Lennard-Jones interaction sites modeling the two methyl groups. For propane, the rigid version of the force ﬁeld proposed by Martin and Siepmann31 was used. It consists of three Lennard-Jones interactions sites modeling the methylene group and the two methyl groups. For ethanol, the force ﬁeld proposed by Schnabel et al32was 8 used, which consists of three Lennard-Jones interaction sites modeling the the methylene and methyl group as well as the oxygen atom. It moreover consists of three point charges modeling the polarity of the molecule – in particular for capturing the ability for hydrogen bonding. For methane, the model proposed by Vrabec et al.33was used, which consists of a single Lennard-Jones interaction site. Therefore, the Lennard-Jones ﬂuid simulation results were used for modeling methane by applying the reduced units principle. Initial Guess Values and Setup of Simulations Molecular-Based Equation of State A molecular-based EOS is used for computing initial guess values for the four characteristic curves as well as the VLE. For molecular-based EOS, the component-speciﬁc parameters are usually ﬁtted to VLE data, i.e. saturated densities and the vapor pressure. Due to the strong physical backbone of molecular-based EOS, these models often allow an extrap- olation to extreme conditions.3,34Examples are the SPHCT EOS35and the SAFT-VR Mie EOS.36Equations of state used within the computational procedure proposed in this work are required to produce physically correct characteristic curves, which is not the case for all EOS.3,37Moreover, the EOS should have a physically realistic behavior in the VLE two-phase region as the Boyle curve originates on the spinodal. Hence, the EOS has to be carefully selected. Moreover, depending on the accuracy of the EOS, the width of the stimulation state point grid has to be adapted. Here, the SAFT-VR Mie EOS36was used for demonstration purposes. The pure com- ponent parameters used in this work for modeling methane, ethane, propane, ethanol, and toluene are reported in the Supporting Information. For ethane, propane, and toluene, the chain term was used. For ethanol, additionally, the SAFT association term38was used. In the case that no pure component parameters for a reliable molecular-based EOS are avail- able to determine initial guess values, such can simply be obtained using the corresponding states principle. Hence, using EOS results for the characteristic curves for a simple ﬂuid, e.g. 9 the Lennard-Jones ﬂuid,3and scaling those with the critical parameters of the studied ﬂuid. This approach is demonstrated here (see Supporting Information) for the Lennard-Jones truncated and shifted (LJTS) ﬂuid, where no EOS is presently available that extrapolates well to high pressure. Nevertheless, upon applying this corresponding states principle route, the simulation state point grid should be adapted. Characteristic Curve Zero-Density Limit from Second Virial Coeﬃcient The characteristic curve initial guess values at low densities are sensitive for the compu- tational procedure. If the zero-density limit obtained from the EOS deviates signiﬁcantly from the corresponding force ﬁeld values, the simulations at low densities become erroneous. Therefore, the force ﬁeld second virial coeﬃcient is taken into account already for the initial guess values. This is done by scaling the EOS results in the vicinity of the zero-density limit. Thereby, the corrected initial guess values converge to force ﬁeld characteristic curve results at→0. The scaling is applied for all four characteristic curves using the corresponding charac- teristic second virial coeﬃcient values, cf. Fig. 1. Therefore, the second virial coeﬃcient is determined in a wide temperature range by sampling Mayer’s f-function from the inter- molecular interactions at diﬀerent conﬁgurations as39,40 B=−2/integral.disp∞ 0/uni27E8.alt2exp/parenleft.alt4−uij(rij;!i;!j) kBT/parenright.alt4−1/uni27E9.alt2 !i;!jr2 ijdr ; (6) where uijindicates the intermolecular potential, rijindicates the center of mass distance between two molecules iandj, and the /uni27E8:::/uni27E9brackets indicate the averaging over diﬀerent ori- entation !iand!j. Monte Carlo molecular simulations are carried out for evaluating Eq. (6) at a given temperature. Thereby, randomly generated distances and orientations are evalu- ated. For large ﬂexible molecules (not considered in this work), also diﬀerent conformations are to be considered. 10 For simple spherical model potentials, Eq. (6) simpliﬁes to B=−2/integral.disp∞ 0/parenleft.alt2exp/parenleft.alt4−uij(rij) kBT/parenright.alt4−1/parenright.alt2r2 ijdr : (7) From either Eq. (6) or (7), the second virial coeﬃcient is determined as a function of the temperature B(T). From these results, the three characteristic temperatures are determined that correspond to the zero-density limit state points of the four characteristic curves – the Boyle and the Zeno curve have a common zero-density limit (cf. Fig. 1). The three characteristic virial coeﬃcient points are:3,5B=0,B=Bmax, and @B/slash.left@T=B/slash.leftT, which correspond to the four characteristic temperatures Tlimit ;Z=Tlimit ;B,Tlimit ;A, and Tlimit ;C, respectively (cf. Fig. 1). For applying Eq. (6), MC simulations are carried out in a wide temperature range and the characteristic virial coeﬃcient points are determined by interpolation. For applying Eq. (7), a Newton solver is applied to directly solve the conditions for the characteristic virial coeﬃcient points. Hence, while small statistical uncertainties apply to the simulation route using Eq. (6), practically no uncertainties apply to the direct potential evaluation route using Eq. (7). For the latter, the zero-density limit characteristic curve state point results are exact (within the signiﬁcant digits used for the computer precision calculations). In this work, Eq. (6) was applied for the evaluation of the toluene, ethane, propane, and ethanol force ﬁeld using MC simulations in ms2 (see Ref.24for details), whereas Eq. (7) was applied for the evaluation of the LJ and LJTS potential. The zero-density limit characteristic curve temperatures Tlimit ;^iwith ^i=Z; A; B; C ob- tained from the virial route are used to correct the EOS initial guess values. Thereby, the virial-corrected characteristic curve temperatures Tcorr ^iare calculated as Tcorr ^i=TEOS ^i+0:5/parenleft.alt21+tanh/parenleft.alt1K^i(TEOS ^i−0:95Tvirial limit ;^i)/parenright.alt1/parenright.alt2(Tvirial limit ;^i−TEOS limit ;^i); (8) where TEOS ^iindicates the vector of temperature values predicted by the EOS for a given 11 characteristic curve ^i=Z; A; B; C andTcorr ^iindicating the virial-corrected vector of tem- perature values. The temperatures TEOS limit ;^iandTvirial limit ;^iindicate the zero-density limit values obtained from the EOS and the virial route, respectively. Based on Eq. (8), the temperature of a given characteristic curve is shifted by Tvirial limit ;^i−TEOS limit ;^istarting at 95% of the zero-density limit temperature. This shift is only applied in the vicinity of the zero-density limit using a tanh function. Thereby the virial-corrected initial guess values are continuous and exhibit the force ﬁeld zero-density limit. The parameter K^iin Eq. (8) is calculated by K^i=2:647 TEOS limit ;^i−0:95Tvirial limit ;^i(9) to scale the tangent hyperbolic function. Setup of the Simulations Based on the virial coeﬃcient-corrected EOS initial guess data, the simulation state point grid is generated. In general, the simulations can be carried out in diﬀerent ensembles. Here, the NVT and NpT ensemble are considered. Based on the ensemble used in the simulations, the grid has to be speciﬁed such that the simulations from a set have the same thermodynamic variable X. Fig. 3 schematically shows the simulation grid conﬁgurations considered in this work. The simulations used for determining one characteristic curve state point are thereby speciﬁed either at isothermal, isochoric, or isobaric conditions. Moreover, a method is tested where the sets of simulation state points are approximately orthogonal to the respective characteristic curve in either the density-temperature or pressure-temperature projection. This "gradient" method can only be used to evaluate conditions which are not linked to an iso-condition, e.g. for the compressibility factor condition Z=1. Most simulation grid types depicted in Fig. 3 can be combined with the diﬀerent ther- modynamic characteristic curve conditions (cf. Table 1). However, not all combinations of simulation grid and thermodynamic condition are feasible, e.g. the evaluation of a derivative 12 at isochoric conditions requires an isochoric arrangement of the simulation set. Table 2 sum- marizes the simulation settings considered in this work for all four characteristic curves. For a given characteristic curve, the initial guess value was the same for all studied simulation settings, cf. Table 2. The simulation grid is generated based on a heuristic function depending on both the temperature and the distance to speciﬁc points, e.g. the zero-density limit and the critical point (cf. Fig. 3). Details on the grid generator are given in the Supporting Information. Simulations in the unstable region (as predicted by the EOS) are circumvented. Also, the large diﬀerence of the compressibility of the ﬂuid in the diﬀerent state regions is considered by the grid generator such that the characteristic curve state points can be determined with high accuracy. The grid generator covers three orders of magnitude in pressure and two orders of magnitude in temperature and, moreover, covers the metastable VLE and SLE region. The grid generator has some method-speciﬁc parameters such as the number of state points Nsimused in a simulation set, the number of simulation sets Nset(corresponding to the number of characteristic curve state points to be determined), and the parameter Lcharacterizing the scattering width of the data points of a simulation set around the initial guess value (see Supporting Information for details). The simulation parameters are recommended to be in the range Nsim=5:::10,Nset=10:::30, and LNVT=0:02:::0:03for NVT andLNpT=0:05:::0:08for NpT simulations. Molecular Simulations Simulation Settings Molecular dynamics simulations are carried out at each state point on the simulation grid. The simulations are carried out either in the NVT or NpT ensemble (cf. Table 2). For the barostat, the piston mass parameter Mpiston is to be speciﬁed carefully such that a stable simulation is carried out at each state point. This is not trivial as a large pressure and 13 temperature range is covered by the simulations. The piston mass parameter is given in the dimension of kg/slash.leftm4. Based on preliminary simulations, a simple heuristic function was developed that describes the barostat piston mass parameter Mpiston as a function of the temperature and pressure as Mpiston =10−7/parenleft.alt2101 3log2(p 10pc)−1 4log2(T 2:5Tc)/parenright.alt2kg m4; (10) where Tcandpcindicate the critical temperature and pressure, respectively. For the sampling of the Helmholtz energy derivatives, the Lustig formalism26–28is applied. Thereby, the Helmholtz energy derivatives up to second order are directly sampled in the NVT simulation from the ﬂuctuation theorem. Based on this data, the Helmholtz energy derivative conditions of the characteristic curves are applied, cf. Table 1. Fig. 4 shows the relations of the thermodynamic properties sampled in ms2 that are used in this work for the diﬀerent characteristic curve conditions, cf. Table 1. The left column indicates the primary variables sampled in the simulation, i.e. the velocity, potential energy, and coordinates of the particles as well as the box volume. The blue highlighted variables in Fig. 4 indicate the properties used for evaluating the diﬀerent thermodynamic conditions (cf. Table 1). As the sampling of these properties in the simulation goes hand in hand with diﬀerent degrees of noise, i.e. yield in general diﬀerent accuracy and precision, it is a main goal of this work to determine the most beneﬁcial sampling route. In this work, the system size of the MD simulations was 2,500 particles. The time step was =0:001/radical.alt1 m/slash.left"for the Lennard-Jones model ﬂuid simulations. The equilibration was carried out for 50,000 time steps in the NVT ensemble and additional 50,000 time steps in the NpT ensemble if production was carried out in the NpT ensemble. The number of production steps was 750,000. The cut-oﬀ radius was 5 for the full Lennard-Jones ﬂuid and in the case of the LJTS ﬂuid, the cut-oﬀ was 2.5 .21For the real substance ﬂuids, the time step was set to  =2:3 fsfor toluene, ethane, and propane and  =1:16 fs for 14 ethanol. The cut-oﬀ radius was rc=15Å for the real substance simulations. For the cut-oﬀ radius, the center of mass method was used. For the electrostatic long-range interactions, the reaction ﬁeld method was used. Classical velocity scaling was applied for thermostating and the Anderson barostat for prescribing the pressure. Determining Characteristic Curve State Points Based on the respective thermodynamic condition (cf. Table 1), the results from each simu- lation set are evaluated to determine a characteristic curve state point. Therefore, a second- order polynomial function is ﬁtted to the simulation results, cf. Fig. 5-left. Then, the function is evaluated such as the state point is determined that fulﬁlls the thermodynamic condition for a given characteristic curve. For determining the characteristic curve state points from the respective thermodynamic condition, one of two procedures is applied: (1) a thermodynamic property Wdirectly sampled in the simulation is to be equal to zero or unity – depending on the condition; (2) the derivative of a thermodynamic property Wwith respect to a state variable Yis to be equal to zero. The procedure (1) applies to the four Helmholtz energy derivative conditions (green labeled line in Table 1), the condition Z=1for the Zeno curve, and the two derivatives (@u/slash.left@v)T;ms 2and(@h/slash.left@p)T;ms 2(directly sampled by ms2). The procedure (2) applies to all remaining thermodynamic conditions considered in this work (cf. Table 1). For procedure (1) cases, the zero (or unity) crossing point of the ﬁtted curve determines the characteristic curve state point regarding Y. For procedure (2) cases, the minimum of the ﬁtted curve determines Y. Examples for the ﬁts to the simulation results for all considered routes (cf. Table 1) are presented in the Supporting Information. From the evaluation of the simulation set results, each characteristic curve state point is suﬃciently speciﬁed by two independent thermodynamic state variables XandY. One state variable Xwas prescribed as boundary condition and one variable Yis determined from the ﬁtting procedure. The remaining variable Z(out of pvT) is – as an intermediate step – estimated by a cubic interpolation 15 with not-a-knot end conditions between the results of the simulation set. The full set of the three state variables pvTat a given characteristic curve state point is determined in an independent simulation, see below. Statistical Uncertainties The results for the characteristic curve state points are subject to uncertainties. Such arise from the molecular simulations themselves, the quality of the ﬁtting procedure described above, and the quality of the initial guess values. Here, a new method for estimating the statistical uncertainties is proposed. Thereby, both the uncertainties from the simulations as well as those from the ﬁtting procedure are considered. The uncertainty estimation method – termed error ellipse method – is schematically shown in Fig. 5. Therein, the error ellipse method is demonstrated on procedure (2) type data. Hence, the derivative of a given property Wwith respect to a given state property Yis to fulﬁll the condition @W/slash.left@Y=0. Both properties WandYare subject to statistical uncertainties due to the numerical nature of the molecular simulation. The statistical uncertainties WmnandYmnare obtained from the individual simulations with m=1:::N setandn=1:::N sim. In this work, the statistical uncertainties WmnandYmnwere computed in ms2 using the classical block averaging method.41 To determine the statistical uncertainty of the characteristic curve state point Ym, an ensemble N tof worst case ﬁts is evaluated. Therefore, an ellipse for each simulation result data point (Wmn;Ymn)is computed with the semi-major and semi-minor axes being the statistical uncertainty of the data point WmnandYmn, cf. Fig. 5-left. By randomly selecting a set of points on the ellipsis of the data points of a simulation set, a pseudo data point set is generated (green data points in Fig. 5-left). For the pseudo data point set, the ﬁtting procedure (as described above) is applied again. From that ﬁt, a new pseudo characteristic curve state point Y′ mois determined with o=1:::N t. The standard deviation 16 is then computed as Ym=/roottop /rootmod/rootmod/rootbot1/slash.leftN tN t /summation.disp o=1/parenleft.alt2Y′mo−/parenleft.alt41/slash.leftN tN t /summation.disp o=1Y′mo/parenright.alt4/parenright.alt22 : (11) This procedure is repeated N t-times until the standard deviation Ymconverges, cf. Fig. 5-right. As a criterion for convergence, the ﬂuctuation of Ymis to decay below 1% over the previous 50 evaluation steps. Simulations at Characteristic Curve State Points For each determined characteristic curve state point, an additional NVT MD simulation is carried out. This is done for two reasons: (a) to fully determine the thermodynamic charac- teristic curve state point pvT(one of the three variables is only known by spline interpolation of the sample points of the simulation set) and (b) to apply the Lustig formalism26–28such that ~aklwith k;l=1;2are determined at the characteristic curve state points. Thereby, the characteristic curve state points can be assessed regarding their thermodynamic consistency using the Helmholtz energy derivative conditions, cf. Table 1. Moreover, a large number of thermodynamic properties at the characteristic curve state points becomes available. The simulation settings are equal to the NVT simulation settings described above. Compiling Results and Assessment of Thermodynamic Consistency In a ﬁnal step, the results are compiled and assessed regarding the thermodynamic consis- tency. For the data compilation, the zero-density characteristic curve state points obtained from the virial route are added to the molecular simulation results obtained from the grid evaluation. For the latter, the statistical uncertainties were determined from the error ellipse method, whereas practically no uncertainties apply for the zero-density characteristic curve state points. Three consistency tests are applied to the ﬁnal data set: (1) for each characteristic curve, 17 the data points obtained from molecular simulations at low densities are compared to the zero-density limit data points obtained from the virial coeﬃcient route. The latter has – especially in the case of simple model ﬂuids – a signiﬁcantly higher precision and accuracy. Hence, the molecular simulation results at low densities are to converge into the zero-density limit data to be thermodynamically consistent. (2) The Helmholtz energy derivative data ~aklsampled at a given characteristic curve state point is used to validate the thermodynamic condition for the characteristic curve (cf. green colored line in Table 1). (3) It is veriﬁed that the molecular simulation characteristic curve results comply with the regularities postulated byBrown .1 Implementation The computational procedure outlined above for determining Brown’s characteristic curves is available in the electronic Supporting Information as a script-based implementation of the algorithm. These workﬂow manager scripts are written in Matlab and shell. The presented implementation contains functions for pre-processing, post-processing, and the evaluation of the three simulation steps. The implementation of the algorithm is depicted in Fig. 6. The implementation enables a highly automated execution of the entire algorithm (cf. Fig. 2). Thereby, the three simulation steps are subsequently carried out: (1) the calculation of the virial coeﬃcient, (2) the determination of the characteristic curve state points, and (3) a ﬁ- nal simulation at the estimated characteristic curve state points. In each of these simulation steps, the workﬂow manager can be used for four main tasks: (1) to create the input ﬁles for the simulations, (2) to launch the simulations on a given cluster machine, (3) to carry out the post-processing, and (4) to perform the physical evaluation of the results. Accord- ingly, the main part of the workﬂow manager consists of four functions: PreProcessing.m , CreateParFiles.sh ,PostProcessing.m , and Evaluation.m . Since these four tasks are similar for the three simulations, the workﬂow manager was eﬃciently implemented in four functions that are used in each simulation step. The simulation step is speciﬁed by the user 18 via the argument MODE , cf. Fig. 6. Both, the SI units system and the Lennard-Jones units system are implemented in the workﬂow manager. The SI units system can be favorably used for real substances, e.g. toluene, whereas the reduced units system can be favorably used for model substances such as the Lennard-Jones ﬂuid14or the 2CLJQ model system.42,43The user can specify the desired units system via the UNITS argument in the function call. The implementation was used in this work in conjunction with ms2,15but it can be straightforwardly adapted to other simulation engines. The force ﬁeld input ﬁles can be directly taken from the MolMod database.16Thereby, the implementation presented here can be ﬂexibly integrated in diﬀerent user workﬂows. The workﬂow has some simulation parameters, e.g. for specifying the number of tem- peratures to be evaluated for the second virial coeﬃcient, the number of simulations in a simulation set Nsim, and the number of sets Nsetfor a given characteristic curve. Details on the simulation parameters of the workﬂow manager are given in the Supporting Information. Comparison of Sampling Methods For each of the four characteristic curves, multiple sampling approaches were tested in this work (all available in the implementation). Table 2 gives an overview of all tested approaches, which were applied to the Lennard-Jones ﬂuid and toluene. For the Zeno and the Boyle curve, six approaches were compared; seven approaches for the Amagat; and eight approaches for the Charles curve. Each approach consists of a thermodynamic condition (cf. Table 1) and a simulation grid type (cf. Fig. 3). To identify the most favorable approach for each curve, the following criteria were used in the comparison: (i) the scattering of the determined characteristic curve state points, (ii) the magnitude of the statistical uncertainty Y, (iii) the quality of the ﬁt curves, (iv) the behavior of the results in the low-density region and their consistency with the zero-density 19 limit data point, (v) the behavior of the results at low temperatures in the vicinity of the VLE binodal and in the VLE metastable region, (vi) the behavior of the results in the high- density region in the vicinity of the solid-ﬂuid equilibrium, (vii) the complexity of the grid, (viii) the computational complexity of the molecular simulation method, (ix) the accessibility of the required properties, and (x) the required simulation time. Based on this systematic evaluation, for each characteristic curve, the most beneﬁcial approach was identiﬁed. Based on the comprehensive study on the Lennard-Jones ﬂuid and toluene, the most fa- vorable sampling route was identiﬁed among several ones that were considered, cf. Table 2. The results for the Lennard-Jones ﬂuid are discussed in detail in the main body of this work. The results for toluene are in most parts presented in the Supporting Information. The ﬁnd- ings reported in the Supporting Information for toluene support the ﬁndings reported here for the Lennard-Jones ﬂuid regarding the beneﬁts and drawbacks of the diﬀerent sampling routes. Finally, the most favorable sampling route was applied to determine the character- istic curves of ethane, propane, and ethanol. The ﬁnal results of these are presented in the main body of this work and conﬁrm the robustness of the determined sampling approach. Details are reported in the Supporting Information. For the graphical representation of the results, a consistent color and symbol code was used. The symbol color thereby indicates the thermodynamic condition applied in the data set (cf. Table 1) and the symbol type indicates the grid type (cf. Fig. 3). For all four characteristic curves, the Lennard-Jones ﬂuid results from all tested sampling routes are presented and discussed in both the pressure-temperature and density-temperature projection – including the corresponding deviation plots. The corresponding toluene results are given in the Supporting Information. Also, the (uncorrected) initial guess equation of state data are shown. The shown error bars are obtained from the error ellipse method. 20 Zeno Curve Fig. 7 shows the results for the Zeno curve of the Lennard-Jones ﬂuid. Overall, the results from the diﬀerent methods are in very good agreement. The results from this work deviate from the results from Deiters & Neumaier6,14– in particular in the vicinity of the pressure maximum of the Zeno curve. However, no uncertainties were reported for this data set.6,14 Also, only the pressure-temperature data is available for this data set.6,14 The data sets from all approaches obtained in this work show very small scattering in a wide temperature range. Nevertheless, in the vicinity of the zero-density limit, the results from all approaches show signiﬁcant scattering. All methods have similarly small error bars, which however increase upon approaching the zero-density limit. The NpT simulations yield slightly smaller error bars in the low-temperature region. The quality of the results of the ﬁt curves (presented in the Supporting Information) is overall similar for the diﬀerent ap- proaches. In the vicinity of the zero-density limit, all data sets show an increasing scattering. In the high-density region, the data points from the " Z=1NVT gradient" data set exhibit larger scattering compared to the results from other approaches. Regarding the complexity of the grid, placing the simulation points perpendicular to the curve predicted by the EOS does not provide signiﬁcant improvements. Therefore, the simple isothermal arrangement is preferred. Regarding the computational complexity, only the NpT simulations are slightly more demanding due to the barostat and the additional simulation settings required for the piston mass parameter. Also, the NpT simulations required about 7% more computer time compared to the NVT simulation (due to the additional NpT equilibration and the baro- stating). Overall, the approach " Z=1NVT isothermal" is found to be most beneﬁcial for determining the Zeno curve. The Zeno curve exhibits the peculiarity that it exhibits a linear trend in the density- temperature projection. This has been discovered by Batschinski44and more recently studied byApfelbaum et al.45,46for diﬀerent substances. The applicability of this linearity-law was tested here based on the simulation data from this work. The results presented in the 21 Supporting Information show that the computer experiment results are in good agreement with the linearity-law, but not exact. Amagat Curve The results for the Amagat curve are shown in Fig. 8. The results from the diﬀerent ap- proaches do not agree as well as for the Zeno curve. However, for the Amagat curve, the results from Deiters & Neumaier6,14agree with the results from this work within the scat- tering and the uncertainties. The ﬁrst three data points from Deiters & Neumaier6are probably located in the metastable solid-ﬂuid region.3 The results obtained from the ~a11=0route and the (@Z/slash.left@T)v=0route show some scat- tering, whereas the results obtained from the (@u/slash.left@v)T=0route are relatively smooth. Also, the results from the routes ~a11=0and(@Z/slash.left@T)v=0yield larger statistical uncertainties, which is due to the fact that the simulation results for a given set of state points exhibit some scattering and the ﬁtted polynomial curves cannot describe the data perfectly (see Supporting Information for details). Both (@u/slash.left@v)T=0data sets are in excellent agreement with the zero-density limit data point, whereas the other data sets show important scatter- ing in the vicinity of the zero-density limit. Comparing the two (@u/slash.left@v)T=0data sets in the low-temperature region, freezing was observed in the NVT simulations for the last data points (not shown in Fig. 8), whereas for the NpT simulations, freezing in the metastable region was only observed for the lowest studied temperature. Hence, the NpT simulations are more beneﬁcial for studying the solid-liquid metastable region of the Amagat curve. For the Amagat curve, the NVT ensemble simulations are computationally less expensive, i.e. they require on average 8% less computer time than the NpT simulations. The simulation engine used in this work directly samples (@u/slash.left@v)in NVT simulations. This data was also used for computing the characteristic curve state points. The results are only presented in the Supporting Information as they exhibit signiﬁcantly larger scat- tering and statistical uncertainties compared to the results obtained form the evaluation of 22 (@u/slash.left@v)T=0taking uandvsampled in the simulations. Overall, for determining the Amagat curve, the method " (@u/slash.left@v)T=0NpT isotherm" is found to be most beneﬁcial. Boyle Curve Fig. 9 shows the results for the Boyle curve. The results from this work agree reasonably well with the results from Deiters & Neumaier .6,14The results obtained from the Helmholtz energy route ~a01+~a02=0show large scattering compared to the other approaches. The four data sets from the compressibility factor routes (two from the (@Z/slash.left@p)T=0route and two from the (@Z/slash.left@v)T=0route) are in excellent agreement in a wide temperature range. The zero-density limit is captured very well by all four compressibility factor route data sets. Nevertheless, in the vicinity of the vapor-liquid equilibrium, important deviations are observed: interestingly, the results obtained from the (@Z/slash.left@p)T=0route signiﬁcantly deviate – especially in the -Tprojection – from the (@Z/slash.left@v)T=0results and the EOS initial guess value (which is known to give a realistic description of that ﬂuid region47). The results of the ﬁt curves for the " (@Z/slash.left@p)T=0NVT isotherm" data set (see Supporting Information) reveal that the second order polynomial does not describe the simulation results well at low temperatures. The results of the ﬁt curves for the " (@Z/slash.left@p)T=0NpT isotherm" data set indicate that several simulations already contained a vapor phase, which was conﬁrmed by an inspection of the simulation trajectories. The same was observed for the " (@Z/slash.left@v)T=0 NpT isotherm" simulations at low temperatures, but the eﬀect has only a minor inﬂuence on the determined state point. In the low-temperature region, only the " (@Z/slash.left@v)T=0NVT isotherm" data set yields reliable results with low scattering. The " (@Z/slash.left@v)T=0NVT isotherm" route is also most beneﬁcial regarding the computational complexity, the grid complexity, and the computational eﬀort. Hence, the " (@Z/slash.left@v)T=0NVT isotherm" route is overall found to be the most favorable route for determining the Boyle curve. 23 Charles Curve Fig. 10 shows the results for the Charles curve. Opposite to the Zeno, Amagat, and Boyle curve, several data sets are available in the literature for the Lennard-Jones ﬂuid Charles curve.6–8,10,48–50The results from the literature overall exhibit signiﬁcant scattering (dis- cussed in detail in the Supporting Information). The results from this work agree with the literature results within this scattering range. For both isobaric approaches, no results near the pressure maximum can be determined, which is simply due to the fact that the minimum of (@Z/slash.left@T)pcannot be reliably determined from the seven simulation points. Moreover, the uncertainties of the results from the isobaric approaches increase upon approaching the pressure maximum. Also, the results from the isobaric approaches exhibit signiﬁcant scattering in the vicinity of the zero-density limit. The results from the Helmholtz energy derivative route show important scattering in the entire temperature range, which also goes hand in hand with large statistical uncertainties. The quality of the ﬁt curves (see Supporting Information) is signiﬁcantly better for the two(@h/slash.left@p)T=0data sets compared to the results from the other approaches. Moreover, the results from the (@h/slash.left@p)T=0route are in excellent agreement with the zero-density limit characteristic curve state point and show an overall very smooth trend in the entire temperature range. The quality of the two data sets " (@h/slash.left@p)T=0NVT isotherm" and "(@h/slash.left@p)T=0NpT isotherm" is similar. Since the NVT simulations are computationally less complex and less expensive, the approach " (@h/slash.left@p)T=0NVT isotherm" is identiﬁed to be the most favorable route for determining the Charles curve. The simulation engine used in this work directly samples (@h/slash.left@p)in NpT simulations. These data were also used for computing the characteristic curve state points. The results are only presented in the Supporting Information as they exhibit signiﬁcantly larger scattering and statistical uncertainties compared to the results obtained form the evaluation of (@h/slash.left@p) taking handpsampled in the simulations. 24 Assessment of Thermodynamic Consistency The thermodynamic consistency of the results from the most favorable approaches for deter- mining the four characteristic curves was assessed. Therefore, an additional NVT simulation was carried out at the determined state point. Fig. 11 shows the results for the Helmholtz energy condition sampled in the ﬁnal MD simulations. Most data points satisfy the respec- tive criterion within their statistical uncertainties. Fig. 12 shows the ﬁnal results of the four characteristic curves (for both the Lennard-Jones ﬂuid and toluene). For the Lennard-Jones ﬂuid results, the criteria postulated by Brown are well satisﬁed, e.g. regarding the curva- ture, maxima, intersection points etc. Moreover, the ﬁnal MD results are consistent with the zero-density limit data points obtained from the virial route. Also, the toluene results are conform with Brown’s criteria and the zero-density limit data points. More details on the toluene results are given in the Supporting Information. Determining the Characteristic Curves of Methane, Ethane, Propane, and Ethanol The robustness and applicability of the computational procedure was further tested by de- termining the characteristic curves of methane, ethane, propane, and ethanol. For methane, the Lennard-Jones ﬂuid results were used. For ethane, propane, and ethanol, only the most favorable sampling route was applied. The results for the three n-alkanes are depicted in Fig. 13. The results for ethanol are depicted in Fig. 14. The results from the grid simulations are in excellent agreement with the zero-density limit data point obtained from the virial route – for all studied real substances. The results for the Helmholtz energy consistency test are presented in the Supporting Information. The vast majority of data points is in reason- able accordance with this consistency test. Moreover, the simulation results for all studied real substances are conform with Brown’s criteria. Interestingly, the EOS initial guess values of ethanol (cf. Fig. 14) obtained from the SAFT-VR Mie model show some deﬁciencies, i.e. an undulation for the Charles and Amagat curve. This is not the case for the three studied alkanes. From a model perspective, the main diﬀerence between the alkanes and ethanol is 25 that the association SAFT term was used in the latter. Yet, the computational procedure is suﬃciently robust such that characteristic curve state points were reliably determined in the entire state range. Based on the results for the three studied n-alkanes (cf. Fig. 13), the inﬂuence of the chain length on Brown’s characteristic curves can be discussed: the overall shape does not change with increasing alkane chain length, which is captured by both the simulation results and the EOS. With increasing chain length, the zero-density limit data point is shifted to higher temperatures. This is in line with the behavior of the critical point. The pressure maximum value is not aﬀected as strongly by the chain length. Conclusions In this work, a method for determining Brown’s characteristic curves by molecular simulation was developed and tested. The method is generalized such that it can be applied to an arbitrary force ﬁeld of a molecular ﬂuid. Initial guess values are taken from a molecular- based equation of state. At the core of the computational procedure lies a three-stage molecular simulation sequence: ﬁrst, the second virial coeﬃcient is sampled in a Monte Carlo run. Second, simulations are carried out on a grid of state points (derived from the virial-corrected EOS initial guess data). From the results, characteristic curve state points are determined. Third, a ﬁnal simulation at the determined state point is carried out. Finally, the results are assessed regarding the thermodynamic consistency. For each characteristic curve, diﬀerent simulation routes were systematically tested and compared. The tests were carried out for two ﬂuids, the Lennard-Jones ﬂuid and toluene. From this systematic comparison, for each characteristic curve, the most beneﬁcial approach was identiﬁed that yields reliable results in the entire temperature range. The four most beneﬁcial approaches are: •Z=1in NVT simulations for the Zeno curve, 26 •(@u/slash.left@v)T=0in NpT simulations for the Amagat curve, •(@Z/slash.left@v)T=0in NVT simulations for the Boyle curve, and •(@h/slash.left@p)T=0in NVT simulations for the Charles curve. Moreover, the method was further tested on ethane, propane, ethanol, and the Lennard- Jones truncated & shifted ﬂuid. The results support the robustness and general applicability of the approach. An implementation of the developed algorithm is provided as an automated workﬂow manager tool. Despite the fact that the tool was used in combination with the simulation engine ms2 in this work, the implementation can readily be used with other simulation engines. Based on the computational procedure and tool developed in this work, reliable predic- tions for Brown’s characteristic curves can be obtained for a given force ﬁeld model, which can then be used for testing the extrapolation behavior of equations of state. In this work, rigid united-atom force ﬁelds for relatively small molecules were used. For future work, it would be interesting to compare diﬀerent force ﬁeld types, e.g. all-atom, coarse grain, and reactive force ﬁelds. Moreover, classical Lennard-Jones force ﬁelds were used here, which have the repulsive exponent conveniently chosen as 12, which is known to be problematic – especially at elevated pressure. Hence, it would be interesting to study Brown’s character- istic curves of the Mie potential, which is known to provide a more accurate description of the repulsive interactions of real ﬂuids. Conﬂicts of Interest There are no conﬂicts of interest to declare. 27 Acknowledgement The authors gratefully acknowledge ﬁnancial support by the Deutsche Forschungsgemein- schaft (DFG) within the IRTG 2057 Physical Modelling for Virtual Manufacturing Systems and Processes (project number 252408385), by the BMBF within the WindHPC project, and by the KSB foundation. The simulations were carried out on the ELWE supercomputer at Regional University Computing Center Kaiserslautern (RHRK) under the grant TUK-MTD. Supporting Information Available 28 Tables & ﬁgures Table 1: Thermodynamic conditions of the four characteristic curves. The color bar indicates thermodynamic conditions considered in this work. The color bar corresponds to the color code used in the plots. The approach identiﬁed as the most favorable one is colored in gray. Zeno (Z) Amagat (A) Boyle (B) Charles (C) ~a01=0 ~ a11=0 ~ a01+~a02=0 ~ a01+~a02+~a11=0 Z=1 /parenleft.alt1@Z @T/parenright.alt1v=0 /parenleft.alt1@Z @p/parenright.alt1T=0 /parenleft.alt1@Z @T/parenright.alt1p=0 /parenleft.alt1@u @v/parenright.alt1T=0 /parenleft.alt1@Z @v/parenright.alt1T=0 /parenleft.alt1@Z @v/parenright.alt1p=0 /parenleft.alt1@h @p/parenright.alt1T=0 (A) /parenleft.alt1@u @v/parenright.alt1T;ms 2=0 /parenleft.alt1@h @p/parenright.alt1T;ms 2=0 /parenleft.alt1@p @T/parenright.alt1v=p T/parenleft.alt1@p @v/parenright.alt1T=−p v/parenleft.alt1@v @T/parenright.alt1p=v T /parenleft.alt1@T @p/parenright.alt1h=0 (A)Derivatives sampled directly by the molecular simulation engine ms2. 29 Table 2: Considered computational approaches for determining the characteristic curves: Zeno (Z), Amagat (A), Boyle (B), and Charles (C). All combinations of thermodynamic condition (cf. Table 1) and simulation grid type (cf. Fig. 3) considered in this work. The approach identiﬁed as the most favorable one is colored in gray. condition grid condition grid Z~a01=0 NVT isotherm B~a01+~a02=0 NVT isotherm ~a01=0 NVT gradient ~a01+~a02=0 NVT gradient Z=1 NVT isotherm (@Z @p)T=0 NVT isotherm Z=1 NVT gradient (@Z @p)T=0 NpT isotherm Z=1 NpT isotherm (@Z @v)T=0 NVT isotherm Z=1 NpT gradient (@Z @v)T=0 NpT isotherm A~a11=0 NVT isotherm Ca01+~a02+~a11=0NVT isotherm ~a11=0 NVT gradient ~a01+~a02+~a11=0NVT gradient (@Z @T)v=0 NVT isochor (@Z @T)p=0 NpT isobar (@u @v)T=0 NVT isotherm (@Z @v)p=0 NpT isobar (@u @v)T=0 NpT isotherm (@h @p)T=0 NVT isotherm (@u @v)T;ms 2=0NVT isotherm (@h @p)T=0 NpT isotherm (@u @v)T;ms 2=0NVT gradient (@h @p)T;ms 2=0 NpT isotherm (@h @p)T;ms 2=0 NpT gradient 30 Amagat (A) Charles (C) Zeno (Z) BoyleCZ BC(B)AZ BZB C A Tp BZ Tlimit,ATlimit,CTlimit, limit,Z B=TFigure 1: Scheme of Brown’s characteristic curves in the double logarithmic pressure- temperature projection (top) and the second virial coeﬃcient (bottom). Red dashed lines indi- cate the link between the zero-density limit of the characteristic curves and the second virial co- eﬃcient. Top plot: thick black line indicates vapor-liquid equilibrium binodal, dashed black lines indicate vapor-liquid equilibrium spinodals, and thin black lines indicate the four characteristic curves. Red symbols indicate characteristic points; the black star indicates the critical point. 31 BCC estimate from molecular-based EOS converged? yesnomolecular simulation of Mayer function (complex force ﬁelds) direct integration of Mayer function (model ﬂuids) compute virial coeﬃcient compute zero-density limit BCC state points from virial coeﬃcient adjust zero-density limit of BCC EOS estimate initialise state point grid perform molecular simulations of state point sets of grid derive BCC state points from molecular simulation resultscompute error bar ellipsis select set of random points on ellipsis derive BCC state point from selected ellipsis points compute standard deviation for BCC /c100Ym state point m compute statistical uncertainty /c100Ymfor given BCC state point mperform molecular simulations found BCCat state points add zero-density limit BCC state points to ﬁnal results assessment of thermodynamic consistency of ﬁnal BCC resultsFigure 2: Flowchart of the proposed method for determining Brown’s characteristic curves (BCC). Purple indicates the computation of an initial guess of the characteristic curves – includ- ing the exact zero-density BCC state points from the virial route. Green indicates the central molecular simulations for determining the characteristic curve state points. Red indicates the cal- culation of the statistical uncertainty. Blue indicates the data compilation and assessment of the thermodynamic consistency. 32 /s49 /s50 /s51 /s114/s102/s97/s118/s46 /s102/s97/s118/s46/s112/s114/s112/s112/s114/s90/s101/s110/s111 /s65/s109 /s97/s103/s97/s116 /s66/s111/s121/s108/s101 /s67/s104/s97/s114/s108/s101/s115 /s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109 /s78/s112/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109 /s78/s86/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116 /s78/s112/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s78/s86/s84/s32/s105/s115/s111/s99/s104/s111/s114 /s78/s112/s84/s32/s105/s115/s111/s98/s97/s114 /s84 /s84/s84/s84/s102/s97/s118/s46 /s102/s97/s118/s46Figure 3: Grid types for simulation sets tested for the four characteristic curves. Gray dots in- dicate the simulation state points. Black lines and triangles indicate the virial-corrected EOS initial guess. The most favorable grid type identiﬁed in this work is indicated by ’fav. ’ . The black symbol depicted in the ﬁrst column is consistently used in the following for the respective grid type. primary data sampled in simulationVboxvparticles ρEkinh Zu v akl~Txyzparticles virialEpot uparticles u p Figure 4: Schematic representation of the relation of the sampled properties (exemplaricly shown for ms2 here). The potential energy of the particles uparticles , the position of the parti- clesxyzparticles , the velocity of the particles vparticles , and the volume of the simulation box Vbox are the primary variables sampled in the simulation. The other thermodynamic properties are derived from these: the potential energy of the system Epot, the virial of the system, the total kinetic energy of the system Ekin, the internal energy u, its derivative with respect to volume @u/slash.left@v, the pressure p, the enthalpy h, the compressibility factor Z, and the Helmholtz energy derivatives ~akl. The properties indicated in blue were used in this work for determining the char- acteristic curves. 33 /s87 /s109 /s89 /s109/s88 /s109/s32/s61/s32/s99/s111/s110/s115/s116 /s46 /s48 /s49/s48/s48 /s50/s48/s48 /s51/s48/s48/s100/s89 /s109/s32/s47/s32/s49/s48/s45/s52 /s78 /s102/s105/s116Figure 5: Schematic representation of the ﬁtting method (left) and the convergence of the error ellipse method (right). The ﬁtting method is used to determine the characteristic curve state point m(red point) with respect to the state point variable Y, where Yindicates one of ( pvT). The second-order polynomial curve (gray line) was ﬁtted to the results from the set of simulation points with n=1:::N sim(gray points), which were at constant X, where Xindicates one other state variable out of ( pvT). For determining the uncertainty Ymof the characteristic curve state point m, ellipsis were computed from the error bars WmnandYmn. Randomly, points on the diﬀerent ellipsis (green points) were used to re-ﬁt a polynomial curve (green curve) and derive Y′ m. This was repeated o=1:::N t-times and the standard deviation was computed (cf. Eq. (11)). Once Y(N t)was converged (right plot), the statistical uncertainty Ym(red error bar) of the state point variable Ymwas obtained. EOS BCC estimate molecular simulation #1PreProcessing.m (MODE = 1) PostProcessing.m (MODE = 1) Evaluation.m (MODE = 1)CreateParFiles.sh molecular simulation #2PreProcessing.m (MODE = 2) PostProcessing.m (MODE = 2) Evaluation.m (MODE = 2)CreateParFiles.sh molecular simulation #3PreProcessing.m (MODE = 3) PostProcessing.m (MODE = 3)CreateParFiles.sh Figure 6: Flowchart of the implementation of the algorithm proposed in this work. The under- lying source code ﬁles are presented in the Supporting Information. In the ﬁrst column, molecu- lar simulations are set up, carried out, and evaluated for determining the second virial coeﬃcient and its characteristic points (purple section in Fig. 2). In the second column, molecular simula- tions are set up (based on the initial guess data from the EOS and the virial coeﬃcient), carried out, and evaluated for determining the characteristic curve state points with respect to two state point variables X; Y (two out of pvT) as well as the statistical uncertainties (green and red sec- tions in Fig. 2). In the third column, molecular simulations are set up, carried out, and evaluated at the determined characteristic state point. Based on that data, the thermodynamic consistency of the results is evaluated (blue section in Fig. 2). The diﬀerent thermodynamic routes (cf. Fig. 3 and Table 1) can be applied and speciﬁed in the implementation by the user. 34 /s48/s46/s48/s48/s46/s50/s48/s46/s52/s48/s46/s54/s48/s46/s56/s49/s46/s48 /s97 /s48/s49/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s97 /s48/s49/s32/s61/s32/s48/s32/s78/s86/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s32 /s90 /s32/s61/s32/s49/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s90 /s32/s61/s32/s49/s32/s78/s86/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s32 /s32/s90 /s32/s61/s32/s49/s32/s78/s112/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s90 /s32/s61/s32/s49/s32/s78/s112/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s32 /s108/s105/s109 /s114 /s48/s32/s118/s105/s114/s105/s97/s108/s32/s114/s111/s117/s116/s101/s32 /s108/s105/s109 /s114 /s48/s32/s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s32 /s108/s105/s116/s101/s114/s97/s116/s117/s114/s101/s32 /s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s112/s32 /s47/s32 /s101/s115/s45/s51 /s48/s46/s48/s48/s46/s50/s48/s46/s52/s48/s46/s54/s48/s46/s56/s114/s32 /s47/s32 /s115/s45/s51 /s49 /s50 /s51/s48/s46/s48/s50/s48/s46/s48/s48/s48/s46/s48/s50/s112 /s45/s112 /s69/s79/s83/s32/s47/s32 /s101/s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66/s45/s48/s46/s48/s52 /s49 /s50 /s51/s48/s46/s48/s48/s53/s48/s46/s48/s48/s48/s48/s46/s48/s48/s53/s48/s46/s48/s49/s48/s114 /s32/s45 /s114 /s69/s79/s83/s32/s47/s32 /s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66Figure 7: Results for the Zeno curve of the Lennard-Jones ﬂuid. Left: Pressure-temperature projection; Right: density-temperature projection. Red and green symbols indicate the simu- lation results from this work; black crosses the results from Deiters & Neumaier6taken from the Lennard-Jones database.14Thick solid line and ﬁlled black triangle indicate the EOS ini- tial guess value for the Zeno curve. The open triangle indicates the zero-density limit Zeno curve state point obtained from the virial route. The bottom row shows the deviation plots using the EOS results as a base line. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS36 results for the VLE binodal, spinodal, and critical point, respectively. The gray shaded area indi- cates the most favorable approach. 35 /s48/s52/s56/s49/s50 /s97 /s49/s49/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s97 /s49/s49/s32/s61/s32/s48/s32/s78/s86/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s32 /s40 /s90 /s47 /s84 /s41 /s118/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s99/s104/s111/s114/s32 /s40 /s117 /s47 /s118 /s41 /s84/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s40 /s117 /s47 /s118 /s41 /s84/s32/s61/s32/s48/s32/s78/s112/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s108/s105/s109 /s114 /s48/s32/s118/s105/s114/s105/s97/s108/s32/s114/s111/s117/s116/s101/s32 /s108/s105/s109 /s114 /s48/s32/s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s32 /s32/s108/s105/s116/s101/s114/s97/s116/s117/s114/s101 /s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s112/s32 /s47/s32 /s101/s115/s45/s51 /s48/s46/s48/s48/s46/s50/s48/s46/s52/s48/s46/s54/s48/s46/s56/s49/s46/s48/s114/s32 /s47/s32 /s115/s45/s51 /s53 /s49/s48 /s49/s53 /s50/s48 /s50/s53/s49/s48/s49/s112 /s45/s112 /s69/s79/s83/s32/s47/s32 /s101/s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66/s53 /s49/s48 /s49/s53 /s50/s48 /s50/s53/s48/s46/s48/s50/s48/s46/s48/s48/s48/s46/s48/s50/s114 /s32/s45 /s114 /s69/s79/s83/s32/s47/s32 /s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66Figure 8: Results for the Amagat curve of the Lennard-Jones ﬂuid. Left: pressure-temperature projection; Right: density-temperature projection. Red, green, and blue symbols indicate the simulation results from this work; black crosses the results from Deiters & Neumaier6taken from the Lennard-Jones database.14Thick solid line and ﬁlled black triangle indicate the EOS data for the Amagat curve. The open triangle indicates the zero-density limit Zeno curve state point obtained from the virial route. The bottom row shows the deviation plots using the EOS results as a base line. In the deviation plots, for data points located at higher temperatures than the zero-density limit of the EOS, a linear extrapolation of the EOS results was applied. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS36results for the VLE binodal, spinodal, and critical point, respectively. The gray shaded area indicates the most favorable approach. 36 /s48/s46/s48/s48/s46/s50/s48/s46/s52 /s97 /s48/s49/s43/s32 /s97 /s48/s50/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s97 /s48/s49/s43/s32 /s97 /s48/s50/s32/s61/s32/s48/s32/s78/s86/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s32 /s40 /s90 /s47 /s112 /s41 /s84/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109 /s40 /s90 /s47 /s112 /s41 /s84/s32/s61/s32/s48/s32/s78/s112/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s40 /s90 /s47 /s118 /s41 /s84/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s40 /s90 /s47 /s118 /s41 /s84/s32/s61/s32/s48/s32/s78/s112/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109 /s108/s105/s109 /s114 /s48/s32/s118/s105/s114/s105/s97/s108/s32/s114/s111/s117/s116/s101/s32 /s108/s105/s109 /s114 /s48/s32/s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s32 /s108/s105/s116/s101/s114/s97/s116/s117/s114/s101 /s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s112/s32 /s47/s32 /s101/s115/s45/s51 /s48/s46/s48/s48/s46/s50/s48/s46/s52/s48/s46/s54/s114/s32 /s47/s32 /s115/s45/s51 /s50 /s51/s48/s46/s48/s52/s48/s46/s48/s50/s48/s46/s48/s48/s48/s46/s48/s50/s48/s46/s48/s52 /s32/s112 /s45/s112 /s69/s79/s83/s32/s47/s32 /s101/s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66/s50 /s51/s48/s46/s48/s52/s48/s46/s48/s50/s48/s46/s48/s48/s48/s46/s48/s50/s48/s46/s48/s52/s114 /s45/s114 /s69/s79/s83/s32/s47/s32 /s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66Figure 9: Results for the Boyle curve of the Lennard-Jones ﬂuid. Left: pressure-temperature projection; Right: density-temperature projection. Red, green, and blue symbols indicate the simulation results from this work; black crosses the results from Deiters & Neumaier6taken from the Lennard-Jones database.14Thick solid line and ﬁlled black triangle indicate the EOS data for the Boyle curve. The open triangle indicates the zero-density limit Zeno curve state point obtained from the virial route. The bottom row shows the deviation plots using the EOS results as a base line. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS36results for the VLE binodal, spinodal, and critical point, respectively. The gray shaded area indicates the most favorable approach. 37 /s48/s46/s48/s48/s46/s53/s49/s46/s48/s49/s46/s53 /s97 /s48/s49/s43/s32 /s97 /s48/s50/s43/s32 /s97 /s49/s49/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s97 /s48/s49/s43/s32 /s97 /s48/s50/s43/s32 /s97 /s49/s49/s32/s61/s32/s48/s32/s78/s86/s84/s32/s103/s114/s97/s100/s105/s101/s110/s116/s32 /s40 /s90 /s47 /s84 /s41 /s112/s32/s61/s32/s48/s32/s78/s112/s84/s32/s105/s115/s111/s98/s97/s114 /s40 /s90 /s47 /s118 /s41 /s112/s32/s61/s32/s48/s32/s78/s112/s84/s32/s105/s115/s111/s98/s97/s114/s32 /s40 /s104 /s47 /s112 /s41 /s84/s32/s61/s32/s48/s32/s78/s86/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109/s32 /s40 /s104 /s47 /s112 /s41 /s84/s32/s61/s32/s48/s32/s78/s112/s84/s32/s105/s115/s111/s116/s104/s101/s114/s109 /s108/s105/s109 /s114 /s48/s32/s118/s105/s114/s105/s97/s108/s32/s114/s111/s117/s116/s101/s32 /s108/s105/s109 /s114 /s48/s32/s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s32 /s108/s105/s116/s101/s114/s97/s116/s117/s114/s101 /s83/s65/s70/s84/s45/s86/s82/s32/s77/s105/s101/s32/s69/s79/s83/s112/s32 /s47/s32 /s101/s115/s45/s51 /s48/s46/s48/s48/s46/s50/s48/s46/s52/s48/s46/s54/s114/s32 /s47/s32 /s115/s45/s51 /s50 /s52 /s54/s48/s46/s50/s48/s46/s49/s48/s46/s48/s48/s46/s49/s48/s46/s50/s112 /s45/s112 /s69/s79/s83/s32/s47/s32 /s101/s115/s45/s51 /s84 /s32/s47/s32 /s101 /s107/s45/s49 /s66/s50 /s52 /s54/s48/s46/s48/s50/s48/s46/s48/s48/s48/s46/s48/s50/s114 /s45/s114 /s69/s79/s83/s32/s47/s32 /s115/s45/s51 /s84 /s32/s47/s32 /s101 /s107/s45/s49 /s66Figure 10: Results for the Charles curve of the Lennard-Jones ﬂuid. Left: pressure-temperature projection; Right: density-temperature projection. Red, green, orange, and blue symbols indicate the simulation results from this work; black crosses results from the literature6–8,10,48–50taken from the Lennard-Jones database.14Thick solid line and ﬁlled black triangle indicate the EOS data for the Charles curve. The open triangle indicates the zero-density limit Zeno curve state point obtained from the virial route. The bottom row shows the deviation plots using the EOS results as a base line. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS36re- sults for the VLE binodal, spinodal, and critical point, respectively. The gray shaded area indi- cates the most favorable approach. 38 /s49 /s50 /s51/s48/s46/s48/s48/s50/s48/s48/s46/s48/s48/s50/s97 /s48/s49/s32/s47/s32 /s101/s78 /s65/s90/s101/s110/s111 /s65/s109/s97/s103/s97/s116 /s67/s104/s97/s114/s108/s101/s115 /s66/s111/s121 /s108/s101/s53 /s49/s48 /s49/s53 /s50/s48 /s50/s53/s48/s46/s49/s48/s48/s46/s49/s97 /s49/s49/s32/s47/s32 /s101/s78 /s65 /s49 /s50 /s51/s48/s46/s50/s48/s48/s46/s50/s97 /s48/s49/s43/s32 /s97 /s48/s50/s32/s47/s32 /s101/s78 /s65 /s84 /s32/s47/s32 /s101 /s107/s45/s49 /s66/s50 /s52 /s54/s48/s46/s53/s48/s48/s46/s53/s97 /s48/s49/s43/s32 /s97 /s48/s50/s43/s32 /s97 /s49/s49/s32/s47/s32 /s101/s78 /s65 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66Figure 11: Consistency test for the Helmholtz energy criteria (cf. Table 1) using the Lennard- Jones ﬂuid results from the ﬁnal simulations at the determined characteristic curve state points (blue section in Fig. 2). Only results shown for the most favorable approach for each curve. /s48/s46/s49 /s49 /s49/s48/s48/s46/s48/s49/s48/s46/s49/s49/s49/s48/s112/s32 /s47/s32 /s101/s115/s45/s51 /s84/s32 /s47/s32 /s101 /s107/s45/s49 /s66/s49/s48/s48/s48 /s49/s48/s48/s48/s48/s49/s49/s48/s49/s48/s48/s49/s48/s48/s48 /s84 /s32/s47/s32/s75/s112/s32 /s47/s32/s77/s80/s97 Figure 12: Characteristic curves for the Lennard-Jones ﬂuid (left) and toluene (right) in the double logarithmic p-Tprojection. Only results shown for the most favorable approach for each curve. Symbols indicate the ﬁnal results; thick lines the (un-corrected) EOS initial guess value. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS results for the VLE binodal, spinodal, and critical point, respectively. 39 /s49/s48/s48 /s49/s48/s48/s48 /s49/s48/s48/s48/s48/s49/s49/s48/s49/s48/s48/s49/s48/s48/s48 /s84 /s32/s47/s32/s75/s112/s32 /s47/s32/s77/s80/s97/s109/s101/s116/s104/s97/s110/s101/s101/s116/s104/s97/s110/s101/s112/s114/s111/s112/s97/s110/s101 /s49/s48/s48 /s49/s48/s48/s48 /s49/s48/s48/s48/s48/s49/s49/s48/s49/s48/s48/s49/s48/s48/s48/s112/s32 /s47/s32/s77/s80/s97/s49/s48/s48 /s49/s48/s48/s48 /s49/s48/s48/s48/s48/s49/s49/s48/s49/s48/s48/s49/s48/s48/s48/s112/s32 /s47/s32/s77/s80/s97Figure 13: Characteristic curves for three alkanes in the double logarithmic p-Tprojection: methane (bottom), ethane (middle), and propane (top). Symbols indicate the ﬁnal results; thick lines the (un-corrected) EOS initial guess value. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS results for the VLE binodal, spinodal, and critical point, respectively. 40 /s49/s48/s48/s48 /s49/s48/s48/s48/s48/s49/s49/s48/s49/s48/s48/s49/s48/s48/s48 /s84 /s32/s47/s32/s75/s112/s32 /s47/s32/s77/s80/s97Figure 14: Characteristic curves for ethanol in the double logarithmic p-Tprojection. Sym- bols indicate the ﬁnal results; thick lines the (un-corrected) EOS initial guess value. Thin dashed lines, thin solid lines, and ﬁlled star indicate the EOS results for the VLE binodal, spinodal, and critical point, respectively. 41