**Part 5**

**Dynamics of Nanomachines** 

20 Will-be-set-by-IN-TECH

338 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Thomas, J., McGaughey, A. & Kuter-Arnebeck, O. (2010). Pressure-driven water flow through

Werder, T., Walther, J., Jaffe, R., Halicioglu, T. & Koumoutsakos, P. (2003). On the water-carbon

Whitby, M., Cagnon, L., Thanou, M. & Quirke, N. (2008). Enhanced fluid flow through

Ybert, C., Barentin, C., Cottin-Bizonne, C., Joseph, P. & Bocquet, L. (2007). Achieving large

Zhang, H., Zhang, Z. & Ye, H. (2011). Molecular dynamics-based prediction of boundary slip

Zhu, F., Tajkhorshid, E. & Schulten, K. (2004). Collective diffusion model for water permeation

nanotubes, *The Journal of Physical Chemistry B* 107(6): 1345–1352.

of fluids in nanochannels, *Microfluidics and Nanofluidics* pp. 1–9.

through microscopic channels, *Physical review letters* 93(22): 224501.

nanoscale carbon pipes, *Nano letters* 8(9): 2632–2637.

*Journal of Thermal Sciences* 49(2): 281–289.

*fluids* 19: 123601.

carbon nanotubes: Insights from molecular dynamics simulation, *International*

interaction for use in molecular dynamics simulations of graphite and carbon

slip with superhydrophobic surfaces: Scaling laws for generic geometries, *Physics of*

## **Analysis of the Atomization Process by Molecular Dynamics Simulation**

Yeh Chun-Lang

*Department of Aeronautical Engineering, National Formosa University Taiwan, ROC* 

### **1. Introduction**

Investigation of the atomizer flows has received a lot of attentions due to their broad applications in many fields, e.g. industrial processes, agriculture, meteorology, medicine, etc. Historically, less attention has been devoted to the internal flow development in atomizers because of the small size of practical atomizers, which makes measurements difficult. However, internal flows in atomizers are of interest due to their potential effect on the subsequent atomization process that takes place external to the atomizers.

Owing to the difficulties inherent in the simulation of the atomizer internal and external flows, no systematic numerical researches have been carried out on this subject. One way to study internal flows in atomizers is to conduct experimental investigations by the use of large-scale models. Another approach is to do detailed numerical simulations of the flow. An advantage of the numerical simulation is that practical-scale atomizers can be studied as easily as large scale ones. The methods of analyzing atomizer flows can be divided into two categories, namely, the macroscopic analysis and the microscopic analysis, each with its own advantages and disadvantages. The former approach is usually adopted as an aided tool for atomizer design due to its relatively lower computational cost. However, when atomization occurs, the liquid breaks into ultra-fine nano-scale droplets. Under such situation, the Navier-Stokes equation based on continuum concept is not suitable and the microscopic analysis should be used for subsequent simulation of the atomization process. In the past decade, the author has carried out a series of researches on the macroscopic analysis of the atomizer internal and external flows (Yeh, 2002, 2003, 2004, 2005, 2007). The macroscopic behavior and characteristics of the atomizer flows, such as discharge coefficient, spray angle, film thickness, atomizer geometry, etc. were analyzed. However, as mentioned above, the macroscopic analysis is not suitable when atomization occurs. The author's previous studies were therefore initially confined to the investigation of the early stage of the atomization process, during which the liquid had not yet ruptured to ultra-fine nano-scale droplets. In order to improve this inherent deficiency of the macroscopic analysis, the author has employed molecular dynamics (MD) simulation to investigate the microscopic evolution of the atomizer internal and external flows. In the past few years, the author has studied the atomization processes of a nano-scale liquid thread and a nanojet. Liquid thread and nanojet are two of the most fundamental and important phenomena during atomization process. Although the size of a nano-scale atomizer is much smaller than that of a practical-scale atomizer, analysis of the nano-scale atomizer flow can provide a preliminary understanding of the detailed process of a practical-scale atomizer flow, especially when the available computing resources can not directly facilitate the microscopic analysis of a practical-scale atomizer flow.

### **2. Analysis by the molecular dynamics simulation**

From one of the author's previous studies (Yeh, 2005) for atomizer flow by macroscopic analysis, it can be seen that the liquid evolves into threads after leaving the atomizer. In the author's recent papers (Yeh, 2009a, 2009b, 2010), the atomization processes of the nano-scale liquid thread and nanojet have been discussed in detail by molecular dynamics simulation.

#### **2.1 Molecular dynamics analysis of the instability for a nano-scale liquid thread**

In the scope of fluid flow researches, the instability has received much attention due to its potential influence on the flow development. Previous instability analyses were focused mainly on large scale flow fields. The applicabilities of these theories to micro- or nano-scale flow fields are still uncertain. Lord Rayleigh (1879) studied the instability of cylindrical thin films. He analyzed an inviscid liquid cylinder and a viscous one. Later on, Weber (1931) and Tomotika (1935) considered more realistic cases of liquid threads in unbounded domains. Goren (1962) studied the annular liquid films in contact with a solid, i.e., supported on a wire or lining the interior walls of a capillary. He used linear stability analysis to determine the fastest growing mode when either inertia or viscous forces are negligible. Koplik and Banavar (1993) studied the Rayleigh's instability of a cylindrical liquid thread in vacuum by three-dimensional MD simulation. The maximum number of molecules they used consists of 8,192 liquid argon Lennard-Jones molecules for a cylindrical liquid thread of nondimensionalized radius of 7.5 in a box of non-dimensionalized length of 54.7. For this simulation condition, only one liquid particle was formed. If a smaller computational domain was used instead, no liquid particle was found from their study. Kawano (1998) applied 10,278 Lennard-Jones molecules of liquid and vapor coexisting argon in three dimensions to analyze the interfacial motion of a cylindrical liquid thread of nondimensionalized radii of 2.0 to 4.0 in a box of non-dimensionalized length up to 120. For this condition of larger computational domain, a maximum number of 8 to 9 liquid particles were observed. Min and Wong (2006) studied the Rayleigh's instability of nanometer scale Lennard-Jones liquid threads by MD simulation and concluded that Rayleigh's continuum prediction holds down to the molecular scale. Kim, Lee, Han and Park (2006) applied MD simulation to investigate the thermodynamic properties and stability characteristics of the nano-scale liquid thread. They found that the overall trends of the simulation results agree with the classical stability theory. However, the classical theory overpredicts the region of stable domain compared to the MD results as the radius decreases.

In the following discussion, the instability of a liquid thread is investigated by MD simulation. The influences of liquid thread radius, fundamental cell length, and temperature will be discussed. Snapshots of molecules, number of liquid particles formed, and density field are analyzed. Two linear stability criteria, namely Rayleigh's stability criterion and Kim's stability criterion, are accessed for their validity in molecular scale. This can provide insights into the mechanism and prediction of the atomization process.

#### **2.1.1 Molecilar dynamics simulation method**

342 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

that of a practical-scale atomizer, analysis of the nano-scale atomizer flow can provide a preliminary understanding of the detailed process of a practical-scale atomizer flow, especially when the available computing resources can not directly facilitate the microscopic

From one of the author's previous studies (Yeh, 2005) for atomizer flow by macroscopic analysis, it can be seen that the liquid evolves into threads after leaving the atomizer. In the author's recent papers (Yeh, 2009a, 2009b, 2010), the atomization processes of the nano-scale liquid thread and nanojet have been discussed in detail by molecular dynamics simulation.

In the scope of fluid flow researches, the instability has received much attention due to its potential influence on the flow development. Previous instability analyses were focused mainly on large scale flow fields. The applicabilities of these theories to micro- or nano-scale flow fields are still uncertain. Lord Rayleigh (1879) studied the instability of cylindrical thin films. He analyzed an inviscid liquid cylinder and a viscous one. Later on, Weber (1931) and Tomotika (1935) considered more realistic cases of liquid threads in unbounded domains. Goren (1962) studied the annular liquid films in contact with a solid, i.e., supported on a wire or lining the interior walls of a capillary. He used linear stability analysis to determine the fastest growing mode when either inertia or viscous forces are negligible. Koplik and Banavar (1993) studied the Rayleigh's instability of a cylindrical liquid thread in vacuum by three-dimensional MD simulation. The maximum number of molecules they used consists of 8,192 liquid argon Lennard-Jones molecules for a cylindrical liquid thread of nondimensionalized radius of 7.5 in a box of non-dimensionalized length of 54.7. For this simulation condition, only one liquid particle was formed. If a smaller computational domain was used instead, no liquid particle was found from their study. Kawano (1998) applied 10,278 Lennard-Jones molecules of liquid and vapor coexisting argon in three dimensions to analyze the interfacial motion of a cylindrical liquid thread of nondimensionalized radii of 2.0 to 4.0 in a box of non-dimensionalized length up to 120. For this condition of larger computational domain, a maximum number of 8 to 9 liquid particles were observed. Min and Wong (2006) studied the Rayleigh's instability of nanometer scale Lennard-Jones liquid threads by MD simulation and concluded that Rayleigh's continuum prediction holds down to the molecular scale. Kim, Lee, Han and Park (2006) applied MD simulation to investigate the thermodynamic properties and stability characteristics of the nano-scale liquid thread. They found that the overall trends of the simulation results agree with the classical stability theory. However, the classical theory overpredicts the region of

**2.1 Molecular dynamics analysis of the instability for a nano-scale liquid thread** 

analysis of a practical-scale atomizer flow.

**2. Analysis by the molecular dynamics simulation** 

stable domain compared to the MD results as the radius decreases.

insights into the mechanism and prediction of the atomization process.

In the following discussion, the instability of a liquid thread is investigated by MD simulation. The influences of liquid thread radius, fundamental cell length, and temperature will be discussed. Snapshots of molecules, number of liquid particles formed, and density field are analyzed. Two linear stability criteria, namely Rayleigh's stability criterion and Kim's stability criterion, are accessed for their validity in molecular scale. This can provide In this study, the vaporization process of a liquid thread is investigated by MD simulation. The inter-atomic potential is one of the most important parts of MD simulation. Many possible potential models exist, such as hard sphere, soft sphere, square well, etc (Haile, 1992). In this research, the Lennard-Jones 12-6 potential model, which is widely used, is adopted for calculation. It is

$$\phi(\mathbf{r}) = 4\varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right] \tag{1}$$

where *r* denotes the distance between two molecules, *ε* and *σ* are the representative scales of energy and length, respectively. The Lennard-Jones fluid in this research is taken to be argon for its ease of physical understanding. The parameters for argon are as follows (Kawano, 1998) : the length parameter *σ*=0.354 nm, the energy parameter *ε/k*B=93.3K, and the molecular weight *m*=6.64×10-26 kg, where *k*B=1.38×10-23 J/K denotes the Boltzmann constant. The cut-off radius *rc* beyond which the intermolecular interaction is neglected is 5.0*σ*.

The simulation domain is schematically shown in Fig.1, with periodic boundary conditions applied in all three directions. Simulation domain dimensions, temperatures and number of molecules, are listed in Table 1, together with some simulation results. The time integration of motion is performed by Gear's fifth predictor-corrector method (Haile, 1992) with a time step of ᇞ*t*\*=0.001(i.e. 2.5 fs). Note that all quantities with an asterisk in this paper, such as L\*, *R*\*, ρ\*, ᇞ*t*\*, etc., are non-dimensionalized in terms of *σ*, *ε*, and *m*, i.e. L\*=L/*σ*, *R*\*=*R*/*σ*, ρ\*=*Nσ*3/V, ᇞ*t*\*=ᇞ*t* (*ε*/*m*)1/2/*σ*, *T*\*= *k*B*T*/*ε*.

Fig. 1. Illustration of the computational domain for the simulation of a liquid thread


Table 1. Simulation domain dimensions, temperatures, number of molecules, and simulation results

In this research, a cylindrical liquid thread of length L\* and radius *R*\* is placed at the center of the computational domain and the remaining space is vacuum. The initial density of the liquid argon is ρL\*=0.819. The system temperature is kept at *T*\*=0.75 or 1.0. These dimensionless values correspond to ρL=1223 kg/m3 for argon and *T*=70K or 93.3K. Note that these two temperatures are below the critical temperature(150K)of argon.

The procedure for MD simulation includes three stages : initialization, equilibration and production. Initially, equilibration is performed for liquid argon molecules in a rectangular parallelepiped with length and width equal to the liquid thread diameter (*D*\*=2*R*\*=4, 6, 8) and with height equal to the side length of the computational domain (L\*=10, 16, 20, 24, 30, 60, 120, 240, 480). The initial velocities of molecules are decided by normal random numbers.

Table 1. Simulation domain dimensions, temperatures, number of molecules, and simulation

In this research, a cylindrical liquid thread of length L\* and radius *R*\* is placed at the center of the computational domain and the remaining space is vacuum. The initial density of the liquid argon is ρL\*=0.819. The system temperature is kept at *T*\*=0.75 or 1.0. These dimensionless values correspond to ρL=1223 kg/m3 for argon and *T*=70K or 93.3K. Note that

The procedure for MD simulation includes three stages : initialization, equilibration and production. Initially, equilibration is performed for liquid argon molecules in a rectangular parallelepiped with length and width equal to the liquid thread diameter (*D*\*=2*R*\*=4, 6, 8) and with height equal to the side length of the computational domain (L\*=10, 16, 20, 24, 30, 60, 120, 240, 480). The initial velocities of molecules are decided by normal random numbers.

these two temperatures are below the critical temperature(150K)of argon.

results

Np

Case No. L\* *R*\* T\* Nmol *f*

Velocity rescaling is performed at each time step by Eq.(2) to make sure that the molecules are at the desired temperature *T*\* :

$$
\upsilon\_i^{new} = \upsilon\_i^{old} \sqrt{\frac{T\_D}{T\_A}} \tag{2}
$$

where *vi new* and *vi old* are the velocities of molecule *i* after and before correction, respectively, and *TD* and *TA* are the desired and the actual molecular temperatures, respectively. The liquid molecules are equilibrated for 106 time steps at the desired temperature *T*\*. The achievement of equilibrium state is confirmed by obtaining the radial distribution function. After the liquid molecules are equilibrated, the rectangular parallelepiped for the liquid molecules is truncated to the desired cylindrical liquid thread by removing unwanted regions. The cylindrical liquid thread then is put into the computational domain and the production stage proceeds. A minimum image method and the Verlet neighbor list scheme (Haile, 1992) to keep track of which molecules are actually interacting at a given time interval of 0.005 are used in the equilibration and the production stages.

#### **2.1.2 Liquid thread vaporization process**

In the following discussion, a cylindrical liquid thread of length L\* and radius *R*\* is placed at the center of the simulation domain and the remaining space is in vacuum, as illustrated in Fig.1. Simulation domain dimensions, temperatures and number of molecules, are listed in Table 1.

Figures 2~4 show the vaporization processes of liquid threads at T\*=0.75(cases 1~23 in Table 1). Note that L\*=120 and *R*\*=3 correspond to L=42.5nm and *R*=1.06nm, respectively. It is found that when *R*\*=2, the liquid thread breaks up into drops for all fundamental cell lengths. For L\*=10, 20, 30 and 60, the liquid thread ruptures only from its two ends, i.e. the top and bottom surfaces of the fundamental cell, and gets shorter due to the contraction motion in its axial direction. Only one liquid particle is formed. On the other hand, for L\*=120, 240 and 480, the thread ruptures not only from the top and bottom surfaces of the fundamental cell but also from its interior section. The number of liquid particles produced for L\*=120, 240 and 480 is 2, 5 and 10, respectively. During the vaporization process, collision and coalescence of the liquid particles occur. The two liquid particles produced for L\*=120 coalesce into one liquid particle. The number of liquid particles at t\*=1000 for L\*=240 and 480 is 3 and 5, respectively.

**R\*=2, L\*=30, T\*=0.75**

**R\*=2, L\*=60, T\*=0.75**



X Y

Z

X Y

Z



0

10

**t\*=500**

**R\*=2, L\*=30, T\*=0.75**

**t\*=500**

**R\*=2, L\*=60, T\*=0.75**

**R\*=2, L\*=20, T\*=0.75**

**t\*=600**

10


0




0

10

X Y

Z



X Y

Z

X Y

Z



0

10

**R\*=2, L\*=20, T\*=0.75**

**t\*=100**

**R\*=2, L\*=30, T\*=0.75**

**t\*=100**

**R\*=2, L\*=60, T\*=0.75**

**t\*=100**

10





0

0

10

X Y

Z

Fig. 2. Vaporization processes of liquid threads of *R*\*=2 at T\*=0.75

**R\*=3, L\*=10, T\*=0.75**

**R\*=3, L\*=16, T\*=0.75**




0

0

5



0

X Y





0

15

15

**R\*=3, L\*=30, T\*=0.75**

0

15

X Y

Z

X Y

Z

0

10

10

**R\*=3, L\*=20, T\*=0.75**

**t\*=500**

X Y

Z


Z

5

**R\*=3, L\*=10, T\*=0.75**

**t\*=500**

**R\*=3, L\*=16, T\*=0.75**

**t\*=500**

**t\*=500**

0





0

0

0

10

0

5







0 15

0

0

10

0

X Y


X Y

Z





0 15

0

15

X Y

Z

X Y

Z

0

10

10

**R\*=3, L\*=20, T\*=0.75**

**t\*=100**

**R\*=3, L\*=30, T\*=0.75**

Z

5

**R\*=3, L\*=10, T\*=0.75**

**t\*=100**

**R\*=3, L\*=16, T\*=0.75**

**t\*=100**

**t\*=100**

0

0

5

**R\*=3, L\*=60, T\*=0.75**

**R\*=3, L\*=120, T\*=0.75**

**R\*=3, L\*=240, T\*=0.75**

**R\*=3, L\*=60, T\*=0.75**

**R\*=3, L\*=120, T\*=0.75**

**R\*=3, L\*=240, T\*=0.75**

Fig. 3. Vaporization processes of liquid threads of *R*\*=3 at T\*=0.75

When *R*\*=4, the liquid thread remains intact for L\*=10 and 20. For L\*=24, 30, 60, 120 and 240, the liquid thread ruptures merely from its two ends and only one liquid particle is formed. For L\*=480, the thread ruptures not only from the top and bottom surfaces of the fundamental cell but also from its interior section. Three liquid particles are produced. The liquid particles eventually coalesce into one liquid particle.

**t\*=100**

**R\*=4, L\*=24, T\*=0.75**

**t\*=100**

**R\*=4, L\*=30, T\*=0.75**

**t\*=100**

**R\*=4, L\*=60, T\*=0.75**

**t\*=100**

**R\*=4, L\*=120, T\*=0.75**



Z

0


0


0 50

50

Z


0 15

0

15

Z

Z

0 10






0 50

0 15

0 10

350 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

When *R*\*=4, the liquid thread remains intact for L\*=10 and 20. For L\*=24, 30, 60, 120 and 240, the liquid thread ruptures merely from its two ends and only one liquid particle is formed. For L\*=480, the thread ruptures not only from the top and bottom surfaces of the fundamental cell but also from its interior section. Three liquid particles are produced. The





0

0

10

0

0 200 -200

0




X Y

Z

0

10

10

**R\*=4, L\*=20, T\*=0.75**

0

5

**t\*=900**

**t\*=900**

5

**R\*=4, L\*=10, T\*=0.75**

0

5

X Y

Z

0 200

**t\*=900**

**R\*=3, L\*=480, T\*=0.75**

200

X Y

Z


0


0 200

Fig. 3. Vaporization processes of liquid threads of *R*\*=3 at T\*=0.75



X Y

Z

0 10

0

5

**t\*=100**

**R\*=4, L\*=20, T\*=0.75**

**t\*=100**

5

**R\*=4, L\*=10, T\*=0.75**




0 10

0

10

0

0

5

X Y

Z

liquid particles eventually coalesce into one liquid particle.

**t\*=400**

**R\*=3, L\*=480, T\*=0.75**


0 200 200

X Y

Z



0

50

Z

Z

Fig. 4. Vaporization processes of liquid threads of *R*\*=4 at T\*=0.75

Figure 5 shows the vaporization processes for liquid threads of *R*\*=4 at T\*=1.0(cases 24~30 in Table 1). Unlike the liquid threads of the same radius *R*\*=4 at a lower temperature T\*=0.75(cases 16~23 in Table 1), the liquid thread remains intact only for L\*=10. For L\*=20, 30, 60 and 120, the liquid thread ruptures from its two ends, i.e. the top and bottom surfaces of the fundamental cell, and only one liquid particle is formed. For L\*=240 and 480, the thread ruptures not only from the top and bottom surfaces of the fundamental cell but also from its interior section. The number of liquid particles produced for L\*=240 and 480 is 2 and 4, respectively. During the vaporization process, collision and coalescence of the liquid particles occur. The liquid particles eventually coalesce into one liquid particle.





Figure 5 shows the vaporization processes for liquid threads of *R*\*=4 at T\*=1.0(cases 24~30 in Table 1). Unlike the liquid threads of the same radius *R*\*=4 at a lower temperature T\*=0.75(cases 16~23 in Table 1), the liquid thread remains intact only for L\*=10. For L\*=20, 30, 60 and 120, the liquid thread ruptures from its two ends, i.e. the top and bottom surfaces of the fundamental cell, and only one liquid particle is formed. For L\*=240 and 480, the thread ruptures not only from the top and bottom surfaces of the fundamental cell but also from its interior section. The number of liquid particles produced for L\*=240 and 480 is 2 and 4, respectively. During the vaporization process, collision and coalescence of the liquid particles occur. The liquid particles eventually

0


0 200

200

0 200

0 100 -100




X Y

Z

0 200

X Y

Z

0 200

0 100

**t\*=500**

**R\*=4, L\*=240, T\*=0.75**

**t\*=200**

**R\*=4, L\*=480, T\*=0.75**

**t\*=900**

**R\*=4, L\*=480, T\*=0.75**

0

100

X Y

Z




X Y

Z


0 200

Fig. 4. Vaporization processes of liquid threads of *R*\*=4 at T\*=0.75

0 100

0

100

X Y

Z


0 100

**t\*=100**

**R\*=4, L\*=240, T\*=0.75**

**t\*=900**

**R\*=4, L\*=240, T\*=0.75**

**t\*=500**

**R\*=4, L\*=480, T\*=0.75**




0


coalesce into one liquid particle.

0 200

200

0 100

0 100 0

100

X Y

Z

**R\*=4, L\*=60, T\*=1.0**

**R\*=4, L\*=10, T\*=1.0**

**R\*=4, L\*=20, T\*=1.0**

**R\*=4, L\*=30, T\*=1.0**

**R\*=4, L\*=60, T\*=1.0**

Fig. 5. Vaporization processes of liquid threads of *R*\*=4 at T\*=1.0

From the above discussion for Figs.2~5, the following observations can be obtained.

First, if the fundamental cell length is long enough, more than one liquid particle may be produced in the cell. Under such situation, the liquid thread ruptures not only from its two ends but also from the interior of the fundamental cell. On the other hand, if the fundamental cell length is small, the liquid thread may remain intact or produce only one liquid particle in the cell. If the thread breaks up, it ruptures from its two ends only, i.e. the top and bottom surfaces of the fundamental cell, but not from its interior.

Second, a liquid thread with a longer fundamental cell length may produce more liquid particles in the cell.

Third, a thinner liquid thread may produce more liquid particles in the cell.

Fourth, a liquid thread at a higher temperature may produce more liquid particles.

#### **2.1.3 Stability analysis**

354 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules




0 200



0 200

**t\*=900**

0

200

X Y

Z

**R\*=4, L\*=480, T\*=1.0**

0 100 -50 0 50


0




0 200

0

200

X Y

Z

0 100

100

X Y

Z

X Y

Z

**t\*=900**

**R\*=4, L\*=120, T\*=1.0**

**t\*=900**

**R\*=4, L\*=240, T\*=1.0**

**t\*=500**





0





Fig. 5. Vaporization processes of liquid threads of *R*\*=4 at T\*=1.0

0 200

**R\*=4, L\*=480, T\*=1.0**

0 200

0

200

X Y

Z

0 100

100

X Y

Z

0 50

**t\*=100**

**R\*=4, L\*=120, T\*=1.0**

**t\*=100**

**R\*=4, L\*=240, T\*=1.0**

**t\*=180**




0 200

**R\*=4, L\*=480, T\*=1.0**

0 100

0 50 0

50

X Y

Z

According to the classical theory by Rayleigh (1879), the radius of the liquid thread, *R*, is related to the critical wavelength of perturbation, *λc* , as

$$
\lambda\_c = 2\pi R \tag{3}
$$

A liquid thread will break up into drops if the axial wavelength of the surface perturbation *L*>*λc*. If *L*<*λc*, the thread is stable and will remain intact. Owing to the periodic boundary conditions in this study, the fundamental cell size *L*\* can be regarded as the longest wavelength of the perturbation.

For *R*\*=2, *λ<sup>c</sup>* \* is 12.6. According to Rayleigh's stability criterion, a liquid thread of radius *R*\*=2 and length *L*\*=10 is expected to be stable and remain intact. However, from Fig.21 or Table 1 (Case 1) , it is seen that the liquid thread of *R*\*=2 and *L*\*=10 is unstable, which implies that Rayleigh's stability criterion overpredicts the stable domain as compared to the MD simulation results. Rayleigh's stability criterion holds for the remaining liquid threads with *L*\*>*λ<sup>c</sup>* \* (*L*\*=20~480) , which are unstable and break up into drop(s).

For *R*\*=3, *λ<sup>c</sup>* \* is 18.9. From Fig.3 or Table 1 (Cases 8 and 9) , the liquid thread of *R*\*=3 and *L*\*=10 remains intact while the liquid thread of *R*\*=3 and *L*\*=16 is unstable. Rayleigh's stability criterion holds for the former case but violates the MD simulation results for the latter case. The remaining liquid threads with *L*\*>*λ<sup>c</sup>* \* (*L*\*=20~480) are unstable and break up into drop(s). Rayleigh's stability criterion holds for these cases.

For *R*\*=4, *λ<sup>c</sup>* \* is 25.2. From Fig.4 or Table 1 (Cases 16~23) , the liquid threads of *L*\*=10 and 20 remain intact while the liquid thread of *L*\*=24 is unstable. Rayleigh's stability criterion holds for the first two cases but violates the MD simulation results for the third case. The remaining liquid threads with *L*\*>*λ<sup>c</sup>* \* (*L*\*=30~480) are unstable and break up into drop(s). Rayleigh's stability criterion holds for these cases.

If the temperature of the liquid threads of *R*\*=4 is increased to T\*=1.0, from Fig.5 or Table 1 (Cases 24~30) , only the liquid thread of *L*\*=10 remains intact while, in contrast to the situation of T\*=0.75, the liquid thread of *L*\*=20 is unstable. Rayleigh's stability criterion holds for the former case but violates the MD simulation results for the latter case. The remaining liquid threads with *L*\*>*λ<sup>c</sup>* \* (*L*\*=30~480) are unstable and break up into drop(s). Rayleigh's stability criterion holds for these cases.

Recently, Kim, Lee, Han and Park (2006) proposed a new linear relation from their MD simulation results. From their study, Rayleigh's linear relation can be modified as

$$
\lambda\_c = 2.44\pi (\text{R-1.23}) \tag{4}
$$

From Eq.(4), *λ<sup>c</sup>* \* for *R*\* of 2, 3 and 4 are 5.9, 13.6 and 21.2, respectively. Then, from Table 1 (Cases 1~23) , it is seen that Kim's stability criterion holds for the liquid threads at T\*=0.75. On the other hand, at T\*=1.0 (Cases 24~30) , the liquid thread of *R*\*=4 and *L*\*=10 remains intact while the liquid thread of *R*\*=4 and *L*\*=20 is unstable. Kim's stability criterion holds for the former case but violates the MD simulation results for the latter case. The remaining liquid threads with *L*\*>*λ<sup>c</sup>* \* (*L*\*=30~480) are unstable and break up into drop(s). Kim's stability criterion holds for these cases.

From the above discussion, the following observations can be made.

First, a liquid thread with a longer fundamental cell length is more unstable.

Second, a thinner liquid thread is more unstable.

Third, a liquid thread at a higher temperature is more unstable.

Fourth, the trends of linear stability theories agree with MD simulation results. However, Rayleigh's stability criterion overpredicts stable domain as compared to the MD simulation results. Kim's stability criterion gives more accurate predictions. However, it overpredicts the stable domain at a higher temperature as compared to the MD simulation results.

#### **2.1.4 Density distribution**

In practical applications, e.g. combustor or printer, faster vaporization is usually desirable. Criteria have to be made to quantify the discussion regarding the vaporization process of a liquid thread. In this research, a liquid thread is considered to vaporize faster if the distribution of molecules reaches uniform state quicker during the vaporization process. This criterion essentially concerns with the evolution of the density distribution. The density at a specified point in the fundamental cell can be defined as

$$\rho = \lim\_{\delta V \to 0} \frac{\delta N}{\delta V} \tag{5}$$

where *δV* is a small volume surrounding the point considered and *δN* is the number of molecules inside the volume *δV*. The density defined by Eq.(5) is actually an averaged density of a small volume surrounding the point considered. The value will approach the density of a specified point if the volume *δV* shrinks to that point. However, for a meaningful density field, the volume *δV* can not be too small because when *δV* becomes too small, it is difficult to obtain a definite value for *δN/δV*. In this study, the volume *δV* is taken to be a sphere with non-dimensionalized radius *R*\*=2 and with its center located at the point considered. This is an optimal choice after numerical test.

holds for the former case but violates the MD simulation results for the latter case. The

Recently, Kim, Lee, Han and Park (2006) proposed a new linear relation from their MD

(Cases 1~23) , it is seen that Kim's stability criterion holds for the liquid threads at T\*=0.75. On the other hand, at T\*=1.0 (Cases 24~30) , the liquid thread of *R*\*=4 and *L*\*=10 remains intact while the liquid thread of *R*\*=4 and *L*\*=20 is unstable. Kim's stability criterion holds for the former case but violates the MD simulation results for the latter case. The remaining

Fourth, the trends of linear stability theories agree with MD simulation results. However, Rayleigh's stability criterion overpredicts stable domain as compared to the MD simulation results. Kim's stability criterion gives more accurate predictions. However, it overpredicts

In practical applications, e.g. combustor or printer, faster vaporization is usually desirable. Criteria have to be made to quantify the discussion regarding the vaporization process of a liquid thread. In this research, a liquid thread is considered to vaporize faster if the distribution of molecules reaches uniform state quicker during the vaporization process. This criterion essentially concerns with the evolution of the density distribution. The density

> 0 lim *V*

 

where *δV* is a small volume surrounding the point considered and *δN* is the number of molecules inside the volume *δV*. The density defined by Eq.(5) is actually an averaged density of a small volume surrounding the point considered. The value will approach the density of a specified point if the volume *δV* shrinks to that point. However, for a meaningful density field, the volume *δV* can not be too small because when *δV* becomes too small, it is difficult to obtain a definite value for *δN/δV*. In this study, the volume *δV* is taken to be a sphere with non-dimensionalized radius *R*\*=2 and with its center located at the point

 *V* 

 *N*

(5)

the stable domain at a higher temperature as compared to the MD simulation results.

\* for *R*\* of 2, 3 and 4 are 5.9, 13.6 and 21.2, respectively. Then, from Table 1

simulation results. From their study, Rayleigh's linear relation can be modified as

From the above discussion, the following observations can be made.

Third, a liquid thread at a higher temperature is more unstable.

at a specified point in the fundamental cell can be defined as

considered. This is an optimal choice after numerical test.

First, a liquid thread with a longer fundamental cell length is more unstable.

\* (*L*\*=30~480) are unstable and break up into drop(s).

*λc=*2.44*π*(*R*-1.23) (4)

\* (*L*\*=30~480) are unstable and break up into drop(s). Kim's

remaining liquid threads with *L*\*>*λ<sup>c</sup>*

From Eq.(4), *λ<sup>c</sup>*

liquid threads with *L*\*>*λ<sup>c</sup>*

**2.1.4 Density distribution** 

stability criterion holds for these cases.

Second, a thinner liquid thread is more unstable.

Rayleigh's stability criterion holds for these cases.

Fig. 6. Comparison of time averaged density uniformity factor, *f* , for liquid threads of different radius, length and temperature

The time averaged density uniformity factor, *f* , in a time interval of *t\**=0 to 1000, as shown in Fig.6 or listed in Table 1, can be used to indicate the vaporization speed of a liquid thread. The density uniformity factor, *fρ* , is defined as

$$f\_{\rho} = \frac{\sum\_{N} \left(\rho^\* - \rho^\* \_{\text{eq}}\right)\_{t^\*} \Delta V}{\sum\_{N} \left(\rho^\* - \rho^\* \_{\text{eq}}\right)\_{t^\*=0} \Delta V} \tag{6}$$

where *N* is the total number of molecules in the fundamental cell, i.e. Nmol in Table 1, ρ\* and ᇞ*V* are the density and volume of molecule *i*, respectively, as defined by Eq.(5), and ρ\*eq is the density value when the molecules are uniformly distributed, i.e. ρ\*eq≡Nmol / Vol, where Vol is the volume of the fundamental cell. The density uniformity factor, *fρ* , as defined by Eq.(6) represents the deviation from uniform state. From Fig.6 or Table 1, the following observations can be drawn.

First, a liquid thread with a shorter fundamental cell length evaporates quicker. Take the liquid threads of *R*\*=2 (Cases 1~7) as an example. The liquid thread of L\*=10 (Case 1) produces one liquid particle while the liquid thread of L\*=480 (Case 7) produces ten liquid particles. The length of the latter liquid thread is forty-eight times of the former one. This implies that if forty-eight liquid threads of L\*=10 are connected in series, there will be fortyeight liquid particles. It is known that molecular interaction plays an important role in the vaporization process. More liquid particles can provide more molecular interactions and this is conducive to vaporization.

Second, a liquid thread evaporates quicker at a higher temperature.

Third, a liquid thread with a higher *f* is more unstable and produces more liquid particles in the fundamental cell.

#### **2.2 Molecular dynamics analysis of the vaporization process for two nano-scale liquid threads coexisting in a periodic fundamental cell**

As mentioned in section 2.1, previous studies of nano-scale liquid threads have mostly been devoted to the investigation of a single liquid thread in a periodic fundamental cell. Because of the interaction between the two nano-scale liquid threads coexisting in a periodic fundamental cell, the vaporization process is different from that of a single liquid thread in a periodic fundamental cell. This section discusses the influences of the liquid thread radius, fundamental cell length, and relative position of the two threads. Snapshots of molecules, the number of liquid particles formed, and density field are analyzed. Two linear stability criteria, namely, Rayleigh's stability criterion and Kim's stability criterion, are accessed for their validity in molecular scale. This approach will be helpful for the understanding and prediction of the atomization process.

In this study, two cylindrical liquid threads of length L\* and radius *R*\* are placed in the computational domain and the remaining space is a vacuum. The relative position of the two threads is controlled by L1\*, L2\* and L4\*, as shown in Fig.7. L3\* is set to be zero, i.e., the two liquid threads are kept at *x*=0 initially and are shifted only in the *y* direction. The initial density of the liquid argon is ρL\*=0.819 and the system temperature is kept at *T*\*=0.75. These dimensionless values correspond to ρL=1223 kg/m3 for argon and *T*=70K, which is below the critical temperature (150K) of argon. Simulation domain dimensions, the relative position of the two threads and number of molecules are listed in Table 2, together with some simulation results.

**(The origin of the coordinates is located at the center of the fundamental cell.)**

Fig. 7. Illustration of the computational domain for the simulation of two nano-scale liquid threads coexisting in a periodic fundamental cell

#### **2.2.1 Liquid thread vaporization process**

358 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

**2.2 Molecular dynamics analysis of the vaporization process for two nano-scale liquid** 

As mentioned in section 2.1, previous studies of nano-scale liquid threads have mostly been devoted to the investigation of a single liquid thread in a periodic fundamental cell. Because of the interaction between the two nano-scale liquid threads coexisting in a periodic fundamental cell, the vaporization process is different from that of a single liquid thread in a periodic fundamental cell. This section discusses the influences of the liquid thread radius, fundamental cell length, and relative position of the two threads. Snapshots of molecules, the number of liquid particles formed, and density field are analyzed. Two linear stability criteria, namely, Rayleigh's stability criterion and Kim's stability criterion, are accessed for their validity in molecular scale. This approach will be helpful for the understanding and

In this study, two cylindrical liquid threads of length L\* and radius *R*\* are placed in the computational domain and the remaining space is a vacuum. The relative position of the two threads is controlled by L1\*, L2\* and L4\*, as shown in Fig.7. L3\* is set to be zero, i.e., the two liquid threads are kept at *x*=0 initially and are shifted only in the *y* direction. The initial density of the liquid argon is ρL\*=0.819 and the system temperature is kept at *T*\*=0.75. These dimensionless values correspond to ρL=1223 kg/m3 for argon and *T*=70K, which is below the critical temperature (150K) of argon. Simulation domain dimensions, the relative position of the two threads and number of molecules are listed in Table 2, together with

(vacuum)

**(The origin of the coordinates is located at the center of the fundamental cell.)**

Fig. 7. Illustration of the computational domain for the simulation of two nano-scale liquid

*y*

*x*

*z*

liquid thread #1 liquid thread #2

L4\*

*R\**

L2 L \* 1\*

L\*

*R\** L3\*

**threads coexisting in a periodic fundamental cell** 

prediction of the atomization process.

L\*

threads coexisting in a periodic fundamental cell

some simulation results.

L\*

Figure 8 shows the vaporization process of two liquid threads of *R*\*=4, L\*=30, L1\*=10, L2\*=10 and L4\*=15 (case 1 in Table 2) . Note that L\*=30 and *R*\*=4 correspond to L=10.6nm and *R*=1.41nm, respectively. It can be seen that the two liquid threads quickly coalesce into a single thread and remain intact. No liquid particles are formed during the vaporization process. Figure 9 shows the vaporization process of two liquid threads of smaller radius *R*\*=2 but with the same length and relative position : L\*=30, L1\*=10, L2\*=10 and L4\*=15 (case 2 in Table 2) . For this case in which the liquid threads are thinner, it is seen that the two threads rupture from their two ends, i.e., the top and bottom surfaces of the fundamental cell, and get shorter due to the contraction motion in their axial directions. Two liquid particles are formed and they quickly coalesce into a single particle by collision and coalescence. The coalesced liquid particle prevails during the subsequent vaporization process.


\* The number in the parenthesis in the Np column denotes the number of liquid particles formed at time t\*=1000, while the number outside the parenthesis in the Np column denotes the maximum number of liquid particles formed during the vaporization process. If no parenthesis is denoted, it means that the maximum number of liquid particles formed prevails over the vaporization process.

Table 2. Simulation domain dimensions, relative position of the two threads, number of molecules, and simulation results (L3\*=0)

From Table 2, it is seen that the value of *f* is smaller for the thinner liquid threads. As stated in the previous section, this implies that the thinner liquid threads evaporate more quickly. In addition, it can be seen from Table 2 that the thinner liquid threads produce more liquid particles. In our previous study (Yeh, 2009b), we found that a single liquid thread is more unstable and produces more liquid particles when it is thinner. The present results for two liquid threads coexisting in a fundamental cell corroborate the findings of the previous study of a single liquid thread.

Fig. 8. Vaporization process of two liquid threads of *R*\*=4, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15

Fig. 9. Vaporization process of two liquid threads of *R*\*=2, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15

360 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

is smaller for the thinner liquid threads. As

**t\*=800**

**R\*=4, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15**

**t\*=900**

**R\*=2, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15**

0

0

X Y

Z

0

X Y

Z

0

**t\*=200**

**R\*=4, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15**

Fig. 8. Vaporization process of two liquid threads of *R*\*=4, L\*=30, L1\*=10, L2\*=10, L3\*=0,

**t\*=20**

**R\*=2, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15**

Fig. 9. Vaporization process of two liquid threads of *R*\*=2, L\*=30, L1\*=10, L2\*=10, L3\*=0,

stated in the previous section, this implies that the thinner liquid threads evaporate more quickly. In addition, it can be seen from Table 2 that the thinner liquid threads produce more liquid particles. In our previous study (Yeh, 2009b), we found that a single liquid thread is more unstable and produces more liquid particles when it is thinner. The present results for two liquid threads coexisting in a fundamental cell corroborate the findings of the

0

0

X Y


Z

0

X Y


Z

0

From Table 2, it is seen that the value of *f*

previous study of a single liquid thread.

0

0

X Y


Z

0

X Y


Z

0


L4\*=15


L4\*=15

**t\*=1**

**R\*=4, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15**

**t\*=1**

**R\*=2, L\*=30, L1\*=10, L2\*=10, L3\*=0, L4\*=15**

Fig. 10. Vaporization process of two liquid threads of *R*\*=4, L\*=60, L1\*=30, L2\*=15, L3\*=0, L4\*=30

Fig. 11. Vaporization process of two liquid threads of *R*\*=4, L\*=60, L1\*=10, L2\*=25, L3\*=0, L4\*=30

Figure 10 shows the vaporization process of two liquid threads of *R*\*=4, L\*=60, L1\*=30, L2\*=15 and L4\*=30 (case 3 in Table 2) . It is observed that the two threads rupture from their two ends, i.e., the top and bottom surfaces of the fundamental cell, and get shorter due to the contraction motion in their axial directions. Two liquid particles are formed and prevail during the subsequent vaporization process. In contrast to cases 1 and 2 for shorter liquid threads, the two liquid particles do not coalesce. If the two liquid threads get closer to each other, e.g., *R*\*=4, L\*=60, L1\*=10, L2\*=25 and L4\*=30 (case 4 in Table 2) , the vaporization process, which is shown in Fig.11, would be very different from case 3. From Fig.11, it is seen that because the two threads are very close to each other, they quickly coalesce into a single thread. The coalesced liquid thread gets shorter due to the contraction motion in its axial direction and finally evolves into a single liquid particle.

From Table 2, it can be seen that the value of *f* is larger for case 4. This implies that vaporization is slower when the two liquid threads are close to each other. In addition, Table 2 shows that more liquid particles are formed when the separation of the two threads

is larger (case 3 in Table 2) . In our previous study (Yeh, 2009b), we found that molecular interaction plays an important role in the vaporization process. More liquid particles can provide more molecular interactions, which is conducive to vaporization. Comparison of cases 3, 4 and 5 also reveals this tendency. Case 5, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as in cases 3 and 4. From Table 2, it is observed that only one liquid particle is formed for case 5 and the vaporization speed for case 5 is slower than case 3, for which two liquid particles are formed. The present results for two liquid threads coexisting in a fundamental cell corroborate the findings of our previous study for a single liquid thread.

Fig. 12. Vaporization process of two liquid threads of *R*\*=2, L\*=60, L1\*=30, L2\*=15, L3\*=0, L4\*=30

Fig. 13. Vaporization process of two liquid threads of *R*\*=2, L\*=60, L1\*=10, L2\*=25, L3\*=0, L4\*=30

Figure 12 shows the vaporization process of two liquid threads of *R*\*=2, L\*=60, L1\*=30, L2\*=15 and L4\*=30 (case 6 in Table 2) . Similar to case 3, the two threads rupture from their two ends and get shorter due to the contraction motion in their axial directions. Two liquid

is larger (case 3 in Table 2) . In our previous study (Yeh, 2009b), we found that molecular interaction plays an important role in the vaporization process. More liquid particles can provide more molecular interactions, which is conducive to vaporization. Comparison of cases 3, 4 and 5 also reveals this tendency. Case 5, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as in cases 3 and 4. From Table 2, it is observed that only one liquid particle is formed for case 5 and the vaporization speed for case 5 is slower than case 3, for which two liquid particles are formed. The present results for two liquid threads coexisting in a fundamental

0

0

X Y


Z

0

X Y


X Y

Z

0

0

X Y

Z

0

**t\*=900**

**R\*=2, L\*=60, L1\*=30, L2\*=15, L3\*=0, L4\*=30**

**t\*=900**

**R\*=2, L\*=60, L1\*=10, L2\*=25, L3\*=0, L4\*=30**

Z

0

cell corroborate the findings of our previous study for a single liquid thread.

**t\*=200**

**R\*=2, L\*=60, L1\*=30, L2\*=15, L3\*=0, L4\*=30**

Fig. 12. Vaporization process of two liquid threads of *R*\*=2, L\*=60, L1\*=30, L2\*=15, L3\*=0,

**t\*=60**

**R\*=2, L\*=60, L1\*=10, L2\*=25, L3\*=0, L4\*=30**

Fig. 13. Vaporization process of two liquid threads of *R*\*=2, L\*=60, L1\*=10, L2\*=25, L3\*=0,

Figure 12 shows the vaporization process of two liquid threads of *R*\*=2, L\*=60, L1\*=30, L2\*=15 and L4\*=30 (case 6 in Table 2) . Similar to case 3, the two threads rupture from their two ends and get shorter due to the contraction motion in their axial directions. Two liquid

0

0

X Y


Z

0

X Y


Z

0


L4\*=30


L4\*=30

**t\*=1**

**R\*=2, L\*=60, L1\*=30, L2\*=15, L3\*=0, L4\*=30**

**t\*=1**

**R\*=2, L\*=60, L1\*=10, L2\*=25, L3\*=0, L4\*=30**

particles are formed and prevail during the subsequent vaporization process. For two closer liquid threads, e.g., *R*\*=2, L\*=60, L1\*=10, L2\*=25 and L4\*=30 (case 7 in Table 2) , the vaporization process, which is shown in Fig.13, is different from case 6. The two threads also rupture from their two ends and get shorter due to the contraction motion in their axial directions. Two liquid particles are formed but they quickly coalesce into a single particle and prevail during the subsequent vaporization process. This vaporization process is similar to but slightly different from that for thicker liquid threads (case 4) . The thinner liquid threads (case 7) evolve into two liquid particles and then quickly coalesce into a single particle while the thicker liquid threads (case 4) quickly coalesce into a single liquid thread and then evolve into a single liquid particle.

From Table 2, it is seen that the value of *f* is larger for case 7, i.e., vaporization is slower when the two liquid threads are close to each other. In addition, Table 2 shows that more liquid particles are produced when the separation of the two threads is larger (case 6) . Case 8, which was investigated in our previous study (2009b), denotes the situation for a single liquid thread of the same radius and length as those in cases 6 and 7. From Table 2, it is seen that only one liquid particle is formed for case 8 and the vaporization speed for case 8 is slower than case 6, for which two liquid particles are formed. This corroborates the previous results for thicker liquid threads (cases 3, 4 and 5) . Furthermore, comparing cases 6~8 for thinner liquid threads (*R*\*=2) with cases 3~5 for thicker liquid threads (*R*\*=4) , it can be seen that the value of *f* is smaller for thinner liquid threads, which implies that thinner liquid threads evaporate more quickly. This also corroborates the previous results for shorter liquid threads (cases 1 and 2).

Figure 14 shows the vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=60, L2\*=30 and L4\*=60 (case 9 in Table 2) . Similar to case 3 for shorter liquid threads, the two threads rupture from their two ends and get shorter due to the contraction motion in their axial directions. Two liquid particles are formed and prevail during the subsequent vaporization process. For two closer liquid threads, e.g., *R*\*=4, L\*=120, L1\*=40, L2\*=40 and L4\*=60 (case 10 in Table 2) , the vaporization process, which is shown in Fig.15, is similar to case 9. Two liquid particles are formed and prevail during the subsequent vaporization process. If the two liquid threads get even closer, e.g., *R*\*=4, L\*=120, L1\*=10, L2\*=55 and L4\*=60 (case 11 in Table 2) , the vaporization process, which is shown in Fig.16, is quite different from cases 9 and 10. Similar to case 4 for shorter liquid threads, it is observed that because the two threads in case 11 are very close to each other, they quickly coalesce into a single thread before they separately evolve into liquid particles. The coalesced liquid thread then gets shorter due to the contraction motion in its axial direction and finally evolves into a single liquid particle.

Similar to cases 4 and 7, Table 2 shows that the value of *f* is larger for case 11, i.e., vaporization is slower when the two liquid threads are very close to each other. Similarly, it is seen that more liquid particles are produced when the separation of the two threads is larger (cases 9 and 10) . Case 12, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as those in cases 9~11. From Table 2, it is seen that only one liquid particle is formed for case 12 and the vaporization speed for case 12 is slower than in cases 9 and 10, for which two liquid particles are formed. This corroborates the previous results for shorter liquid threads (cases 3~8) .

Fig. 14. Vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=60, L2\*=30, L3\*=0, L4\*=60

Fig. 15. Vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=40, L2\*=40, L3\*=0, L4\*=60

Fig. 16. Vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=10, L2\*=55, L3\*=0, L4\*=60

Figure 17 shows the vaporization process of two liquid threads of *R*\*=2, L\*=120, L1\*=60, L2\*=30 and L4\*=60 (case 13 in Table 2) . For this case of thinner liquid threads, it is seen that the two threads rupture from their two ends and quickly evolve into ten liquid particles. During the vaporization process, collision and coalescence of the liquid particles occur and finally two liquid particles are formed. For two closer liquid threads, e.g., *R*\*=2, L\*=120, L1\*=40, L2\*=40 and L4\*=60 (case 14 in Table 2) , the vaporization process, which is shown in Fig.18, is similar to but slightly different from that of case 13. Nine liquid particles are formed initially and they finally evolve into two liquid particles by collision and coalescence. If the two liquid threads get even closer, e.g., *R*\*=2, L\*=120, L1\*=10, L2\*=55 and L4\*=60 (case 15 in Table 2) , the vaporization process, which is shown in Fig.19, is different

Fig. 14. Vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=60, L2\*=30, L3\*=0,

Fig. 15. Vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=40, L2\*=40, L3\*=0,

Fig. 16. Vaporization process of two liquid threads of *R*\*=4, L\*=120, L1\*=10, L2\*=55, L3\*=0,

Figure 17 shows the vaporization process of two liquid threads of *R*\*=2, L\*=120, L1\*=60, L2\*=30 and L4\*=60 (case 13 in Table 2) . For this case of thinner liquid threads, it is seen that the two threads rupture from their two ends and quickly evolve into ten liquid particles. During the vaporization process, collision and coalescence of the liquid particles occur and finally two liquid particles are formed. For two closer liquid threads, e.g., *R*\*=2, L\*=120, L1\*=40, L2\*=40 and L4\*=60 (case 14 in Table 2) , the vaporization process, which is shown in Fig.18, is similar to but slightly different from that of case 13. Nine liquid particles are formed initially and they finally evolve into two liquid particles by collision and coalescence. If the two liquid threads get even closer, e.g., *R*\*=2, L\*=120, L1\*=10, L2\*=55 and L4\*=60 (case 15 in Table 2) , the vaporization process, which is shown in Fig.19, is different

L4\*=60

L4\*=60

L4\*=60

from those in cases 13 and 14. Nine liquid particles are formed, but they finally evolve into a single particle by collision and coalescence.

Fig. 17. Vaporization process of two liquid threads of *R*\*=2, L\*=120, L1\*=60, L2\*=30, L3\*=0, L4\*=60

Fig. 18. Vaporization process of two liquid threads of *R*\*=2, L\*=120, L1\*=40, L2\*=40, L3\*=0, L4\*=60

Similar to cases 4, 7 and 11, Table 2 shows that the value of *f* is larger for case 15, i.e., vaporization is slower when the two liquid threads are very close to each other. Similarly, it can also be seen that more liquid particles are formed when the separation of the two threads is larger (cases 13 and 14) . Case 16, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as those in cases 13~15. From Table 2, it is seen that less liquid particles are formed for case 16 and the vaporization speed for case 16 is slower than in cases 13 and 14, in which more liquid particles are formed. This corroborates the previous results (cases 3~12) . Furthermore, comparing cases 13~16 for thinner liquid threads (*R*\*=2) with cases 9~12 for thicker liquid threads (*R*\*=4) , it is seen that the value of *f* is smaller for thinner liquid threads, which implies that thinner liquid threads evaporate more quickly. This also corroborates the previous results for shorter liquid threads (cases 1~8).

Figure 20 shows the vaporization process of two liquid threads of *R*\*=4, L\*=240, L1\*=120, L2\*=60 and L4\*=120 (case 17 in Table 2) . Similar to cases 3 and 9 for shorter liquid threads, the two threads rupture from their two ends and get shorter due to the contraction motion in their axial directions. Two liquid particles are formed and prevail during the subsequent vaporization process. For two closer liquid threads, e.g., *R*\*=4, L\*=240, L1\*=40, L2\*=100 and L4\*=120 (case 18 in Table 2) , the vaporization process, which is shown in Fig.21, is similar to that in case 17. If the two liquid threads get even closer, e.g., *R*\*=4, L\*=240, L1\*=10, L2\*=115 and L4\*=120 (case 19 in Table 2) , the vaporization process, which is shown in Fig.22, is quite different from those in cases 17 and 18. Similar to cases 4 and 11 for shorter liquid threads, it is seen that because the two threads in case 19 are very close to each other, they quickly coalesce into a single thread before they separately evolve into liquid particles. The coalesced liquid thread then becomes shorter due to the contraction motion in its axial direction and finally evolves into a single liquid particle.

Fig. 19. Vaporization process of two liquid threads of *R*\*=2, L\*=120, L1\*=10, L2\*=55, L3\*=0, L4\*=60

From Table 2, it is observed that the value of *f* is larger for case 19, i.e., vaporization is slower when the two liquid threads are very close to each other. Similarly, it can be seen that more liquid particles are produced when the separation of the two threads is larger (cases 17 and 18) . Case 20, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as those in cases 17~19. From Table 2, it is seen that only one liquid particle is formed for case 20 and the vaporization speed for case 20 is slower than those in cases 17 and 18, for which two liquid particles are formed. This corroborates the previous results for shorter liquid threads (cases 3~16) .

Fig. 20. Vaporization process of two liquid threads of *R*\*=4, L\*=240, L1\*=120, L2\*=60, L3\*=0, L4\*=120

If the length of the liquid thread is further increased, e.g. *R*\*=4 and L\*=480, more liquid particles are formed. Figure 23 shows the vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=160, L2\*=160 and L4\*=240 (case 21 in Table 2) . Eight liquid particles are formed and they evolve into four liquid particles by collision and coalescence. For two

that in case 17. If the two liquid threads get even closer, e.g., *R*\*=4, L\*=240, L1\*=10, L2\*=115 and L4\*=120 (case 19 in Table 2) , the vaporization process, which is shown in Fig.22, is quite different from those in cases 17 and 18. Similar to cases 4 and 11 for shorter liquid threads, it is seen that because the two threads in case 19 are very close to each other, they quickly coalesce into a single thread before they separately evolve into liquid particles. The coalesced liquid thread then becomes shorter due to the contraction motion in its axial

Fig. 19. Vaporization process of two liquid threads of *R*\*=2, L\*=120, L1\*=10, L2\*=55, L3\*=0,

slower when the two liquid threads are very close to each other. Similarly, it can be seen that more liquid particles are produced when the separation of the two threads is larger (cases 17 and 18) . Case 20, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as those in cases 17~19. From Table 2, it is seen that only one liquid particle is formed for case 20 and the vaporization speed for case 20 is slower than those in cases 17 and 18, for which two liquid particles are formed. This corroborates the previous results for shorter liquid threads (cases

Fig. 20. Vaporization process of two liquid threads of *R*\*=4, L\*=240, L1\*=120, L2\*=60, L3\*=0,

If the length of the liquid thread is further increased, e.g. *R*\*=4 and L\*=480, more liquid particles are formed. Figure 23 shows the vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=160, L2\*=160 and L4\*=240 (case 21 in Table 2) . Eight liquid particles are formed and they evolve into four liquid particles by collision and coalescence. For two

is larger for case 19, i.e., vaporization is

direction and finally evolves into a single liquid particle.

From Table 2, it is observed that the value of *f*

L4\*=60

3~16) .

L4\*=120

closer liquid threads, e.g., *R*\*=4, L\*=480, L1\*=40, L2\*=220 and L4\*=240 (case 22 in Table 2) , the vaporization process, which is shown in Fig.24, is similar to but slightly different from that in case 21. Six liquid particles are formed initially and they finally evolve into two liquid particles by collision and coalescence. If the two liquid threads get even closer, e.g., *R*\*=4, L\*=480, L1\*=10, L2\*=235 and L4\*=240 (case 23 in Table 2) , the vaporization process, which is shown in Fig.25, is very different from those in cases 21 and 22. Because the two threads in case 23 are very close to each other, they quickly coalesce into a single thread before they separately evolve into liquid particles. The coalesced liquid thread ruptures from its ends and interior and gets shorter due to the contraction motion in its axial direction and finally evolves into a single liquid particle by collision and coalescence. The vaporization process for case 23 (Fig.25) is similar to but slightly different from that in case 19 (Fig.22) . For both cases, the two threads quickly coalesce into a single thread because they are very close to each other. However, the coalesced liquid thread in case 19 evolves into a single liquid particle without rupturing in its interior while the coalesced liquid thread for case 23 ruptures in its interior before it finally evolves into a single liquid particle.

Fig. 21. Vaporization process of two liquid threads of *R*\*=4, L\*=240, L1\*=40, L2\*=100, L3\*=0, L4\*=120

Fig. 22. Vaporization process of two liquid threads of *R*\*=4, L\*=240, L1\*=10, L2\*=115, L3\*=0, L4\*=120

From Table 2, it is observed that the value of *f* is larger for case 23, i.e., vaporization is slower when the two liquid threads are very close to each other. Similarly, it can be seen that more liquid particles are produced when the separation of the two threads is larger (cases 21 and 22) . Case 24, which was investigated in our previous study (Yeh, 2009b), denotes the situation for a single liquid thread of the same radius and length as those in cases 21~23. From Table 2, it can be seen that less liquid particles are formed for case 24 and the vaporization speed for case 24 is slower than those in cases 21 and 22, in which more liquid particles are formed. This corroborates the previous results for shorter liquid threads (cases 3~20).

Fig. 23. Vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=160, L2\*=160, L3\*=0, L4\*=240

Fig. 24. Vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=40, L2\*=220, L3\*=0, L4\*=240

Fig. 25. Vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=10, L2\*=235, L3\*=0, L4\*=240

From the above discussion, the following observations can be made.

First, if the fundamental cell length is small, the liquid threads may remain intact or evolve into only one liquid particle in the cell. If the threads break up in this case, they rupture from their ends only, i.e. the top and bottom surfaces of the fundamental cell, but not from their interiors. On the other hand, if the fundamental cell length is larger, more than one liquid particle may be produced in the cell and the liquid threads may rupture not only from their ends but also from their interiors.

Second, thinner liquid threads may produce more liquid particles in the cell and evaporate more quickly.

Third, more liquid particles are formed when the separation of the two threads is larger. Moreover, vaporization is slower when the two liquid threads are close to each other.

Fourth, on the basis of identical liquid thread radius and length, liquid threads that produce more liquid particles evaporate more quickly.

### **2.2.2 Stability analysis**

368 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

vaporization speed for case 24 is slower than those in cases 21 and 22, in which more liquid particles are formed. This corroborates the previous results for shorter liquid threads (cases

Fig. 23. Vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=160, L2\*=160, L3\*=0,

Fig. 24. Vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=40, L2\*=220, L3\*=0,

Fig. 25. Vaporization process of two liquid threads of *R*\*=4, L\*=480, L1\*=10, L2\*=235, L3\*=0,

First, if the fundamental cell length is small, the liquid threads may remain intact or evolve into only one liquid particle in the cell. If the threads break up in this case, they rupture from their ends only, i.e. the top and bottom surfaces of the fundamental cell, but not from

From the above discussion, the following observations can be made.

3~20).

L4\*=240

L4\*=240

L4\*=240

In the following discussion, two linear stability criteria, namely Rayleigh's stability criterion (1879) and Kim's stability criterion (2006), are accessed for their validity in molecular scale. This is the first study to investigate the instability of two nano-scale liquid threads coexisting in a periodic fundamental cell by MD simulation.

According to classical theory by Rayleigh (1879), the radius of a liquid thread, *R*, is related to the critical wavelength of perturbation, *λc* , by Eq.(3) in section 2.1.3. A liquid thread will break up into drops if the axial wavelength of the surface perturbation *L*>*λc*. If *L*<*λc*, the thread is stable and will remain intact. Owing to the periodic boundary conditions in this study, the fundamental cell size *L*\* can be regarded as the longest wavelength of the perturbation.

For *R*\*=2, *λ<sup>c</sup>* \* is 12.6. According to Rayleigh's stability criterion, liquid threads of radius *R*\*=2 and length *L*\*=30, 60 and 120 are expected to be unstable and will break up into drop(s). From Figs.9, 12, 13, 17, 19 or Table 2 (cases 2, 6~8, 13~16) , Rayleigh's stability criterion holds for these cases.

For *R*\*=4, *λ<sup>c</sup>* \* is 25.2. According to Rayleigh's stability criterion, liquid threads of radius *R*\*=4 and length *L*\*=30, 60, 120, 240 and 480 are expected to be unstable and will break up into drop(s). However, from Figs.8, 10, 11, 14~16, 20~25 or Table 2 (cases 1, 3~5, 9~12, 17~24) , the liquid threads of *R*\*=4 and *L*\*=30 (case 1) remain intact while the remaining liquid threads are unstable and break up into drop(s). Rayleigh's stability criterion violates the MD simulation result for case 1 but holds for the remaining cases. The invalidity of Rayleigh's stability criterion for case 1 can be interpreted according to the vaporization process shown in Fig.8. It is seen that the two liquid threads quickly coalesce into a single thread and remain intact. This is somewhat like combining the two threads to form a thread of a larger radius, which can increase the critical wavelength of perturbation, *λc*, according to Eq.(3) and, therefore, broaden the stable domain.

Recently, Kim, Lee, Han and Park (2006) proposed a new linear relation from their MD simulation results. From their study, Rayleigh's linear relation can be modified by Eq.(4) in section 2.1.3. From Eq.(4), *λ<sup>c</sup>* \* for *R*\* of 2 and 4 are 5.9 and 21.2, respectively. Then, from Table 2, it is seen that Kim's stability criterion violates the MD simulation result for case 1 while it holds for the remaining cases. The validity of Kim's stability criterion is similar to that of Rayleigh's stability criterion.

From the above discussion, it is found that the trends of Rayleigh's stability criterion and Kim's stability criterion agree with MD simulation results. However, when the two liquid threads coalesce into a single thread and remain intact, the critical wavelength of perturbation may be increased and the stable domain is broadened. In such a situation, Rayleigh's stability criterion and Kim's stability criterion underpredict the stable domain.

### **2.3 Molecular dynamics simulation for the atomization process of a nanojet**

In this section, liquid argon nanojets made of 44000 Lennard-Jones molecules are investigated under various simulation parameters to examine their influence on the nanojet atomization process. Snapshots of the molecules, evolution of the density field, and evolution of the intermolecular force are analyzed. This can provide insight into the fundamental mechanism of the atomization process and will be helpful for the design of nanojet devices such as nano-printer or nano-sprayer.

The nano-atomizer is schematically shown in Fig.26. The simulation domain comprises a cubical box of side length 3600, with periodic boundary conditions applied in all three directions. The nano-atomizer is placed at the center of the box. Simulation parameters are listed in Table 3, which include nano-atomizer dimensions, temperatures, number of molecules and simulation results.

The argon molecules inside the nano-atomizer are liquid and the nano-atomizer is made of rigid argon molecules. A push panel composed of 600 argon molecules is constructed with a downward velocity of 120m/s. In this study, the interactions among liquid argon molecules, nano-atomizer and push panel are taken into account.

Fig. 26. Illustration of the nano-atomizer configuration and dimensions

From the above discussion, it is found that the trends of Rayleigh's stability criterion and Kim's stability criterion agree with MD simulation results. However, when the two liquid threads coalesce into a single thread and remain intact, the critical wavelength of perturbation may be increased and the stable domain is broadened. In such a situation, Rayleigh's stability criterion and Kim's stability criterion underpredict the stable domain.

In this section, liquid argon nanojets made of 44000 Lennard-Jones molecules are investigated under various simulation parameters to examine their influence on the nanojet atomization process. Snapshots of the molecules, evolution of the density field, and evolution of the intermolecular force are analyzed. This can provide insight into the fundamental mechanism of the atomization process and will be helpful for the design of

The nano-atomizer is schematically shown in Fig.26. The simulation domain comprises a cubical box of side length 3600, with periodic boundary conditions applied in all three directions. The nano-atomizer is placed at the center of the box. Simulation parameters are listed in Table 3, which include nano-atomizer dimensions, temperatures, number of

The argon molecules inside the nano-atomizer are liquid and the nano-atomizer is made of rigid argon molecules. A push panel composed of 600 argon molecules is constructed with a downward velocity of 120m/s. In this study, the interactions among liquid argon molecules,

**push panel**

**D/2**

/2

**surroundings (vacuum)**

**L3**

Fig. 26. Illustration of the nano-atomizer configuration and dimensions

**2.3 Molecular dynamics simulation for the atomization process of a nanojet** 

nanojet devices such as nano-printer or nano-sprayer.

nano-atomizer and push panel are taken into account.

**L2**

**L4**

**L5**

molecules and simulation results.


Table 3. Nano-atomizer dimensions, temperatures, number of molecules and simulation results

In the following discussion, a liquid argon nanojet of length L5\* and diameter 2L3\*+D\* is pushed by a panel into vacuum through a nano-nozzle of orifice diameter D\*, as illustrated in Fig.26. Simulation conditions are listed in Table 3.

#### **2.3.1 Nanojet atomization process**

Figure 27 shows the atomization process for a nanojet of L2\*=L3\*=L4 \*=5.73, L5\*=76.7, D\*/2=8.81 and T\*=0.75, which corresponds to a nanojet of length 26.2 nm and diameter 10 nm, and a nano-nozzle of orifice length 2nm, diameter 6 nm, as well as an actual temperature of 70 K. The dot in Fig.27 indicates the center of the molecule. From the figure it is found that the nanojet does not break up. Owing to the low temperature, the molecular kinetic energies are so low that the molecules congregate near the orifice exit. Very few liquid molecules evaporate at this low temperature. At a higher temperature T\*=1.5 (140 K) , as depicted in Fig.28, the molecules leave the orifice exit earlier than at T\*=0.75, due to their higher molecular kinetic energies. More liquid molecules evaporate at this higher temperature. However, like at T\*=0.75, the nanojet has not broken up before t\*=80. If the temperature is further increased to T\*=2.0 (187 K) , as shown in Fig.29, evident evaporation is observed. Many evaporated molecules are produced and the non-evaporated liquid molecules concentrate within the central region. Figure 30 depicts the snapshots for temperature T\*=3.0 (278 K) . It is observed that breakup of the nanojet occurs and its spray angle is larger than at T\*=2.0. The spurted molecules from the nano-atomizer are more evenly distributed at this temperature. If the temperature is further increased to T\*=4.5 (420 K) , as shown in Fig.31, the spray angle is even larger than at T\*=3.0 and the spurted molecules from the nano-atomizer are much more uniformly distributed as compared to the lower temperature cases. Comparison of Figs.27, 28, 29, 30 and 31 reveals that the liquid nanojet evaporates quicker at higher temperatures. This will be further illustrated in later sections discussing the density distribution and the intermolecular force.

Fig. 27. Atomization process for case 1 in Table 3 (L2 \* =5.73, L3 \* =5.73, L4 \* =5.73, L5 \* =76.7, D\*/2=8.81, TD\*=0.75)

Fig. 27. Atomization process for case 1 in Table 3 (L2\* =5.73, L3\* =5.73, L4

D\*/2=8.81, TD\*=0.75)

\* =5.73, L5

\* =76.7,

Fig. 28. Atomization process for case 2 in Table 3 (L2\* =5.73, L3 \* =5.73, L4 \* =5.73, L5\* =76.7, D\*/2=8.81, TD\*=1.5)

To investigate the influence of nozzle geometry on the nanojet atomization, comparison of the snapshots at t\*=80 for four different nozzle geometries and T\*=2.0 or 3.0 is shown in Figs.32 and 33. Note that Figs.32(a), (b), (c) and (d) correspond to cases 6, 3, 7 and 8, respectively, in Table 3; while Figs.33(a), (b), (c) and (d) correspond to cases 9, 4, 10 and 11, respectively, in Table 3. In Figs.32(a), (b), (c) and 9(a), (b), (c), the nozzle orifice diameters are equal (6 nm) but the nozzle orifice lengths are varied; while in Figs.32(d) and 33(d), the nozzle orifice length is the same as for Figs.32(b) and 33(b) (2 nm) but the nozzle orifice diameter is smaller (4 nm) . By a careful comparison of Figs.32(a), (b) and (c), it can be observed that, on the basis of identical nozzle orifice diameter, a nanojet from a nozzle with a shorter orifice length (L2) moves farther. On the other hand, from Figs.32(b) and (d), on the basis of identical nozzle orifice length, a nanojet from a nozzle with a larger orifice diameter moves farther. Figure 52 reveals similar tendency. Note that the nano-atomizer in this research is basically a plain-orifice atomizer. As pointed out by Lefebvre (1989) , in a practical plain-orifice atomizer, resistance increases with nozzle orifice length/diameter ratio. Therefore, a nanojet from a nozzle with a smaller orifice length/diameter ratio moves farther due to its smaller resistance. This will be further illustrated in later sections discussing the density distribution and the intermolecular force.

Fig. 29. Atomization process for case 3 in Table 3 (L2\* =5.73, L3\* =5.73, L4 \* =5.73, L5\* =76.7, D\*/2=8.81, TD\*=2.0)

#### **2.3.2 Density distribution**

It is important that the system be in equilibrium state before statistical values of the local properties can be taken. However, owing to the computational capacity limitations, the MD simulation can not proceed to a macroscopically long period. Nevertheless, the purpose of this paper is not to discuss statistical values of the local properties but to investigate the atomization process of a nanojet, which is important and conducive to the understanding of the fundamental mechanism of the atomization process. Criteria have to be made to quantify the discussion regarding the nanojet atomization process. Unfortunately, such criteria are still arbitrary in the literature. Because the system temperature in this study is kept at the desired temperature, a constant temperature criterion is not suitable for the discussion of the atomization process. In this research, a nanojet is considered to vaporize faster if the distribution of molecules reaches a uniform state quicker during the atomization process. This criterion essentially concerns with the evolution of the density distribution. The density at a specified point in the fundamental cell can be defined as Eq.(5). In this study, the volume δV is taken to be a sphere with nondimensionalized radius R\*=2 and with its center located at the point considered. This is an optimal choice after numerical test.

Fig. 29. Atomization process for case 3 in Table 3 (L2\* =5.73, L3\* =5.73, L4

It is important that the system be in equilibrium state before statistical values of the local properties can be taken. However, owing to the computational capacity limitations, the MD simulation can not proceed to a macroscopically long period. Nevertheless, the purpose of this paper is not to discuss statistical values of the local properties but to investigate the atomization process of a nanojet, which is important and conducive to the understanding of the fundamental mechanism of the atomization process. Criteria have to be made to quantify the discussion regarding the nanojet atomization process. Unfortunately, such criteria are still arbitrary in the literature. Because the system temperature in this study is kept at the desired temperature, a constant temperature criterion is not suitable for the discussion of the atomization process. In this research, a nanojet is considered to vaporize faster if the distribution of molecules reaches a uniform state quicker during the atomization process. This criterion essentially concerns with the evolution of the density distribution. The density at a specified point in the fundamental cell can be defined as Eq.(5). In this study, the volume δV is taken to be a sphere with nondimensionalized radius R\*=2 and with its center located at the point considered. This is an

D\*/2=8.81, TD\*=2.0)

**2.3.2 Density distribution** 

optimal choice after numerical test.

\* =5.73, L5\* =76.7,

Fig. 30. Atomization process for case 4 in Table 3 (L2\* =5.73, L3 \* =5.73, L4 \* =5.73, L5\* =76.7, D\*/2=8.81, TD\*=3.0)

Figure 34 shows the evolution of density uniformity factor for nanojets at different temperatures and the conditions of L2\*=L3\*=L4\*=5.73, L5\*=76.7, D\*/2=8.81 ( cases 1~5 in Table 3 ). The density uniformity factor is defined as Eq.(6). From Fig.34 it is observed that a higher temperature nanojet evaporates faster than a lower temperature one and this corroborates the results of Figs.27~31 as discussed in section 2.3.1. The time averaged value of the density uniformity factor, *f* , in a time interval of t\*=0 to 80, as listed in Table 3, also reveals this observation. In Fig.34, it is noted that at lower temperatures (T\*=0.75 and 1.5) , the density uniformity factor increases first and then decreases. For a lower temperature nanojet, the momenta of the liquid molecules away from the push panel in the nano-atomizer are low while the molecules near the push panel have relatively higher momenta due to the action of the push panel. This results in a compression effect that leads to the increase of the density uniformity factor at the earlier stage of the atomization process; while at a later stage, the density uniformity factor drops because of the ejection of the molecules.

Fig. 31. Atomization process for case 5 in Table 3 (L2\* =5.73, L3\* =5.73, L4 \* =5.73, L5\* =76.7, D\*/2=8.81, TD\*=4.5)

(T\*=0.75 and 1.5) , the density uniformity factor increases first and then decreases. For a lower temperature nanojet, the momenta of the liquid molecules away from the push panel in the nano-atomizer are low while the molecules near the push panel have relatively higher momenta due to the action of the push panel. This results in a compression effect that leads to the increase of the density uniformity factor at the earlier stage of the atomization process; while at a later stage, the density uniformity factor drops

Fig. 31. Atomization process for case 5 in Table 3 (L2\* =5.73, L3\* =5.73, L4

D\*/2=8.81, TD\*=4.5)

\* =5.73, L5\* =76.7,

because of the ejection of the molecules.

Fig. 32. Comparison of the snapshots at t\*=80 and T\*=2.0 for four different nozzle geometries

Fig. 33. Comparison of the snapshots at t\*=80 and T\*=3.0 for four different nozzle geometries

Fig. 33. Comparison of the snapshots at t\*=80 and T\*=3.0 for four different nozzle geometries

Fig. 34. Evolution of the density uniformity factor for different temperatures (cases 1~5 in Table 3)

Fig. 35. Evolution of the density uniformity for different orifice lengths at ( cases 3, 6 and 7 in Table 3 )

Fig. 36. Evolution of the density uniformity factor for different orifice lengths at T\*=3.0 (cases 4, 9 and 10 in Table 3)

Fig. 37. Evolution of the density uniformity factor for different orifice radii at T\*=2.0 (cases 3 and 8 in Table 3)

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> 0.4

Fig. 36. Evolution of the density uniformity factor for different orifice lengths at T\*=3.0

*orifice radius 3nm orifice radius 2nm*

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> 0.75

Fig. 37. Evolution of the density uniformity factor for different orifice radii at T\*=2.0 (cases 3

**t\***

*orifice length 3nm orifice length 2nm orifice length 1nm*

**t\***

0.5

1

**N**

and 8 in Table 3)

**\***

**v**

/

**N**

**\*t\*=0**

**v**

1.25

1.5

**N**

(cases 4, 9 and 10 in Table 3)

**\***

**v**

/

**N**

**\*t\*=0**

**v**

0.6

0.7

0.8

0.9

1

Fig. 38. Evolution of the density uniformity factor for different orifice radii at T\*=3.0 (cases 4 and 11 in Table 3)

Figure 35 shows the evolution of density uniformity factor for different orifice lengths on the basis of identical nozzle orifice diameter (6 nm) and the conditions of L3 \*=5.73, L5 \*=76.7, T\*=2.0 ( cases 3, 6 and 7 in Table 3 ). It is observed that a nanojet with a shorter orifice length evaporates quicker. The time averaged value of the density uniformity factor, *f* , in a time interval of *t\**=0 to 80, as listed in Table 3, also reveals this observation. This corroborates the results of Figs.32(a), (b) and (c) as discussed in section 2.3.1. In addition, it is also observed that, as time elapsed, the influence of the orifice length mitigates. This is because as time elapsed, more and more molecules spurt from the atomizer and hence the interaction between the liquid molecules and the rigid atomizer molecules mitigates due to the decrease of number of molecules inside the atomizer. Figure 36 also reveals this tendency. However, as can be observed from Fig.36, the influence of the orifice length becomes less pronounced at a higher temperature because of the higher molecular kinetic energy to overcome the resistance caused by the orifice. Figure 37 shows the evolution of density uniformity factor for different orifice diameters on the basis of identical nozzle orifice length (2 nm) and the conditions of L4\*=5.73, L5\*=76.7, T\*=2.0 (cases 3 and 8 in Table 3). It is observed that a nanojet with a larger orifice diameter evaporates quicker. The time averaged value of the density uniformity factor, *f* , in a time interval of *t\**=0 to 80, as listed in Table 3, also reveals this observation. This corroborates the results of Figs.32(a) and (d) and also reveals previous observation that a nanojet from a nozzle with a smaller orifice length/diameter ratio evaporates quicker. Similar tendency is obtained from Fig.38 for a higher temperature nanojet.

#### **2.3.3 Intermolecular force**

The intermolecular force is an indication of the surface tension experienced by the liquid particles and has a great effect upon the atomization process. Lefebvre (1989), Chigier (1999) and Hiroyasu (2000) pointed out that surface tension and interfacial force are the major controlling mechanisms for atomization. Owing to the vacuum environment, the aerodynamic effect on the atomization process is negligible in this study. Thus, the surface tension becomes the major controlling mechanism for the atomization process. Figure 39 shows the evolution of averaged non-dimensionalized intermolecular force for nanojets with different temperatures and the conditions of L2 \*=L3 \*=L4 \*=5.73, L5 \*=76.7, D\*/2=8.81 ( cases 1~5 in Table 3 ). The averaged non-dimensionalized intermolecular force at time t\* is defined as

$$\frac{\mathbf{F}\_{i^\*}^\*}{\mathbf{F}\_{i^\*}^\*} = \frac{\sum\_{i=1}^N \mathbf{F}\_{i,t^\*}^\*}{\mathbf{N}} \tag{7}$$

where N is the total number of molecules in the fundamental cell and F\* i,t\* is the resultant force of the non-dimensionalized intermolecular force vector acting on molecule i at time t\* , i.e. F\* i,t\* = ( F\* x,i,t\*2 + F\* y,i,t\*2 + F\* z,i,t\*2 )1/2 , where F\* x,i,t\*, F\* y,i,t\* and F\* z,i,t\* are the components of the intermolecular force vector at the x, y and z directions, respectively, acting on molecule i at time t\*. Note that in the above definition of \* F , N is the total number of molecules in the t\* fundamental cell, which includes liquid, vapor and solid molecules (atomizer and push panel) ; while in the definition of density uniformity factor, Eq.(6), N is only the initial number of liquid molecules in the fundamental cell, i.e. the solid molecules are excluded. The intermolecular force diminishes with time because of the increase of distances between molecules as the nanojet vaporizes. From Fig.39, it is observed that a higher temperature nanojet evaporates faster than a lower temperature one. This corroborates the results of Figs.27~31 discussed in section 2.3.1 and Fig.34 in section 2.3.2. In Fig.39, it is also noted that although a higher temperature nanojet has a larger intermolecular force at the earlier stage of the atomization process due to its higher momentum, it evaporates faster and therefore the intermolecular force decays quicker. Figure 40 shows the time averaged value of the averaged non-dimensionalized intermolecular force for different orifice lengths on the basis of identical nozzle orifice diameter (6 nm) and the conditions of L3 \*=5.73, L5 \*=76.7, T\*=2.0 ( cases 3, 6 and 7 in Table 3). It is observed that a nanojet with a shorter orifice length evaporates quicker. This corroborates the results of Figs.32(a), (b) and (c) as discussed in section 2.3.1 and Fig.35 as discussed in section 2.3.2. In addition, it is also observed that, as time elapsed, the influence of the orifice length mitigates. As explained in section 2.3.2, more and more molecules spurt from the atomizer as time elapsed. This causes the interaction between the liquid molecules and the rigid atomizer molecules to mitigate due to the decrease of number of molecules inside the atomizer. Figure 41 also reveals this tendency. However, as can be observed from Fig.41, the influence of the orifice length becomes less pronounced at a higher temperature because of the higher molecular kinetic energy to overcome the resistance caused by the orifice. Figure 42 shows the evolution of averaged non-dimensionalized intermolecular force for different orifice diameters on the basis of

The intermolecular force is an indication of the surface tension experienced by the liquid particles and has a great effect upon the atomization process. Lefebvre (1989), Chigier (1999) and Hiroyasu (2000) pointed out that surface tension and interfacial force are the major controlling mechanisms for atomization. Owing to the vacuum environment, the aerodynamic effect on the atomization process is negligible in this study. Thus, the surface tension becomes the major controlling mechanism for the atomization process. Figure 39 shows the evolution of averaged non-dimensionalized intermolecular force for nanojets

cases 1~5 in Table 3 ). The averaged non-dimensionalized intermolecular force at time t\* is

<sup>N</sup> \* , \*

F

N *i t*

x,i,t\*, F\*

\* i 1 t\*

 

force of the non-dimensionalized intermolecular force vector acting on molecule i at time t\* ,

the intermolecular force vector at the x, y and z directions, respectively, acting on molecule i at time t\*. Note that in the above definition of \* F , N is the total number of molecules in the t\* fundamental cell, which includes liquid, vapor and solid molecules (atomizer and push panel) ; while in the definition of density uniformity factor, Eq.(6), N is only the initial number of liquid molecules in the fundamental cell, i.e. the solid molecules are excluded. The intermolecular force diminishes with time because of the increase of distances between molecules as the nanojet vaporizes. From Fig.39, it is observed that a higher temperature nanojet evaporates faster than a lower temperature one. This corroborates the results of Figs.27~31 discussed in section 2.3.1 and Fig.34 in section 2.3.2. In Fig.39, it is also noted that although a higher temperature nanojet has a larger intermolecular force at the earlier stage of the atomization process due to its higher momentum, it evaporates faster and therefore the intermolecular force decays quicker. Figure 40 shows the time averaged value of the averaged non-dimensionalized intermolecular force for different orifice lengths on the basis of identical nozzle orifice diameter (6 nm) and the conditions of L3\*=5.73, L5\*=76.7, T\*=2.0 ( cases 3, 6 and 7 in Table 3). It is observed that a nanojet with a shorter orifice length evaporates quicker. This corroborates the results of Figs.32(a), (b) and (c) as discussed in section 2.3.1 and Fig.35 as discussed in section 2.3.2. In addition, it is also observed that, as time elapsed, the influence of the orifice length mitigates. As explained in section 2.3.2, more and more molecules spurt from the atomizer as time elapsed. This causes the interaction between the liquid molecules and the rigid atomizer molecules to mitigate due to the decrease of number of molecules inside the atomizer. Figure 41 also reveals this tendency. However, as can be observed from Fig.41, the influence of the orifice length becomes less pronounced at a higher temperature because of the higher molecular kinetic energy to overcome the resistance caused by the orifice. Figure 42 shows the evolution of averaged non-dimensionalized intermolecular force for different orifice diameters on the basis of

F

where N is the total number of molecules in the fundamental cell and F\*

z,i,t\*2 )1/2 , where F\*

\*=L3 \*=L4

y,i,t\* and F\*

\*=5.73, L5

\*=76.7, D\*/2=8.81 (

i,t\* is the resultant

z,i,t\* are the components of

(7)

**2.3.3 Intermolecular force** 

defined as

i.e. F\*

i,t\* = ( F\*

x,i,t\*2 + F\*

y,i,t\*2 + F\*

with different temperatures and the conditions of L2

identical nozzle orifice length (2 nm) and the conditions of L4\*=5.73, L5\*=76.7, T\*=2.0 ( cases 3 and 8 in Table 3 ). It is observed that a nanojet with a larger orifice diameter evaporates quicker. This corroborates the results of Figs.51(a) and (d) as discussed in section 2.3.1 and Fig.37 as discussed in section 2.3.2 and also reveals that a nanojet from a nozzle with a smaller orifice length/diameter ratio evaporates quicker. Similar tendency can be observed from Fig.43 for a higher temperature nanojet.

Fig. 39. Evolution of the averaged non-dimensionalized intermolecular force for different temperatures (cases 1~5 in Table 3)

Fig. 40. Evolution of the averaged non-dimensionalized intermolecular force for different orifice lengths at T\*=2.0 ( cases 3, 6 and 7 in Table 3 )

Fig. 41. Evolution of the averaged non-intermolecular force for different orifice lengthsdimensionalized at T\*=3.0 (cases 4, 9 and 10 in Table 3)

Fig. 42. Evolution of the averaged non-dimens ionalized intermolecular force for different orifice radii at T\*=2.0 (cases 3 and 8 in Table 3)

Fig. 43. Evolution of the averaged non-dimensionalized intermolecular force for different orifice radii at T\*=3.0 ( cases 4 and 11 in Table 3 )

### **3. Conclusions**

384 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

*orifice length 3nm orifice length 2nm orifice length 1nm*

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>15</sup>

<sup>0</sup> <sup>20</sup> <sup>40</sup> <sup>60</sup> <sup>80</sup> <sup>10</sup>

Fig. 42. Evolution of the averaged non-dimens ionalized intermolecular force for different

*orifice radius 3nm orifice radius 2nm*

**t\***

Fig. 41. Evolution of the averaged non-intermolecular force for different orifice

lengthsdimensionalized at T\*=3.0 (cases 4, 9 and 10 in Table 3)

**t\***

20

20

orifice radii at T\*=2.0 (cases 3 and 8 in Table 3)

30

40

F\*i,t\*

/

N

**N** 50

60

70

25

30

35

F\*i,t\*

/

N

**N** 40

45

50

The purpose of this study is to provide a thorough investigation and understanding of the microscopic evolution and detailed mechanism of the flow inside and outside the atomizers. The main findings of this study are listed below.


liquid particles are formed when the separation of the two threads is larger. Moreover, vaporization is slower when the two liquid threads are close to each other. On the basis of identical liquid thread radius and length, liquid threads that produce more liquid particles evaporate more quickly. Finally, the trends of Rayleigh's stability criterion and Kim's stability criterion agree with MD simulation results. However, when the two threads coalesce into a single one and remain intact, the critical wavelength of perturbation may be increased and the stable domain is broadened. In such a situation, Rayleigh's stability criterion and Kim's stability criterion underpredict the stable domain.

3. For the atomization process of a nanojet, it is found that a liquid nanojet evaporates faster at a higher temperature. On the basis of identical nozzle orifice diameter, a nanojet from a nozzle with a shorter orifice length evaporates quicker. However, the influence of the orifice length mitigates as time elapsed. In addition, the influence of the orifice length becomes less pronounced at a higher temperature. On the other hand, on the basis of identical nozzle orifice length, a nanojet from a nozzle with a larger orifice diameter evaporates quicker. The present simulation results reveal that a nozzle with a smaller orifice length/diameter ratio produces better atomization. This corroborates the results from conventional macroscopic analysis.

### **4. Acknowledgment**

The author gratefully acknowledges the grant support from the National Science Council, Taiwan, R.O.C., under the contract NSC100-2221-E-150-047.

### **5. References**


3. For the atomization process of a nanojet, it is found that a liquid nanojet evaporates faster at a higher temperature. On the basis of identical nozzle orifice diameter, a nanojet from a nozzle with a shorter orifice length evaporates quicker. However, the influence of the orifice length mitigates as time elapsed. In addition, the influence of the orifice length becomes less pronounced at a higher temperature. On the other hand, on the basis of identical nozzle orifice length, a nanojet from a nozzle with a larger orifice diameter evaporates quicker. The present simulation results reveal that a nozzle with a smaller orifice length/diameter ratio produces better atomization. This corroborates the

The author gratefully acknowledges the grant support from the National Science Council,

Ashgriz, N. & Poo, J. Y. (1991). FLAIR:Flux Line-Segment Model for Advection and Interface Reconstruction. *Journal of Computational Physics*, Vol.93, pp.449-468. Brackbill, J. U., Kothe, D. B. & Zemach, C. (1992). A Continuum Method for Modeling Surface Tension. *Journal of Computational Physics*, Vol.100, pp.335-354. Chen, S., Johnson, D. B. & Raad, P. E. (1995). Velocity Boundary Conditions for the

Chen, S. K. & Lefebvre, A. H. (1994). Discharge Coefficients for Plain-Orifice Effervescent

Goren, S. (1962). The Instability of an Annular Thread of Fluid. *Journal of Fluid Mechanics*,

Incompressible Flow of Fluid with Free Surface. *The Physics of Fluids*, Vol.8,

Haile, J. M. (1992). *Molecular Dynamics Simulation*, John Wiley & Sons, New York, Chap.5. Harlow, F. H. & Welch, J. E. (1965). Numerical Calculation of Time-Dependent Viscous

Crank, J. (1984). *Free and Moving Boundary Problems*, Oxford University Press, New York. Floryan, J. M. & Rasmussen, H. (1989). Numerical Methods for Viscous Flows with Moving

Atomizers. *Atomization and Sprays*, Vol.4, pp.275-290. Chigier, N. (1999). Breakup of Liquid Sheets and Jet. *AIAA paper* 99-3640.

Boundaries. *Appl. Mech. Rev.*, Vol.42, pp.323-340.

Simulation of Free Surface Fluid Flow. *Journal of Computational Physics*, Vol.116,

results from conventional macroscopic analysis.

Taiwan, R.O.C., under the contract NSC100-2221-E-150-047.

domain.

**4. Acknowledgment**

pp.262-276.

Vol.12, pp.309-319.

pp.2182-2189.

**5. References**

liquid particles are formed when the separation of the two threads is larger. Moreover, vaporization is slower when the two liquid threads are close to each other. On the basis of identical liquid thread radius and length, liquid threads that produce more liquid particles evaporate more quickly. Finally, the trends of Rayleigh's stability criterion and Kim's stability criterion agree with MD simulation results. However, when the two threads coalesce into a single one and remain intact, the critical wavelength of perturbation may be increased and the stable domain is broadened. In such a situation, Rayleigh's stability criterion and Kim's stability criterion underpredict the stable


## **Molecular Dynamics Simulation of Nanoscale Machining**

Akinjide Oluwajobi

*Obafemi Awolowo University, Ile-Ife Nigeria* 

### **1. Introduction**

388 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Tomotika, S. (1935). On the Instability of a Cylindrical Thread of a Viscous Liquid

Weber, C. (1931). Zum zerfall eines ussigkeitsstrahles. *Zeitschrift fur Angewandte Mathematik* 

Yeh, Chun-Lang (2010). Molecular Dynamics Analysis of the Vaporization Process for Two

Yeh, Chun-Lang (2009b). Molecular Dynamics Analysis of the Instability for a Nano-Scale

Yeh, Chun-Lang (2007). Numerical Simulation of Turbulent Liquid Jet Emanating from

Yeh, Chun-Lang (2005). Turbulent Flow Investigation inside and outside Plain-Orifice

Yeh, Chun-Lang (2004). Numerical Investigation of Liquid Jet Emanating from Plain-Orifice

Yeh, Chun-Lang (2003). Effect of Inlet Turbulence Intensity on Discharge Coefficients for

Yeh, Chun-Lang (2002). Numerical Study of Inlet and Geometry Effects on Discharge

*Aeronautics Astronautics and Aviation*, Vol.35, No.3, pp.299-306.

*Computer Modeling in Engineering & Sciences*, Vol.67, No.3, pp.175-209. Yeh, Chun-Lang (2009a). Molecular Dynamics Simulation for the Atomization Process of a

*Series A, Mathematical and Physical Sciences*, Vol.150, No.870, pp.322-337. Van Doormaal, J. P. & Raithby, G. D. (1984). Enhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows. *Numerical Heat Transfer*, Vol.7, pp.147-163. Vicelli, J. A. (1969). A Method for Including Arbitrary External Boundaries in the MAC

pp.543-551.

200.

pp.253-283.

810-823.

*A*, Vol.51, pp.1187-1212.

*Series B*, Vol.47, No.1, pp.37-47.

*Mechanics, Series A*, Vol.18, No.3, pp.153-161.

*und Mechanik*, Vol.11, No.2, pp.136-154.

Surrounded by Another Viscous Liquid. *Proceedings of the Royal Society of London,* 

Incompressible Fluid Computing Technique. *Journal of Computational Physics*, Vol.4,

Nano-Scale Liquid Threads Coexisting in a Periodic Fundamental Cell. *CMES:* 

Nanojet. *CMES: Computer Modeling in Engineering & Sciences*, Vol.39, No.2, pp.179-

Liquid Thread. *CMES: Computer Modeling in Engineering & Sciences*, Vol.50, No.3,

Plain-Orifice Atomizer and Pressure-Swirl Atomizer. *Numerical Heat Transfer, Part* 

Atomizers with Rounded Orifice Inlets. *Heat and Mass Transfer*, Vol.41, No.9, pp.

Atomizers with Chamfered or Rounded Orifice Inlets. *JSME International Journal,* 

Liquid Jet Emanating from a Plain-Orifice Atomizer:a Numerical Study. *Journal of* 

Coefficients for Liquid Jet Emanating from a Plain-Orifice Atomizer. *Journal of* 

Product miniaturization is a major motivation for the development of ultra-precision technologies and processes which can achieve high form and excellent surface finish. Of all the available manufacturing approaches, mechanical nanometric machining is still a good option for machining complex 3D devices in a controllable way (Jackson, 2008). As the dimension goes down to the nanoscale, the machining phenomena take place in a limited region of tool-workpiece interface. At this length scale and interface, the material removal mechanisms are not fully understood, so more insight is needed, which on the long run will help to achieve high precision manufacturing with predictability, repeatability and productivity (Luo, 2004). At present, it is very difficult to observe the diverse microscopic physical phenomena occurring through experiments at the nanoscale (Rentsch, 2008). Subsequently, the other alternative is to explore available simulation techniques. Continuum mechanics approach is not adequate, as the point of interest/interface cannot be assumed to be homogeneous, but rather discrete, so, atomistic simulation methods are the suitable techniques for modelling at the nanoscale.

The aim of this chapter is to review the use of classical molecular dynamics (MD) method for nanoscale machining. The steps on how the MD method can be applied to model nanometric machining and the importance of choosing appropriate interatomic potentials for the simulation are considered. Also, some examples of MD simulation of nanoscale machining are presented.

## **2. The MD method**

Molecular dynamics (MD) is a computer simulation technique used in the study of the motions of a set of particles – atoms and molecules (Allen & Tildesley, 1988; Haile, 1997; Field, 1999; Frenkel & Smit, 2001; Leach, 2001; Schlick, 2002; Rapaport, 2004). The technique works by following the time evolution of a set of interacting atoms while integrating the equations of their motion. It is employed to study the classical motion of a system and then to extract experimental observables from the dynamics (Turkerman & Martyna, 2000). The MD is deterministic; once the positions, velocities and accelerations of the particles are known, the state of the system can be predicted. The method is also based on statistical mechanics – a way to obtain a set of configurations distributed according to some statistical distribution functions (Ercolessi, 1997, Hernandez, 2008).

Fig. 1. Schematic of the MD Simulation of Nanometric Cutting (2D) (Komanduri & Raff, 2001)

The MD method was initiated in the late 1950s at Lawrence Radiation Laboratory in the US by Alder and Wainwright in the study of interactions of hard spheres (Alder & Wainwright, 1957, 1959). Gibson et al., 1961, carried out probably the first MD calculations with a continuous potential based on a finite difference time integration method. In 1964, Rahman carried out the simulation of liquid argon by using a realistic potential (Rahman, 1964). Later on, the first MD simulation of a realistic system was conducted by (Stillinger & Rahman, 1974) in 1974. Since then, the use of the simulation method has spread from Physics to Materials Science and now to Mechanical Engineering and other disciplines.

In the field of nanometric cutting, Belak pioneered work on the study of cutting copper with a diamond tool (Belak & Stowers, 1990). Initially, the method was used extensively to model indentation and cutting. Figure 1 shows the schematic of the MD simulation of nanoscale cutting. In 1991, Belak and Stowers first applied the MD to abrasive processes (Belak & Stowers, 1991) and Rentsch and Inasaki's study later presented the first results of simulations targeted on the pile-up phenomenon in abrasive machining (Rentsch & Inasaki, 1994). The MD method has been used in a variety of other application areas, namely; indentation (Kenny & Mulliah, 2005, Smith et al., 2003), wear and friction (Cheng et al., 2003, Shimizu et al., 1998), void growth (Potirniche et al., 2006, Zhao et al., 2009), etcetera.

The MD simulation is based on Newton's second law of motion. It consists of the numerical step-by-step integration of the classical equations of motion. For a set of N atoms,

$$F\_i = m\_i a\_i \tag{1}$$

Where *mi* is the mass of atom i, 2 2 *i i d r a dt* is the acceleration of the atom i and *Fi* is the resultant force acting on atom i. These forces should be balanced by the potential energy between atoms, which are usually presented as the gradient of a potential energy function.

### **2.1 Thermodynamic ensembles**

390 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Fig. 1. Schematic of the MD Simulation of Nanometric Cutting (2D) (Komanduri &

The MD method was initiated in the late 1950s at Lawrence Radiation Laboratory in the US by Alder and Wainwright in the study of interactions of hard spheres (Alder & Wainwright, 1957, 1959). Gibson et al., 1961, carried out probably the first MD calculations with a continuous potential based on a finite difference time integration method. In 1964, Rahman carried out the simulation of liquid argon by using a realistic potential (Rahman, 1964). Later on, the first MD simulation of a realistic system was conducted by (Stillinger & Rahman, 1974) in 1974. Since then, the use of the simulation method has spread from Physics to Materials Science and now to Mechanical Engineering and other disciplines.

In the field of nanometric cutting, Belak pioneered work on the study of cutting copper with a diamond tool (Belak & Stowers, 1990). Initially, the method was used extensively to model indentation and cutting. Figure 1 shows the schematic of the MD simulation of nanoscale cutting. In 1991, Belak and Stowers first applied the MD to abrasive processes (Belak & Stowers, 1991) and Rentsch and Inasaki's study later presented the first results of simulations targeted on the pile-up phenomenon in abrasive machining (Rentsch & Inasaki, 1994). The MD method has been used in a variety of other application areas, namely; indentation (Kenny & Mulliah, 2005, Smith et al., 2003), wear and friction (Cheng et al., 2003,

Shimizu et al., 1998), void growth (Potirniche et al., 2006, Zhao et al., 2009), etcetera.

step-by-step integration of the classical equations of motion. For a set of N atoms,

2 2 *i*

*d r*

*i*

*a*

Where *mi* is the mass of atom i,

function.

The MD simulation is based on Newton's second law of motion. It consists of the numerical

resultant force acting on atom i. These forces should be balanced by the potential energy between atoms, which are usually presented as the gradient of a potential energy

*F ma i ii* (1)

*dt* is the acceleration of the atom i and *Fi* is the

Raff, 2001)

In MD simulations, various thermodynamic ensembles are employed. An ensemble is a large group of atoms or systems which are in different microscopic states, but have the same macroscopic or thermodynamic states. If a system of N atoms in a given macrostate is defined in terms of thermodynamic quantities, number of atoms, N; pressure, P; temperature, T; entropy, S; volume, V etc, there are many configurations at the atomic scale, which will lead to the same macrostate. The microstate of a system, defined by the atomistic positions and momenta cannot be known in a deterministic manner, because of the uncertainty principle of quantum mechanics. To avoid this problem, statistical mechanics approach is used for the atomic description.

The commonly used ensembles are the following, namely;

*Microcanonical Ensemble (NVE)* 

This is an isolated system, with N atoms, which occupies a constant volume and the overall energy, E of the system is constant.

*Canonical Ensemble (NVT)* 

This is a system in a temperature bath, with N atoms; and the volume, V and the temperature, T of the system are kept constant.

*Isobaric Isothermal Ensemble (NPT)* 

This is a system in a temperature and pressure bath, with N atoms; and the pressure, P and the temperature, T of the system are kept constant.

*Grand Canonical Ensemble ( VT)* 

This is a system with a constant chemical potential and temperature T.

### **2.2 Steps in MD simulation**

The outline of the MD simulation is as follows:


### **2.2.1 Select the model**

The model should be selected correctly to reproduce what is expected (whether 2D/3D). The level of details and accuracy needed should guide this decision.

### **2.2.2 Choose an appropriate interatomic potential**

Then, the next thing is a major task in MD simulation, which is the selection of the potential function for the atomic interactions. This is a function, ( ,......... ) *Vr r i N* of the position of the nuclei which represents the potential energy of the system. Where ( ,......... ) *r r i N* represents the complete co-ordinate position of the atoms. Then forces are derived from it as,

$$F\_i = -\frac{\partial}{\partial r\_i} V(r\_i, \dots, \dots, r\_N) \tag{2}$$

Strictly speaking, the problem of modelling a material is that of finding a potential function for that material. If the potential function doesn't model the atomic behaviour correctly, the results produced from the simulation would be wrong. Consider the energy of N interacting particles, which can be written as (Tersoff, 1988);

$$E = \sum\_{i} V\_1(r\_i) + \sum\_{i$$

Where , , *<sup>i</sup> <sup>j</sup> <sup>k</sup> rrr* are the positions of the particles and the functions 123 *VVV V* , , ,.... *<sup>N</sup>* are the mbody potentials.

From equation (3), the second term, 2 (,) *<sup>i</sup> <sup>j</sup> i j V rr* is the two-body or pair potential and the third

term is the three-body potential and so on. The pair potentials are the simplest interatomic potentials used for the interaction of a set of particles. The most popularly used is the Lennard-Jones potential; others are Morse potential, Born-Mayer potential and et cetera. Apart from the pair potentials, there are multi-body potentials, like Tersoff and Embedded-Atom Method (EAM) potentials. The following is a consideration of some widely used interatomic potentials.

*Lennard-Jones Potential (Lennard-Jones, 1924)* 

$$V\_{\
u} = 4\varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right] \tag{4}$$

Where and are constants which are dependent on the interacting particles. The LJ potential is ideal for rare gases, where the interactions between the non- bonded and uncharged atoms are due to weak van der Waal forces.

*Morse Potential (Morse, 1929)* 

$$V\_{\vec{\eta}} = D\{\exp[-2\alpha(r\_{\vec{\eta}} - r\_e)] - 2\exp[-\alpha(r\_{\vec{\eta}} - r\_e)]\}\tag{5}$$

Where, *ij r* and *er* are instantaneous and equilibrium distances between atoms *i* and *j*  respectively; and D are constants determined on the basis of the physical properties of the material. The Morse potential is suitable for cubic metals and they can be used to model the interactions between an atom and a surface.

#### *Born-Mayer Potential*

Born and Mayer suggested that the repulsion between the atoms would have a roughly exponential dependence on distance (Born & Mayer, 1932).

$$V\_{i\bar{j}} = A\{\exp[-\frac{r}{a\_{\text{BM}}}]\}\tag{6}$$

Where A and *aBM* are constants dependent on the material. This potential is used for metals, Group III-V semiconductors and ceramics.

#### *Tersoff Potential*

Tersoff modelled the total energy of the system as a sum of pair like interactions and as a function of the atomic coordinates, given as equation (7). The potential is based on the concept that the strength of a bond between two atoms is not a constant, but depends on the local environment (Tersoff, 1988a ,1988b).

$$E = \sum\_{i} E\_{i} = \frac{1}{2} \sum\_{i} \sum\_{i \neq j} V\_{ij} \tag{7}$$

and

392 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

nuclei which represents the potential energy of the system. Where ( ,......... ) *r r i N* represents the

( ,.......... ) *i iN i F Vr r r*

Strictly speaking, the problem of modelling a material is that of finding a potential function for that material. If the potential function doesn't model the atomic behaviour correctly, the results produced from the simulation would be wrong. Consider the energy of N interacting

*E Vr V rr Vrrr V rrr r* 

Where , , *<sup>i</sup> <sup>j</sup> <sup>k</sup> rrr* are the positions of the particles and the functions 123 *VVV V* , , ,.... *<sup>N</sup>* are the m-

term is the three-body potential and so on. The pair potentials are the simplest interatomic potentials used for the interaction of a set of particles. The most popularly used is the Lennard-Jones potential; others are Morse potential, Born-Mayer potential and et cetera. Apart from the pair potentials, there are multi-body potentials, like Tersoff and Embedded-Atom Method (EAM) potentials. The following is a consideration of some widely used interatomic potentials.

> <sup>4</sup> *<sup>V</sup> ij r r*

 

potential is ideal for rare gases, where the interactions between the non- bonded and

{exp[ 2 ( )] 2exp[ ( ]} *VD rr rr ij*

Where, *ij r* and *er* are instantaneous and equilibrium distances between atoms *i* and *j* 

the material. The Morse potential is suitable for cubic metals and they can be used to model

Born and Mayer suggested that the repulsion between the atoms would have a roughly

{exp[ )]} *ij*

*<sup>r</sup> V A*

*BM*

*ij e*

( ) ( , ) ( , , ) ........ ( , , ...... ) *i i <sup>j</sup> <sup>i</sup> <sup>j</sup> <sup>k</sup> N i <sup>j</sup> k N*

(3)

12 6

 

are constants which are dependent on the interacting particles. The LJ

*<sup>a</sup>* (6)

and D are constants determined on the basis of the physical properties of

(2)

.......

is the two-body or pair potential and the third

(4)

*ij e* (5)

complete co-ordinate position of the atoms. Then forces are derived from it as,

*i ij ijk ijk N*

*V rr*

*i j*

particles, which can be written as (Tersoff, 1988);

From equation (3), the second term, 2 (,) *<sup>i</sup> <sup>j</sup>*

*Lennard-Jones Potential (Lennard-Jones, 1924)* 

uncharged atoms are due to weak van der Waal forces.

the interactions between an atom and a surface.

exponential dependence on distance (Born & Mayer, 1932).

body potentials.

Where

 and 

respectively;

*Born-Mayer Potential* 

*Morse Potential (Morse, 1929)* 

12 3

$$V\_{i\dagger} = f\_{\mathcal{C}}(r\_{i\dagger})[a\_{i\dagger}f\_R(r\_{i\dagger}) + b\_{i\dagger}f\_A(r\_{i\dagger})] \tag{8}$$

$$\begin{split} & \text{where} \\ & f\_R(r) = A \exp(-\mathcal{A}\_1 r), \\ & f\_A(r) = -B \exp(-\mathcal{A}\_2 r), \\ & f\_C(r) = \left[ \frac{1}{2} - \frac{1}{2} \sin\left[\frac{\pi}{2} (r - R) / D \right], R - D < r < R + D \right], \\ & \text{or} \quad r > R + D \\ & b\_{\vec{\eta}} = \left( 1 + \beta^n \zeta\_{\vec{\eta}}^{\vec{\eta}} \right)^{-1/2 \cdot n}, \\ & \zeta\_{\vec{\eta}} = \sum\_{k \left( \vec{\pi}\_{\vec{\eta}} \right)} f\_C(r\_{\vec{\eta}}) g(\theta\_{\vec{\eta}k}) \exp\left[ \lambda\_3^3 (r\_{\vec{\eta}} - r\_{\vec{\eta}})^3 \right], \\ & g(\theta) = 1 + \frac{p^2}{q^2} - \frac{p^2}{\left[ q^2 + \left( l - \cos \theta \right)^2 \right]}, \\ & a\_{\vec{\eta}} = \left( 1 + \alpha^n \eta\_{\vec{\eta}}^n \right)^{-1/2 \cdot n}, \\ & \eta\_{\vec{\eta}} = \sum\_{k \left( \vec{\eta} \right)} f\_C(r\_{\vec{\eta}}) \exp\left[ \lambda\_3^3 \left( r\_{\vec{\eta}} - r\_{\vec{\eta}} \right)^3 \right] \end{split}$$

Where R and D are cutoff parameters; A,B, , , , , ,n,p,q,h <sup>123</sup> are fitting parameters of the Tersoff potential.

The Tersoff potential is used for covalently bonded materials like silicon and carbon atoms. *Embedded-Atom Method (EAM) Potential* 

For the EAM potential, the total energy of the system can be written as (Foiles ,1985, Foiles et al., 1986).

$$E\_{\text{tot}} = \sum\_{i} G\_{i}(\rho\_{h,i}) + \frac{1}{2} \sum\_{i,j} V\_{ij}(r\_{ij}) \tag{9}$$

*h i*, is the total electron density at atom i due to the rest of the atoms in the system.

*Gi* is the embedding energy for placing an atom into the electron density

*Vi*, *<sup>j</sup>* is the short range pair interaction representing the core-core repulsion

*ij r* is the separation of atoms i and j

The EAM potential was developed for a wide range of metals. It incorporates an approximation for the many-atom interactions neglected by the pair-potential scheme.

*Modified Embedded-Atom Method (MEAM) Potential* 

The MEAM is theoretically an extension of the EAM potential (Baskes, 1992) with modifications to include the directionality of the bonding. The bond-angle was explicitly handled so as to accommodate covalent systems. The total energy is given by (Equation 9). The MEAM is suitable for modelling metals and alloys with fcc, bcc, hcp and cubic structures, and also for covalent materials such as silicon and carbon (Oluwajobi & Chen, 2010a).

The list of the potentials is not exhaustive; and there are other potentials which are modified forms of the ones already discussed. To obtain physically meaningful results from atomistic simulations, it is essential that reliable interatomic potentials are used. A reliable potential would reproduce various physical properties of the relevant elements or alloys, including the elastic, structural, defect, surface and thermal properties etc (Lee et al., 2005). For example, Tersoff potential was designed for the description of covalent materials like silicon, germanium, carbon, silicon carbide etc. and it cannot adequately model metals. Also, the EAM potential was designed for metals, as it describes the bonding in metals more satisfactorily, but the MEAM potential can be used for the modelling of both metallic systems and covalently bonded materials. Table 1 shows some interatomic potentials and their suitability (Oluwajobi & Chen, 2010).

It should be stated that EAM and MEAM potentials should be used for metals, where appropriate; Tersoff and MEAM potentials should be used for covalent materials, and for the interface of materials where suitable potentials are not available, appropriate existing potentials like LJ and Morse can be used with caution.

### **2.2.3 Algorithms for the Integration of the Equations of Motion**

The next step, after choosing the interatomic potential/s is to select an appropriate algorithm for the integration of the Newton's equations of motions. This is to numerically solve for positions and velocities at a time, t and at a later time, t+t. (Generally, the algorithm is already chosen in the MD software that one uses for the simulation). MD simulations use time steps from a few femto seconds ( <sup>15</sup> 10 s) (Shimada et al., 1993). If a large time step is used, the computational time is reduced, but when the time step is too large, this can cause instability and inaccuracy in the solutions.

For the EAM potential, the total energy of the system can be written as (Foiles ,1985, Foiles

,

*h i*, is the total electron density at atom i due to the rest of the atoms in the system.

*Gi* is the embedding energy for placing an atom into the electron density *Vi*, *<sup>j</sup>* is the short range pair interaction representing the core-core repulsion

<sup>1</sup> ( ) () <sup>2</sup> *tot i h i ij ij i i j E G Vr* 

The EAM potential was developed for a wide range of metals. It incorporates an approximation for the many-atom interactions neglected by the pair-potential scheme.

The MEAM is theoretically an extension of the EAM potential (Baskes, 1992) with modifications to include the directionality of the bonding. The bond-angle was explicitly handled so as to accommodate covalent systems. The total energy is given by (Equation 9). The MEAM is suitable for modelling metals and alloys with fcc, bcc, hcp and cubic structures, and also for covalent materials such as silicon and carbon (Oluwajobi & Chen,

The list of the potentials is not exhaustive; and there are other potentials which are modified forms of the ones already discussed. To obtain physically meaningful results from atomistic simulations, it is essential that reliable interatomic potentials are used. A reliable potential would reproduce various physical properties of the relevant elements or alloys, including the elastic, structural, defect, surface and thermal properties etc (Lee et al., 2005). For example, Tersoff potential was designed for the description of covalent materials like silicon, germanium, carbon, silicon carbide etc. and it cannot adequately model metals. Also, the EAM potential was designed for metals, as it describes the bonding in metals more satisfactorily, but the MEAM potential can be used for the modelling of both metallic systems and covalently bonded materials. Table 1 shows some interatomic potentials and

It should be stated that EAM and MEAM potentials should be used for metals, where appropriate; Tersoff and MEAM potentials should be used for covalent materials, and for the interface of materials where suitable potentials are not available, appropriate existing

The next step, after choosing the interatomic potential/s is to select an appropriate algorithm for the integration of the Newton's equations of motions. This is to numerically solve for positions and velocities at a time, t and at a later time, t+t. (Generally, the algorithm is already chosen in the MD software that one uses for the simulation). MD simulations use time steps from a few femto seconds ( <sup>15</sup> 10 s) (Shimada et al., 1993). If a large time step is used, the computational time is reduced, but when the time step is too

,

(9)

et al., 1986).

*ij r* is the separation of atoms i and j

*Modified Embedded-Atom Method (MEAM) Potential* 

their suitability (Oluwajobi & Chen, 2010).

potentials like LJ and Morse can be used with caution.

**2.2.3 Algorithms for the Integration of the Equations of Motion** 

large, this can cause instability and inaccuracy in the solutions.

2010a).


Table 1. Comparison of the Interatomic Potentials (Oluwajobi & Chen, 2010a)

Many numerical schemes have been developed for the integration of the equations of motions. Some of these are the Verlet algorithm (Verlet, 1967), the predictor-corrector algorithm (Rahman, 1964, Gear & Tu, 1974, Gear & Watanabe, 1974) and the Beeman's algorithm (Beeman, 1976). All the above integration schemes make the assumption, that the positions, velocities and accelerations can be approximated using a Taylor series expansion:

$$\mathbf{r}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{r}(\mathbf{t}) + \mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{1}{2}\mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^2 + \dots$$

$$\mathbf{v}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{\bar{r}}(\mathbf{t}) + \mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{1}{2}\mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^2 + \dots \tag{10}$$

$$\mathbf{a}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{\bar{r}}(\mathbf{t}) + \mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{1}{2}\mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^2 + \dots$$

Where t is a finite time step, r is the position, v and . r are the velocity, a and .. r are the acceleration; ... <sup>r</sup> and .... r are the third and the fourth derivative of position with time.

*The Basic Verlet Algorithm (Verlet, 1967)*

Using this method, the next step of position can be predicted as follows;

$$\mathbf{r}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{r}(\mathbf{t}) + \mathbf{\dot{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{1}{2}\mathbf{\dot{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^2 + \dots \tag{11}$$

In the same way, the previous step;

$$\mathbf{r}(\mathbf{t} - \omega \mathbf{t}) = \mathbf{r}(\mathbf{t}) - \mathbf{\dot{r}}(\mathbf{t})\omega \mathbf{t} + \frac{1}{2}\mathbf{\dot{r}}(\mathbf{t})\omega \mathbf{t}^2 - \dots \tag{12}$$

Summing equations (11) and (12), we have

$$\mathbf{r(t+\Delta t) = 2r(t) - r(t-\Delta t) + \bar{r(t)}\Delta t^2} \tag{13}$$

The method uses no explicit velocities, it is quite straightforward, easy to implement and its storage requirements are modest, but it is not too accurate.

#### *The Verlet Leapfrog Algorithm*

In this algorithm, the velocities are calculated by taking the average value halfway between position steps. The equations are as follows;

$$\mathbf{r}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{r}(\mathbf{t}) + \mathbf{\dot{r}}(\mathbf{t} + \frac{1}{2}\mathcal{A}\mathbf{t})\mathcal{A}\mathbf{t} \tag{14}$$

$$\mathbf{v}\left(\mathbf{t} + \frac{1}{2}\mathcal{A}\mathbf{t}\right) = \mathbf{\dot{r}}\left(\mathbf{t} - \frac{1}{2}\mathcal{A}\mathbf{t}\right) + \mathbf{\ddot{r}}\left(\mathbf{t}\right)\mathcal{A}\mathbf{t} \tag{15}$$

In this method, the velocities leap over the positions and then, the positions leap over the velocities. Consequently, the positions and the velocities are not known simultaneously, but the velocities are calculated explicitly.

#### *The Velocity Verlet Algorithm*

Using equations deduced from equations (10) and (13), and ignoring infinitesimals, this algorithm calculates new positions, velocities and accelerations using their values at time, t.

algorithm (Rahman, 1964, Gear & Tu, 1974, Gear & Watanabe, 1974) and the Beeman's algorithm (Beeman, 1976). All the above integration schemes make the assumption, that the positions, velocities and accelerations can be approximated using a Taylor series expansion:

<sup>1</sup> r(t t) r(t) r(t) t r(t) t ... <sup>2</sup>

Using this method, the next step of position can be predicted as follows;

 

storage requirements are modest, but it is not too accurate.

Where t is a finite time step, r is the position, v and .

acceleration; ...

<sup>r</sup> and ....

*The Basic Verlet Algorithm (Verlet, 1967)*

In the same way, the previous step;

*The Verlet Leapfrog Algorithm* 

Summing equations (11) and (12), we have

position steps. The equations are as follows;

the velocities are calculated explicitly.

*The Velocity Verlet Algorithm* 

<sup>1</sup> v(t t) r(t) r(t) t r(t) t ... <sup>2</sup>

<sup>1</sup> a(t t) r(t) r(t) t r(t) t ... <sup>2</sup>

. .. <sup>1</sup> <sup>2</sup> r(t t) r(t) r(t) t r(t) t ... <sup>2</sup>

. .. <sup>1</sup> <sup>2</sup> r(t t) r(t) r(t) t r(t) t ... <sup>2</sup>

.. <sup>2</sup> r(t t) 2r(t) r(t t) r(t) t

The method uses no explicit velocities, it is quite straightforward, easy to implement and its

In this algorithm, the velocities are calculated by taking the average value halfway between

. <sup>1</sup> r(t t) r(t) r(t t) t <sup>2</sup> 

. .. 1 1 v(t t) r(t t) r(t) t 2 2 

In this method, the velocities leap over the positions and then, the positions leap over the velocities. Consequently, the positions and the velocities are not known simultaneously, but

Using equations deduced from equations (10) and (13), and ignoring infinitesimals, this algorithm calculates new positions, velocities and accelerations using their values at time, t.

. .. <sup>2</sup>

(10)

(13)

(14)

(15)

r are the

r are the velocity, a and ..

(11)

(12)

 

> 

 

. .. ... <sup>2</sup>

.. ... .... <sup>2</sup>

r are the third and the fourth derivative of position with time.

$$\mathbf{r}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{r}(\mathbf{t}) + \dot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{1}{2}\ddot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^2 \tag{16}$$

$$\mathbf{v}\left(\mathbf{t} + \mathcal{A}\mathbf{t}\right) = \mathbf{\dot{r}}\left(\mathbf{t}\right) + \frac{1}{2}\mathbf{\ddot{r}}\left(\mathbf{\dot{r}}\left(\mathbf{t}\right) + \mathbf{\ddot{r}}\left(\mathbf{t} + \mathcal{A}\mathbf{t}\right)\right)\mathcal{A}\mathbf{t} \tag{17}$$

The velocity Verlet algorithm requires low memory.

*The Predictor-Corrector Algorithm (Rahman, 1964; Gear & Tu, 1975; Gear & Watanabe, 1974)* 

Velocities at time t+t are predicted and the forces calculated, and then the corrected forms of the velocities are later calculated. Combining the *r(t +*t*)* and *a(t +*t*)* in equation 10, the position can be expressed as:

$$\mathbf{r}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{r}(\mathbf{t}) + \dot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{2}{3}\ddot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^{2} - \frac{1}{6}\ddot{\mathbf{r}}(\mathbf{t} - \mathcal{A}\mathbf{t})\mathcal{A}\mathbf{t}^{2} + \mathbf{Q}(\mathcal{A}\mathbf{t}^{4}) \tag{18}$$

The velocities at time t = t+t are then calculated from the positions. Combining the *r(t +* t*)* and *a(t +*t*)* in equation (18), the position can be expressed as:

$$\text{(Predicted)}\quad \mathbf{v}(\mathbf{t}+\mathcal{A}\mathbf{t}) = \mathbf{\bar{r}}(\mathbf{t}) + \frac{3}{2}\mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t} - \frac{1}{2}\mathbf{\bar{r}}(\mathbf{t}-\mathcal{A}\mathbf{t})\mathcal{A}\mathbf{t} \tag{19}$$

The acceleration at t = t+t are calculated from the positions and the predicted velocities.

$$\text{(Corrected)}\ \mathrm{v}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{\dot{r}}(\mathbf{t}) + \frac{1}{3}\mathbf{\dot{r}}(\mathbf{t} + \mathcal{A}\mathbf{t})\mathbf{A}\mathbf{t} + \frac{\mathbf{5}}{6}\mathbf{\dot{r}}(\mathbf{t})\mathbf{A}\mathbf{t} - \frac{1}{6}\mathbf{\dot{r}}(\mathbf{t} - \mathcal{A}\mathbf{t})\mathbf{A}\mathbf{t} \tag{20}$$

This algorithm has the advantage that, by comparing the predicted and the corrected values of parameters, it is possible to perform a self-check on the algorithm for accuracy.

#### *The Beeman's Algorithm*

The Beeman's algorithm (Beeman, 1976) is based on equations (21) and (22), which can be deduced from the equations (18) and (20). The algorithm is more complicated and requires more memory than the velocity Verlet, but it provides more accurate expression for the velocities and better energy conservation.

$$\mathbf{r}(\mathbf{t} + \mathcal{A}\mathbf{t}) = \mathbf{r}(\mathbf{t}) + \dot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{2}{3}\mathbf{\bar{r}}(\mathbf{t})\mathcal{A}\mathbf{t}^2 - \frac{1}{6}\mathbf{\bar{r}}(\mathbf{t} - \mathcal{A}\mathbf{t})\mathcal{A}\mathbf{t}^2 \tag{21}$$

$$\mathbf{r}\,\mathbf{v}(\mathbf{t}+\mathcal{A}\mathbf{t}) = \dot{\mathbf{r}}(\mathbf{t}) + \dot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t} + \frac{\nabla}{6}\ddot{\mathbf{r}}(\mathbf{t})\mathcal{A}\mathbf{t} - \frac{1}{6}\ddot{\mathbf{r}}(\mathbf{t}-\mathcal{A}\mathbf{t})\mathcal{A}\mathbf{t} \tag{22}$$

When choosing an integration algorithm following factors need to be considered:


### **2.2.4 Initialize the model**

To initialize the simulation, the MD 'box' (the control volume) must be defined; then initial positions and velocities of the atoms must be assigned – this is a kind of initial randomization. Positions of the atoms can be defined, by assuming certain crystal structure and the initial velocities can be randomly assigned.

### **2.2.5 Relax the model from its initial state to its dynamically equilibrium condition**

The model has to be relaxed from its artificially assigned initial conditions to its natural, dynamically equilibrium condition, to mimic real materials. This can be achieved by running the MD program under constant temperature for a pre- determined time steps, so that the velocities of atoms that are initially assigned randomly or based on a normal distribution, will gradually reach equilibrium at the specified temperature of the simulation (Cheong et al., 2001). This relaxation time steps will vary, depending on the time needed for the model to reach its natural, dynamically equilibrium state, that is consistent with the environmental temperature.

### **2.2.6 Run the simulation and analyze the results**

The simulation is then run and the results analyzed. Using the MD, the effect of such variables as edge effect, depth of cut, velocity and other machining parameters can be defined and the simulations conducted accordingly. It is also easy to vary the properties of the work materials and the cutting tools in MD simulations (Komanduri & Raff, 2001).

From the results of the simulation, there are various quantities can that be derived for nano machining predictions. Some of which are cutting forces, cutting temperature, potential and kinetic energies, pressure etc.

### **2.3 Cutting forces**

In nanometric machining, the cutting forces are the interatomic forces (Luo et al., 2003). These forces are the superposition of the interaction forces between the cutting tool and the workpiece atoms. A low cutting force is as a result of fine cutting conditions, which will in turn decrease the vibration of the cutting system and then result in better surface roughness (Jackson, 2008).

Fig. 2. Atomistic Interaction in Nanometric Machining (Ikawa et al., 1991, Promyoo et al., 2008)

To initialize the simulation, the MD 'box' (the control volume) must be defined; then initial positions and velocities of the atoms must be assigned – this is a kind of initial randomization. Positions of the atoms can be defined, by assuming certain crystal structure

**2.2.5 Relax the model from its initial state to its dynamically equilibrium condition** 

The model has to be relaxed from its artificially assigned initial conditions to its natural, dynamically equilibrium condition, to mimic real materials. This can be achieved by running the MD program under constant temperature for a pre- determined time steps, so that the velocities of atoms that are initially assigned randomly or based on a normal distribution, will gradually reach equilibrium at the specified temperature of the simulation (Cheong et al., 2001). This relaxation time steps will vary, depending on the time needed for the model to reach its natural, dynamically equilibrium state, that is consistent with the

The simulation is then run and the results analyzed. Using the MD, the effect of such variables as edge effect, depth of cut, velocity and other machining parameters can be defined and the simulations conducted accordingly. It is also easy to vary the properties of the work materials and the cutting tools in MD simulations (Komanduri & Raff, 2001).

From the results of the simulation, there are various quantities can that be derived for nano machining predictions. Some of which are cutting forces, cutting temperature, potential and

In nanometric machining, the cutting forces are the interatomic forces (Luo et al., 2003). These forces are the superposition of the interaction forces between the cutting tool and the workpiece atoms. A low cutting force is as a result of fine cutting conditions, which will in turn decrease the vibration of the cutting system and then result in better surface roughness

Fig. 2. Atomistic Interaction in Nanometric Machining (Ikawa et al., 1991, Promyoo et al., 2008)

**2.2.4 Initialize the model** 

environmental temperature.

kinetic energies, pressure etc.

**2.3 Cutting forces** 

(Jackson, 2008).

and the initial velocities can be randomly assigned.

**2.2.6 Run the simulation and analyze the results** 

(a) (b)

Figure 2 shows the time increment *(Δt*), in which every atom changes its position and interacts with its surrounding neighbour atoms in a manner that is determined from the interatomic potential function.

E.g. for atoms described by Lennard Jones potential,*Vij* , the interatomic force between atoms i and j is

$$F(r\_{ij}) = -\frac{dV\_{ij}}{dr\_{ij}}\tag{23}$$

The force acting on the i-th workpiece atom is thus a summation of the interaction with the surrounding atoms (Luo et al., 2003);

$$F\_{wi} = \sum\_{j \neq 1}^{N\_t} F\_{wtij} + \sum\_{j \neq 1}^{N\_w} F\_{wwij} = \sum\_{j \neq 1}^{N\_t} -\frac{dV(r\_{wtij})}{dr\_{wtij}} + \sum\_{j \neq 1}^{N\_w} -\frac{dV(r\_{wwij})}{dr\_{wwij}} \tag{24}$$

Similarly, the force on each tool atom is

$$F\_{ti} = \sum\_{j \neq 1}^{N\_t} F\_{ttij} + \sum\_{j \neq 1}^{N\_w} F\_{twij} = \sum\_{j \neq 1}^{N\_t} -\frac{dV(r\_{ttij})}{dr\_{ttij}} + \sum\_{j \neq 1}^{N\_w} -\frac{dV(r\_{twij})}{dr\_{twij}} \tag{25}$$

#### **2.4 Cutting temperature**

In MD simulations, it is assumed that the cutting energy is totally transferred into cutting heat and this results in an increase of the cutting temperature and kinetic energy of the system.

#### **2.5 Minimum depth cut**

The Minimum Depth of Cut (MDC) is defined as the minimum undeformed chip thickness that can be removed stably from a work surface at a cutting edge under perfect performance of a machine tool (Ikawa et al., 1992). The concept of MDC is that the depth of cut must be over a certain critical thickness before any chip is formed. This phenomenon of MDC leads to a rising of slipping forces, burr formation and surface roughness (Ducobu et al., 2009). Conventionally, the tool - workpiece material interface has been considered to be homogeneous and continuum mechanics are used in the analysis of the MDC. In nanomachining, analysis is based on discrete atoms whose interactions are governed by appropriate interatomic potentials. The understanding and the accurate prediction of the MDC is very crucial in improving the ultra-precision metal removal technologies, as this would assist in the selection of appropriate machining parameters and optimal geometry design (Oluwajobi & Chen, 2010b).

#### **2.6 Nanometric cutting versus conventional cutting mechanics**

Nanometric cutting is based on a very small region of the tool-workpiece interface, which contains few atoms and so discrete mechanics apply. On the other hand, conventional cutting is based on continuum mechanics. The cutting forces are based on the interatomic potentials in nanometric cutting, but they are dependent on the plastic deformation in conventional cutting. Table 2 shows the comparison of nanometric cutting and conventional


cutting mechanics, and it should be noted that the basic physics of the cutting mechanism on the nanoscale is quite different from that on the macro/conventional scale.

Table 2. Comparison of Nanometric Cutting and Conventional Cutting Mechanics (Adapted from Luo et al., 2003)

### **3. The simulation set-up and procedure**

There are many issues to address for the MD simulation set-up, namely; the software and the hardware. For the hardware, suitable computer system configuration should be utilized depending on whether serial or parallel processing is required. For the software, there are many open source MD simulation software that are available. Two of the most popular software are DLPOLY and LAMMPS. (Another option is to develop one's own MD software from the scratch.) There are also open source software for the visualization of the simulation like the VMD. A typical structure can be seen in Figure 3, which comprises the preprocessing, the main MD processing and the post-processing parts.

Fig. 3. Software Methodology Flowchart

### **3.1 Pre-processing**

A suitable software has to be used to set-up the simulation model. A way to do this is to place the atoms on a perfect crystal lattice structure that represents the atomic structure of the material to be studied. The positions of atoms from an earlier simulation can also be used to initialize a new simulation (Cheong et al., 2001).

### **3.2 MD processing**

400 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

cutting mechanics, and it should be noted that the basic physics of the cutting mechanism on

**Nanometric Cutting Conventional Cutting** 

Discrete Molecular Mechanics Continuum Mechanics

Shear/Friction Power

(Grain Boundary Void)

Inter Crystal Deformation

the nanoscale is quite different from that on the macro/conventional scale.

**Energy Consideration** Interatomic Potential

**3. The simulation set-up and procedure** 

Fig. 3. Software Methodology Flowchart

used to initialize a new simulation (Cheong et al., 2001).

**3.1 Pre-processing** 

Adequate Software

**Chip Formation** Inner Crystal Deformation

processing, the main MD processing and the post-processing parts.

Pre-processing Main MD Processing Post-processing

A suitable software has to be used to set-up the simulation model. A way to do this is to place the atoms on a perfect crystal lattice structure that represents the atomic structure of the material to be studied. The positions of atoms from an earlier simulation can also be

LAMMPS/

DLPOLY VMD

**Workpiece Material** Heterogeneous Homogeneous **Cutting Physics** Atomic Cluster Model Shear Plane Model

Functional

**Deformation and Stress** Discontinuous Continuous **Cutting Tool Edge Radius** Significant Ignored **Cutting Tool Wear** Cutting Face and Cutting Edge Rake Face

**Cutting Force** Interatomic Forces Plastic Deformation

(Point Defects or Dislocation)

Table 2. Comparison of Nanometric Cutting and Conventional Cutting Mechanics (Adapted

There are many issues to address for the MD simulation set-up, namely; the software and the hardware. For the hardware, suitable computer system configuration should be utilized depending on whether serial or parallel processing is required. For the software, there are many open source MD simulation software that are available. Two of the most popular software are DLPOLY and LAMMPS. (Another option is to develop one's own MD software from the scratch.) There are also open source software for the visualization of the simulation like the VMD. A typical structure can be seen in Figure 3, which comprises the pre-

**Fundamental Cutting** 

from Luo et al., 2003)

**Principles** 

The main MD software accepts the input data of the atomic configuration and then integrates the equation of motion for the interacting atoms.

### *LAMMPS– Large-Scale Atomic/Molecular Massively Parallel Simulator*

This is a classical MD software that models an ensemble of atoms using a variety of empirical potentials and boundary conditions. It runs on single and parallel processor computers. It is an open-source code and it is distributed under the terms of the GNU public licence. The LAMMPS computes the Newton's equation of motion for a system of interacting atoms. It requires as its input; the types of atoms and the list of their initial coordinates, molecular topology information and the empirical potentials assigned to all the atoms.

### *DL\_POLY*

This is a general purpose MD software for the simulation of a wide range of atomic and molecular systems. It affords the use of standard and user defined interatomic potentials for the selected system. It also accommodates many boundary conditions, namely; cubic periodic, slab, orthorhombic periodic, parallelepiped periodic etc. (Smith & Todorov, 2005).

### **3.3 Post-processing**

#### *VMD – Visual molecular dynamics*

The VMD is a molecular visualization software for displaying and analyzing atomic systems by using 3-D graphics technology. The VMD can be employed for the visual display of the LAMMPS/DL\_POLY MD simulation results.

### **3.4 The simulation configuration**

Various configurations are used for the simulation of nanoscale machining, to optimize computational resources. All the configuration share the same features in that they consist of the workpiece and the tool. Generally, the workpiece can be divided into boundary, Newtonian and thermostat atoms (See Figure 4).


Table 3. A Typical MD Simulation Parameters

Fig. 4. A Simulation Model

Table 3 shows the simulation conditions that can be applied in a study. The workpiece is crystalline copper and the tool is also crystalline diamond. The configuration has a total of 54232 atoms. The workpiece consists of 43240 copper atoms with FCC lattice It includes 3 kinds of atoms namely; boundary atoms, thermostat atoms and Newtonian atoms. The boundary atoms are kept fixed to reduce edge effects. The thermostat atoms conduct the heat generated during the cutting process out of the cutting region. This is achieved by the velocity scaling of the thermostat atoms, (with the conversion between the kinetic energy (KE) and temperature) via equation (26),

$$\sum\_{i} \frac{1}{2} m\_i v\_i^2 = \frac{3}{2} N k\_B T\_i \tag{26}$$

Where *mi* is the mass of the ith atom, *<sup>i</sup> v* is the resultant velocity of the ith atom, *N* is the number of the thermostat atoms, *Ti* is the temperature of the ith atom and *Bk* is the Boltzmann constant (1.3806504 x10-23 JK-1)

Whenever the temperature of the thermostat atoms exceeds the preset bulk temperature of 293K, their velocities are scaled by using equation (27),

$$\upsilon\_{i,new} = \upsilon\_i \sqrt{\frac{T\_{desired}}{T\_{current}}} \tag{27}$$

Where *Tcurrent* is the current temperature that is calculated from the KE and the *Tdesired* is the desired temperature.

Table 3 shows the simulation conditions that can be applied in a study. The workpiece is crystalline copper and the tool is also crystalline diamond. The configuration has a total of 54232 atoms. The workpiece consists of 43240 copper atoms with FCC lattice It includes 3 kinds of atoms namely; boundary atoms, thermostat atoms and Newtonian atoms. The boundary atoms are kept fixed to reduce edge effects. The thermostat atoms conduct the heat generated during the cutting process out of the cutting region. This is achieved by the velocity scaling of the thermostat atoms, (with the conversion between the kinetic energy

> 1 3 <sup>2</sup> 2 2 *ii Bi*

Where *mi* is the mass of the ith atom, *<sup>i</sup> v* is the resultant velocity of the ith atom, *N* is the number of the thermostat atoms, *Ti* is the temperature of the ith atom and *Bk* is the

Whenever the temperature of the thermostat atoms exceeds the preset bulk temperature of

, *desired i new i*

Where *Tcurrent* is the current temperature that is calculated from the KE and the *Tdesired* is the

*v v*

*T*

*current*

*m v Nk T* (26)

*<sup>T</sup>* (27)

*i*

Fig. 4. A Simulation Model

(KE) and temperature) via equation (26),

Boltzmann constant (1.3806504 x10-23 JK-1)

desired temperature.

293K, their velocities are scaled by using equation (27),

The Newtonian atoms obey the Newton's equation of motion. The cutting tool consists of 10992 carbon atoms with diamond lattice structure. In this case, the cutting tool is pointed shaped and it can be modelled as a rigid or as a non-rigid body.

The atomic interactions in the simulation are the following, namely;


Suitable interatomic potentials are required for all the above atomic interactions.

### **4. Some examples of MD simulation of nanomachining**

Belak and Stowers pioneered work on the MD simulation of nanometric machining (See Figure 5). Copper was used as the workpiece with 103 – 106 atoms at room temperature. The simulations were in 2D and the focus was on the chip formation process, mechanisms of plastic deformation, flow of material and the flow of energy into the chip and workpiece. They modelled the workpiece by using the EAM potential and they used a shifted and truncated LJ potential for the tool-workpiece interface. Shimada et al 1992 also carried out MD machining simulations using copper as the workpiece (Figure 6). They used the Morse potential for the atomic interactions and showed the chip formation for the simulation model of 6076 atoms at the cutting speed of 200m/s.

(Rentsch & Inasaki, 1994) modelled a copper work material and a diamond tool for their study. They used the Lennard-Jones potential function for the copper atom interactions, but kept the boundaries and the tool stiff. A total of 11476 atoms in 13 horizontal (1,1,1) – layers of fcc-lattice were used for the copper, and the tool was shaped from a diamond lattice block by clearing on the four (1,1,1) –planes. Using a cutting speed of 100m/s, they observed a pronounced build-up phenomenon after 25000 time steps (see Figure 7).

Fig. 5. MD Simulation of the Orthogonal Cutting of Copper; (a) with cutting speed of 2500m/s and (b) with cutting speed of 500m/s (Belak & Stowers, 1990)

Fig. 6. MD Simulation of Cutting of Copper (Shimada et al., 1992)

Fig. 7. Advanced MD Simulation with Straight Aligned Tool (Rentsch & Inasaki, 1994)

The MD simulation of nanometric cutting was carried out by Komanduri team with a range of negative-rake tools to simulate the Ultra-Precision Grinding (UPG) process (Komanduri et al., 1999). A copper work material and an infinitely hard tool (tungsten) were used in the simulations. A pairwise sum of Morse potential was used for the study, in which they concluded that simulation tests can facilitate a better understanding of the process without the need for expensive and time consuming machining or grinding experimental work.

Fig. 6. MD Simulation of Cutting of Copper (Shimada et al., 1992)

experimental work.

Fig. 7. Advanced MD Simulation with Straight Aligned Tool (Rentsch & Inasaki, 1994)

The MD simulation of nanometric cutting was carried out by Komanduri team with a range of negative-rake tools to simulate the Ultra-Precision Grinding (UPG) process (Komanduri et al., 1999). A copper work material and an infinitely hard tool (tungsten) were used in the simulations. A pairwise sum of Morse potential was used for the study, in which they concluded that simulation tests can facilitate a better understanding of the process without the need for expensive and time consuming machining or grinding The investigation of the fundamental atomistic processes in chemical mechanical polishing of copper was carried out by Ye et al (2002). They simulated the nanoscale polishing of a copper surface by a single abrasive particle, using the embedded-atom potential. The temperature was controlled by maintaining 1.2nm of the substrate at 300K and the rescaling of atom velocities was performed whenever the temperature deviated more than 10K from the specified value. This allowed the transfer of heat from the machined region to the bulk of the work-material. They focused on the mechanical abrasion aspect of material removal and found that dislocations and atomically rough planarized surface were formed. They also studied the nature of the material removal, chip formation, materials defects and frictional forces as functions of the cutting speed, depth of cut and abrasive geometry. They established that the material removal rate scales linearly with the depth of planarization and is directly proportional to the velocity of cut.

To be more realistic, (Rentsch & Inasaki, 2006) extended the MD material modelling to consider fluids, like coolants. They considered the impact of such fluids on the surface generation and the tribological contact conditions. The fluid-fluid interactions were calculated on the basis of the Lennard-Jones potential function and the embedded-atom potential function was used for the inner workpiece reactions (the internal tool dynamics were ignored). They observed an intensive self-diffusion of the fluid atoms, and these filled the whole free space above the workpiece. No impact on the stress distribution was observed, but the whole fluid-tool/workpiece contact was heated up in a narrow range (See Figures 8 (a) and (b).

Fig. 8. MD Machining Simulation (a) without fluid (vacuum) after 230.7ps (b) with fluid after 230ps , (Rentsch & Inasaki, 2006)

Shimizu et al (2006) carried out MD simulations on the effect of vibration, velocity and acceleration in vibration assisted cutting. They used aluminium as the workpiece and diamond as the tool. Using Morse potentials for the interactions, they obtained that the effect of vibration on the plastic flow and cutting forces is more than the effect due to acceleration. See some snapshots of the simulation in Figure 9.

To extend the capability of MD studies, Pei et al., (2009), performed large scale MD simulations with model size of up to 10million atoms (Figures 10 (a), (b) and (c). Their results showed that as the cutting depth decreases, the tangential cutting force decreases and that size effects exists in nanometric cutting.

Oluwajobi & Chen, (2010a, 2011a, 2010b, 2011b, 2012) have conducted several MD simulations of nanometric machining. To show the effect of interatomic potentials on nanomachining, three popular potentials namely; EAM, Morse and the Lennand-Jones, were employed to model nanometric machining. The following shows the results;

### *Modelling with LJ Potential*

The LJ parameters used for the atomic interactions are = 2.2277 Angstroms and = 0.415eV (Hwang et al 2004), which apply to both the Cu-Cu and the Cu-C interactions. The cutting forces are shown in Figure 11.

Fig. 9. Snapshots of Atomic Arrays in Vibration-assisted Cutting Process and Travelling Distance from Initial Arrays. Cutting speed, *V ms <sup>c</sup>* 50 / , Amplitude, A=4nm and Frequency, f =4GHz (Shimizu et al., 2006).

Shimizu et al (2006) carried out MD simulations on the effect of vibration, velocity and acceleration in vibration assisted cutting. They used aluminium as the workpiece and diamond as the tool. Using Morse potentials for the interactions, they obtained that the effect of vibration on the plastic flow and cutting forces is more than the effect due to

To extend the capability of MD studies, Pei et al., (2009), performed large scale MD simulations with model size of up to 10million atoms (Figures 10 (a), (b) and (c). Their results showed that as the cutting depth decreases, the tangential cutting force decreases

Oluwajobi & Chen, (2010a, 2011a, 2010b, 2011b, 2012) have conducted several MD simulations of nanometric machining. To show the effect of interatomic potentials on nanomachining, three popular potentials namely; EAM, Morse and the Lennand-Jones, were

= 0.415eV (Hwang et al 2004), which apply to both the Cu-Cu and the Cu-C interactions.

Fig. 9. Snapshots of Atomic Arrays in Vibration-assisted Cutting Process and Travelling Distance from Initial Arrays. Cutting speed, *V ms <sup>c</sup>* 50 / , Amplitude, A=4nm and

= 2.2277 Angstroms and

employed to model nanometric machining. The following shows the results;

acceleration. See some snapshots of the simulation in Figure 9.

The LJ parameters used for the atomic interactions are

and that size effects exists in nanometric cutting.

The cutting forces are shown in Figure 11.

Frequency, f =4GHz (Shimizu et al., 2006).

*Modelling with LJ Potential* 

Fig. 10. MD Simulation Models with Number of the Atoms in the Workpiece around (a) 2 Million, (b) 4 Million and (c) 10 Million (Pei et al., 2009)

Fig. 11. Cutting Forces for LJ Potential (Oluwajobi & Chen, 2011a)

*Modelling with Morse Potential* 

*For Cu-Cu interactions:* (Girifalco & Weizer, 1959, Pei et al., 2006)

<sup>1</sup> 0.3429 , 0.13588( ) , 0.2866 *D eV nm r nm <sup>e</sup>* 

*For Cu-C interactions:* (Hwang et al., 2004)

<sup>1</sup> 0.087 , 0.17( ) , 0.22 *D eV nm r nm e* 

The cut-off distance chosen was 6.4 Angstroms (that is, the interactions between atoms separated by more than this distance are neglected). The cutting forces are shown in Figure 12.

*Modelling with EAM Potential* 

For the EAM potential, parameters used can be found in (Oluwajobi & Chen 2011a) and the cutting forces are shown in Figure 13.

Fig. 12. Cutting Forces for Morse Potential (Oluwajobi & Chen, 2011a)

Fig. 13. Cutting Forces for EAM Potential (Oluwajobi & Chen, 2011a)

From Figures 11 – 13, it can be observed that the cutting forces variation is lowest for the EAM potential. The results of the EAM model are comparable with those of (Pei et al., 2006, Promyoo et al., 2008), and the EAM potential best describes the metallic bonding in the copper atoms. In contrast, the pair potentials (both LJ and Morse potentials), do not incorporate the many-body effects; they have no environmental dependence and they do not account for the directional nature of bonding in metals (Li et al., 2008). Through this investigation, it was identified that the EAM potential is the most appropriate of the 3 potentials commonly used for the modelling of nanomachining of copper with diamond tool. This is because the EAM potential provides the best description of the metallic bonding in the workpiece, also, the cutting forces variation is smallest; the potential and total energies are most stable for the depth of cut considered. The study recommended that the EAM potential should be used, rather than LJ and Morse potentials for the modelling of copper and other fcc metals in MD simulations of nanomachining (Oluwajobi & Chen, 2011a).

**Cutting Forces for Morse Potential**

0 20000 40000 60000 80000 100000 120000

Fx Fy Fz

Fx Fy Fz

**No of Steps**

**Cutting Forces for EAM Potential**

0 20000 40000 60000 80000 100000 120000

**No of Steps**

From Figures 11 – 13, it can be observed that the cutting forces variation is lowest for the EAM potential. The results of the EAM model are comparable with those of (Pei et al., 2006, Promyoo et al., 2008), and the EAM potential best describes the metallic bonding in the copper atoms. In contrast, the pair potentials (both LJ and Morse potentials), do not incorporate the many-body effects; they have no environmental dependence and they do not account for the directional nature of bonding in metals (Li et al., 2008). Through this investigation, it was identified that the EAM potential is the most appropriate of the 3 potentials commonly used for the modelling of nanomachining of copper with diamond tool. This is because the EAM potential provides the best description of the metallic bonding in the workpiece, also, the cutting forces variation is smallest; the potential and total energies are most stable for the depth of cut considered. The study recommended that the EAM potential should be used, rather than LJ and Morse potentials for the modelling of copper and other fcc metals in MD

Fig. 12. Cutting Forces for Morse Potential (Oluwajobi & Chen, 2011a)

Fig. 13. Cutting Forces for EAM Potential (Oluwajobi & Chen, 2011a)

simulations of nanomachining (Oluwajobi & Chen, 2011a).

**Cutting Forces (eV/Angs)**

**Cutting Forces (eV/Angs)**

Fig. 14. Groove Profiles for Different Cutting Directions (a) 0deg (b) 30deg (c) 45deg (d) 60deg (e) 90deg (Zhang et al., 2009)

Many existing MD simulation studies on nanometric cutting have been limited to single pass or simple line-type groove. As an extension of the single pass studies, Zhang et al., (2009) modelled folder- line grooves for AFM-based nanometric cutting process of copper workpiece with diamond tool (Figure 14). They used the EAM potential for the coppercopper interactions and the Morse potential for the copper-diamond interactions. They treated the diamond tool as rigid and concluded that the normal, lateral and the resultant forces were almost symmetric with respect to the critical folder angle of 45. Shi et al., (2011) investigated the multi-groove simulation of single-point tuning of copper workpiece with diamond tool. They used two diamond tools, offset by a fixed distance to simulate a two-groove cutting and modelled the copper-copper and the copper-diamond interactions by using the Morse potential (Figure 15). They also treated the tool as a rigid body and observed that the tool forces increase with increase in feed rate and depth of cut. In practice, most machining processes involve the use of multiple passes to create new surface patterns and the diamond tool is deformable and subjected to wear. To address the above problem, novel multi-pass nanometric machining MD simulations have been carried out, which show the consecutive passes of cut. See Figures 16 and 17 (Oluwajobi & Chen, 2011b).

Fig. 15. MD Simulation of Multi-groove Cutting with Two Diamond Tools: (a) 3D perspective view (b) Stress Distribution (Shi et al., 2011)

Many existing MD simulation studies on nanometric cutting have been limited to single pass or simple line-type groove. As an extension of the single pass studies, Zhang et al., (2009) modelled folder- line grooves for AFM-based nanometric cutting process of copper workpiece with diamond tool (Figure 14). They used the EAM potential for the coppercopper interactions and the Morse potential for the copper-diamond interactions. They treated the diamond tool as rigid and concluded that the normal, lateral and the resultant forces were almost symmetric with respect to the critical folder angle of 45. Shi et al., (2011) investigated the multi-groove simulation of single-point tuning of copper workpiece with diamond tool. They used two diamond tools, offset by a fixed distance to simulate a two-groove cutting and modelled the copper-copper and the copper-diamond interactions by using the Morse potential (Figure 15). They also treated the tool as a rigid body and observed that the tool forces increase with increase in feed rate and depth of cut. In practice, most machining processes involve the use of multiple passes to create new surface patterns and the diamond tool is deformable and subjected to wear. To address the above problem, novel multi-pass nanometric machining MD simulations have been carried out, which show the consecutive passes of cut. See Figures 16 and 17

Fig. 15. MD Simulation of Multi-groove Cutting with Two Diamond Tools: (a) 3D

perspective view (b) Stress Distribution (Shi et al., 2011)

(Oluwajobi & Chen, 2011b).

Fig. 16. a. Cross Section of the Machined Grooves with Passes 1-3 (direction of cut is perpendicular to the paper face) 16b: Tool Tip Dimensions (Oluwajobi & Chen, 2011b)

Fig. 17. MD Simulation of Multipass Nanometric Machining (a) pass 1 (b) pass 2 (c) pass 3 (Oluwajobi & Chen, 2011b)

From the simulations, it was observed that the magnitude of the tangential and the normal cutting force components increase with the increase in the depth of cut and they decrease with the consecutive passes. Also, the lateral force components are influenced by atomic vibrations and the cutting cross sectional area (Oluwajobi & Chen, 2011b, 2012).

## **5. Conclusions**

The MD simulation of nanoscale machining has come a long way from inception; from the use of 'primitive' two dimensional systems to two and a half dimensional system and now to 'realistic' three dimensional systems. More developments have occurred from the first simulations where not many or no quantitative analyses were carried out, to now, where rigorous quantitative analyses for the evaluation of the machine process parameters are the norm. Also, the overall configuration have increased from a few thousands to millions atoms at present. The commonly used potentials for MD simulations of metals include LJ, Morse, EAM and MEAM. In many cases, the rationale or justification for their use are not given except for the reference to their use in earlier studies. Ideally, multi-body potentials (EAM and MEAM) should be used for metals, rather than pairwise potentials (LJ and Morse); because they more correctly model the metallic bonding. It is important to use appropriate interatomic potentials for MD simulations so as to obtain useful results. It is now feasible to carry out multipass nanometric machining MD simulations. From the results of such simulations, minimum surface roughness can be effectively evaluated, which is very critical in the industry.

At present, MD methods can only handle small length and time scales. This is far from what obtains during experimental studies. To tackle these problems, multiscale methods are been used. Spatial and temporal multiscale techniques would considerably extend the length and time scales in simulations respectively (Medyanik, 2007, Medyanik & Liu, 2008).

## **6. Nomenclature**


From the simulations, it was observed that the magnitude of the tangential and the normal cutting force components increase with the increase in the depth of cut and they decrease with the consecutive passes. Also, the lateral force components are influenced by atomic

The MD simulation of nanoscale machining has come a long way from inception; from the use of 'primitive' two dimensional systems to two and a half dimensional system and now to 'realistic' three dimensional systems. More developments have occurred from the first simulations where not many or no quantitative analyses were carried out, to now, where rigorous quantitative analyses for the evaluation of the machine process parameters are the norm. Also, the overall configuration have increased from a few thousands to millions atoms at present. The commonly used potentials for MD simulations of metals include LJ, Morse, EAM and MEAM. In many cases, the rationale or justification for their use are not given except for the reference to their use in earlier studies. Ideally, multi-body potentials (EAM and MEAM) should be used for metals, rather than pairwise potentials (LJ and Morse); because they more correctly model the metallic bonding. It is important to use appropriate interatomic potentials for MD simulations so as to obtain useful results. It is now feasible to carry out multipass nanometric machining MD simulations. From the results of such simulations, minimum surface roughness can be effectively evaluated, which is very

At present, MD methods can only handle small length and time scales. This is far from what obtains during experimental studies. To tackle these problems, multiscale methods are been used. Spatial and temporal multiscale techniques would considerably extend the length and

3 3 *d r dt ij b A parameter that provides information for the direction and the length of bond* 

> 4 4 *d r dt*

time scales in simulations respectively (Medyanik, 2007, Medyanik & Liu, 2008).

2 2 *<sup>i</sup> d r dt ij a Function that limit the range of interaction in Tersoff's potential* 

vibrations and the cutting cross sectional area (Oluwajobi & Chen, 2011b, 2012).

**5. Conclusions** 

critical in the industry.

**6. Nomenclature** 

*<sup>i</sup> a Acceleration of atom i,* 

*A Material dependent constant* 

*b Third derivative of position,* 

*c Fourth derivative of position,* 

*D Material dependent constant* 

*E Energy of interacting particles Etot Total embedding energy f Material dependent constant* 

 *(Bond strength)* 

*e Elementary Charge* 



### **7. References**


ttij r *The displacement vector of the i-th cutting tool atom and the interacting cutting* 

V (r) <sup>e</sup> *Electronic energy calculated by considering the two nuclei as fixed in space a* 

 *Effective co-ordination parameter (The count of the number of other bonds to atom* 

twij r *The displacement vector of the i-th cutting tool atom and the interacting* 

U (r) ij *Short range pair interaction representing the core-core repulsion* 

*V Short range pair interaction representing the core-core repulsion i*, *<sup>j</sup>*

*r Atomic electron density of atom j at the distance ij r from the nucleus* 

*h i*, *Total electron density at atom i due to the rest of the atoms in the system* 

(,) *r t Wave function (The probability of finding any particle at position r, at time t)* 

Alder B.J.and T.E. Wainwright (1957), Phase Transition for Hard Sphere System*, Journal of* 

Alder B.J.and T.E. Wainwright (1959), Studies in Molecular Dynamics. I. General Method,

Allen M.P. and D.J. Tildesley (1988), *Computer Simulation of Liquids*, (Oxford University

Beeman D. (1976), Some Multistep Methods for use in Molecular Dynamics Calculations,

Belak, J. and Stowers I. F. (1990), A Molecular Dynamics Model of the Orthogonal

Press) Baskes M.I. (1992), Modified Embedded-Atom Potentials for Cubic Materials

Cutting Process, *Proceedings of the American Society of Precision Engineering*, pp 76-

Va *Potential due to attractive forces between atoms i and j* 

*Vr Potential due to repulsive forces between atoms i and j* 

 *workpiece atom j* 

V ,V(r ) ij ij *Interatomic Pair Potential* 

 *tool atom j* 

v *Velocity* 

*ij* 

( ) *a j* 

**7. References** 

79.

 *distance r apart* 

*W Allowed energy levels* 

 *i,besides the ij bond).* 

 *Average electron density* 

( )*r Electron density* 

*t Finite time step* 

 *Material dependent constant* 

 *Material dependent constant* 

*ijk The bond angle between bonds ij and ik* 

*Chemical Physics,* Vol. 27, pp 1208-1209

*Journal of Chemical Physics*, Vol. 31 pp 459-466

and Impurities, *Physical Review B*, Vol. 46, pp 2727-2742

*Journal of Computational Physics* Vol. 20 pp 130-139


Hernandez E.R. (2008), Molecular Dynamics from Basic Techniques to Applications: A

Hwang H.J., O-K Kwon & J. W. Kang (2004), Copper Nanocluster Diffusion in Carbon

Ikawa N., S. Shimada, H. Tanaka & G. Ohmori (1991*)*, An Atomistic Analysis of Nanometric

Ikawa N., S. Shimada & H. Tanaka (1992), Minimum Thickness of Cut in Micromachining,

Jackson M.J. (2008), Micro and Nanomachining, In J.P. Davin – *Machining: Fundamentals and* 

Kenny S.D., D. Mulliah, C.F. Sanz-Navarro & R. Smith (2005), Molecular Dynamics

Khukhryansky Y. P., E.A. Shunikov and V.V. Emelyanov (2004), Study of Cluster Forming

Komanduri R., and L.M. Raff (2001), A Review on the Molecular Dynamics Simulation of

Kragelsky I.V., Dobychin M.N., & Kombalov V.S. (1982), *Friction and Wear-Calculation* 

Lee B, B.D.Wirth, J. Shim, J. Kwon, S.C. Kwon & J. Hong (2005), Modified Embedded-

Lennard-Jones J.E. (1924), On the Forces between Atoms and Ions, *Proc. Royal Soc*. Vol. 109,

Li J.H., X.D. Dai, S.H. Liang, K.P. Tai, Y. Kong & B.X. Lin (2008), Interatomic Potentials of

Luo X., K. Cheng, X. Guo & R. Holt (2003), An Investigation on the Mechanics of

Luo X. (2004), *'High Precision Surfaces Generation: Modelling, Simulation and Machining* 

Medyanik S.N. (2007), Atomistic Simulation Methods and Multiscale Modelling, Accessed

Medyanik S.N. and W.K. Liu (2008), Multiple Time Scale Method for Atomistic Simulations,

Morse P.M. (1929), Diatomic Molecules according to Wave Mechanics II Vibrational Levels,

Murrel J.N. and R.E. Mottram (1990), Potential Energy Functions for Atomic Solids,

October 2011 http://multiscale.emsl.pnl.gov/docs/atomistic.pdf

LAMMPS Manual, http://lammps.sandia.gov/doc/Manual.pdf, (Accessed in 2009) Leach A. (2001), *Molecular Modelling: Principles and Applications*, Prentice Hall

Nanotube, *Solid State Communications*, Vol. 129, pp 687-690

*Proceedings*, Vol. 1077, Iss.1, pp 95-123

*the CIRP*, Vol. 40, Iss. 1, pp 551-554

*Recent Advances*, pp 271-297, Springer

*Crystal Growth*, Vol. 269, pp 292-297

*Methods*, Pergamon Press, New York

*Physics Reports* ,Vol. 455, pp 1-134

*Physical Review,* Vol. 34, pp 57-64

*Molecular Physics*, Vol. 69, pp 571-588

*Production Research*, Vol.41, No.7, pp 1149-1165

*Computational Mechanics*, Vol. 42, pp 569-577

*Verification'*, PhD Thesis, Leeds Metropolitan University

Vol. 215, Part B, pp 1639-1672

(1-15)

pp 584-597

*Royal Society A*., Vol. 363, No 1833, pp 1949-1959

*Nanotechnology V*ol. 3, pp 6-9

Molecular Dynamics Primer, *Frontiers in Contemporary Physics, AIP Conference* 

Chip Removal as Affected by Tool-Work Interaction in Diamond Turning, *Annals of* 

Simulations of Nanoindentation and Nanotribology, *Philosophical Transactions of the* 

at Growth of A3B5 Semiconductor Compounds from Liquid Phase*, Journal of* 

Machining at the Atomic Scale, *Proceedings of the Institution of Mechanical Engineers*,

Atom Method Interatomic Potential for the Fe-Cu Alloy System and Cascade Simulations on Pure Fe and Fe-Cu Alloys, *Physical Review B*, Vol. 71, pp 184205

the Binary Transition Metal Systems and Some Applications in Materials Physics,

Nanometric Cutting and the Development of its Test-bed*, International Journal of* 


## **High-Throughput Simulations of Protein Dynamics in Molecular Machines: The 'Link' Domain of RNA Polymerase**

Robert O. J. Weinzierl *Imperial College London United Kingdom* 

### **1. Introduction**

418 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

Shimada S., N. Ikawa, H. Tanaka, G. Ohmori and J. Uchikoshi (1993), Feasibility Study on

Shimizu J. H. Eda, M. Yoritsune and E. Ohmura (1998), Molecular Dynamics Simulations of

Shimizu J., L. Zhou & H. Eda (2006), Molecular Dynamics Simulation of Vibration-Assisted

Smith R., D. Christopher & S. Kenny (2003), Defect Generation and Pileup of Atoms

Smith W. and I.T. Todorov (2005), The DL\_POLY Molecular Dynamics Package,*Z.* 

Stillinger F.H. and A. Rahman (1974), Improved Simulation of Liquid Water by Molecular

Stillinger F.H. and T.A. Weber (1985), Computer Simulation of Local Order in Condensed

Sutton A.P. and J. Chen (1990), Long-Range Finnis Sinclair Potentials, *Philosophical Magazine* 

Tersoff J. (1988a), New Empirical Approach for the Structure and Energy of Covalent

Tersoff J. (1988b), Empirical Interatomic Potential for Silicon with Improved Elastic

Tuckerman M.E. and G.L. Martyna (2000), Understanding Modern Molecular Dynamics:

Van Gunsteren W. F. and H.J.C. Berendsen (1990), Computer Simulation of Molecular

Verlet L. (1967), Computer Experiments on Classical Fluids I. Thermodynamics Properties of

Visual Molecular Dynamics (VMD), http://www.ks.uiuc.edu/Research/vmd/ (Accessed in

Wang C., T. Yu, W. Duan & L.Wang (1995), A First Principles Interatomic Potential and

Zhang J.T. , T. Sun, Y. Yan, Y. Liang & S. Dong, Molecular Dynamics Study of Groove

*Physics A (Materials Science and Processing)*, Vol. 94, 2009, pp 593-600 Zhao K.J., C.Q. Chen, Y.P. Shen and T.J. Lu (2009), Molecular Dynamics Study on the Nano-

Lennard-Jones Molecules, *Physical Review D,* Vol.159, pp 98-103

Techniques and Applications, *Journal of Physical Chemistry*, Vol. 104, pp 159-

Dynamics: Methodology, Applications, and Perspectives in Chemistry, *Anew.* 

Application to the GrainBoundary in Ni, *Physics Letters A*, Vol. 197, Iss. 5-6, pp 449-

Fabrication Process using AFM-Based Nanometric Cutting Technique, *Applied* 

void Growth in Face-centered Cubic Single Crystal Copper, *Computational Materials* 

Friction on the Atomic Scale, *Nanotechnology*, Vol. 9, pp 118-123

Dynamics, *Journal of Chemical Physics*, Vol. 60, pp 1545-1557

Phases of Silicon, *Physical Review B*, Vol. 31, pp 5262-5271

Systems, *Physical Review B*, Vol. 37 No 12, pp 6991-7000

Properties, *Physical Review B*, Vol. 38 No 14, pp 9902-9905

*of the CIRP,* Vol. 42, No.1, pp 91-94

*Kristallogra*, Vol. 220, pp 563-566

*Letters*, Vol. 61, Iss. 3, pp 139-146

*Chem. Int. Ed. Engl.,* Vol. 29, pp 992-1023

*Science*, Vol. 46, pp 749-754

10

178

2010)

457

*Nanomanufacturing*, Vol. 1, No. 1, pp 105-116

Ultimate Accuracy in Microcutting using Molecular Dynamics Simulation, *Annals* 

Cutting: Influences of Vibration, Acceleration and Velocity*, International Journal of* 

during Nanoindentation of Fe Single Crystals, *Physical Review B*, Vol. 67, pp 245405-

Complex molecular machines are involved in all information-processing steps in cells. The handling of information-bearing macromolecules (predominantly nucleic acids and proteins) requires a range of catalytic capabilities, such as the sequence-specific synthesis of macromolecular entities from smaller precursors and further processing by breaking/joining of pre-existing molecular components. These catalytic activities have to be spatially and temporally precisely controlled to ensure accuracy and efficiency levels that are commensurate with the associated biological functions (Lyubimov et al., 2011). Unlike metabolic enzymes, which are often medium-sized enzymes with freely accessible active sites, the enzymes associated with replication, transcription, translation, recombination are typically large multi-subunit protein complexes capable of a complex conformational spectrum. This conformational spectrum manifests itself through the dynamic features of individual domains within these molecular machines: many of the domains act explicitly as nanomechanical elements by serving as allosteric sensors, molecular hinges and motor units (Bustamante et al., 2011; Heindl et al., 2011). A deeper understanding of the functions of molecular machines is currently a high-priority topic and ultimate success will strongly depend on a combination of 'wet-lab' experimental techniques (X-ray crystallography, NMR, spectroscopic methods, site-directed mutagenesis and chemical cross-linking) with sophisticated computational simulations. Fully atomistic molecular dynamics (MD) simulations offer the most comprehensive and detailed insights into the behavior of molecular entities and are therefore the method of choice for attempts to unravel the structural basis of molecular mechanisms (reviewed in Karplus & McCammon, 2002; Freddolino et al., 2010; McGeagh et al. 2011; Schlick et al., 2011).

The ultimate goal of MD studies is the complete and comprehensive simulation of molecular machines operating within a relevant time span. Since the functions of molecular machines critically depend on a series of tightly orchestrated conformational changes that typically occur in the micro- and milliseconds range, such a goal is not yet routinely achievable. Currently available computational resources are still insufficient for simulating the molecular dynamics of macromolecular entities for the length required for complex conformational changes to take place. We are therefore focusing on the more practicable goal of studying the dynamic properties of individual components of molecular machines in order to obtain a better understanding of their nanomechanical properties. In any natural or man-made machine, the mechanical properties of individual components determine both their functional capabilities, as well as their interactions with other parts. In many cases a description of these properties allows either a clearcut deduction of the role of such components within the overall mechanism, or offers at least valuable suggestions for further experimental approaches. Most importantly, once the nanomechanical properties of individual components are understood in sufficient detail, their mode of interactions within defined assemblies can be modelled with increased accuracy and confidence. Such approaches result in insights into complex molecular mechanisms generated from the coordinated interplay between multiple protein domains.

Efforts in my laboratory are predominantly focused on understanding the intricate mechanisms of an enzyme playing a key role in gene expression: RNA polymerase (RNAP). Cellular RNAPs are highly conserved in evolution (ranging from bacteria to archaea and eukaryotes) and consist of a number of different subunits (reviewed in Werner, 2008). The large size and intrinsic complexity of these enzymes presents a particular challenge to understand the various mechanisms involved, such as separation and re-annealing of the DNA strands and the template-directed synthesis of an RNA transcript. Although structures of RNAPs in various states/conformations have been determined (e.g. Zhang et al., 1999; Cramer et al., 2001; Gnatt et al., 2001; Vassylyev, 2002), we are currently still far from having a complete catalogue of the most important conformations. Maybe even more importantly, the intermediate structural transitions between stable structural states are likely to be kinetically short-lived and are therefore very unlikely to be captured in the crystallized state required for high-resolution structural studies (Kaplan & Kornberg, 2008). Having a detailed understanding of such short-lived and metastable intermediate structures will, however, be critical for piecing together the conformational changes required for the various biological activities of RNAPs. Within the context of a larger strategy aimed at investigating the nanomechanical mechanisms underlying RNAP function, we are pursuing a two-pronged strategy: a laboratory-based robotic high-throughput approach capable of studying the functional consequences of a large number of site-directed mutants (the "RNAP Factory"; Nottebaum et al., 2008), and a computer-based in-depth investigation of nanomechanical properties of domains (or domain assemblies) using MD methods. The combination of these two complementary methods has already been proven to be highly effective for obtaining further insights into the nanomechanical 'motor' of RNAPs, as represented by the catalytic site and its surrounding protein domains (Fig. 1).

Our most recent studies focused on a prominent component of the catalytic site, the Bridge Helix. This structure is a 35 amino acid long -helix which spans the catalytic site and multiple lines of evidence suggest that the Bridge Helix is, in conjunction with the Trigger Loop, predominantly involved in the translocation of nucleic acid substrates through the catalytic site. High-throughput mutagenesis experiments, combined with systematic molecular dynamics simulations, have demonstrated the presence of two hinges ('BH-HN' and 'BH-HC') which bend during the nucleotide addition cycle to propagate the DNA template strand and the nascent transcript through the active site (Klug, 2001; Bar-Nahum et al., 2005; Tan et al., 2008; Brueckner et al., 2009; Seibold et al., 2010; Weinzierl, 2010a and 2010b; Heindl et al., 2011).

practicable goal of studying the dynamic properties of individual components of molecular machines in order to obtain a better understanding of their nanomechanical properties. In any natural or man-made machine, the mechanical properties of individual components determine both their functional capabilities, as well as their interactions with other parts. In many cases a description of these properties allows either a clearcut deduction of the role of such components within the overall mechanism, or offers at least valuable suggestions for further experimental approaches. Most importantly, once the nanomechanical properties of individual components are understood in sufficient detail, their mode of interactions within defined assemblies can be modelled with increased accuracy and confidence. Such approaches result in insights into complex molecular mechanisms generated from the coordinated interplay between multiple

Efforts in my laboratory are predominantly focused on understanding the intricate mechanisms of an enzyme playing a key role in gene expression: RNA polymerase (RNAP). Cellular RNAPs are highly conserved in evolution (ranging from bacteria to archaea and eukaryotes) and consist of a number of different subunits (reviewed in Werner, 2008). The large size and intrinsic complexity of these enzymes presents a particular challenge to understand the various mechanisms involved, such as separation and re-annealing of the DNA strands and the template-directed synthesis of an RNA transcript. Although structures of RNAPs in various states/conformations have been determined (e.g. Zhang et al., 1999; Cramer et al., 2001; Gnatt et al., 2001; Vassylyev, 2002), we are currently still far from having a complete catalogue of the most important conformations. Maybe even more importantly, the intermediate structural transitions between stable structural states are likely to be kinetically short-lived and are therefore very unlikely to be captured in the crystallized state required for high-resolution structural studies (Kaplan & Kornberg, 2008). Having a detailed understanding of such short-lived and metastable intermediate structures will, however, be critical for piecing together the conformational changes required for the various biological activities of RNAPs. Within the context of a larger strategy aimed at investigating the nanomechanical mechanisms underlying RNAP function, we are pursuing a two-pronged strategy: a laboratory-based robotic high-throughput approach capable of studying the functional consequences of a large number of site-directed mutants (the "RNAP Factory"; Nottebaum et al., 2008), and a computer-based in-depth investigation of nanomechanical properties of domains (or domain assemblies) using MD methods. The combination of these two complementary methods has already been proven to be highly effective for obtaining further insights into the nanomechanical 'motor' of RNAPs, as represented by the catalytic

Our most recent studies focused on a prominent component of the catalytic site, the Bridge Helix. This structure is a 35 amino acid long -helix which spans the catalytic site and multiple lines of evidence suggest that the Bridge Helix is, in conjunction with the Trigger Loop, predominantly involved in the translocation of nucleic acid substrates through the catalytic site. High-throughput mutagenesis experiments, combined with systematic molecular dynamics simulations, have demonstrated the presence of two hinges ('BH-HN' and 'BH-HC') which bend during the nucleotide addition cycle to propagate the DNA template strand and the nascent transcript through the active site (Klug, 2001; Bar-Nahum et al., 2005; Tan et al., 2008; Brueckner et al., 2009; Seibold et al., 2010; Weinzierl, 2010a and

protein domains.

site and its surrounding protein domains (Fig. 1).

2010b; Heindl et al., 2011).

Fig. 1. Part of the nanomechanical 'motor unit' in the catalytic center of RNA polymerase. The Bridge Helix is shown in green (the two hinges BH-HN and BH-HC are red and blue, respectively), the F-Loop in grey and the Link domain in purple. The nucleic acid substrates are shown as stick models (DNA light blue, RNA pink, rNTP orange).

While the C-terminal hinge ('BH-HC') may be actively involved in pushing the nascent RNA out of the active site to create space for the next incoming ribonucleotide triphosphate (rNTP; Klug, 2001; Bar-Nahum et al., 2005; Tan et al., 2008; Brueckner et al., 2009), the role of the more recently discovered N-terminal hinge (BH-HN) is still poorly understood. In contrast to BH-HC, the BH-HN region is not in direct contact with any of the substrates, suggesting that BH-HN interacts with the catalytic site more indirectly via allosteric contacts (Weinzierl, 2010b). Of the domains surrounding BH-HN, the 'Link' domain emerges as the most likely candidate for playing such a communication role. As can be seen in Fig. 1, the Link domain consists of a short -helix which contacts directly the N-terminal portion of the Bridge. One particular residue present in the -helical portion of the Link domain (arginine 766 in the RPB2 subunit of *S. cerevisiae* RNAPII) is evolutionarily invariant and structural studies suggest that it plays a key role in contacting the -phosphate of the rNTP in the active site (Fig. 2).

A simple working hypothesis therefore suggests that the Bridge Helix N-terminal hinge can 'sense' - via the Link domain - whether the catalytic site is loaded with a correctly positioned rNTP before it is joined to the nascent transcript. When the - pyrophosphate group is cleaved off during nucleotide addition and diffuses out of the catalytic site, this particular contact is lost and may result in kinking of BH-HN. This model thus proposes that the Link domain acts as a conformational sensor that communicates the state of the catalytic site (empty, rNTP loaded, nucleotide addition and pyrophosphate release) in order to influence the conformational state of Bridge Helix. In mechanistic terms, completion of the nucleotide addition reaction signals to the Bridge Helix to initiate translocation of the DNA-RNA hybrid to create space for the next rNTP (Fig. 3)

Fig. 2. Arginine 766 (R766) in the RPB2 subunit of *S. cerevisiae* RNAPII contacts the rNTP loaded in the catalytic site. The R766 residue (light brown), rNTP (orange) and the two catalytic metal ions (grey) are shown in space-filling mode. The remainder of the color scheme is as in Fig. 1.

Fig. 3. Role of the Link domain as a conformational sensor for the detection of an rNTP in the catalytic site of RNAP.

Taking into account the possible role of the Link domain in communicating the rNTP loading state of the catalytic site to the Bridge Helix, we decided to investigate the mechanistic properties of this system in more detail. In particular, we are interested in studying the conformational spectrum of the Link domain to gain as much information as possible before initiating an experimental high-throughput mutagenesis approach. As pointed out earlier, it seems highly probable that molecular dynamics simulations are nowadays sufficiently reliable to deduce key nanomechanical properties inherent in a protein domain. We specifically wanted to simulate the Link domain extensively under highly standardized conditions to capture a broad spectrum of conformational changes the domain is capable of. We expected that the data obtained in this way would provide quantitative evidence for the relative stability of different parts of the domain and highlight regions that could potentially be involved in allosteric switching processes.

### **2. Implementation**

422 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

the conformational state of Bridge Helix. In mechanistic terms, completion of the nucleotide addition reaction signals to the Bridge Helix to initiate translocation of the DNA-RNA

Fig. 2. Arginine 766 (R766) in the RPB2 subunit of *S. cerevisiae* RNAPII contacts the rNTP loaded in the catalytic site. The R766 residue (light brown), rNTP (orange) and the two catalytic metal ions (grey) are shown in space-filling mode. The remainder of the color

**R766** 

Fig. 3. Role of the Link domain as a conformational sensor for the detection of an rNTP in

Taking into account the possible role of the Link domain in communicating the rNTP loading state of the catalytic site to the Bridge Helix, we decided to investigate the mechanistic properties of this system in more detail. In particular, we are interested in studying the conformational spectrum of the Link domain to gain as much information as possible before initiating an experimental high-throughput mutagenesis approach. As pointed out earlier, it seems highly probable that molecular dynamics simulations are nowadays sufficiently reliable to deduce key nanomechanical properties inherent in a protein domain. We specifically wanted to simulate the Link domain extensively under highly standardized conditions to capture a broad spectrum of conformational changes the

hybrid to create space for the next rNTP (Fig. 3)

scheme is as in Fig. 1.

the catalytic site of RNAP.

The implementation of the strategy outlined above requires a considerable computational and organisational infrastructure. Here I will discuss the various issues that need to be considered and will provide an overview of the methods we developed to map out an universally applicable approach that can also be used by other researchers to study different systems.

### **2.1 What do we want to achieve?**

Before proceeding further, we need to define the desirable outcome of our MD simulation strategy. In order to 'get to know' the dynamic properties of a protein domain, two key requirements must be met: First and most obvious, we need to simulate a domain for a sufficient amount of time so that we can gather quantitative data regarding the overall features of the domains, such as conformational spectrum and overall stability, local variations in secondary structure RMSD values and other information of interest. The simulation time that is considered to be 'sufficient' has increased over the last few years with the steady increase in computational resources available, but is nowadays generally considered to be in the 10-100 nanosecond range. For many stably folded domains, simulation over such a time-frame will allow data collection sufficient to quantitate, for example, the half-life of secondary structure elements. The second key requirement is the need to run multiple independent simulations. In previous work on the Bridge Helix domain (Weinzierl, 2010b; Heindl et al., 2011) we encountered structural features that behaved stochastically, i.e. not every simulation run displayed the full conformational spectrum. We found that kinking of the Bridge Helix hinges observed in some molecular dynamics production runs resulted in very stable alternative conformations that were essentially irreversible within the tens of nanosecond simulation time. Such 'trapped' conformations play almost certainly an important role in molecular switching processes occurring in the millisecond range, but would result in a biased picture of the conformational spectrum of a domain if only a small number of simulations was carried out. For practical purposes (and based on practical experience) we therefore typically perform approximately 50 independent simulations of 20 nanoseconds each to maximize the chance of detecting and quantitating relevant conformational changes occurring within small protein domains of fewer than 30 amino acids (this results in an overall simulation timeof one microsecond).

### **2.2 Constructing a personal high-performance MD simulation environment**

Although we have used national grid computer resources in the past to carry out molecular dynamics simulations of protein domains (Weinzierl, 2010b; Heindl et al., 2011), we switched for our most recent work to high-performance personal computer systems based on workstations. The availability of multi-core CPUs, powerful graphics cards and affordable high-capacity mass storage media has over the last few years enabled the construction of workstation computing systems that are capable of providing the processing power for executing molecular dynamics simulations on a useful time scale. The custom-built workstation is based on liquid cooled dual Xeon 5690 Westmere CPUs clocked at 3.47 GHz with 24 GByte of DDR-3 RAM and running 64-bit Windows 7 Ultimate as operating system (more technical details are available on request). This 12 core (24 hyperthread) platform is capable of delivering 'raw' computing power consistently in excess of 120 GFlops in Linpack output. The MD simulation software of choice, 'YASARA Structure' (www.yasara.org), enables real-time simulations with standard (e.g. AMBER-03; Duan et al., 2003), or various custom force fields (NOVA, YAMBER, YASARA; Krieger et al., 2002; 2004; 2009). All structures were energyminimized in pre-equilibrated simulation boxes filled with TP3 water, and sodium and chloride ions added to a final concentration of approximately 150 mM. The simulation runs were executed using a 2 femtosecond time step in a fully solvated atomistic production mode without restraints employing an AMBER03 force field (Duan et al., 2003) at 300 K (27°C) at pH 7.0. The solvent density was kept constant at 0.997 g/l. Longrange Coulomb interactions were calculated using the Particle Mesh Ewald algorithm with periodic boundary conditions (cut-off 7.86Å; Essman et al., 1995).

In practical terms, the MD workstation allows the atomistic simulation of complex macromolecular assemblies (e.g. complete RNAPs in fully hydrated simulation cells [~300, 000 atoms] to be carried out at a rate of 1 nanosecond per 28 hours computation (average CPU usage ~75%), which compares favourably with resources typically allocated to average users from supercomputing facilities.

The resulting MD trajectories were analyzed with a series of analysis packages to study the spatial and temporal stability of secondary structures, distance relationships, and correlated/concerted conformational movements using YANACONDA macro scripts (http://www.yasara.org/yanaconda.htm). The structures shown in this text were generated using PyMOL (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC).

### **2.3 High-throughput molecular dynamics simulations of the RNAP link domain**

On high-performance computational platforms the simulation of small domains is relatively inefficient because the multi-core CPU processing power is not used to the fullest extent under such circumstances. On the other hand, running multiple simulations at the same time in parallel applications allows more effective use of processing power, but the reading in and out of data for the different simulations wastes time that otherwise would be available for simulations. We have come up with a simple solution to this problem which allows simulations of multiple domains to run in parallel, but within a single application. Construction of a simulation cell that contains multiple copies of the target structure as a regular threedimensional array allows an increased number of small domains to be simulated as a large structure (Fig. 4). The linear dimensions of simulation cells in any particular direction are often limited by the simulation software used, so we typically use three-dimensional clusters that contain the simulated domains in a crystal-like array or cube structure.

The domains within such an array can either be identical copies, thus allowing the systematic sampling of the conformational properties of a particular domain within a single

enabled the construction of workstation computing systems that are capable of providing the processing power for executing molecular dynamics simulations on a useful time scale. The custom-built workstation is based on liquid cooled dual Xeon 5690 Westmere CPUs clocked at 3.47 GHz with 24 GByte of DDR-3 RAM and running 64-bit Windows 7 Ultimate as operating system (more technical details are available on request). This 12 core (24 hyperthread) platform is capable of delivering 'raw' computing power consistently in excess of 120 GFlops in Linpack output. The MD simulation software of choice, 'YASARA Structure' (www.yasara.org), enables real-time simulations with standard (e.g. AMBER-03; Duan et al., 2003), or various custom force fields (NOVA, YAMBER, YASARA; Krieger et al., 2002; 2004; 2009). All structures were energyminimized in pre-equilibrated simulation boxes filled with TP3 water, and sodium and chloride ions added to a final concentration of approximately 150 mM. The simulation runs were executed using a 2 femtosecond time step in a fully solvated atomistic production mode without restraints employing an AMBER03 force field (Duan et al., 2003) at 300 K (27°C) at pH 7.0. The solvent density was kept constant at 0.997 g/l. Longrange Coulomb interactions were calculated using the Particle Mesh Ewald algorithm

In practical terms, the MD workstation allows the atomistic simulation of complex macromolecular assemblies (e.g. complete RNAPs in fully hydrated simulation cells [~300, 000 atoms] to be carried out at a rate of 1 nanosecond per 28 hours computation (average CPU usage ~75%), which compares favourably with resources typically allocated to average

The resulting MD trajectories were analyzed with a series of analysis packages to study the spatial and temporal stability of secondary structures, distance relationships, and correlated/concerted conformational movements using YANACONDA macro scripts (http://www.yasara.org/yanaconda.htm). The structures shown in this text were generated using PyMOL (The PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger,

**2.3 High-throughput molecular dynamics simulations of the RNAP link domain** 

contain the simulated domains in a crystal-like array or cube structure.

On high-performance computational platforms the simulation of small domains is relatively inefficient because the multi-core CPU processing power is not used to the fullest extent under such circumstances. On the other hand, running multiple simulations at the same time in parallel applications allows more effective use of processing power, but the reading in and out of data for the different simulations wastes time that otherwise would be available for simulations. We have come up with a simple solution to this problem which allows simulations of multiple domains to run in parallel, but within a single application. Construction of a simulation cell that contains multiple copies of the target structure as a regular threedimensional array allows an increased number of small domains to be simulated as a large structure (Fig. 4). The linear dimensions of simulation cells in any particular direction are often limited by the simulation software used, so we typically use three-dimensional clusters that

The domains within such an array can either be identical copies, thus allowing the systematic sampling of the conformational properties of a particular domain within a single

with periodic boundary conditions (cut-off 7.86Å; Essman et al., 1995).

users from supercomputing facilities.

LLC).

run, or could represent different (e.g. containing various mutations) versions. The presence of multiple copies within the same simulation cell automatically guarantees comparability because the same simulation parameters (e.g. ionic concentration, temperature, pressure and force field) will be applied to all domains in a consistent manner. A particular concern applying to such an arrangement is that the individual domains must be spaced sufficiently far apart to avoid artefactual contacts between the different molecules. Since the domains move randomly during MD simulation (diffusion), we automatically reposition the domains to their original position at fixed intervals (typically after every nanosecond simulation time), thus avoiding any 'clashes' within the simulation cell.

Fig. 4. A three-dimensional array of copies of the Link domain arranged for highthroughput MD simulation.

We will illustrate how this procedure works in practice using the Link domain described above. For the experiments described below we ran three independent repeats of arrays containing 18 distinct copies of the Link domain (as illustrated in Fig. 4) from the yeast RNAPII structure (PDB 2E2H; Wang et al., 2006) for 20 nanoseconds simulation time each, thus resulting in data representing just over one microsecond of total simulation data. Snapshots of the simulation cell were stored at ten picosecond intervals and stored for further processing and analysis.

### **2.3.1 The link domain contains a surprisingly stable -helical core**

When subjecting individual domains to molecular dynamics simulations, a common concern is that a domain may not be intrinsically stable, resulting in rapid unfolding and loss of defined structure. We expected this to happen in the case of the isolated Link domain, which consists only of a relatively short -helical region and irregularly structured N- and C-termini (Fig. 5).

Fig. 5. Detailed structure (A) and evolutionary conservation (B) of the Link domain. *H. sapiens* (human) and *S. cerevisae* (yeast) represent eukaryotes, *M. jannaschii* and *S. solfataricus* are archaea and *E.coli*, *T. aquaticus* and *T. thermophilus* signify bacteria. Residues identical to the *M. jannaschii* sequence are shown in red and the position of the invariant arginine is highlighted with a pale green bar.

Surprisingly, this small domain turned out to be exceptionally stable, with an essentially invariant -helical 'core' (amino acid residues spanning serine 764 to alanine 772) remaining in constant -helical conformation throughout most, or even all (asparagine 767 to glutamine 770) of the 20 nanosecond simulation periods (Fig. 6). In contrast, the irregularly structured regions, encompassing proline 757 to glutamine 763 undergoes a variety of conformations, including (random) coil, short-lived and 310 helix and turns (no -sheets were observed at any stage).

From this initial analysis we can therefore conclude that the Link domain consist of a nanomechanically robust portion and a highly flexible part (Fig. 7). A strategically placed proline reside (proline 765) delimits the N-terminal boundary of the -helix (Cochran et al., 2001; Rey et al., 2010; Tendulkar & Wangikar, 2011) and thus determines the position of the sharp kink present in the L-shaped Link domain (Figs. 5 & 7).

It is remarkable that, despite the apparently random movement of the N-terminal portion of the Link domain (proline 757 to glutamine 763), the overall 'L-shape' of the domain is nevertheless maintained in more than 90% of all the simulations carried out so far (Fig. 7). This observation suggests the presence of interacting side-chains that stabilize the intramolecular positioning of the mobile N-terminus relative to the rigid C-terminal -helix. A detailed investigation of the trajectories highlights such interactions involving histidine 761 ('H761'), glutamine 763 ('Q763'), arginine 766 ('R766') and - to a lesser extent - glutamine 770 ('Q770'; Fig. 8).

Fig. 5. Detailed structure (A) and evolutionary conservation (B) of the Link domain. *H. sapiens* (human) and *S. cerevisae* (yeast) represent eukaryotes, *M. jannaschii* and *S. solfataricus* are archaea and *E.coli*, *T. aquaticus* and *T. thermophilus* signify bacteria. Residues identical to the *M. jannaschii* sequence are shown in red and the position of the invariant arginine is

**R766** 

Surprisingly, this small domain turned out to be exceptionally stable, with an essentially invariant -helical 'core' (amino acid residues spanning serine 764 to alanine 772) remaining in constant -helical conformation throughout most, or even all (asparagine 767 to glutamine 770) of the 20 nanosecond simulation periods (Fig. 6). In contrast, the irregularly structured regions, encompassing proline 757 to glutamine 763 undergoes a variety of conformations, including (random) coil, short-lived and 310 helix and turns (no -sheets

From this initial analysis we can therefore conclude that the Link domain consist of a nanomechanically robust portion and a highly flexible part (Fig. 7). A strategically placed proline reside (proline 765) delimits the N-terminal boundary of the -helix (Cochran et al., 2001; Rey et al., 2010; Tendulkar & Wangikar, 2011) and thus determines the position of the

It is remarkable that, despite the apparently random movement of the N-terminal portion of the Link domain (proline 757 to glutamine 763), the overall 'L-shape' of the domain is nevertheless maintained in more than 90% of all the simulations carried out so far (Fig. 7). This observation suggests the presence of interacting side-chains that stabilize the intramolecular positioning of the mobile N-terminus relative to the rigid C-terminal -helix. A detailed investigation of the trajectories highlights such interactions involving histidine 761 ('H761'), glutamine 763 ('Q763'), arginine 766 ('R766') and - to a lesser extent - glutamine

sharp kink present in the L-shaped Link domain (Figs. 5 & 7).

highlighted with a pale green bar.

were observed at any stage).

770 ('Q770'; Fig. 8).

Fig. 6. Secondary structure analysis of the Link domain conformation based on 50 independent 20 nanosecond simulation runs. The proportion of four separate secondary structures adopted is shown for each residue (P757 to Q776).

Fig. 7. Superimposition of trajectory traces of a representative 20 nanosecond MD simulation of the Link domain. The -helical part (serine 764 to alanine 772) is structurally exceptionally stable, whereas the N-terminal portion undergoes significant conformational change.

At the beginning of the simulation (Fig. 8, t=0), the R766 side-chain engages in extensive van der Waals and/or hydrophobic interactions with Q763, which - due to the -helical repeat, occupies a similar orientation. The side-chains appear tightly lined up relative to each other and even move synchronously (Fig. 8, compare t=0 and t=30 ps). R766 can, however, make alternative contacts with Q770, whose side-chain (similar to Q763) occupies the same side of the -helix (Fig. 8, t=120 ps). While H761 is not initially involved in these interactions of R766 with Q763 or Q770, there is a clear tendency at a later stage of the simulation for it to

Fig. 8. Snapshots of the Link domain during 20 nanoseconds of MD simulation. See text for further information.

approach R766 (Fig. 8, t=350 ps). A conformational change moving Q763 out of the way appears necessary for the kink between the N-terminal flexible part and the -helix to become more distinct. This is followed by a tightening of the contacts between H761 and R766 (Fig. 8, t=590 ps). Strong interactions involving side-chains of histidine and arginine have been documented previously (Heyda et al., 2010) and it therefore appears that this observation is based on similar principles. In the final stages we observe additional stabilization of the H761-R766 interaction through the adjacent glutamine side-chains that were already previously identified as R 766 interaction partners (Fig. 8, t=1310 ps). This association, however, appears to be reversible within the nanosecond time frame of the MD simulation since we observe a distinct break-up of the side-chain cluster towards the end of the run (Fig. 8, t=1910 ps and t= 19970 ps).

## **3. Conclusion**

428 Molecular Dynamics – Studies of Synthetic and Biological Macromolecules

**H761**

**Q763**

**R766**

**t= 590 ps**

**t= 1310 ps**

**t= 1910 ps**

**t= 19970 ps**

Fig. 8. Snapshots of the Link domain during 20 nanoseconds of MD simulation. See text for

approach R766 (Fig. 8, t=350 ps). A conformational change moving Q763 out of the way appears necessary for the kink between the N-terminal flexible part and the -helix to become more distinct. This is followed by a tightening of the contacts between H761 and R766 (Fig. 8, t=590 ps). Strong interactions involving side-chains of histidine and arginine have been documented previously (Heyda et al., 2010) and it therefore appears that this observation is based on similar principles. In the final stages we observe additional

**t= 0 (simulation start)**

**Q770**

**t= 30 ps**

**t= 120 ps**

**t= 350 ps**

further information.

The gap between reality (conformational changes occurring in complex molecular machines in the millisecond to second range) and currently existing methodology (atomistic simulation of macromolecular structures in the nano- to microsecond range) is large and prevents us from obtaining a comprehensive understanding of the detailed functions of molecular machines that play fundamental role in key cellular processes. While it is very likely that this gap will narrow with the further development of computing power, viable intermediate solutions have to be found that can bridge the gap in the meantime. The approach outline here is aimed at breaking up complex macromolecular structures into components that provide distinctly identifiable nanomechanical functions which can be extensively studied using MD simulations. By gaining an understanding of the conformational spectrum of the individual components we can obtain two categories of complementary insights: the structural changes that a component is naturally capable of and the conformations that are either only achievable through energetic input or are unlikely to occur at all (this concept is already familiar to many molecular biologists in simpler systems, such as the conformational space allowed to a amino acids as represented by Ramachandran plots).

The work presented here has provided many new insights into the nanomechanical property of a newly identified structure of RNAP, the Link domain. As briefly outlined in the working model presented in Fig. 3, the Link domain appears to be strategically placed to act as an rNTP sensor, capable of transmitting the occupancy state of the catalytic site to the Bridge Helix. The Bridge Helix itself is a conformationally dynamic domain and almost certainly plays a critical role in translocating the nucleic acid substrate through the active site after successful rNTP incorporation (Weinzierl 2010a; 2010b). It thus appears essential for the translocation mechanism (Bridge Helix) to be allosterically linked to a sensor capable of reporting the successful completion of the incorporation of a nucleotide building block.

Although the initial identification of the Link domain was made entirely for structural reasons (i.e. its proximity to the N-terminal hinge region of the Bridge Helix), many of the observations emerging from the extensive molecular dynamic characterizations described here strengthen various aspects of the working model. The MD results show very clearly that the C-terminal -helical portion of the Link domain is very stable and conformationally essentially invariant. A proline residue, P765, separates this rigid portion of the domain from a much more flexible part, which also contacts the N-terminal portion of the Bridge Helix (probably with a key involvement of a portion of the F-Loop (Fig. 1; Miropolskaya et al., 2009; 2010). An arginine residue, R766, appears to play multiple critical roles. Apart from contacting the -phosphate of the rNTP in the catalytic site (Fig. 2), our MD results suggest that R766 may have a decisive influence on the conformation of the Link domain. Strong contacts (Heyda et al., 2010) between H761 and R766 determine the kinking of the L-shaped Link domain (Fig. 8), suggesting that the R766 side-chain makes exclusive contacts with the rNTP or H761, thus providing a plausible model for the transmission of the sensor information to the Bridge Helix (Fig. 3).

The exploration and definition of conformational space accessible to an individual protein domain provides new insights and better-informed working hypotheses for future experimental work. We will shortly initiate a high-throughput mutagenesis screen aimed at replacing each of the residues of the Link domain in the structurally highly orthologous archaeal RNAP (Fig. 5) in order to test the preliminary insights obtained from the computer modelling.

In our MD simulation work we have up to now focused on characterizing the nanomechanical properties of the wildtype domain. The approach illustrated here can also be expanded to include *in silico* studies of mutant variants to aid the interpretation of experimental results. MD simulations of domain assemblies of increasing size and complexity (e.g. Link domain, F-Loop, Bridge Helix and nucleic acid substrates) will also make additional contributions towards achieving a better understanding of inter-domain communication through allosteric signalling mechanisms in the foreseeable future.

### **4. Acknowledgment**

The research presented here was funded by MRC project grant G1100057 to ROJW.

### **5. References**


rNTP or H761, thus providing a plausible model for the transmission of the sensor

The exploration and definition of conformational space accessible to an individual protein domain provides new insights and better-informed working hypotheses for future experimental work. We will shortly initiate a high-throughput mutagenesis screen aimed at replacing each of the residues of the Link domain in the structurally highly orthologous archaeal RNAP (Fig. 5) in order to test the preliminary insights obtained from the computer

In our MD simulation work we have up to now focused on characterizing the nanomechanical properties of the wildtype domain. The approach illustrated here can also be expanded to include *in silico* studies of mutant variants to aid the interpretation of experimental results. MD simulations of domain assemblies of increasing size and complexity (e.g. Link domain, F-Loop, Bridge Helix and nucleic acid substrates) will also make additional contributions towards achieving a better understanding of inter-domain

communication through allosteric signalling mechanisms in the foreseeable future.

The research presented here was funded by MRC project grant G1100057 to ROJW.

120, No. 2, pp. 183-193, ISSN 0092-8674

16, pp. 1999–2012, ISSN 1096-987X

19, pp. 8577-8593, ISSN 0003-6951

Bar-Nahum, G., Epshtein, V., Ruckenstein, A.E., Rafikov, R., Mustaev, A. & Nudler, E.

Brueckner, F., Ortiz, J. & Cramer, P. (2009). A movie of the RNA polymerase nucleotide

Bustamante, C., Cheng, W. & Mejia, Y.X. (2011). Revisiting the central dogma one molecule

Cochran, D.A., Penel, S. & Doig A.J. (2001). Effect of the N1 residue on the stability of the

Cramer, P., Bushnell, D.A. & Kornberg, R.D. (2001). Structural basis of transcription: RNA

Duan, Y., Wu, C., Chowdhury, S., Lee, M.C., Xiong, G., Zhang, W., Yang, R., Cieplak, P.,

Essmann, U., Perera, L., Berkowitz, M.L., Darden, T., Lee, H. & Pedersen, L.G. (1995). A

at a time. *Cell*, Vol. 144, No. 4, pp. 480-497, ISSN 0092-8674

(2005). A ratchet mechanism of transcription elongation and its control. *Cell*, Vol.

addition cycle. *Current Opinions in Structural Biology*, Vol. 19, No. 3, pp. 294-299,

alpha-helix for all 20 amino acids. *Protein Science*, Vol. 10, No. 3, pp. 463-70, ISSN

polymerase II at 2.8 angstrom resolution. *Science*, Vol. 292, No. 5523, pp. 1863-1876,

Luo, R., Lee, T.; Caldwell, J., Wang, J. & Kollman, P. (2003). A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. *Journal of Computational Chemistry*, Vol. 24, No.

smooth particle mesh Ewald method. *The Journal of Chemical Physics*, Vol. 103, No.

information to the Bridge Helix (Fig. 3).

modelling.

**4. Acknowledgment** 

ISSN 0959-440X

1469-896X

ISSN 0036-8075

**5. References** 


Rey, J., Deville, J. & Chabbert, M. (2010). Structural determinants stabilizing helical

Seibold, S.A., Singh, B.N., Zhang, C., Kireeva, M., Domecq, C., Bouchard, A., Nazione, A.M.,

Schlick, T., Collepardo-Guevara, R., Halvorsen, L.A., Jung, S. & Xiao X. (2011). Biomolecular

Tan, L., Wiesler, S., Trzaska, D., Carney, H.C. & Weinzierl, R.O. (2008). Bridge helix and

Tendulkar, A.V. & Wangikar, P.P. (2011). Characterization and sequence prediction of

Vassylyev, D.G., Sekine, S.; Laptenko, O.; Lee, J.; Vassylyeva, M.N., Borukhov, S. &

at 2.6 Å resolution. *Nature*, Vol. 417, No. 6890, pp. 712-719, ISSN 0028-0836 Wang, D., Bushnell, D.A., Westover, K.D., Kaplan, C.D. & Kornberg, R.D. (2006). Structural

Weinzierl, R.O. (2010a). Nanomechanical constraints acting on the catalytic site of cellular

Weinzierl, R.O. (2010b). The nucleotide addition cycle of RNA polymerase is controlled by

Werner, F. (2008). Structural evolution of multisubunit RNA polymerases. *Trends in* 

Zhang, G., Campbell, E.A., Minakhin, L., Richter, C., Severinov, K. & Darst, S.A. (1999).

276, ISSN: 1047-8477

575-587, ISSN 0006-3002

2105

(printed) 0300-5127

1741-7007

44, No. 2, pp. 191-228, ISSN 0033-5835

*Biology*, Vol. 7, No. 10:40, 1741-7007, ISSN 1741-7007

*Cell*, Vol. 127, No. 5, pp. 941-54, ISSN 0092-8674

*Microbiology*, Vol. 16, No. 6, pp. 247-250, ISSN 0966-842X

*Cell*, Vol. 98, No. 6, pp. 811-824, ISSN 0092-8674

distortions related to proline. *Journal of Structural Biology*, Vol. 171, No. 3, pp. 266-

Feig, M., Cukier, R.I., Coulombe, B., Kashlev, M., Hampsey, M. & Burton, Z.F. (2010). Conformational coupling, bridge helix dynamics and active site dehydration in catalysis by RNA polymerase. *Biochimica et Biophysica Acta*, Vol. 1799, No. 8, pp.

modeling and simulation: a field coming of age. *Quarterly Reviews in Biophysics*, Vol.

trigger loop perturbations generate superactive RNA polymerases. *Journal of* 

structural variations in α-helix. *BMC Bioinformatics*, Vol. 12, Suppl 1:S20, ISSN 1471-

Yokoyama, S. (2002). Crystal structure of a bacterial RNA polymerase holoenzyme

basis of transcription: role of the trigger loop in substrate specificity and catalysis.

RNA polymerases. *Biochemical Society Transactions*, Vol. 38, No. 2, pp. 428-432, ISSN

two molecular hinges in the Bridge Helix domain. *BMC Biology*, Vol. 8:134, ISSN

Crystal structure of Thermus aquaticus core RNA polymerase at 3.3 Å resolution.

## *Edited by Lichang Wang*

Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This second book begins with an introduction of molecular dynamics simulations to macromolecules and then illustrates the computer experiments using molecular dynamics simulations in the studies of synthetic and biological macromolecules, plasmas, and nanomachines.

Photo by agsandrew / iStock

Molecular Dynamics - Studies of Synthetic and Biological Macromolecules

Molecular Dynamics

Studies of Synthetic and Biological

Macromolecules

*Edited by Lichang Wang*