econst.xmo

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

XMD (Version 2.5.29 Nov 27 2000)

Not using pthread library.
Not using asm("finit"); patch in thread routines.
#######################################################################
#                                                                     #
#  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              #
#                                                                     #
#                                                                     #
#                                                                     #
#######################################################################
#
#
#  Allocate storage
allocate 2000 40000
WARNING:  Allocate command is now obsolete.
#
#  Read in NiAl potential
#
read ../nial.txt
#
#  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
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
EPOT  -7.136969522e-12

#
#
#
#  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
A0 3.32
EPOT  -6.792916665e-12
A0 3.42
EPOT  -7.058508826e-12
A0 3.52
EPOT  -7.136969522e-12
A0 3.62
EPOT  -7.071784393e-12
A0 3.72
EPOT  -6.898437339e-12

#  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
EPOT  -7.136969522e-12


#
#
#  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
AX0 3.52
EPOT  -7.136969522e-12
AX0 3.62
EPOT  -7.128511241e-12
AX0 3.72
EPOT  -7.102921909e-12

#  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
EPOT  -7.136969522e-12

TIME INFORMATION
   Ending   Time: 03/19/2001   13:37:10
   Elapsed  Time: 0 min 01 sec

ERROR STATISTICS
   Number of Fatal Errors:               0
   Number of Unknown Command Errors:     0
   Number of Misc. Warnings:             1   <=======  !!!