####################################################################### # # # 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