ation will be strained.
The actual straining is done in the “for-loop” beginning at line 38,
by increasing the coordinates of the atoms tagged “Right” by the amount
along the z-direction. On lines 27 and 28, the stress and total length of the
configuration are calculated and appended to two separate lists. These will be
used later for plotting.
pre_strain_hook=strainer in the
the class (technically, the method
__call__) is called before the force and energy
calculation in every MD step.
Computing Young’s modulus¶
To properly calculate the Young’s modulus, it is necessary to take the total volume of the simulation domain into account. The stress reported by QuantumATK therefore has to be multiplied by the total volume of the cell and then divided by the volume of the nanotube.
Several definitions of the volume of a CNT. A common
way of defining it is to say that the wall of the nanotube extends for
half a \(\pi\)-bond in and out of the wall, so that it becomes a
cylindrical shell with thickness equal to that of a \(\pi\)-bond,
approximately 1.35 Å, so that the total cross-sectional area is 1.35 Å
times the circumference of the nanotube. The circumference for a
(5,5)-armchair nanotube is
\(5 \times (3 \times x)\), where \(x\)=1.42 Å is the carbon-carbon
bond length. Finally, for the volume, you need the length of the tube,
which is given in the
.nc file saved by QuantumATK.
At the bottom of the script, delete the line
Add instead the following lines of code in order to invoke the proper stress calculation:
# ------------------------------------------------------------- # Finding Young's Modulus # ------------------------------------------------------------- # Find center of cnt from mean x- and y-coordinates coords = bulk_configuration.cartesianCoordinates() x_coords = coords[:,0] y_coords = coords[:,1] x_0 = numpy.sum(x_coords)/len(x_coords) y_0 = numpy.sum(y_coords)/len(y_coords) # Find radius by mean distance from center to atoms' x- and y-coordinates radius = numpy.sum( ( (x_coords - x_0)**2 + (y_coords - y_0)**2 )**(1/2.0) ) / len(x_coords) # Find cross-sectional area by multiplying circumference with length of a pi-bond pi_bond = 1.35 circumference = radius * 2 * numpy.pi area = circumference * pi_bond # Find volumes of CNT and of whole cell volume_cnt = area * strainer.lengths[-1] u1, u2, u3 = bulk_configuration.bravaisLattice().primitiveVectors() volume_cell = u1 * u2 * u3 # Factor the stresses to the volume of the CNT stresses_factored = numpy.array(strainer.stresses * (volume_cell / volume_cnt) )*eV/Angstrom**3 # Calculate the strain strain = (strainer.length() - strainer.length() ) / strainer.length() # Calculate Young's modulus in TPa print "Young's Modulus: ", (stresses_factored[-1].inUnitsOf(Pa) - stresses_factored.inUnitsOf(Pa) ) / float(strain[-1]) / 1e12, " TPa"
Finally, add a few lines of code that will save different parts of the calculation and plot the results using pylab.
# Save the result to a file. print strainer.length() nlsave("mdtrajectory_cnt.hdf5", md_trajectory, "MD") nlsave("mdtrajectory_cnt.hdf5",strainer.length(), "C length") nlsave("mdtrajectory_cnt.hdf5",strainer.stress(), "C stress") nlsave("mdtrajectory_cnt.hdf5",stresses_factored, "C stresses factored") # Plot stress-strain curve. import pylab pylab.figure() pylab.plot(strain * 100, stresses_factored.inUnitsOf(GPa)) pylab.title('CNT with Stone-Wales defect') pylab.xlabel('Strain (%)') pylab.ylabel('Stress (GPa)') pylab.grid(True) pylab.savefig("mdtrajectory_cnt_fig.png")
The script is now ready to be executed. You may download the script to
make sure you did the manual editing correctly:
Visualize and analyse the results¶
The nanotube lengths and stresses that occur during the MD simulation run
are saved in the
.nc file with different ID’s (“C length”, “C stress”,
and “C stresses factored”). These are used to produce the plot titled
mdtrajectory_cnt_fig.png, which is automatically saved in the
project folder. It shows the stress vs. strain relationship for the nanotube,
which is also shown below. The linear (elastic) behavior for small