# A Carbon Nanoswitch¶

Difficulty: High

## Introduction¶

In recent times the idea of creating machines that work on the nanoscale has captured the imagination of many scientists. The simplest machine imaginable is a nanoswitch; something that is extremely small and can switch an electric current on or off on the timescale of a few picoseconds, see refs. [Zheng2002], [Rivera2003] and [Rivera2005]. In this exercise we will consider such a device, which will be made from carbon nanotubes.

The idea is simple. We shall create a model switch from two nanotubes of carbon, one of which fits inside the other like a telescope drawtube. Our objective will be to determine some important properties of the switch, such as how quickly it can operate, what variables control its function and the characteristics of its dynamical behaviour. We shall also attempt to describe the system with a simple mathematical model and see how well this stands up to a full atomistic molecular dynamics simulation.

Let us start with the mathematical model.

## Theory¶

Consider two smooth tubes \(T_1\) and \(T_2\), of length L and diameters \(D_1\) and \(D_2\) as shown in Fig Two telescopic nanotubes , with \(D_1 > D_2\). The diameters are such that \(T_2\) can slide into \(T_1\) telescopically.

The displacement \(Δx\) of one tube with respect to the other is

Where \(x_1\) and \(x_2\) define the locations of the centre of mass of each tube. Note we have defined \(Δx\) so that it is always positive. It is reasonable to assume that the interaction energy U between the two tubes is proportional to the length of the overlap between them. Thus

where F is a constant, which has the units of force. Notice that when \(Δx=0\), \(U(Δx)=-FL\), which corresponds to the minimum energy. Also when \(Δx=L\), \(U(Δx)=0\), as you would expect, since the tubes are then fully withdrawn from each other. A plot of the system energy is shown in a), which is plotted as a function of \(x_2-x_1\) to show that the inner tube may extend in the negative direction as well as the positive. The forces acting on the tubes are obtained from the (negative) derivative of the energy, from which we have:

which is plotted in , b). So, when the inner tube \(T_2\) is displaced in the positive direction, it experiences a constant force -F pulling in the negative direction, and when it is displaced in the negative direction it experiences a constant +F force pulling it in the positive direction. The opposite is true for the outer tube \(T_1\). This pattern of forces implies an oscillatory system.

## Mathematical Exercise¶

Assuming the system starts from a *stationary* condition where
\(Δx=X_0\), and Newton’s Laws of Motion, derive formulae expressing
the position of each tube as a function of time, neglecting any possible
friction. You may assume that the tubes have masses \(M_1\) and
\(M_2\) and that when \(Δx=0\), then \(x_1=x_2=0\) (i.e.
this corresponds to the system centre of mass). Then use the model to
answer the following queries:

Derive a formula for the period of oscillation of the system in terms of the parameters of the model.

The forces do not represent simple harmonic motion in this case, in what ways is it different?

How can we use the formula derived for the period of oscillation to determine the important parameter F?

Try to work these out before seeing the solution, see subsection [sec:solution]. When you have answered these queries, you will be ready to tackle the molecular dynamics simulation.

## Instructions¶

Copy the Exercise8.tar.xz file into the $DL_POLY/data directory and unpack:

```
cd $DL_POLY/data
wget https://ccp5.gitlab.io/dlpoly/Exercise8.tar.xz
tar -xvf Exercise8.tar.xz
```

Now go to the $DL_POLY/execute directory and start up the GUI:

```
cd $DL_POLY/execute
```

## Task¶

Copy the contents of the subdirectory **Exercise8** into the **execute**
subdirectory. To do this using the GUI, select from the main menu
**Execute > Store/Fetch Data** and in the Data Archive window type
Exercise8 in the Fetch box and click Fetch. You may notice that the only
file copied across is a **CONTROL** file. Unlike other exercises in this
series the **Exercise8** sub-directory does not contain all the required
DL_POLY input files. The **CONTROL** file is however a valid one, so
take care not to delete or damage it as we will need it later. The
**CONFIG** and **FIELD** files will be generated using the GUI in this
exercise.

The **Exercise8** directory contains other useful material we shall need
later, but these will have to be copied across using linix commands,
since the GUI cannot see them. Note also that **Exercise8** is a good
place to save the files you create in this exercise.

To perform the simulation of this system using DL_POLY, we must first
create a **CONFIG** file containing the two carbon nanotubes - with one
tube inside the other. Next we must create a **FIELD** file to describe
the interactions between the atoms in the tubes.

## Creating the **CONFIG** File¶

This will be a somewhat lengthy procedure using the DL_POLY GUI, but it is much easier with the GUI than without!

Now go to the

**$DL_POLY/execute**directory and start up the GUIWhen the GUI appears, move the mouse cursor to the FileMaker menu and, with the mouse button held down, find the CONFIG/Bucky option. Select this option and the Make Fullerene panel will appear.

Without changing any of the options appearing on the panel, click on the Tube button. This will do two things: draw a picture of a carbon nanotube on the screen; and create a file called CFGBUK.0 in the execute directory. This fact will be announced in the Monitor window of the GUI. You have now created the first carbon tube.

To create the second nanotube, you must change the default panel parameter for the tube circumference. Click on the text box on the left of the panel below the specification

*Tube size: rings X,Y*and change the number from 8 to an integer between 16 and 21. Click the Tube button once again and the new nanotube will appear on screen (and a file CFGBUK.1 will be created in the execute directory). Now we must merge the two tubes into one system, and set up a displacement, as in our mathematical model.With the second carbon tube still on screen click the Edt button on the right hand side of the GUI. This will activate the GUI molecular structure editor.

From the newly appeared Editor menu on top of the GUI, select the Frg/Options Search. In the file browser that appears, enter the name CFGBUK.0 and click the Open button.

Now click the Frg button on the bottom left of the GUI to activate the

*fragment insertion*option. A single click on the graphics window will insert the first carbon tube inside the second.Notice that the new tube is highlighted in red, and the word

*SHIFT*has appeared in the top left of the graphics window. This means the inner tube is group selected and you can displace it by dragging the mouse cursor vertically up the screen. You should try to move the tube so that it stays central in the outer tube, otherwise there may be problems later! We suggest you use the vertical line of the MD cell as a guide. Displace the tube about half its length, but no more.Now we need to remove the MD cell surrounding the system. If this is done, DL_POLY will remove any net system rotation, which would complicate our simulation. From the Options menu on top of the GUI, select the Box None option. Now click on the Box button on the bottom right of the GUI and the MD cell will disappear.

Click the Sav button, to save your new structure in a numbered file CFGEDT.0. Now click the Edt button once again to close the editor and prevent you accidentally editing your structure any more! (Note: if you close the editor without saving, the GUI will save the structure for you in an unnumbered file called CFGEDT. This is safe until you open the Editor again!)

## Creating the **FIELD** File¶

In principle constructing a **FIELD** file for such a system is
difficult, but fortunately the GUI has sufficient capability to
construct a sensible (if not perfect) one quickly. We shall use it.

Proceed with the structure you have on screen. If you have already shut down the GUI meantime, start it up again and load up the CFGEDT.0 file again using the New button.

From the FileMaker menu, select the FIELD/Dreiding option. The Dreiding Field Maker panel will appear. The Dreiding field is a simple, universal, all-atom force field suitable for organic molecules. You do not need to change any of the options on the panel, just click the Make button.

After less than a minute the GUI will create a

**FIELD**file named FLDDRE.0. While this is going on, take a look at the Monitor window, where the GUI reports the information it has deduced about the structure.When the FLDDRE.0 file is written, we are ready to proceed to the next stage.

If the construction of the FLDDRE.0 file fails, then it is likely that you have brought the two nanotubes too close together during the setup and confused the GUI’s bond analysis which is based on atom proximity. In which case start again and try harder to keep the two tubes apart.

## Running the Simulation¶

The simulation can be set running in two ways:

From the

**$DL_POLY/execute**directory type the command**DLPOLY.Z**and the program will start running. (But remember firstly to rename the two files FLDDRE.0 and CFGEDT.0 to**FIELD**and**CONFIG**respectively.) Control will be handed back to you when the program has stopped (or crashed!)From the GUI, select the Execute/Run DL_POLY option from the menu. This opens the Run DL_POLY panel. Use the

**FIELD**and**CONFIG**buttons to select the files you want (i.e. those you have prepared above). You will not need to use the**CONTROL**button, since this file is already correctly named. Click the Run button, see , and the job will be start in the background. If you click the Status button at frequient intervals the GUI will eventually tell you that the job has finished. It should take less than 200 seconds.

## Inspecting the Results¶

You can use the GUI to do some simple checking of the results.

Using the Analysis/Display/REVCON menu option you can look at the final structure, to see if it has evolved to a sensible structure.

Using the Analysis/Tools/Statistics menu option to plot the system’s energy, temperature etc as a function of time. Of particular interest is the configuration energy (does this make sense in terms of our model)?

However there is a much more interesting thing you can do using a
program called VMD, which is available from here. Download and install
the program on your machine as instructed. But before you use it you
must firstly convert the DL_POLY **HISTORY** file into an XYZ file that
VMD can read. There is a small program in the
**$DL_POLY/data/Exercise8/files** directory to make this conversion
called *his2xyz.F90*. Copy this into the execute directory and compile
it using gfortran.

```
gfortran -O3 -o his2xyz his2xyz.F90
```

Run the program using the command: **his2xyz** (or on some linux
systems: **./his2xyz**), which will automatically pick up the
**HISTORY** file and make a file called HIS.xyz.

**Note:** You can run his2xyz, and other command line programs (see
below), from within the GUI. You simply need to open the Run DL_POLY
panel from the Execute menu and type the command in the text box on the
panel. Clicking the Run button will set the program running.

When you start VMD you need to pick up the HIS.xyz file (you must direct VMD to seek xyz files). VMD will play this like a film. You should see clearly what the molecular dynamics have produced. Experiment with the graphics of VMD, to see how different renderings may be produced.

## Analysing the Results¶

Having amused yourself looking at movies, it’s time to get quantitative.
We need to determine the period of oscillation of the system so we can
determine the force F in our model. It is hard to do this from the
movie! To get a quantitative answer, copy the *comhis.f* program from
the **$DL_POLY/data/Exercise8/Files** directory into the execute
directory and compile it with gfortran. Use this to calculate the
position of the centre of mass of the inner tube \(T_2\), using data
in the **HISTORY** file and plot it as a function of time. The output
file is called COM.XY.

```
gfortran -O3 -o comhis comhis.F90
```

The command to execute the program is simply: **comhis noOfAtoms**,
where *noOfAtoms* is the number of atoms in the outer tube, \(T_1\).

From such a plot you should be able to estimate the period of
oscillation with fair accuracy. It is also necessary to determine the
sweep of the oscillation – how far the tube moves between stationary
points on the plot. Both are needed to calculate the force F. Another
way to determine the period of oscillation is to Fourier transform the
plot of the centre of mass motion. This will project out of the data the
true frequency of oscillation, which is the reciprocal of the period.
Unfortunately, this is not as trivial as it sounds, so we have provided
a program, *fft.f* in the **$DL_POLY/data/Exercise8/Files** directory.

```
gfortran -O3 -o fft fft.f
```

Use this to obtain the Fourier transform of the file COM.XY and plot the output file FFT.XY. Use the command:

```
fft < COM.XY > FFT.XY
```

Determine the frequency of the oscillation and from it, the period. How does this compare with your previous determination?

Those of you who have a particular interest in using Fourier transforms
for frequency analysis should look at the program *fft.f* and see how
this works. Don’t worry about the FFT subroutine itself; this is
normally a ‘black box’ routine from a numerical library. Look instead
at:

The subtraction of the mean value of the input function before the Fourier transform is done;

The way the

*Dirichlet*boundary conditions are applied to the first and last points of the input data; andThe use made of a window function – in this case a Gaussian.

These *tricks* are often neglected in naive applications of FFTs, and
the result is a confusing, noisy function. Note also that *fft.f*
calculates the modulus of the (complex) Fourier transform, since we are
interested only in the amplitude, and not the phase of the oscillation.

Key question: Using the information gathered above what is your estimate of the force constant F in our mathematical model? If you have time, try to do at least one more simulation, using a different diameter for the outer tube \(T_1\).

## Further Investigations¶

You should already have seen how the system configuration energy compares with the mathematical model, but what about the force? The model says this should be constant either side of the origin of the coordinates, but switch sign as we pass through the origin. How can we verify these? We can use the Verlet algorithm!

The standard Verlet algorithm is:

which is easily rearranged to:

So we see that we may perform a numerical differentiation of the position of the tube with respect to time to obtain the force.

We have provided the program *acceln.f* in the
**$DL_POLY/data/Exercise8/Files** directory, which uses this approach to
calculate the acceleration of a tube, using the data in the file COM.XY.
Compile the program and run it using the commands:

```
gfortran -O3 -o acceln acceln.f
acceln < COM.XY > ACC.XY.
```

Now plot the file ACC.XY. How well does this result compare with the model? In what ways is it different and why?

We have made no mention of friction in this treatment. How could you go about determining the effect of this? Devise an experiment that you could do to test this. Speculate on the effect that temperature might have on the simulation. Would it increase or decrease friction?

Now you’ve done at least one simulation and seen at least one movie, you should think about how the experiment can be improved. Did anything unexpected happen? Would this affect the validity of our simple model? How would you amend the model to make it more consistent with the simulation?

## Some example results¶

A sample of results are presented in

\(T_2\) |
Frequency [\(ps^{-1}\)] |
Period [ps] |
\(X_0\) [Å] |
m [D] |
Force [\(DAps^{-2}\)] |
---|---|---|---|---|---|

16 |
0.0915 |
10.93 |
8.085 |
138.67 |
300.3 |

17 |
0.0831 |
12.04 |
9.390 |
141.44 |
293.2 |

18 |
0.0864 |
11.57 |
9.935 |
144.00 |
342.0 |

19 |
0.0729 |
13.72 |
9.300 |
146.37 |
231.4 |

20 |
0.0661 |
15.13 |
8.900 |
148.57 |
184.8 |

21 |
0.0746 |
13.41 |
9.505 |
150.62 |
254.8 |

Sample graph of center of mass movement in fig Centre of mass movement and a sample video can be seen on youtube at https://youtu.be/S6q2dr-p5DU.

## Solution¶

The potential energy is:

and the derivatives w.r.t. \(x_1\) and \(x_2\) (i.e. the forces) are:

So the equations of motion are:

Integrating w.r.t. time gives:

In which the integration constants are zero, for a stationary start. This however cannot be the whole story. We see clearly from , b) that there is a discontinuity in the force when \(x_2-x_1\) changes sign. This means that the integration over time must be performed piecewise, with each piece being evaluated between the discontinuities. For now however, we shall integrate from the time t=0 to the first instance when \(x_2-x_1\) changes sign. This will be sufficient for our needs. Integrating w.r.t. time a second time (subject to the condition \(x_2 - x_1 > 0\)) gives:

Where, following from the fact that the centre of mass of the system is at the origin, and using the initial displacement \(X_0\), we have:

Concentrating now on the inner tube \(T_2\), we see that when t=0, the position of the tube is given as \(x_2 = C/M_2\), which is positive. At the same time \(x_1\) is negative, so the sign of \(x_2 – x_1\) is positive. Thus as t increases, \(x_2\) must decrease. In fact the fomulae above show that the change in \(x_2\) is quadratic in time and resembles a free falling object under gravity. At the point \(x_2 = x_1\), the force switches sign but retains its magnitude, so the motion thereafter will be the exact opposite of what has happened until now (replace \(x_2\) with \(-x_2\), \(v_2\) with \(–v_2\), F with –F), it will be a deceleration that carries the tube \(T_2\) to the point \(–C/M_2\). The motion from there on retraces the whole sequence so far and an oscillation is therefore established. In the first part of the motion, the time taken for \(T_2\) to ‘fall’ from position \(x_2=C/M_2\) to \(x_2=0\) is given by:

This represent one quarter of a full oscillation, so the period of oscillation \(t_2\) of tube \(T_2\) is given by:

Note that the period of the oscillation depends on the magnitude of C i.e. on the initial displacement of the tubes. This is different from simple harmonic motion, in which the period is independent of the initial displacement. Another way in which the motion differs from SHM, is that it is piecewise quadratic in time, so it does not follow a simple trigonometric function [e.g. \(x_2 = A sin(t/t_2)\)]. To determine the force constant F, we can simply invert the above formula and obtain

We can obtain \(X_0\) conveniently by noting that the sweep of the oscillation (i.e. the distance travelled from one stationary point to another – from one peak of the plot to the next trough) is \(2X_0\).

- Zheng2002
Zheng, J. Z. Liu, and Q. Jiang, Excess van der Waals interaction energy of a multiwalled carbon nanotube with an extruded core and the induced core oscillation, Phys. Rev. B, 65, p. 245409, May 2002, doi: 10.1103/PhysRevB.65.245409

- Rivera2003
Rivera, C. McCabe, and P. T. Cummings, Oscillatory behavior of double-walled nanotubes under extension: a simple nanoscale damped spring, Nano Letters, 3(8), p. 1001, 2003, doi: 10.1021/nl034171o

- Rivera2005
Rivera, C. McCabe, and P. T. Cummings, The oscillatory damped behaviour of incommensurate double-walled carbon nanotubes, Nanotechnology, 16(2), p. 186, 2005, doi: 10.1088/0957-4484/16/2/003