Molecular Dynamics Simulation of Liquid Crystals in the Nematic Phase

Introduction

Molecular dynamics is a continually evolving computer simulation technique that has been used to model the behavior of many different systems of atoms. Given a model potential describing the atomic interactions and initial conditions, it is possible to calculate the position and momenta of all of the atoms (ranging from a couple hundred to hundreds of thousands) in simulation for a given period of time by numerically integrating the Newtonian equations of motion.

The goal of this project is to construct a working molecular dynamics simulation of a liquid crystal in the nematic phase. This has been in the past using an anisotropic, interatomic potential called the Gay-Berne potential, this model potential has been used to study the nematic-isotropic interface and elastic constants of the simulated liquid crystal among other properties.

Theory

LIQUID CRYSTALS

A liquid crystal has properties of both a liquid and a crystal. A liquid crystal can flow and form droplets, but it also has anisotropic properties and differently levels of order in the orientation and position of the constituent molecules. Molecules that make up a liquid crystal tend to have a long axis and a short axis instead of a spherical shape. The three main phases of liquid crystals are isotropic, nematic, and smectic. In an isotropic phase, there exists no order in the orientation or position. The nematic phase is characterized by the molecules having orientational order and aligning along a vector aptly named the director. The smectic phases displays both orientational order and also order in the positions of the center of masses of the particles. The isotropic phase occurs at higher temperatures than the nematic phase, and the nematic phase occurs at higher temperatures than smectic phases. The goal of this project is to successfully simulate a nematic phase using parameters that fit the behavior of para-Azoxyanisole (PAA). It is important to note that while the geometry and well depth are chosen based on PAA, the Gay-Berne potential is not simulating PAA at an atomistic level.

MOLECULAR DYNAMICS

Three items are necessary in order to construct a molecular dynamics simulation. To begin the simulation, initial conditions for the system must be specified. Next, a model potential is needed to describe the interaction of a particle with its neighbors. The final ingredient is a numerical method for integrating the equations of motion for each particle. This integration should use a particle's current position and momentum, along with the calculated forces acting upon it, to advance the particle to a new position by a specified time step.

The Gay-Berne potential is one of the most widely used potentials for molecular dynamics simulations of liquid crystals. As such, it is the one used in the project.

Once the forces and torques have been calculated using the Gay-Berne potential, a leapfrog integration algorithm is used to advance the positions and orientations of the particle by a time step.

In order to initialize the system the particles are arranged in a cubic lattice to prevent large repulsion due to overlapping. Each particle is given a random velocity and a scaling factor is calculated in order to set the correct temperature.

Experimental Methods

The simulation will be performed using 256 particles. The Gay-Berne potential used to simulate the liquid crystal will form a nematic phase at a density of 0.32 and temperatures below 1.7 (both given in reduced units). Using reduced units the mass of the particles, side-to-side length and well depth parameter are set to unity and other quantities are adjusted accordingly.

While the simulation is running, the position and orientation of each particle for each time step is written to a data file. Additionally, macroscopic quantities such as the potential and kinetic energies, the pressure, and the temperature of the system can be recorded. The simulation will be examined for correct behavior by using scientific visualization, calculating a scalar order parameter, and calculating the radial distribution function for the system.

The shape and orientation of the particles will need to be represented in order to determine that the liquid crystal phase is correct. This visualization will be created using a program called gnuplot. Although gnuplot does not have the capability to produce the visualization of the ellipsoids in the simulation, it is able to draw the molecular axes in 3-D in a simple manner.

PAIR CORRELATION FUNCTION

The radial distribution function, or pair correlation function, describes the distribution of particles in a system with respect to the random distribution of an ideal gas. Mathematically,

In a liquid, molecules pack into spherical shells. As such, it is expected that there will be peaks in the radial distribution function representing the locations of these shells with respect to a molecule. The radial distribution function can be measured experimentally using x-ray scattering and also calculated theoretically.

ORDER PARAMETER

In order to quantify the amount of order in the system, the scalar order parameter

will be used. Theta is the angle of the molecular axis with respect to another angle. This angle is usually taken to be the director (the average orientation of the molecules). In the simulation S will be calculated with respect to each of the coordinate axes instead due to the difficulty in calculating the director. For a liquid crystal, 0.3<S<0.8. For an isotropic phase S will be zero and for perfect alignment one.

Results

CONFIRMING A LIQUID STATE USING THE PAIR CORRELATION FUNCTION

There are peaks in the radial distribution function indicating a liquid phase. The next step is to check that there is ordering in the system which will confirm that a nematic phase has been successfully simulated.

QUANTIFYING ORDER IN THE SYSTEM

As seen in above the system exhibits orientational ordering but not positional ordering. This is characteristic of a nematic phase. After multiple trials, it was found that the simulation preferred to along either the z-direction or in some direction close to the z-direction. In Figure 9, there appears to be ordering in the system, however, Sx and Sy do not take on values of ~(-0.5) which indicates that the molecules are not completely perpendicular to the xy-plane. The tendency of the system to align in a direction close to the z-axis should not exist as no preference has not been encoded into the simulation. It is likely that this problem stems from the calculation of the torques from the potential. This part of the code is most complex and most prone to error.

Conclusion

Despite experiencing technical difficulties, an order phase of liquid crystals was simulated using molecular dynamics, and the goal of the project was reached to some degree. The calculated radial distribution function matches that of a liquid even with the geometry of the Gay-Berne molecules and the orientational ordering. A scalar order parameter was calculated in order to quantify the orientational ordering in the system, and it was found that it increased as the temperature increased as expected. The next steps in continuing the project would be to fix the problem with the reduced units in order to obtain results that are consistent with the literature.