The calculations described here were simulations at constant
number of atoms, constant volume, and constant temperature
(constant NVT) of 5CB near an
amorphous polyethylene substrate,
and made use of the commercial software package Cerius
.[37]
The thickness of the substrate was about 2.0 nm, which should be
sufficient for representing a surface in an atomistic
molecular dynamics simulation.[38]
The simulations involved 8, 16, or 32 5CB molecules at 300 K,
which is in the nematic range of 5CB.[39]
The length of the simulations ranged from 10 ps to 100 ps
with a time step 
fs.
The positions of the atoms were recorded every 0.1 ps for later analysis.
The temperature was kept constant by scaling the velocities at
each time step.
The initial condition for most of the simulations was formed by positioning an array of 5CB molecules in a planar orientation at the surface of the amorphous polyethylene (Figure 1). The 5CB molecules, which are colored red, lie along the x-axis, the y-axis is perpendicular to the surface, and the z-axis points out of the page. The area of the substrate varied depending on the number of 5CB molecules used. Periodic boundary conditions in all three directions were applied at the boundaries of the cell shown by the dashed blue lines.
The Dreiding force field[40] was used.
The main features of this force field are the energies characterized as follows:
where 
includes the effects of covalent bonding,
includes effects, some of them long range, not directly
related to covalent bonding,
is the energy for stretching of bonds,
is the energy for bending the angles formed by each
pair of consecutive bond directions,
is the energy for rotation about the bond directions,
is the van der Waals interactions, and
is the electrostatic interaction between partial charges.
The bond energy is taken to be of the Hookian form
where 
is the force constant of the 
bond,
is the length of the bond, and
is the equilibrium length of the bond.
The energy for bending the angles formed by each
pair of consecutive bond directions is
where 
is the force constant for bending the angle,
is the angle, and
is the equilibrium value of the angle.
The energy for rotation about the bond directions is
where 
is half the rotational barrier for the bond,
and 
is the periodicity of the rotational potential
and determines the number of minima in the rotational potential.
The energy for van der Waals interactions is
where 
is the depth of the potential well,
is the separation of the and 
atoms, and
is the equilibrium separation of the and atoms.
The function 
is a cutoff function that
goes smoothly from 1 to 0 as goes
from 
to 
.[37]
The values for and are 
nm and
nm.
The energy of the electrostatic interactions is
where 
and 
are the charges on the and atoms,
is the dielectric constant, and
is the separation between the and atoms.
The issue of what to use for is somewhat tricky.
In this case, a distance-dependent dielectric constant
has been used.
A more sophisticated approach exists[41]
but has not been implemented in Cerius.
The Coulomb interaction plays an important role in our simulation
because 5CB is a polar molecule.
The charges used were taken from the work of Wilson and
Allen[34] on simulations of the bulk properties of 5CB.
It is necessary to define a molecular axis direction,
, for each molecule
in order to compute the order parameter of the 5CB molecules.
In models that represent the molecules as elongated rigid objects with
rotational symmetry about the elongated direction,
there is no problem in identifying the
axis of a molecule; it is simply the direction of elongation.
In fully atomistic molecular dynamics simulations such as those done
here, the molecules are not rigid and so there is
ambiguity in specifiying the direction of each molecule.
The method chosen was
simply to use the bond direction of the bond connecting the two
benzene rings in 5CB as the direction of the molecule.
The tensor order parameter S can be defined[39] by its components
through the equation
where 
is a unit vector along the axis of a molecule,
and 
and 
are unit vectors along the space fixed axes,
where i,j=x,y,z.
The order parameter matrix S is diagonalizable.
Its largest eigenvalue is also refered to as the order parameter,
and the eigenvector corresponding to the largest eigenvalue is the director.
The brackets 
denote an average over a probability
density for .
Because directions are equivalent to the points on the unit sphere,
the probability density, f, for can be written as a
function of the polar angle 
and the azimuthal angle 
.
The probability density, 
, for can be interpreted
as either a time average for a single molecule or as an ensemble
average.
For our purposes the ensemble average is more suitable for simulations
run for a limited time, and so we make the approximation
The accuracy of this expression is limited by the modest size of the
greatest N that can feasibly be modeled.
It is thus important to be able to distinguish between an order
parameter that represents appreciable nematic ordering and
random fluctuations.
This question has been addressed by Eppenga and Frenkel,[42]
who derived an approximate expression for the average value of the
largest eigenvalue of S in an ensemble of N molecules of random
orientation.
Their conclusion was that the average value of the order parameter
should be 
.
We have refined this result by generating repeated ensembles of randomly
oriented molecules, and find a more accurate expression to be
.
This is illustrated in Figure 2, which shows the average values
of the largest eigenvalue of S and their standard deviation as
a function of N.
The dashed line is the postulated expression, , and
provides the baseline below which no significance can be attached to the
result of any simulation.