Printer Friendly

Molecular modelling--an enabling technology for chemical engineers.

This article briefly describes the basic concepts involved in the two most commonly used molecular modelling methods-molecular dynamics (MD) and Monte Carlo (MC). The methods are particularly useful for studying structures at the length scale of nanometre. Two examples (both are on the study of the miscibility of polyolefin blends) are used to illustrate the techniques. It is demonstrated that it is the nano-scaled structures formed by the segments of the constituent polyolefins that prevent them from mixing with each other. The examples also show that selection of specific method (MD or MC) depends on the nature of the problem in hand. In general, MC is more efficient than MD in terms of generating equilibrated structure while MD can provide information about the dynamics of a system. This is simply because MD requires the solution of equations of motion (a set of second order differential equations) while MC does not. Nonetheless, both methods need a reasonably accurate force field.

Cet article decrit brievement les concepts fondamentaux des deux methodes de modelisation les plus communement utilisees, soient la dynamique moleculaire (MD) et la technique Monte Carlo (MC). Ces methodes sont particulierement utiles pour l'etude des structures a l'echelle de longueur du nanometre. Deux exemples (les deux sur la miscibilite de melanges de polyolefines) permettent d'illustrer ces techniques. Il est demontre que ce sont les structures d'echelle nanometriques formees par les segments des polyolefines constituantes qui les empechent de se melanger entre elles. Ces exemples montrent egalement que le choix d'une methode specifique (MD ou MC) depend de la nature du probleme en main. En general, la MC est plus efficace que la MD pour ce qui est de produire une structure equilibree, tandis que la MD peut fournir de l'information sur la dynamique d'un systeme. Ceci est simplement du au fait que la MD necessite la resolution des equations de deplacement (une serie d'equations differentielles de second ordre) contrairement a la MC. Toutefois, les deux methodes requierent un champ de force raisonnablement precis.


With ever-changing needs for innovative materials exhibiting unique functional properties, engineering researchers have become increasingly involved in designing materials at the molecular level, which is usually considered to be in the domain of chemistry. This is especially relevant for engineers working in the emerging field of nanotechnology in which manipulation of atoms and/or molecules to create structures in the length scale of nanometre is required. To do so, understanding of the correlations between the chemical composition of the molecules that are used to construct the nanostructures and properties of such structures is most desired. Study of substances at the molecular level (i.e., nanoscience) is not new but utilizing the knowledge about the interactions between atoms and molecules to assemble them in a specific nano-structure that would exhibit the desired functional properties (i.e., nanotechnology) is new. It is believed that nanotechnology will benefit humankind tremendously even though there exist only a few such commercial products at the present time. And since nanotechnology is still in its infancy, most of these products have been developed from an ad hoc rather than a rigorous systematic scientific approach. In this regard, computational chemistry is most suited for such purposes. In fact, with the high costs associated with synthesizing and/or studying nano-structures, computational chemistry definitely offers a relatively less expensive alternative for the materials design purpose. The popularity of using the computational approach to investigate structure-property relationships for both macroscopic and/or nano-scaled systems is reflected by the increasing number of computational chemistry software vendors at various chemical engineering conferences worldwide, as well as of presentations that appear in one of the most attended chemical engineering conference-American Institute of Chemical Engineers (AIChE) annual meetings.


As shown in Figure 1, computational chemistry encompasses phenomena taking place in the time and length scales from femto second to microsecond and from picometre to micrometre (Hung et al., 2004). At different time and length scales, different computational chemistry techniques should be used. For example, if one is interested in the properties of a molecular system that depends on its electronic distribution (e.g. polarizability), one has to carry out quantum mechanical (i.e., ab initio) calculations. Nevertheless, in the domain of nanotechnology, atomistic molecular dynamics (MD) and Monte Carlo (MC) simulation would be most suitable as structures and properties of interest are usually not determined by the electronic motion. Rather, they are simply a function of the nuclear positions. In addition, most of these systems would contain a relatively large number of molecules (a few hundred to a few thousand), quantum mechanical methods are not feasible as the required computational resources are substantial. The purpose of this article is to give the reader a brief introduction to both MD and MC techniques along with two examples, one for each technique, to illustrate the concepts presented here. The examples were drawn from the author's own molecular modelling experience. To learn more about the simulation techniques, which have already developed to a fairly sophisticated state, and the types of molecular systems to which MD and MC methods could be applied, the reader is encouraged to consult molecular modelling textbooks such as those shown in Allen and Tildesley, 1990, Leach, 2001 and Cramer, 2004.


Figures 2 and 3 basically depict the essence of both MD and MC methods. Here, we use monatomic molecules with a number density much lower than a real liquid but higher than a real gas for illustrative purposes only. Obviously, the idea is applicable to polyatomic molecules at real density. As shown in Figure 2, at a finite temperature, the kinetic energy possessed by the molecules would exhibit translational motion in the given volume in such a manner that all the positions and velocities accessed by the molecules, termed as the phase space, would be consistent with the thermodynamic constraints imposed (i.e., the volume and temperature). It is worth noting that for polyatomic molecules, there are extra degrees of freedom for the rotational and vibrational motions. Each combination of positions and velocities represents a microstate. The equilibrium structure or macroscopic properties of the system (e.g. pressure) are simply a manifestation of the average of all accessible microstates. In MD simulation, such microstates are simply generated by solving the Newtonian equations of motion (F = ma) numerically. To do so, forces experienced by the molecules as a result of the presence of other molecules and the container's walls as well as the initial positions and velocities of the molecules are needed. As one can imagine, so long as the mathematical expressions and the associated parameters that are used to describe the forces, so-called the force field, are correct, accurate MD trajectories (i.e., structural and thermodynamic properties) can be obtained. In other words, the selection of a force field is of paramount importance in any MD simulation. There exist many different force fields and they can be made to be from very generic to very specific for high accuracy in predicting various molecular properties. The reader should consult chapter 2 of Cramer, 2004 for a complete list of force fields available. The generic force fields are used to describe a variety of chemical compounds including inorganic compounds, organic substances, metals and transition metals. The mathematical expressions used in such force fields tend to be rather simple to reduce the computational requirements at the expense of accuracy. For example, the intra-molecular energy terms of the force field that we used for the first example discussed below are known not to be so accurate.


The force fields that are developed for a specific class of materials (e.g. the AMBER force field was developed for biological molecules) usually consist of a set of much more complicated mathematical expressions for describing various energy terms. In particular, off-diagonal cross-coupling terms (i.e., terms used to describe the interaction between independent energy terms shown in generic force fields) and high-order energy terms are usually used. Once a force field has been selected and the initial positions and velocities have been assigned, an MD simulation can be started. Here, assignments of initial positions are somewhat arbitrary while initial velocities are done using the Maxwell Boltzmann distribution at the temperature of interest. So long as the simulation time is long enough to equilibrate the system, the choice of the initial positions is not too critical. However, if the simulation time is not long enough, initial positions have to be chosen with caution as the microstates generated by the MD simulation may not correspond to the equilibrium ones. For example, in the simulation of polymer molecules, the rotational isomeric state (RIS) model is commonly used to construct the initial conformations of the molecules since long MD simulations of such systems are practically infeasible (Flory, 1988).

The major drawback of the MD method, as alluded above, is that relatively short trajectories can be produced as computational requirements for solving the Newtonian equations of motion (i.e., a set of second order differential equations, one for each atom) are fairly substantial. In general, for a system that contains a few thousand atoms, a few days of CPU time is required to generate an MD trajectory on the order of one nanosecond on an 8 node parallel workstation cluster. In this regard, MC is a better method to obtain equilibrium structures as only positions of the molecules are required to compute. There are no equations of motion to solve. The reason is that the Hamiltonian (i.e., total energy) is separable (i.e., the kinetic and potential energies are independent of each other). Obviously, if one is interested in the evolution of an equilibrium structure or the dynamics of a molecular process, MD is the method of choice. In practice, the MC method generates equilibrium structures (configurations) by perturbing the current positions of the molecules by moving them into different positions as shown in Figure 3. If the energy of the new configuration is lower than that of the current configuration, the new one is accepted. However, if the energy is higher, the Boltzmann factor (i.e., [e.sup.-[DELTA]Ui[right arrow]i+1/kT]) will be calculated and compared with a random number, generated by the computer, in the range of 0 and 1 to determine whether the new configuration could be accepted or not. This criterion is termed as Metropolis criterion (Metropolis et al., 1953). The method is called Monte Carlo simply because random numbers are being generated throughout the whole simulation to determine the acceptance and rejection of configurations. In general, the perturbation to the configuration is carried out in such a manner that the acceptance rate is in the range of 20 to 50%. But in highly condensed systems such as polymer liquids, low acceptance rates, say 10%, may not be uncommon. Since energy is required to direct the generation of structures, force field is also required in MC.


In the following, two examples (both are on the study of the phase segregation of polyolefin blends) will be used to illustrate the techniques described above. It is interesting to note that in both cases, it is the structure in the length scale of nanometre that prevents the constituent polymer molecules from mingling with each other. The first example is on the MD study of the miscibility of polyethylene (PE)/isotactic polypropylene (iPP) blend in its liquid state. Since both PE and iPP are simple hydrocarbons interacting with each other through non-directional van der Waals interactions, one would expect that these two polymers would be miscible in the liquid state. However, small angle neutron scattering experiments clearly showed that the blend is immiscible at a temperature significantly above their melting temperatures. In order to understand the molecular origin of the observed phase separation, we applied the MD method to a model PE/iPP blend that consisted of only two molecules, one of each initially entangled together as shown in Figure 4. The MD simulation was carried out at 500 K (~ 230[degrees]C) using the Nose MD canonical scheme along with the Verlet leapfrog numerical algorithm and a generic force field Drieing 2.21 (please refer to the original article for the simulation details and the references for the algorithms used (Choi et al., 1995)). As can be seen in Figure 4, the molecules phase separate within a period of slightly less than 1 nanosecond (850 picoseconds). In particular, it was observed that the PE molecule prefers to form local order in the length scale of roughly 3 nanometres that prevents it from mingling with iPP. In fact, the local order revealed by the MD simulation was also observed experimentally for pure PE by X-ray diffraction experiment as well (Pieper and Kilian, 1993). It is speculated that this local order in the liquid state of PE is responsible for PE not mixing with most of the other polymers. The final configuration shown in Figure 4 corresponds to a similar MD simulation that was carried out over a period of 1 nanosecond at the same temperature but with an addition of an ethylene-propylene block copolymer. It is clearly seen in the figure that all three molecules mingled with each other with no sign of phase segregation-at least over the simulation time used. This is consistent with the experimental practice that adding copolymers to immiscible polymer blends is one of the methods for enhancing the adhesion between two immiscible phases.


In the second example, the MC method was used. The results show that phase segregation could also occur to a polymer blend composed of polypropylenes with different tacticity at low temperatures. In particular, blends of isotactic and atactic polypropylenes (iPP and aPP) were studied. Since the tacticity effect is rather subtle and the simulation temperature is low (125[degrees]C in this particular case, a temperature well below the melting temperature of iPP but not aPP), the MD method is not suitable as it can not equilibrate the system using a reasonable computational time, therefore, the MC method was used. In fact, with the use of the MC method, 25 million structures were created (it would take almost a year for MD to create the same number of structures) and the average structure, as characterized by various intermolecular pair correlation functions, shows that it is the nano-scaled structure formed by aPP that is responsible for the phase separation at low temperatures before the crystallization of iPP takes place, which was observed by Keith and Padden almost four decades ago (Choi and Mattice, 2004). Our simulation results show that it is the syndiotactic sequences along the backbone of aPP that cause the formation of such nano-scaled structures. In this simulation, the blend was modelled using 18 PP chains of which 9 were isotactic and the other 9 atactic. They were assembled in simulation boxes with different sizes (i.e., different pressures) subjected to periodic boundary conditions. The periodic boundary conditions were simply applied to mimic the bulk liquid state of the materials. Rather than assigning the initial positions of the molecules in the continuum space as was done in the MD simulation of the PE/iPP blend, they are randomly assigned to the lattice sites of a high coordination diamond lattice with a coordination number of 12. It is worth noting that performing MC simulation on a lattice can further reduce the computational requirements. Please refer to the original article for the simulation details including the selection of the force field (Choi and Mattice, 2004). Since the simulation was carried out on a lattice, both bond length and bond angles are fixed and the only intramolecular energy term considered was the torsional energy. Obviously, van der Waals interaction between the molecules was also included.


In summary, we have presented the basic concepts involved in the MD and MC methods and have demonstrated, using some of the author's previous molecular modelling results, that the methods are fairly powerful for revealing nano-scaled structures that are formed by molecules with different chemical compositions. Also illustrated, the selection of specific method (either MD or MC) depends on the problem in hand.


The author would like to thank the editor, Kumar Nandakumar, for the opportunity to compose this article. Thanks are also extended to Henk Blom (Baxter Healthcare Corporation), Tom A. Kavassalis (Xerox Research Centre of Canada) and Alfred Rudin (University of Waterloo) for calling attention to performing the MD simulation of the PE/iPP blend with the presence of an ethylene-propylene copolymer molecule, which has never been published before.


Allen, M. P. and D. J. Tildesley, "Computer Simulation of Liquids," Oxford University Press, New York, NY (1990).

Choi, P., H. P. Blom, T. A. Kavassalis and A. Rudin, "The Immiscibility of Poly(ethylene) and Poly(propylene): A Molecular Dynamics Study," Macromolecules 28, 8247-8250 (1995).

Choi, P. and W. L. Mattice, "Molecular Origin of Demixing, Prior to Crystallization, of Atactic Polypropylene/Isotactic Polypropylene Blends upon Cooling from the Melt," J. Chem. Phys. 121, 8647-8651 (2004).

Cramer, C. J., "Essentials of Computational Chemistry," 2nd ed., Wiley, Hoboken, NJ (2004).

Flory, P. J., "Statistical Mechanics of Chain Molecules," Oxford University Press, New York, NY (1988).

Hung, F. R., S. Franzen and K. E. Gubbins, "A Graduate Course on Multi-Scale Modelling of Soft Matter," Chem. Eng. Ed. 38(4), 242-249 (2004).

Leach, A. R., "Molecular Modelling: Principles and Applications," 2nd ed., Prentice Hall, New York, NY (2001).

Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, "Equation of State Calculations by Fast Computing Machines," J. Chem. Phys. 21, 1087-1092 (1953).

Pieper, T. and H. G. Kilian, "Packing of Chain Segments--A Method for Describing X-Ray Patterns of Crystalline, Liquid-Crystalline and Noncrystalline Polymers," Adv. Polym. Sci. 108, 49-89 (1993).

Manuscript received March 2, 2006; revised manuscript received May 4, 2006; accepted for publication May 4, 2006.

Phillip Choi

Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB, Canada T6G 2G6

* Author to whom correspondence may be addressed. E-mail address:
COPYRIGHT 2006 Chemical Institute of Canada
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2006 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Choi, Phillip
Publication:Canadian Journal of Chemical Engineering
Geographic Code:1CANA
Date:Jun 1, 2006
Previous Article:Erratum.
Next Article:Experimental investigations of regimes of bubble formation on submerged orifices under constant flow condition.

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters |