碳纳米管以其异常的机械性质而著称；被测得有1TPa量级的杨氏模量[KDE+98]，高于所有已知材料。对于无缺陷的碳纳米管杨氏模量的理论计算相当简单，但是当有缺陷存在时就另当别论了。对此，分子动力学（MD）是计算机械性能的合适工具。 在本实例中你将构造具有Stone–Wales缺陷的碳纳米管（CNT），并使用ATK-Classical来计算杨氏模量。

打开VNL Builder, 使用 File ‣ Add ‣ From Files将之前保存的CNT构象导入Stash。注意到这个结构是一个中间区域有Stone–Wales缺陷的器件构象。

接下来，选中器件的Stash项（显示为高亮），点击按钮将其转化为块体构象（只是将器件电极移除），并使用Selection Tools ‣ Tags工具将纳米管两端的原子添加标签：

• 使用鼠标将最左端的碳原子选中（显示为高亮）并将它们贴上标签“Left”。
• 对最右端的碳原子做同样的操作，如下图所示。

这些标签将会被用来在分子动力学（MD）模拟中对结构施加约束。

将带有Stone–Wales缺陷的CNT块体构象送到Script Generator并添加 New Calculator和Optimization ‣ MolecularDynamics 模块。然后编辑各个模块来对计算进行设置：

• 首先，将默认输出文件改为mdtrajectory_cnt.nc。

• 在 New Calculation模块中选择ATK-Classical 计算器，并确保使用 Tersoff_C_2010 经典势（为模拟石墨烯和纳米管晶格动力学而设计的势）。

• 在MolecularDynamics模块中进行以下修改：

Molecular Dynamics选项：

• 设置 “Type” 为 NVT Berendsen。
• 设置“Steps” 为 20000。
• 设置“Log interval”为 200。

Initial Velocity选项：

• 设置“Type” 为Maxwell-Boltzmann。
• 设置“Temperature” 为300 K。
• 点掉 Remove center-of-mass momentum 勾选项。

NVT Berendsen选项：

• 设置“Thermal coupling”为 1 fs。

如上图所示，保持其他设置为默认值。

将脚本保存为mdtrajectory_cnt.py 并使用按钮将其送到 Editor。

在ATK中可以定义“挂钩”；即在每个MD步之前或之后被触发的Python函数或类。这些挂钩可被用来提供自定义的实时分析或者监测系统，或者-如本实例-在模拟的过程中动态地修改结构。

你现在将定义一个预处理挂钩用以在每一步分子动力学模拟中给纳米管一个小的拉伸：

• 在Editor中打开Python脚本，在脚本底部找到定义md_trajectory的模块。之前定义的“Left”和 “Right” 标签现在可以被用来对系统添加约束。改变语句

constraints=[]

为

constraints=bulk_configuration.indicesFromTags(["Left","Right"])

• 在 md_trajectory 模块添加语句

pre_step_hook=strainer,

并在它紧前面的语句添加一个逗号(method=method,)。

• 你也应当定义什么是“strainer”：在 md_trajectory 模块紧前面添加以下语句：

strainer = StressStrainUtility()

• 脚本底部现在应该如这样：

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------
 
initial_velocity = MaxwellBoltzmannDistribution(
    temperature=300.0*Kelvin,
    remove_center_of_mass_momentum=False
)
 
method = NVTBerendsen(
    time_step=1*femtoSecond,
    reservoir_temperature=300*Kelvin,
    thermal_coupling_constant=1*femtoSecond,
    initial_velocity=initial_velocity
)
 
strainer = StressStrainUtility()
 
md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=bulk_configuration.indicesFromTags(["Left","Right"]),
    trajectory_filename='trajectory.nc',
    steps=20000,
    log_interval=200,
    method=method,
    pre_step_hook=strainer,
)
 
bulk_configuration = md_trajectory.lastImage()

• 上图的17行初始化了类 StressStrainUtility的一个实例，现在你将定义这个Python类。在脚本最顶部嵌入以下语句：

# 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=1000):
        """
        @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):
        """
        Strainer call.
        """
        # Get the lattice vectors.
        u1, u2, u3 = configuration.bravaisLattice().primitiveVectors()
 
#        Don't apply strain in the earliest MD steps
#        if step < self.__skip:
#            return
 
        # Calculate the stress.
        stress = Stress(configuration)
 
        # Store the history
        self.stresses.append( stress.evaluate()[2,2].inUnitsOf(eV/Angstrom**3) )
        self.lengths.append( u3[2].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 constrainted on the right side.
        constraints = configuration.indicesFromTags("Right")
        for c in constraints:
            xyz = numpy.array(configuration.cartesianCoordinates()[c]) + numpy.array(self.__delta)
            configuration._changeAtoms(indices=[c], positions=[xyz]*Angstrom)
 
    def length(self):
        return numpy.array(self.lengths)*Angstrom
 
    def stress(self):
        return numpy.array(self.stresses)*eV/Angstrom**3

这个类定义了一个可以对单胞沿z方向拉伸的方法。你将使用它来研究一个有缺陷的CNT如何对应变做出响应。 正如你在脚本中所见，这个类有两个输出，分别是delta 和skip：

delta 定义了每一步MD模拟中u3矢量（z方向）扩展了多少。

skip 定义了在构象被拉之前经过多少步。

实际的拉伸是通过41行起的“for-loop” 完成的，沿着z方向以delta 值增加被标记为“Right”的原子的坐标。在30和31行，计算应力和构象总长度

并分开地列着。这将在之后的绘图中用到。