terval=200, method=method, pre_step_hook=strainer, )
Line 17 above initializes an instance of the class
StressStrainUtility, which is a Python class that you now will define. Insert the following at the very top of the script:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
# Define a prestep function which can expand the system class StressStrainUtility: """ Perform a strain of a configuration over time. """ def __init__(self, delta=0.0001, skip=0): """ @param delta The strain to apply in each step @param first steps to "skip" strain """ self.__delta = [0.0,0.0,delta]*Angstrom self.__skip = skip self.stresses =  self.lengths =  def __call__(self, step, time, configuration, forces, stress): """ Strainer call. """ # Get the lattice vectors. u1, u2, u3 = configuration.primitiveVectors() # Don't apply strain in the earliest MD steps if step < self.__skip: return # Store the history self.stresses.append( Stress(configuration).evaluate()[2,2].inUnitsOf(eV/Angstrom**3) ) self.lengths.append( u3.inUnitsOf(Angstrom) ) # Increase the lattice vectors. u3 = u3 + self.__delta # Set the new lattice. configuration.setBravaisLattice(UnitCell(u1,u2,u3)) # Move all the atoms that are constrained on the right side. constraints = configuration.indicesFromTags("Right") for c in constraints: xyz = numpy.array(configuration.cartesianCoordinates()[c]) + numpy.array(self.__delta) configuration.setCartesianCoordinates(indices=[c], cartesian_coordinates=[xyz]*Angstrom) def length(self): return numpy.array(self.lengths)*Angstrom def stress(self): return numpy.array(self.stresses)*eV/Angstrom**3
This class defines a method to enable straining of a unit cell along the z-direction. You will use it to study how a CNT with a defect responds to strain.
As you can see from the script, the class has two inputs,
- defines how much the u3 vector (the z-direction) is expanded on every step in the MD calculation.
- defines how many steps should pass before the configuration 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)*