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