From f343a9f4a077514801edb493a170f6b247559e16 Mon Sep 17 00:00:00 2001 From: athomps Date: Mon, 5 Oct 2015 18:22:14 +0000 Subject: [PATCH] Added elastic constant example at finite temperature git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14097 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- examples/ELASTIC/in.elastic | 6 - examples/ELASTIC_T/Si.sw | 18 +++ examples/ELASTIC_T/displace.mod | 140 +++++++++++++++++++++++ examples/ELASTIC_T/in.elastic | 184 +++++++++++++++++++++++++++++++ examples/ELASTIC_T/init.mod | 37 +++++++ examples/ELASTIC_T/potential.mod | 27 +++++ examples/README | 6 +- 7 files changed, 411 insertions(+), 7 deletions(-) create mode 100644 examples/ELASTIC_T/Si.sw create mode 100644 examples/ELASTIC_T/displace.mod create mode 100644 examples/ELASTIC_T/in.elastic create mode 100644 examples/ELASTIC_T/init.mod create mode 100644 examples/ELASTIC_T/potential.mod diff --git a/examples/ELASTIC/in.elastic b/examples/ELASTIC/in.elastic index ccd2b7367d..90f2bf38ca 100644 --- a/examples/ELASTIC/in.elastic +++ b/examples/ELASTIC/in.elastic @@ -42,12 +42,6 @@ # that minimization is not required, you can set maxiter = 0.0 in # init.mod. # -# There are two alternate versions of displace.mod provided. -# They are displace_restart.mod and displace_reverse.mod. -# The former resets the box using a restart file while -# the latter reverses the deformation. Copy whichever -# one you like best to displace.mod. -# include init.mod include potential.mod diff --git a/examples/ELASTIC_T/Si.sw b/examples/ELASTIC_T/Si.sw new file mode 100644 index 0000000000..db4be100ef --- /dev/null +++ b/examples/ELASTIC_T/Si.sw @@ -0,0 +1,18 @@ +# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985) +# Stillinger-Weber parameters for various elements and mixtures +# multiple entries can be added to this file, LAMMPS reads the ones it needs +# these entries are in LAMMPS "metal" units: +# epsilon = eV; sigma = Angstroms +# other quantities are unitless + +# format of a single entry (one or more lines): +# element 1, element 2, element 3, +# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol + +# Here are the original parameters in metal units, for Silicon from: +# +# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985) +# + +Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333 + 7.049556277 0.6022245584 4.0 0.0 0.0 diff --git a/examples/ELASTIC_T/displace.mod b/examples/ELASTIC_T/displace.mod new file mode 100644 index 0000000000..b121f2f694 --- /dev/null +++ b/examples/ELASTIC_T/displace.mod @@ -0,0 +1,140 @@ +# NOTE: This script should not need to be +# modified. See in.elastic for more info. +# +# Find which reference length to use + +if "${dir} == 1" then & + "variable len0 equal ${lx0}" +if "${dir} == 2" then & + "variable len0 equal ${ly0}" +if "${dir} == 3" then & + "variable len0 equal ${lz0}" +if "${dir} == 4" then & + "variable len0 equal ${lz0}" +if "${dir} == 5" then & + "variable len0 equal ${lz0}" +if "${dir} == 6" then & + "variable len0 equal ${ly0}" + +# Reset box and simulation parameters + +clear +box tilt large +read_restart restart.equil +include potential.mod + +# Negative deformation + +variable delta equal -${up}*${len0} +variable deltaxy equal -${up}*xy +variable deltaxz equal -${up}*xz +variable deltayz equal -${up}*yz +if "${dir} == 1" then & + "change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} remap units box" +if "${dir} == 2" then & + "change_box all y delta 0 ${delta} yz delta ${deltayz} remap units box" +if "${dir} == 3" then & + "change_box all z delta 0 ${delta} remap units box" +if "${dir} == 4" then & + "change_box all yz delta ${delta} remap units box" +if "${dir} == 5" then & + "change_box all xz delta ${delta} remap units box" +if "${dir} == 6" then & + "change_box all xy delta ${delta} remap units box" + +# Run MD + +run ${nequil} +include potential.mod +run ${nrun} + +# Obtain new stress tensor + +variable pxx1 equal f_avp[1] +variable pyy1 equal f_avp[2] +variable pzz1 equal f_avp[3] +variable pxy1 equal f_avp[4] +variable pxz1 equal f_avp[5] +variable pyz1 equal f_avp[6] + +# Compute elastic constant from pressure tensor + +variable C1neg equal ${d1} +variable C2neg equal ${d2} +variable C3neg equal ${d3} +variable C4neg equal ${d4} +variable C5neg equal ${d5} +variable C6neg equal ${d6} + +# Reset box and simulation parameters + +clear +box tilt large +read_restart restart.equil +include potential.mod + +# Positive deformation + +variable delta equal ${up}*${len0} +variable deltaxy equal ${up}*xy +variable deltaxz equal ${up}*xz +variable deltayz equal ${up}*yz +if "${dir} == 1" then & + "change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} remap units box" +if "${dir} == 2" then & + "change_box all y delta 0 ${delta} yz delta ${deltayz} remap units box" +if "${dir} == 3" then & + "change_box all z delta 0 ${delta} remap units box" +if "${dir} == 4" then & + "change_box all yz delta ${delta} remap units box" +if "${dir} == 5" then & + "change_box all xz delta ${delta} remap units box" +if "${dir} == 6" then & + "change_box all xy delta ${delta} remap units box" + +# Run MD + +run ${nequil} +include potential.mod +run ${nrun} + +# Obtain new stress tensor + +variable tmp equal pe +variable e1 equal ${tmp} +variable tmp equal press +variable p1 equal ${tmp} +variable tmp equal pxx +variable pxx1 equal ${tmp} +variable tmp equal pyy +variable pyy1 equal ${tmp} +variable tmp equal pzz +variable pzz1 equal ${tmp} +variable tmp equal pxy +variable pxy1 equal ${tmp} +variable tmp equal pxz +variable pxz1 equal ${tmp} +variable tmp equal pyz +variable pyz1 equal ${tmp} + +# Compute elastic constant from pressure tensor + +variable C1pos equal ${d1} +variable C2pos equal ${d2} +variable C3pos equal ${d3} +variable C4pos equal ${d4} +variable C5pos equal ${d5} +variable C6pos equal ${d6} + +# Combine positive and negative + +variable C1${dir} equal 0.5*(${C1neg}+${C1pos}) +variable C2${dir} equal 0.5*(${C2neg}+${C2pos}) +variable C3${dir} equal 0.5*(${C3neg}+${C3pos}) +variable C4${dir} equal 0.5*(${C4neg}+${C4pos}) +variable C5${dir} equal 0.5*(${C5neg}+${C5pos}) +variable C6${dir} equal 0.5*(${C6neg}+${C6pos}) + +# Delete dir to make sure it is not reused + +variable dir delete diff --git a/examples/ELASTIC_T/in.elastic b/examples/ELASTIC_T/in.elastic new file mode 100644 index 0000000000..8427192a82 --- /dev/null +++ b/examples/ELASTIC_T/in.elastic @@ -0,0 +1,184 @@ +# Compute elastic constant tensor for a crystal at finite temperature +# +# Written by Aidan Thompson (Sandia, athomps@sandia.gov) +# +# This script uses the following three include files. +# +# init.mod (must be modified for different crystal structures) +# Define units, MD parameters, deformation parameters, +# and initial configuration of the atoms and simulation cell. +# +# +# potential.mod (must be modified for different pair styles) +# Define pair style and other attributes +# not stored in restart file +# +# +# displace.mod (displace.mod should not need to be modified) +# Perform positive and negative box displacements +# in direction ${dir} and size ${up}. +# It uses the resultant changes +# in stress to compute one +# row of the elastic stiffness tensor +# +# Inputs variables: +# dir = the Voigt deformation component +# (1,2,3,4,5,6) +# Global constants: +# up = the deformation magnitude (strain units) +# cfac = conversion from LAMMPS pressure units to +# output units for elastic constants +# +# +# To run this on a different system, it should only be necessary to +# modify the files init.mod and potential.mod. In order to calculate +# the elastic constants correctly, care must be taken to specify +# the correct units in init.mod (units, cfac and cunits). It is also +# important to verify that the MD sampling of stress components +# is generating accurate statistical averages. +# One indication of this is that the elastic constants are insensitive +# to the choice of the variable ${up} in init.mod. Another is to +# check for finite size effects. +# + +include init.mod + +# Compute initial state + +include potential.mod +run ${nequil} +include potential.mod +run ${nrun} + +variable pxx0 equal f_avp[1] +variable pyy0 equal f_avp[2] +variable pzz0 equal f_avp[3] +variable pxy0 equal f_avp[4] +variable pxz0 equal f_avp[5] +variable pyz0 equal f_avp[6] + +variable tmp equal lx +variable lx0 equal ${tmp} +variable tmp equal ly +variable ly0 equal ${tmp} +variable tmp equal lz +variable lz0 equal ${tmp} + +# These formulas define the derivatives w.r.t. strain components +# Constants uses $, variables use v_ +variable d1 equal -(v_pxx1-${pxx0})/(v_delta/v_len0)*${cfac} +variable d2 equal -(v_pyy1-${pyy0})/(v_delta/v_len0)*${cfac} +variable d3 equal -(v_pzz1-${pzz0})/(v_delta/v_len0)*${cfac} +variable d4 equal -(v_pyz1-${pyz0})/(v_delta/v_len0)*${cfac} +variable d5 equal -(v_pxz1-${pxz0})/(v_delta/v_len0)*${cfac} +variable d6 equal -(v_pxy1-${pxy0})/(v_delta/v_len0)*${cfac} + +# Write restart +write_restart restart.equil + +# uxx Perturbation + +variable dir equal 1 +include displace.mod + +# uyy Perturbation + +variable dir equal 2 +include displace.mod + +# uzz Perturbation + +variable dir equal 3 +include displace.mod + +# uyz Perturbation + +variable dir equal 4 +include displace.mod + +# uxz Perturbation + +variable dir equal 5 +include displace.mod + +# uxy Perturbation + +variable dir equal 6 +include displace.mod + +# Output final values + +variable C11all equal ${C11} +variable C22all equal ${C22} +variable C33all equal ${C33} + +variable C12all equal 0.5*(${C12}+${C21}) +variable C13all equal 0.5*(${C13}+${C31}) +variable C23all equal 0.5*(${C23}+${C32}) + +variable C44all equal ${C44} +variable C55all equal ${C55} +variable C66all equal ${C66} + +variable C14all equal 0.5*(${C14}+${C41}) +variable C15all equal 0.5*(${C15}+${C51}) +variable C16all equal 0.5*(${C16}+${C61}) + +variable C24all equal 0.5*(${C24}+${C42}) +variable C25all equal 0.5*(${C25}+${C52}) +variable C26all equal 0.5*(${C26}+${C62}) + +variable C34all equal 0.5*(${C34}+${C43}) +variable C35all equal 0.5*(${C35}+${C53}) +variable C36all equal 0.5*(${C36}+${C63}) + +variable C45all equal 0.5*(${C45}+${C54}) +variable C46all equal 0.5*(${C46}+${C64}) +variable C56all equal 0.5*(${C56}+${C65}) + +variable C11cubic equal (${C11all}+${C22all}+${C33all})/3.0 +variable C12cubic equal (${C12all}+${C13all}+${C23all})/3.0 +variable C44cubic equal (${C44all}+${C55all}+${C66all})/3.0 +variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0 +variable shearmodulus1 equal ${C44cubic} +variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0 +variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic}) + +# For Stillinger-Weber silicon, the analytical results +# are known to be (E. R. Cowley, 1988): +# C11 = 151.4 GPa +# C12 = 76.4 GPa +# C44 = 56.4 GPa + +print "Elastic Constant C11all = ${C11all} ${cunits}" +print "Elastic Constant C22all = ${C22all} ${cunits}" +print "Elastic Constant C33all = ${C33all} ${cunits}" + +print "Elastic Constant C12all = ${C12all} ${cunits}" +print "Elastic Constant C13all = ${C13all} ${cunits}" +print "Elastic Constant C23all = ${C23all} ${cunits}" + +print "Elastic Constant C44all = ${C44all} ${cunits}" +print "Elastic Constant C55all = ${C55all} ${cunits}" +print "Elastic Constant C66all = ${C66all} ${cunits}" + +print "Elastic Constant C14all = ${C14all} ${cunits}" +print "Elastic Constant C15all = ${C15all} ${cunits}" +print "Elastic Constant C16all = ${C16all} ${cunits}" + +print "Elastic Constant C24all = ${C24all} ${cunits}" +print "Elastic Constant C25all = ${C25all} ${cunits}" +print "Elastic Constant C26all = ${C26all} ${cunits}" + +print "Elastic Constant C34all = ${C34all} ${cunits}" +print "Elastic Constant C35all = ${C35all} ${cunits}" +print "Elastic Constant C36all = ${C36all} ${cunits}" + +print "Elastic Constant C45all = ${C45all} ${cunits}" +print "Elastic Constant C46all = ${C46all} ${cunits}" +print "Elastic Constant C56all = ${C56all} ${cunits}" + +print "Bulk Modulus = ${bulkmodulus} ${cunits}" +print "Shear Modulus 1 = ${shearmodulus1} ${cunits}" +print "Shear Modulus 2 = ${shearmodulus2} ${cunits}" +print "Poisson Ratio = ${poissonratio}" diff --git a/examples/ELASTIC_T/init.mod b/examples/ELASTIC_T/init.mod new file mode 100644 index 0000000000..fb13db9981 --- /dev/null +++ b/examples/ELASTIC_T/init.mod @@ -0,0 +1,37 @@ +# NOTE: This script can be modified for different atomic structures, +# units, etc. See in.elastic for more info. +# + +# Define the finite deformation size. Try several values of this +# variable to verify that results do not depend on it. +variable up equal 1.0e-1 + +# metal units, elastic constants in GPa +units metal +variable cfac equal 1.0e-4 +variable cunits string GPa + +# Define MD parameters +variable nevery equal 10 +variable nrepeat equal 10 +variable nfreq equal ${nevery}*${nrepeat} +variable nthermo equal ${nfreq} +variable nequil equal 10*${nthermo} +variable nrun equal 3*${nthermo} +variable temp equal 2000.0 +variable timestep equal 0.001 +variable mass1 equal 28.06 +variable tdamp equal 0.01 +variable seed equal 123457 + +# generate the box and atom positions using a diamond lattice +variable a equal 5.431 + +boundary p p p + +lattice diamond $a +region box prism 0 3.0 0 3.0 0 3.0 0.0 0.0 0.0 +create_box 1 box +create_atoms 1 box +mass 1 ${mass1} + diff --git a/examples/ELASTIC_T/potential.mod b/examples/ELASTIC_T/potential.mod new file mode 100644 index 0000000000..9477e7285a --- /dev/null +++ b/examples/ELASTIC_T/potential.mod @@ -0,0 +1,27 @@ +# NOTE: This script can be modified for different pair styles +# See in.elastic for more info. + +reset_timestep 0 + +# Choose potential +pair_style sw +pair_coeff * * Si.sw Si + +# Setup neighbor style +neighbor 1.0 nsq +neigh_modify once no every 1 delay 0 check yes + +# Setup output + +fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} c_thermo_press mode vector +thermo ${nthermo} +thermo_style custom step temp pe press f_avp[1] f_avp[2] f_avp[3] f_avp[4] f_avp[5] f_avp[6] +thermo_modify norm no + +# Setup MD + +timestep ${timestep} +fix 4 all nve +fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed} + + diff --git a/examples/README b/examples/README index 60fc4bebc1..c2a9023108 100644 --- a/examples/README +++ b/examples/README @@ -144,9 +144,13 @@ either by itself or in tandem with another code or library. See the COUPLE/README file to get started. The ELASTIC directory has an example script for computing elastic -constants, using a zero temperature Si example. See the +constants at zero temperature, using an Si example. See the ELASTIC/in.elastic file for more info. +The ELASTIC_T directory has an example script for computing elastic +constants at finite temperature, using an Si example. See the +ELASTIC_T/in.elastic file for more info. + The KAPPA directory has example scripts for computing the thermal conductivity (kappa) of a LJ liquid using 4 different methods. See the KAPPA/README file for more info.