stress.xm

一个很好的分子动力学程序

#######################################################################
#                                                                     #
#  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 <stress.str>.  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

# 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