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

\[Δx=|x_2-x_1|\]

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

\[U(∆x)=F(∆x-L)\]

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:

\[\begin{split}\begin{split} F_1 = &+F\; and\; F_2 = -F\; when\; x_2-x_1 > 0 \\ F_1 = &-F\; and\; F_2 = +F\; when\; x_2-x_1 < 0 \end{split}\end{split}\]

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 GUI
When 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; and
The 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:

\[x_{n+1}=x_{n-1} + Δt^2 F_n/M\]

which is easily rearranged to:

\[F_n = M(x_{n+1} - 2x_n + x_{n-1}) / Δt^2\]

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

Sample results¶
$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:

\[U=F(∆x-L)\]

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

\[\begin{split}\begin{split} -\frac{dU}{dx_1} &= sgn(x_2 - x_1) F \\ -\frac{dU}{dx_2} &= - sgn(x_2 - x_1) F. \end{split}\end{split}\]

So the equations of motion are:

\[\begin{split}\begin{split} M_1\frac{dv_1}{dt} &= sgn(x_2 - x_1) F \\ M_2\frac{dv_2}{dt} &= - sgn(x_2 - x_1) F \end{split}\end{split}\]

Integrating w.r.t. time gives:

\[\begin{split}\begin{split} M_1v_1 &= sgn(x_2 - x_1) Ft \\ M_2v_2 &= - sgn(x_2 - x_1) Ft \end{split}\end{split}\]

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:

\[\begin{split}\begin{split} M_1x_1 &= sgn(x_2 - x_1) Ft^2/2 - C \\ M_2x_2 &= - sgn(x_2 - x_1) Ft^2/2 + C \end{split}\end{split}\]

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:

\[C = \frac{M_1M_2X_0}{M_1+M_2} \; or\; C = mX_0.\]

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:

\[0 = -Ft^2/2 + C \; i.e. \; t = \sqrt{\frac{2C}{F}}\]

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

\[t_2 = \sqrt{4\frac{2C}{F}}\]

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

\[F = \frac{32C}{t_2^2}\; or\; F = 32 m \frac{X_0}{t_2^2}\]

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]

1. 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]

1. 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