**2. Molecular dynamics (MD) simulation**

Molecular dynamics (MD) simulations technique has been used in this research to perform the atomistic calculations, of frequency dependent electrical conductivity in Nax(Ti8-xCrx)O16, (x = 1.7). MD simulation is generally carried out to compute the motions of individual molecules in models of solids, liquids and gases. The key idea is motion, which describes how positions, velocities, and orientation change with time (Haile, 1992). The MD simulation method is carried out in such a way that atoms are represented by point particles and the classical (Newton) equations of motion, "force equals mass times acceleration or F = ma" are integrated numerically. The motions of these large numbers of atoms are governed by their mutual interatomic interaction. MD simulations are limited largely by the speed and storage constraints of available computers. Hence, simulations are usually done on system containing 100-1000 particles, with a time step of 1x10-15 s. This technique for simulating the motions of a system of particles when applied to biological macromolecules gives the fluctuations in the relative positions of the atoms in a protein or in DNA as a function of time. Knowledge of these motions provides insights into biological phenomena. MD is also being used to determine protein structure from NMR, to refine protein X-ray crystal structures faster from poorer starting models, and to calculate the free energy changes resulting from mutation in proteins (Karpus & Petsko, 1990).

The construction of the molecular dynamics model involved mainly model development and the use of the Molecular Dynamics simulation technique to solve the equation of motion iteratively. In model development, firstly the interaction between ions and the interaction between ion and environment have to be defined. The interaction between ions comprises three components: Lennard-Jones potential, Coulomb potential and Van der Waals' attraction. The rigid lattice approximation is being used in the simulation. This describes the interaction between the ions and the environment. In the rigid lattice approximation only the tunnel ions are allowed to displace from their equilibrium position. The environment ions still give an interaction with the tunnel ions, i.e. they produce a constant potential surface for the tunnel ions. The tunnel ions interact with one another and are allowed to displace. After this has been done, the equations of motion are developed. A program, Gretep (LMGP-Suite Suite of Programs for the interpretation of X-ray Experiments) is used to generate the positions of all the ions by keying in parameters such as space group, unit cell dimensions, atomic parameters and etc. Then titanium ions are replaced by chromium ions randomly placed in unit cells along the tunnel according to the relative proportions of titanium and chromium. A tunnel ion, which is the sodium ion, is added to the model structure for every chromium ion, in order to preserve charge neutrality.

#### **2.1 Modeling of ion-ion interactions**

The molecular interactions are based on the intermolecular potential energy function. The total potential energy between two ions is the sum of the Lennard-Jones potential and Coulomb potential (West, 1988),

$$V(r) = \frac{\mathcal{L}}{r^p} \pm \frac{e^2}{4\pi\varepsilon\_0 r} \tag{1}$$

Molecular Dynamics Simulation and Conductivity

position of the chromium ion.

sodium ions.

Mechanism in Fast Ionic Crystals Based on Hollandite NaxCrxTi8-xO16 375

where |Aion1 – Aion2| is the absolute value of the displacement in a-axis for ion1 and ion2 and r is the distance between ion1 and ion2. The expressions for the component of forces in the b and c-axes are equivalent in form. The tunnel wall is made up by alternate stacking of two different oxygen-square layers which consist of four oxygen atoms denoted O1 at z = 0 (square-plane) and four oxygen atoms denoted O2 at z = 0.5 (cavity). The number and the positions of the basic tunnel ions are obtained from the x-ray analysis (Michiue & Watanabe, 1995a). There will be six sodium ions in the seven cavities, in other words in the 14 oxygensquare layers. This is shown in figure 2. The position of the sodium ion will depend on the

Fig. 2. Schematic representation of a probable local arrangement for sodium ions in the tunnel of NaxCrxTi8-xO16. In the 7 x Cavity (indicated by the rectangular box), there are six

From figure 2, it is clearly shown that there are three possible equilibrium positions for the sodium ions, Na1 (0.72,0.14,0.5), Na2 (0,0,0.2) and Na3 (0,0,0) (Michiue & Watanabe, 1995a). The three positions are very close to each other and therefore it is impossible for Na ions to occupy these positions simultaneously. The positions for the sodium ion are dependent on the position of the chromium ion. When a titanium ion in the square-plane is substituted by a chromium ion, a sodium ion (either Na2 or Na3) will be placed in the hollandite model as shown in the figure 3. The site of the chromium ion is chosen randomly from the titanium ions in the four corners, A, B, C or D. The Na2 and Na3 differ slightly in their position along the caxis. The atomic parameters for Na2 and Na3 in c-axis or z-axis shown in figure 3 are 0.2 and 0 respectively. When a titanium ion in the cavity is substituted by a chromium ion, there will be four possible positions for the sodium ion (Na1). The chromium ion is chosen randomly from the titanium ions in the four corners, A, B, C or D. The sodium ion preferentially resides at an interstitial site within the same unit cell that contains a chromium ion (Michiue & Watanabe, 1996) shown by the arrow in figure 3. In this work, 24 sodium ions have been considered. Therefore, four sets of the six sodium ions shown in figure 2 are stacked up together to form a

where p is about 10 and the values for are different for different ion-pairs, e is the electron charge and r is the distance between two ions. The values for lambda for each ion-pair have to be found for the calculations of the equations of force. The method used is presented elsewhere (Khoo, 2003; Dissado and Khoo, 2006; Khoo et al., 2004, 2007). The relationship between the tunnel system and its surroundings is defined by using boundary conditions. These describe the interactions between the molecules with their surroundings. Rigid lattice approximation has been used in the simulation to simplify the simulation process and to reduce the simulation executable time. Rigid lattice approximation is done in such a way that only the tunnel ions, which are the sodium ions, are free to displace. All other ions in the lattice sites would remain static. However, there will still be interaction between sodium ions and the ions in the lattice sites.

Reflective boundary conditions are used at the two ends of the tunnel. Rebound is a special type of collision involving a direction change; the result of the direction change is a large velocity change. Collisions in which particles rebound with the same speed are known as elastic collisions. Thus, the velocity of the ion that bounds back from the boundary should have the opposite sign and equal in magnitude with that velocity of ion which hits the boundary. The angle of the velocity hitting the boundary is the same as the angle leaving the boundary. Force is the derivative of potential energy. Therefore, the component of the force along the a-axis is given by equation 2.

$$F\_a = -\frac{dV(r)}{da} = \frac{dr}{da} \frac{\lambda p}{r^{p+1}} \mp \frac{dr}{da} \frac{e^2}{4\pi\varepsilon\_0 r^2} \tag{2}$$

$$r = \sqrt{\left(\text{Aion1} - \text{Aion2}\right)^2 + \left(\text{Bion1} - \text{Bion2}\right)^2 + \left(\text{Cion1} - \text{Cion2}\right)^2} \tag{3}$$

$$\frac{dr}{da} = \frac{|Aion1 - Aion2|}{r} \tag{4}$$

ions still give an interaction with the tunnel ions, i.e. they produce a constant potential surface for the tunnel ions. The tunnel ions interact with one another and are allowed to displace. After this has been done, the equations of motion are developed. A program, Gretep (LMGP-Suite Suite of Programs for the interpretation of X-ray Experiments) is used to generate the positions of all the ions by keying in parameters such as space group, unit cell dimensions, atomic parameters and etc. Then titanium ions are replaced by chromium ions randomly placed in unit cells along the tunnel according to the relative proportions of titanium and chromium. A tunnel ion, which is the sodium ion, is added to the model

The molecular interactions are based on the intermolecular potential energy function. The total potential energy between two ions is the sum of the Lennard-Jones potential and

> ( ) <sup>4</sup> *<sup>p</sup> <sup>e</sup> V r*

where p is about 10 and the values for are different for different ion-pairs, e is the electron charge and r is the distance between two ions. The values for lambda for each ion-pair have to be found for the calculations of the equations of force. The method used is presented elsewhere (Khoo, 2003; Dissado and Khoo, 2006; Khoo et al., 2004, 2007). The relationship between the tunnel system and its surroundings is defined by using boundary conditions. These describe the interactions between the molecules with their surroundings. Rigid lattice approximation has been used in the simulation to simplify the simulation process and to reduce the simulation executable time. Rigid lattice approximation is done in such a way that only the tunnel ions, which are the sodium ions, are free to displace. All other ions in the lattice sites would remain static. However, there will still be interaction between sodium

Reflective boundary conditions are used at the two ends of the tunnel. Rebound is a special type of collision involving a direction change; the result of the direction change is a large velocity change. Collisions in which particles rebound with the same speed are known as elastic collisions. Thus, the velocity of the ion that bounds back from the boundary should have the opposite sign and equal in magnitude with that velocity of ion which hits the boundary. The angle of the velocity hitting the boundary is the same as the angle leaving the boundary. Force is the derivative of potential energy. Therefore, the component of the force

( )

4 *<sup>a</sup> <sup>p</sup> dV r dr dr e <sup>p</sup> <sup>F</sup>*

*dr Aion Aion* | 1 2|

*da r*

*da da da r r* 

2

(4)

(2)

1 2 0

<sup>222</sup> *r Aion Aion Bion Bion Cion Cion* ( 1 2) ( 1 2) ( 1 2) (3)

2

*r r*

0

(1)

structure for every chromium ion, in order to preserve charge neutrality.

**2.1 Modeling of ion-ion interactions** 

Coulomb potential (West, 1988),

ions and the ions in the lattice sites.

along the a-axis is given by equation 2.

where |Aion1 – Aion2| is the absolute value of the displacement in a-axis for ion1 and ion2 and r is the distance between ion1 and ion2. The expressions for the component of forces in the b and c-axes are equivalent in form. The tunnel wall is made up by alternate stacking of two different oxygen-square layers which consist of four oxygen atoms denoted O1 at z = 0 (square-plane) and four oxygen atoms denoted O2 at z = 0.5 (cavity). The number and the positions of the basic tunnel ions are obtained from the x-ray analysis (Michiue & Watanabe, 1995a). There will be six sodium ions in the seven cavities, in other words in the 14 oxygensquare layers. This is shown in figure 2. The position of the sodium ion will depend on the position of the chromium ion.

Fig. 2. Schematic representation of a probable local arrangement for sodium ions in the tunnel of NaxCrxTi8-xO16. In the 7 x Cavity (indicated by the rectangular box), there are six sodium ions.

From figure 2, it is clearly shown that there are three possible equilibrium positions for the sodium ions, Na1 (0.72,0.14,0.5), Na2 (0,0,0.2) and Na3 (0,0,0) (Michiue & Watanabe, 1995a). The three positions are very close to each other and therefore it is impossible for Na ions to occupy these positions simultaneously. The positions for the sodium ion are dependent on the position of the chromium ion. When a titanium ion in the square-plane is substituted by a chromium ion, a sodium ion (either Na2 or Na3) will be placed in the hollandite model as shown in the figure 3. The site of the chromium ion is chosen randomly from the titanium ions in the four corners, A, B, C or D. The Na2 and Na3 differ slightly in their position along the caxis. The atomic parameters for Na2 and Na3 in c-axis or z-axis shown in figure 3 are 0.2 and 0 respectively. When a titanium ion in the cavity is substituted by a chromium ion, there will be four possible positions for the sodium ion (Na1). The chromium ion is chosen randomly from the titanium ions in the four corners, A, B, C or D. The sodium ion preferentially resides at an interstitial site within the same unit cell that contains a chromium ion (Michiue & Watanabe, 1996) shown by the arrow in figure 3. In this work, 24 sodium ions have been considered. Therefore, four sets of the six sodium ions shown in figure 2 are stacked up together to form a

Molecular Dynamics Simulation and Conductivity

**3.1 Position of sodium ions in c-axis without applied field** 

**3. Results** 

it belongs to.

Mechanism in Fast Ionic Crystals Based on Hollandite NaxCrxTi8-xO16 377

After ten runs have been carried out with different initial velocities for the individual sodium ions with their average value being set by the defined temperature, the average positions for each of the sodium ions in the c-axis are then calculated. Figure 5a shows the average position for the first six sodium ions along the c-axis with the initial conditions of 273K, 100000 intervals and time step=10-15s. This gives a good comparison of the positions of the sodium ions throughout the 100000 intervals. The sodium ions only vibrate in their equilibrium positions depending upon where their equilibrium positions were. For example, the sixth sodium ion (pink) only vibrates around the cavity

Fig. 5. (a) The average position along the c-axis for the first six sodium ions (b) The arrangement for the first six sodium ions in the tunnel. The arrows show the alignment of

An electric field in the range of 7.43MV/m to 74.3GV/m was applied along the c-axis to the hollandite model at the 5001th time interval. The initial conditions for the results shown below were temperature=273K, time step=10-15s, 100,000 intervals and electric field=743MV/m. Figure 6 shows the positions of the first six sodium ions as a function of time. Over the initial 5000 intervals, which was in the absence of the electric field, the sodium ions just vibrate around their equilibrium positions. Starting from 5001th intervals, the positions of the sodium ions change dramatically. For example, the sixth sodium ion (pink) moved to the next cavity at some points and then back to the original cavity and then

the cavities in the tunnel and the graphs.

to the next cavity again.

**3.2 Position of sodium ions in c-axis with applied field** 

longer tunnel. C++ is used as a tool to carry out the molecular dynamics simulation, to perform the complicated calculations and to store the required results in the specified text files. A total of five C++ programs have been written and the flow chart is shown in figure 4. Details of the calculations can be found in reference (Khoo, 2003).

Fig. 3. Hollandite model projected slightly off the c-axis to give a clearer view of the 3 dimensional structure. When TiA4+ is replaced by Cr3+, the preferable location for the Na1A+ ions is shown by the arrow, a similar situation is found for TiB4+, TiA4+ & TiD4+ as shown by Na1B+, Na1C+ and Na1D+.

Fig. 4. Flow chart showing the five program codes written to perform the MD simulation, to carry out the complicated calculations and to generate the data required in text files.
