Simulation Variables and Properties

Molecular simulations are designed to sample from a particular ensemble. At the simplest level, this means that certain variables are fixed as the 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 simulations 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 one of the independent, specified variables because experimental data are reported at particular temperatures not energies. In this case NVT simulations can be performed that sample the canonical ensemble. In this case, the number of molecules, the system volume, and the temperature are all fixed at specific values, but the energy fluctuates in time. Again, a time average of the instantaneous energy 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, all other macroscopic properties can theoretically be 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, and particulars are given in Chapter xx. The model variables required of course depends upon the complexity of the model chosen. For example, only a value for the molecular diameter, , 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 model potential which parameterizes a spherical intermolecular interaction potential in terms of a molecular size parameter, , and an attractive potential well depth constant, . With increasing 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. That is, what is the most appropriate value of these potential constants to match simulated and experimental macroscopic properties. However, even in this case where the model variables now become parameters, the model variables remain independent variables for each of the several simulations that must be performed as part of the regression analysis in finding their "best" value(s).

Symbol	Definition	Value
1. Constants
k	Boltzmann's constant	1.3806×10-23 joule/(molecK)
N0	Avagadro's number	6.022×1023
2. Simulation Variables
N	Number of molecules	~ 103
V	Simulation cell volume	~ 10-24 m3
m	Molecular mass	~ 10-25 kg/molec
	Number density	~ 1027 molec/m3
E	Energy (total)	~ 10-20 joule/molec
t	time	~ 10-12 s
3. Model Variables
	Size variable	~ 5×10-10 m
	Energy variable	~ 10-21 joule/molec
rb	Bond distance	~ 10-10 m
k	Vibrational spring constant	~ 103 joule/m2

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. Problems associated with overflow, underflow, and numerical truncation can be avoided by performing computations with dimensionless variables. However, an even more important advantage arises from the use of dimensionless variables in a simulation. 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. Section xx below will discuss this more fully.

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, , then the dimensionless volume, V* (defined as V* = V/), 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 * = ³ ~ (10²⁷ molec/m³)(5x10^-10 m)³ = 0.125. The same principle can be applied to make the other simulation variables dimensionless, but we are limited in the number of scale factors that can be applied to the number of independent units.

Unit systems use two types of units: fundamental units which are independent and derived units that are obtained from the fundamental units and are therefore dependent upon them. Fundamental units include length, mass, temperature, and time. Other units are then derived as combinations of these. In the SI system of units, the fundamental units are m, kg, K, and s, respectively. Energy in this system of units is a derived unit: 1 J = 1 kgm²/s². The derived energy unit is given by

energy = (mass)(length/time)²

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

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

Scaling the simulation variables as we did above for volume, can be alternatively thought of as expressing the variable in a new set of units. Thus, V* = 8000 represents the volume as 8000 ³ units. As with other systems of units, only four independent units can be chosen as independent; all other units are then derived from the four independent units. The volume unit used here, ³ is obviously derived from the length unit, . The four independent units chosen from Table 1 as the independent units are generally for length, m for mass, for energy, and /k for temperature. This allows all variables of length, mass, energy and temperature to be put in dimensionless form by dividing by these scaling variables. Thus, dimensionless distance, r, is given by r* = r/; dimensionless energy, E, is given by E* = E/, and dimensionless temperature is given by T* = kT/. These dimensionless variables are shown in Table 2.

Once the fundamental units, or fundamental scaling factors, have been selected, all other derived units can be derived 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, , that has units of time that is a multiplicative combination of the fundamental units. Thus,

t* = t/ ; = (m)^a ()^b ()^c (/k)^d

where the powers a, b, c, and d must be determined so as to make dimensionally consistent with time. To calculate these powers, we solve four dimensional equations such that the above equation is dimensionally consistent in each of the four fundamental units. In this case, 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

= (m/)^0.5

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.

Table 2.

Symbol	Meaning	Definition
r*	dimensionless distance	r/
E*	dimensionless energy	E/
T*	dimensionless temperature	kT/
U*	dimensionless internal energy	U/
t*	dimensionless time	t/[(m/)0.5]
v*	dimensionless velocity	v/(/m)0.5
F*	dimensionless force	F/
P*	dimensionless pressure	P3/
D*	dimensionless self diffusion coefficient	D/[(/m)0.5]

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. 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.

Table 3

Symbol	Order-of-magnitude value
V*	104
m*	100
*	10-1
E*	101
t*	10-1
P*	10-2
T*	101
rb*	10-1
kv*	105
D*	10-2

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).

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 when computed with an MD simulation. By Ergodicity, the probability densities of phase space points will be the same when sampling with MC simulations.

The second, and perhaps most important, reason for using dimensionless variables in simulations is the fact that a simulation only need be done once for a given set of dimensionless conditions. If the values of dimensionless properties are generated at a given dimensionless state point, those values can then be used to calculate the actual property for any conformal fluid. For example, hosts of simulations have been performed for LJ fluids over a wide range of * 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 *. Generally modified Benedict-Webb-Rubin equations of state have been used to correlate the results in the form P* = P*(T*, *). 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 = 0.4418 nm and /k = 230 K at T = 270 K and = 3.116x10^-5 moles/m³, we first find the dimensionless conditions T* = kT/ = 1.174 and * = ³ = 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 corresponds to the actual pressure P = P*/³ = 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

Mixture simulations are 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 and 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 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.

To understand this concept more fully, consider a binary mixture of LJ fluids characterized by the pure-component parameters ₁, ₁ and ₂, ₂. Let us choose ₁ and ₁ as the characteristic scaling variables for the simulation. The dimensionless properties are then defined in terms of these scaling constants as, T* = kT/₁, * = ₁³, P*₁³/. However, these dimensionless properties lose their generality for fluid mixtures because a specific dimensionless state, say T* = 1.0, * = 0.8 and x₁ = 0.5, does not make the resultant properties independent of the fluids. This is obvious since their is no information in these groups for 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 ₂* = ₂/₁ and ₂* = ₂/₁. 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 smaller molecule which 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 ₁₂ and ₁₂, or ₁₂* = ₁₂/₁ and ₁₂* = ₁₂/₁. While simulators will sometimes assume a model combining rule for ₁₂ and ₁₂ 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.