#######################################################################
# #
# Example: Calculate elastic constant C44 for cubic lattice #
# #
#######################################################################
# #
# #
# METHOD #
# #
# Elastic constants can generally be obtained by applying a #
# strain to the lattice and calculating the second #
# derviative of the energy. With XMD, only strains which #
# conform to the repeating boundary conditions can be #
# applied - shears, for example, will not work. Strain are #
# typically applied using the SCALE command to contract or #
# expand the lattice and box together. #
# #
# This example shows how to obtain the shear elastic constant #
# C44. The method used is the same as for the elastic constant #
# C', except that the lattice is oriented within the simulation #
# box as 110, 1-10, 001 (instead of 100, 010, 001 as for the #
# case of C').
# #
# #
# SHEAR CONSTANT C44 #
# #
# #
# Here we subject the lattice to a dilation in one 110 #
# direction and a contraction along the other, while #
# maintaining a constant volume. We orient the fcc cubic #
# lattice within the simulation box so the the lattice 110 #
# direction lies in the x direction, and the 1-10 lies in #
# the y direction. The z direction is 001. #
# #
# The dilation and contraction is implemented with the #
# scale command, such as #
# #
# SCALE 1.01 1/1.01 1.00 #
# #
# which applies the dilation and contraction, but maintains #
# a constant volume. #
# #
# We apply this scale several times, writing the energy each #
# time, in order to get a sample of the energy as a function #
# of strain. The shear constant is then given by the formula #
# #
# #
# 2 2 #
# a0 d E #
# C44 = ---- --- #
# 2 #
# 4 V d a #
# #
# #
# where a is the lattice constant in either the long or short #
# direction. You can use either as long as you #
# remain consistent during the course of the #
# calculation. #
# a0 is the value of a (lattice constant) at equilibrium #
# V is the volume per atom #
# E is the energy per atom as a function of a #
# #
# #
# #
#######################################################################
#
# Read in NiAl potential
#
read ../nial.txt
#
# Use old neighbor search algorithm to accomodate small atomic array
nsearch sort
#
# Set the size of the simulation box
# NX,NY,NZ are the number of unit cells in the X, Y and Z directions
#
calc NX=6
calc NY=6
calc NZ=4
#
#
# Create an FCC lattice with unit cell length of 1 angstrom
# and with orientation 110 1-10 001. This is done by
# first creating a cubic BCC lattice, and then applying a
# "Bain" strain to obtain an FCC lattice in the desired orientation.
#
box NX NY NZ
particle 2
1 0.25 0.25 0.25
1 0.75 0.75 0.75
dup (NX-1) 1 0 0
dup (NY-1) 0 1 0
dup (NZ-1) 0 0 1
#
# Apply "Bain" strain
#
scale 1/sqrt(2) 1/sqrt(2) 1
#
# Scale lattice to unit cell length of Ni, 3.52 angstroms
#
calc A0=3.52
scale A0
# Save energy of original lattice for later comparison
write energy
#
#
# SHEAR MODULUS C44
#
#
# DEL is amount to change lattice constant
# AX0 is old lattice constant in X direction
# AY0 is old lattice constant in Y direction
# AX is new lattice constant in X direction
# AY is new lattice constant in Y direction
#
calc DEL=0.1
calc AX0=A0
calc AY0=A0
repeat 3
# Write energy for current lattice
WRITE AX0
write energy
# Scale X and Y directions separately
calc AX=AX0+DEL
calc AY=AY0-DEL
scale AX/AX0 AY/AY0 1.0
# Save new lattice sizes
calc AX0 = AX
calc AY0 = AY
end
# revert to equilibrium lattice size
calc AX=3.52
calc AY=3.52
scale AX/AX0 AY/AY0 1.0
# Write energy to test that lattice is original one
write energy