#######################################################################
# #
# Example: Calculate elastic constants B and C' #
# #
#######################################################################
# #
# #
# 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 bulk modulus B and #
# the shear constant C'. The elastic constants C11 and C12 #
# are related to B and C' as follows. #
# #
# B = (C11 + 2*C12) / 3 #
# #
# C' = (C11 - C12) / 2 #
# #
# #
# In order to obtain the last of the three cubic lattice #
# elastic constants, C44, we must use another orientation #
# of the lattice. That will wait for another example. #
# #
# BULK MODULUS B #
# #
# We subject the lattice to a uniform dilation, using the #
# command SCALE 1.01. Applying this several times in #
# sucession gives a series of energies from which the #
# second derivative of the energy may be found. The bulk #
# modulus is then obtained from the following formula. #
# #
# 2 2 #
# a0 d E #
# B = ---- --- #
# 2 #
# 9 V d a #
# #
# where a is the lattice constant #
# 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 #
# #
# #
# SHEAR CONSTANT C' #
# #
# Here we subject the lattice to a dilation in one #
# direction (say x for example) a contraction in the second #
# direction (say y) and no change in the remaining #
# direction, while maintaining a constant volume, using a #
# command such as #
# #
# SCALE 1.01 1/1.01 1.00 #
# #
# Again as before we apply this several times, writing the #
# energy as as before, and use these energies to obtain the #
# energy second derivative. The shear constant is then given #
# by the formula #
# #
# #
# 2 2 #
# a0 d E #
# C' = ---- --- #
# 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=4
calc NY=4
calc NZ=4
#
#
# Create an FCC lattice with unit cell length of 1 angstrom
#
box NX NY NZ
particle 4
1 0.25 0.25 0.25
1 0.25 0.75 0.75
1 0.75 0.25 0.75
1 0.75 0.75 0.25
dup (NX-1) 1 0 0
dup (NY-1) 0 1 0
dup (NZ-1) 0 0 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
#
#
#
# BULK MODULUS
#
# Shrink the lattice a little to start so that we sample energies on
# both sides of equlibrium
#
# DEL is amount to change lattice constant
# A0 is old lattice constant
# ANEW is new lattice constant
#
calc DEL=0.1
calc ANEW=A0-2*DEL
scale ANEW/A0
calc A0=ANEW
repeat 5
# Write energy for current lattice size
WRITE A0
write energy
# Scale to larger lattice size
calc ANEW=A0+DEL
scale ANEW/A0
# Save new lattice size
calc A0=ANEW
end
# revert to equilibrium lattice size
calc ANEW=3.52
scale ANEW/A0
calc A0= ANEW
# Write energy to test that lattice is original one
write energy
#
#
# SHEAR MODULUS
#
#
# 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