<header>

<author>Richard L. Rowley</author>

<affiliation>Brigham Young University</affiliation>

<creation>31 October 1997</creation>

<previous></previous>

<next></next>

</header>

Simulation Variables and Units

<!-- (level) 1 2 3 4-->

Independent and Dependent Variables. Molecular simulations are designed to sample a particular ensemble. At the simplest level, this means that certain variables are fixed as independent variables and the other properties of the system are then obtained as dependent variables. The independent variables are fixed as constants; the dependent variables will fluctuate in time about some average value. The value for a particular property is then taken to be the time average (in MD simulations) or the ensemble average (in MC simulations) of the appropriate dependent variables. For example, the simplest MD simulation might sample the microcanonical or NVE ensemble. The capital letters represent the independent variables that are fixed in the simulation. In this case, the number of molecules, N, the simulation cell volume, V, and the total energy, E, are all held constant. This is a natural ensemble to use in MD simulations because the equations of motion are conservative; i.e., E is a constant of motion during the simulation. In this case, the temperature is not constant, but fluctuates in time. The expected or average temperature can be calculated from such a simulation as a time average of the kinetic energy of the particles. On the other hand, one may prefer to perform simulations in which temperature is a fixed or independent variable because experimental data are reported at specified temperatures not energies. In this case NVT simulations can be performed, sampling the canonical ensemble. The number of molecules, the system volume, and the temperature are all fixed, but the energy fluctuates in time. Again, a time average of the instantaneous energies yields a value for the average or expected system total energy. Methods for changing which variables are independent and which are dependent are discussed in Chapter xx. The most common simulation ensembles are NVE (microcanonical), NVT (canonical), NPT (isothermal-isobaric; where P is pressure), and µVT (grand canonical; where µ is chemical potential), but others can also be devised. Once the independent simulation variables have been fixed and the simulation performed, the other macroscopic properties are obtained as dependent properties from the motion of the system in phase space.

In addition to simulation variables associated with the state of the system, variables associated with the molecular model must be used. Currently molecular models cover a wide range of sophistication and complexity, as discussed in Chapter xx. The model variables required depend upon the complexity of the model chosen. For example, only a value for the molecular diameter, $\sigma$, is required for simulation of a pure fluid if the intermolecular dispersion potential is represented with a simple hard-sphere model. Two or three constants are required if simple attractions are included. A simple model in this category is the Lennard-Jones (LJ) model potential which parameterizes a spherical intermolecular interaction potential in terms of a molecular size parameter, $\sigma$, and an attractive potential well depth constant, $\epsilon$. With increasing model complexity, more constants are required. These may include vibrational force constants, equilibrium bond length values, torsional angle potential constants, asymmetric potential parameters, etc. All of these constants are used as independent variables in the simulation to generate phase space points that can be used in obtaining the desired information about system properties or phenomena. The inverse problem can also be considered in which model parameters are obtained by varying them until the resultant simulated properties agree with experimental values. However, even in this case the model parameters are retained as independent variables for each of the several simulations that must be performed as part of the regression analysis in finding their "best" values.

<!-- (shift) -->

Handling Very Large and Small Numbers. Table 1 shows typical values of commonly used molecular variables. Because simulations involve on the order of 10³ molecules instead of 10²³, values for the mass, energy, and volume are very small numbers. These small numbers can cause problems with underflow especially when multiplying two or more numbers together. One remedy might simply to be careful of the order in which variables are used in order to avoid numerical problems within the code. Another might be to use the appropriate prefix (pico, nano, etc.) before the unit so that the value of the number is manageable. This results in a bookkeeping problem with unit conversions. Yet another solution might be to base everything on molar rather than molecular values, in effect scale Avagadro's number since so many of the numbers in Table 1 are on the order of 10^-20 to 10^-25. Both of these latter ideas involve scaling of the numbers to make them manageable.

Most simulation codes use the concept of scaling in a slightly different form to solve the problem of small numbers. Most codes are written in terms of dimensionless variables. Use of dimensionless variables eliminates problems associated with overflow, underflow, and numerical truncation, but an even more important advantage arises from their use. This is particularly true for simulations involving pure fluids, because the results of the simulations are then applicable, in their scaled form, to all conformal fluids, not just to a particular chemical compound. This is discussed more fully below.

<!-- (shift) -->
Dimensionless Variables. The problem of small numbers can be overcome simply by using simulation variables scaled in terms of the molecular model parameters listed in Table 1. For example, if all distances in the simulation are scaled with the molecular diameter parameter, $/sigma$, then the dimensionless volume, V* (defined as V* = V/$/sigma$³), has a magnitude on the order of V* ~ (10^-24 m³)/(5x10^-10 m)³ = 8000, a number much more manageable. Likewise, the dimensionless number density becomes $/rho$* = $/rho$³ ~ (10²⁷ molec/m³)(5x10^-10 m)³ = 0.125. The same principle can be applied to make the other simulation variables dimensionless, but we can only define as many scale factors as there are independent units in our variables.

There are two types of units: fundamental units are independently defined while derived units are obtained from combinations of the fundamental units. Fundamental units in standard systems of units include length, mass, temperature, and time. In the SI system of units, the fundamental units are m, kg, K, and s, respectively. Energy in the SI system is a derived unit: 1 J = 1 kgm²/s². The derived energy unit is given by

energy = (mass)(length/time)²

An alternative (but not as useful) formulation could use length, mass, temperature, and energy as the fundamental units with time as a derived unit, found from

time = (mass/energy)^0.5/(length).

We can in fact think of scaling the simulation variables, as we did above for volume, as defining a new unit system with which to express the variables. Thus, V* = 8000 represents the volume as 8000 $/sigma$³ units. Only four independent units can be independently chosen; all other units are then derived from these four. The volume unit used here, $/sigma$³ is obviously derived from the length unit, $/sigma$. Let us select from Table 1 $/sigma$ for length, m for mass, $/epsilon$ for energy, and $/epsilon$/k for temperature. This allows all variables of length, mass, energy and temperature to be put in dimensionless form by dividing the variable by these characteristic simulation constants. Thus, dimensionless distance, r, is given by r* = r/$/sigma$; dimensionless energy, E*, is given by E* = E/$/epsilon$; and dimensionless temperature is given by T* = kT/$/epsilon$. These dimensionless variables are shown in Table 2.

<!-- (shift) -->

Dimensional Analysis. Once the fundamental units, or scaling factors, have been selected, all other derived units can be obtained from them. Likewise, all other variables can be made dimensionless in terms of these four fundamental scaling variables. A convenient method to determine what combination of the fundamental variables to use to make a dependent variable or property dimensionless is to use a simple dimensionless analysis tool. For example, we need to make time, t, dimensionless. There is some de-dimensionalizing group, $/tau$, that has units of time that is a multiplicative combination of the fundamental units. Thus,

t* = t/$/tau$; $/tau$ = (m)^a ($/epsilon$)^b ($/sigma$)^c (/k)^d

where the powers a, b, c, and d must be determined so as to make $/tau$ dimensionally consistent with time. To obtain these powers, we solve four dimensional equations to ensure that the above equation is dimensionally consistent in each of the four fundamental units. In this case, $/tau/ has units of time, so the exponent on the right hand side of the equation must be zero for mass, length, and temperature, but the exponent for time must be unity. Thus,

s = [kg]^a [kgm²/s²]^b [m]^c [K]^d

exponent for kg: 0 = a + b

exponent for m: 0 = 2b + c

exponent for s: 1 = -2b

exponent for K: 0 = d

From these four equations, we see that d = 0, b = -1/2, c = - 2b = 1, and a = -b = 1/2. These values can be plugged into the original expression for to obtain

$/tau$ = (m/$/epsion$)^0.5

<!-- (shift) -->

In a similar manner, dimensionless properties can be defined for all other properties in terms of the four selected, independent units or scale factors. Some of the more common dimensionless simulation variables and their definitions are shown in Table 2.

As mentioned above, these dimensionless variables can be viewed as either scaled variables or variables in a new unit system. For example, a dimensionless distance of r* = 5 can also be viewed as an actual distance of r = 5$/sigma$. By using this molecular unit system, we eliminate the use of large and small numbers in the simulation. For example, the simulation and model variables listed in Table 1 all scale to orders of magnitude that can easily be handled by the computer in multiple calculations as shown in Table 3.

<!-- (switch) -->

Scaling. The principle of scaling is widely used in engineering and science practice. Scale-up of processes designed at bench- or pilot scale allows the engineer an opportunity to study the response of the design before expending the enormous capital required for the full-scale process. This requires a well-understood relationship between the model and the final process that the model represents. In some cases, this is as easy as multiplying linear dimensions by a scale factor. In other cases it may involve a more complicated scale-up. Chemical engineers become accustomed to dealing with dimensionless groups such as the Reynold's number, the Prandtl number, the Grasshof number, etc. These dimensionless groups arise when the differential equations describing the process are recast in terms of dimensionless variables. The dimensionless groups then appear as coefficients in the equations. This means that the resultant equations become independent of the particular fluid and must have identical solutions for any fluid when the dimensionless groups have the same value. Thus, the transition from laminar to turbulent flow occurs at a particular Reynold's number regardless of the specific fluid involved. The same principles can be applied to the differential equations describing the motion of molecules in a simulation, and, as we shall see, leads to pure-fluid simulations that are independent of the specific fluid being simulate (as long as the fluids are conformal).

<!-- (switch) -->

Conformal Fluids. Two fluids are said to be conformal if their total potential has the same form. Note that this does not mean that the two fluids have the same potential, only that the potential is of the same form. Thus, two fluids both of whose molecules interact pair-wise with Lennard-Jones potentials are conformal because their potentials become identical in dimensionless form as scaled with their individual and values. This is illustrated in Fig. 1.

When the simulation is performed, Newton's equations of motion are solved. Again, in dimensionles form, Newton's equation of motion may be written as

As shown in Figure 1, the dimensionless intermolecular potential is the same for all Lennard-Jones fluids and therefore so is the intermolecular force. This means that the solution of Newton's equation of motion yields identical dimensionless velocities and trajectories for all LJ fluids; i.e., the trajectory of the system through dimensionless phase space is identical for all conformal fluids at the same dimensionless state point. Equivalently, by Ergodicity, the probability densities of phase space points will be the same in dimensionless coordinates when sampling with MC simulations.

The second, and perhaps most important, reason for using dimensionless variables in simulations is the generality of the results of a simulation for all conformal fluids. For example, hosts of simulations have been performed for LJ fluids over a wide range of $/rho$* and T*, so much so that analytical equations of state have been regressed for the dimensionless pressure obtained from the simulations at each T* and $/rho$*. Generally modified Benedict-Webb-Rubin equations of state have been used to correlate the results in the form P* = P*(T*, $/rho$*). The simulation results, hence the correlation of the results, are identical for all LJ fluids. If for example we want to know the pressure of a LJ fluid for which $/sigma$ = 0.4418 nm and $/epsilon//k = 230 K at T = 270 K and $/rho$ = 3.116x10^-5 moles/m³, we first find the dimensionless conditions T* = kT/$/epsilon$ = 1.174 and $/rho/* = $/rho/$/sigma$³ = 0.600. The dimensionless pressure that results from the simulation or correlated equation of state depends only upon these dimensionless variables, and is P* = 0.146. For this particular fluid, this value corresponds to an actual pressure of P = P*$/epsilon$/$/sigma$³ = 36.8 MPa. Once the simulation has been done in dimensionless variables, it is done at that reduced condition for all conformal fluids.

Dimensionless Mixture Properties. The situation for mixture simulations is complicated by the fact that there is no longer an obvious choice for the characteristic molecular model constants with which to scale the variables. A binary mixture of simple LJ molecules, for example, has two sets of $/epsilon$ and $/sigma$ values, one for each component, from which to choose the characteristic energy and length parameters. It really does not matter which is chosen, as long as everything is consistently scaled to the same selected values. Scaling all variables to the selected characteristic parameters will eliminate the numerical magnitude problems just as it does for pure-component simulations, but much of the conformal fluid advantage of scaling is lost for mixtures because of the unlimited number of size and energy ratios of components comprising the mixture.

<!-- (shift) -->

To understand this concept more fully, consider a binary mixture of LJ fluids characterized by the pure-component parameters $/epsilon$₁, $/sigma$₁ and $/epsilon$₂, $/sigma$₂. Let us choose $/epsilon/₁ and $/sigma$₁ as the characteristic scaling variables for the simulation. The dimensionless properties are then defined in terms of these scaling constants as, T* = kT/₁, $/rho$* = $/rho$$/sigma$₁³, P* = P*$/sigma$₁³/$/epsilon$₁. However, these dimensionless properties lose their generality for fluid mixtures because a specific dimensionless state, say T* = 1.0, $/rho$* = 0.8 and x₁ = 0.5, does not make the resultant properties independent of the fluids. This is obvious because thre is no information in these groups about the second component. All variables must be scaled with the same characteristic constants and so the dimensionless model parameters for component two are given by $/epsilon$₂* = $/epsilon$₂/$/epsilon$₁ and $/sigma$₂* = $/sigma$₂/$/sigma$₁. Thus, the results of the simulation are not just dependent upon the dimensionless state of the system, but also upon the nature of the second component; i.e., upon the ratio of the model energy and size parameters. Generally, one still uses dimensionless variables for mixture simulations, but specifies, model parameter ratios for all of the other components. While there is no requirements on which model parameters are chosen as the characteristic de-dimensionalizing constants, the most common procedure is to choose the parameters for the smallest molecule; this makes the size ratios greater than unity.

There is yet additional mixture-specific information required for mixture simulations that is not removed by scaling. The mixture also includes cross interactions between the two components as described by $/epsilon$₁₂ and $/sigma$₁₂, or $/epsilon$₁₂* = $/epsilon$₁₂/$/epsilon$₁ and $/sigma$₁₂* = $/sigma$₁₂/$/sigma$₁. While simulators will sometimes assume a model combining rule for $/epsilon$₁₂ and $/sigma$₁₂ in terms of the pure-component values, one should always realize that their is no known relationship between the actual like interactions and the cross interactions. In fact, it is usually found that model cross interactions given by a fixed average (either algebraic or geometric) of the like interactions does not adequately represent mixture behavior. This seemingly unique nature of the cross interactions that affects the nonideality of mixtures can not be removed simply by scaling the state variables.