#######################################################################
# #
# Example: Calculate the internal stress #
# Use this to find lattice parameter at a finite #
# temperature. #
# #
# NOTE: This method is obsolete. Instead, see xmd documentaiton #
# about the PRESSURE CLAMP command. (JR 30 Jun 1999). #
# #
#######################################################################
# #
# In order to do simulations at a finite temperature, one needs #
# to know what the lattice parameters are at that temperature. (At #
# this point, we are including all lattice, such as tetragonal and #
# monoclinic, which can have up to 6 lattice parameters. Later on #
# we will restrict the discussion to cubic lattice, which have only #
# one lattice parameter). This can be done by plotting the time #
# averaged internal stress of a lattice versus the lattice #
# parameters. From this plot one can interpolate the value of the #
# lattice parameter that gives zero stress. #
# #
# Here we show how to calculate the internal stress for at #
# lattice at a finite temperature and for one value of the lattice #
# parameter. This example uses a cubic NiAl lattice, so there is #
# only one lattice parameter, which we choose to be 2.9 angstroms. #
# We run this simulation with a CLAMP value of 300K. The stress #
# values are calculated and written to a file using the SSAVE #
# command, which is analogous to the ESAVE command. The command #
# #
# SSAVE 1 stress.str #
# #
# saves the stresses in a text file called . The #
# first column of this file holds the step number, the next six #
# columns holds the Voight stresses, s1 through s6. #
# #
# After the run, each column of stress values should be plotted #
# versus time step. During the early time steps the stresses will #
# be equilibrating. After this initial period the stresses should #
# settle into a period of steady fluctuation. The time averages of #
# the stresses after the initial period should be calculated. #
# These average stresses can now be plotted on a different graph as #
# a function of lattice constant. #
# #
# With a cubic lattice, as we have in our example, symmetry #
# dictates that on average, s1, s2 and s3 will be equal and s4, s5 #
# and s6 will be zero. Thus instead of plotting 6 stress in our #
# second graph, we need only plot one, the average of s1, s2 and #
# s3. The other three can be taken to be zero. #
# #
# We then repeat the process for another value of the lattice #
# parameters, producing a stress versus time step plot from each #
# run, which in turn produces new points to add to the second plot, #
# the average stress versus lattice parameter. #
# #
#######################################################################
# Set cubic lattice parameter
calc A0 = 2.9
# Read potential for nial
read ../nial.txt
# Use old neighbor search algorithm to accomodate small atomic array
nsearch sort
# Make repeating box and lattice (in units of a0)
box 6 6 6
particle 2
1 0.25 0.25 0.25
2 0.75 0.75 0.75
dup 5 1 0 0
dup 5 0 1 0
dup 5 0 0 1
# Scale up to units of angstroms (2.8712 unit cell)
scale A0
# Save stresses from every dynamics step in file "stress.str"
ssave 1 stress.str
# Set particle masses (in atomic mass units)
select type 1
mass 58.71
select type 2
mass 26.982
dtime 3.5e-15
# Set adiabtic simulation at starting temperature of 200K
clamp 300
itemp 300
# Perform dynamics
cmd 100