Merge branch 'Elastic_stress' of github.com:Bibobu/lammps into Elastic_stress

doc/src/Howto_elastic.rst

@@ -18,23 +18,52 @@ At zero temperature, it is easy to estimate these derivatives by
 deforming the simulation box in one of the six directions using the
 :doc:`change_box <change_box>` command and measuring the change in the
 stress tensor. A general-purpose script that does this is given in the
-examples/elastic directory described on the :doc:`Examples <Examples>`
+examples/ELASTIC directory described on the :doc:`Examples <Examples>`
 doc page.
 
 Calculating elastic constants at finite temperature is more
 challenging, because it is necessary to run a simulation that performs
-time averages of differential properties. One way to do this is to
-measure the change in average stress tensor in an NVT simulations when
+time averages of differential properties. There are at least
+3 ways to do this in LAMMPS. The most reliable way to do this is
+by exploiting the relationship between elastic constants, stress
+fluctuations, and the Born matrix, the second derivatives of energy
+w.r.t. strain :ref:`(Ray) <Ray>`.
+The Born matrix calculation has been enabled by
+the :doc:`compute born/matrix <compute_born_matrix>` command,
+which works for any bonded or non-bonded potential in LAMMPS.
+The most expensive part of the calculation is the sampling of
+the stress fluctuations. Several examples of this method are
+provided in the examples/ELASTIC_T/BORN_MATRIX directory
+described on the :doc:`Examples <Examples>` doc page.
+
+A second way is to measure
+the change in average stress tensor in an NVT simulations when
 the cell volume undergoes a finite deformation. In order to balance
 the systematic and statistical errors in this method, the magnitude of
 the deformation must be chosen judiciously, and care must be taken to
 fully equilibrate the deformed cell before sampling the stress
-tensor. Another approach is to sample the triclinic cell fluctuations
+tensor. An example of this method is provided in the
+examples/ELASTIC_T/DEFORMATION directory
+described on the :doc:`Examples <Examples>` doc page.
+
+Another approach is to sample the triclinic cell fluctuations
 that occur in an NPT simulation. This method can also be slow to
-converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`
+converge and requires careful post-processing :ref:`(Shinoda) <Shinoda1>`.
+We do not provide an example of this method.
+
+A nice review of the advantages and disadvantages of all of these methods
+is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier>`.
 
 ----------
 
+.. _Ray:
+
+**(Ray)** J. R. Ray and A. Rahman, J Chem Phys, 80, 4423 (1984).
+
 .. _Shinoda1:
 
 **(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
+
+.. _Clavier:
+
+**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).

doc/src/compute_born_matrix.rst

@@ -80,7 +80,7 @@ in this matrix the indices of :math:`C_{k}` value are the corresponding element
 :math:`k` in the global vector output by this compute. Each term comes from the sum
 of the derivatives of every contribution to the potential energy
 in the system as explained in :ref:`(VanWorkum)
-<VanWorkum>` or :ref:`(Voyiatzis) <Voyiatzis>`.
+<VanWorkum>`.
 
 The output can be accessed using usual Lammps routines:
 
@@ -99,27 +99,32 @@ The output can be accessed using usual Lammps routines:
 In this example, the file *born.out* will contain the information needed to
 compute the first and second terms of the elastic constant matrix in a post
 processing procedure. The other required quantities can be accessed using any
-other *LAMMPS* usual method.
+other *LAMMPS* usual method. Several examples of this method are
+provided in the examples/ELASTIC_T/BORN_MATRIX directory
+described on the :doc:`Examples <Examples>` doc page.
 
-NOTE: In the above :math:`C_{i,j}` computation, the term involving the virial
-stress tensor :math:`\sigma` is the covariance between each elements. In a
-solid the virial stress can have large variations between timesteps and average
-values can be slow to converge. This term is better computed using
-instantaneous values.
+NOTE: In the above :math:`C_{i,j}` computation, the fluctuation
+term involving the virial stress tensor :math:`\sigma` is the
+covariance between each elements. In a
+solid the stress fluctuations can vary rapidly, while average
+fluctuations can be slow to converge.
+A detailed analysis of the convergence rate of all the terms in
+the elastic tensor
+is provided in the paper by Clavier et al. :ref:`(Clavier) <Clavier2>`.
 
-Two different computation methods are implemented in this compute and are
-mutually exclusive.
+Two different computation methods for the Born matrix are implemented in this
+compute and are mutually exclusive.
 
 The first one is a direct computation from the analytical formula from the
-different terms of the potential used for the simulations (see :ref:`(VanWorkum)
-<VanWorkum>`). However, the implementation of such derivations must be done
+different terms of the potential used for the simulations :ref:`(VanWorkum)
+<VanWorkum>`. However, the implementation of such derivations must be done
 for every potential form. This has not been done yet and can be very
 complicated for complex potentials. At the moment a warning message is
 displayed for every term that is not supporting the compute at the moment.
 This method is the default for now.
 
 The second method uses finite differences of energy to numerically approximate
-the second derivatives (see :ref:`(Zhen) <Zhen>`). This is useful when using
+the second derivatives :ref:`(Zhen) <Zhen>`. This is useful when using
 interaction styles for which the analytical second derivatives have not been
 implemented. In this cases, the compute applies linear strain fields of
 magnitude *delta* to all the atoms relative to a point at the center of the
@@ -199,9 +204,9 @@ none
 
 **(Van Workum)** K. Van Workum et al., J. Chem. Phys. 125 144506 (2006)
 
-.. _Voyiatzis:
+.. _Clavier2:
 
-**(Voyiatzis)** E. Voyiatzis, Computer Physics Communications 184(2013)27-33
+**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017).
 
 .. _Zhen:
 

examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw

@@ -0,0 +1 @@
+../../../../potentials/Si.sw

examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in

@@ -0,0 +1,68 @@
+# Average moduli for cubic crystals
+
+variable C11cubic equal (${C11}+${C22}+${C33})/3.0
+variable C12cubic equal (${C12}+${C13}+${C23})/3.0
+variable C44cubic equal (${C44}+${C55}+${C66})/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 "========================================="
+#print "Components of the Elastic Constant Tensor"
+#print "========================================="
+
+print "Elastic Constant C11 = ${C11} ${cunits}"
+print "Elastic Constant C22 = ${C22} ${cunits}"
+print "Elastic Constant C33 = ${C33} ${cunits}"
+
+print "Elastic Constant C12 = ${C12} ${cunits}"
+print "Elastic Constant C13 = ${C13} ${cunits}"
+print "Elastic Constant C23 = ${C23} ${cunits}"
+
+print "Elastic Constant C44 = ${C44} ${cunits}"
+print "Elastic Constant C55 = ${C55} ${cunits}"
+print "Elastic Constant C66 = ${C66} ${cunits}"
+
+print "Elastic Constant C14 = ${C14} ${cunits}"
+print "Elastic Constant C15 = ${C15} ${cunits}"
+print "Elastic Constant C16 = ${C16} ${cunits}"
+
+print "Elastic Constant C24 = ${C24} ${cunits}"
+print "Elastic Constant C25 = ${C25} ${cunits}"
+print "Elastic Constant C26 = ${C26} ${cunits}"
+
+print "Elastic Constant C34 = ${C34} ${cunits}"
+print "Elastic Constant C35 = ${C35} ${cunits}"
+print "Elastic Constant C36 = ${C36} ${cunits}"
+
+print "Elastic Constant C45 = ${C45} ${cunits}"
+print "Elastic Constant C46 = ${C46} ${cunits}"
+print "Elastic Constant C56 = ${C56} ${cunits}"
+
+print "========================================="
+print "Average properties for a cubic crystal"
+print "========================================="
+
+print "Bulk Modulus = ${bulkmodulus} ${cunits}"
+print "Shear Modulus 1 = ${shearmodulus1} ${cunits}"
+print "Shear Modulus 2 = ${shearmodulus2} ${cunits}"
+print "Poisson Ratio = ${poissonratio}"
+
+# summarize sampling protocol
+
+variable tmp equal atoms
+print "Number of atoms = ${tmp}"
+print "Stress sampling interval = ${nevery}"
+variable tmp equal ${nrun}/${nevery}
+print "Stress sample count = ${tmp}"
+print "Born sampling interval = ${neveryborn}"
+variable tmp equal ${nrun}/${neveryborn}
+print "Born sample count = ${tmp}"

examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic

@@ -0,0 +1,25 @@
+# Compute elastic constant tensor for a crystal at finite temperature
+# These settings replicate the 1477~K benchmark of
+# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
+# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
+
+# here: Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
+
+include init.in
+
+# Compute initial state
+
+include potential.in
+thermo_style custom step temp pe press density
+run ${nequil}
+
+# Run dynamics
+
+include potential.in
+include output.in
+
+run ${nrun}
+
+# Output final values
+
+include final_output.in

examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in

@@ -0,0 +1,59 @@
+# NOTE: This script can be modified for different atomic structures, 
+# units, etc. See in.elastic for more info.
+#
+
+# Define MD parameters
+# These can be modified by the user
+# These settings replicate the 1477~K benchmark of
+# Kluge, Ray, and Rahman (1986) that is Ref.[15] in:
+# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265
+
+# select temperature and pressure (lattice constant)
+
+variable temp index 1477.0                 # temperature of initial sample
+variable a index 5.457                     # lattice constant
+
+# select sampling parameters, important for speed/convergence
+
+variable nthermo index 1500               # interval for thermo output
+variable nevery index 10                  # stress sampling interval
+variable neveryborn index 100             # Born sampling interval
+variable timestep index 0.000766          # timestep
+variable nlat index 3                     # number of lattice unit cells
+
+# other settings
+
+variable mass1 index 28.06                # mass
+variable tdamp index 0.01                 # time constant for thermostat
+variable seed index 123457                # seed for thermostat
+variable thermostat index 1               # 0 if NVE, 1 if NVT
+variable delta index 1.0e-6               # Born numdiff strain magnitude 
+
+# hard-coded rules-of-thumb for run length, etc.
+
+variable nfreq equal ${nthermo}           # interval for averaging output
+variable nrepeat equal floor(${nfreq}/${nevery}) # number of samples
+variable nrepeatborn equal floor(${nfreq}/${neveryborn}) # number of samples
+variable nequil equal 10*${nthermo}       # length of equilibration run
+variable nrun equal 100*${nthermo}        # length of equilibrated run
+
+# generate the box and atom positions using a diamond lattice
+
+units		metal
+
+boundary	p p p
+
+# this generates a standard 8-atom cubic cell
+
+lattice         diamond $a
+region		box prism 0 1 0 1 0 1 0 0 0
+
+# this generates a 2-atom triclinic cell
+#include tri.in
+
+create_box	1 box
+create_atoms	1 box
+mass 1 ${mass1}
+replicate ${nlat} ${nlat} ${nlat}
+velocity	all create ${temp} 87287
+

examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in

@@ -0,0 +1,140 @@
+# Setup output
+
+# Stress fluctuation term F
+
+compute stress all pressure thermo_temp
+variable s1 equal c_stress[1]
+variable s2 equal c_stress[2]
+variable s3 equal c_stress[3]
+variable s4 equal c_stress[6]
+variable s5 equal c_stress[5]
+variable s6 equal c_stress[4]
+
+variable s11 equal v_s1*v_s1
+variable s22 equal v_s2*v_s2
+variable s33 equal v_s3*v_s3
+variable s44 equal v_s4*v_s4
+variable s55 equal v_s5*v_s5
+variable s66 equal v_s6*v_s6
+variable s33 equal v_s3*v_s3
+variable s12 equal v_s1*v_s2
+variable s13 equal v_s1*v_s3
+variable s14 equal v_s1*v_s4
+variable s15 equal v_s1*v_s5
+variable s16 equal v_s1*v_s6
+variable s23 equal v_s2*v_s3
+variable s24 equal v_s2*v_s4
+variable s25 equal v_s2*v_s5
+variable s26 equal v_s2*v_s6
+variable s34 equal v_s3*v_s4
+variable s35 equal v_s3*v_s5
+variable s36 equal v_s3*v_s6
+variable s45 equal v_s4*v_s5
+variable s46 equal v_s4*v_s6
+variable s56 equal v_s5*v_s6
+
+variable mytemp equal temp
+variable mypress equal press
+variable mype equal pe/atoms
+fix avt all ave/time ${nevery} ${nrepeat} ${nfreq} v_mytemp ave running
+fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} v_mypress ave running
+fix avpe all ave/time ${nevery} ${nrepeat} ${nfreq} v_mype ave running
+fix avs all ave/time ${nevery} ${nrepeat} ${nfreq} v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 ave running
+fix avssq all ave/time ${nevery} ${nrepeat} ${nfreq} &
+v_s11 v_s22 v_s33 v_s44 v_s55 v_s66 &
+v_s12 v_s13 v_s14 v_s15 v_s16 &
+v_s23 v_s24 v_s25 v_s26 &
+v_s34 v_s35 v_s36 &
+v_s45 v_s46 &
+v_s56 &
+ave running
+
+# bar to GPa
+variable pconv equal 1.0e5/1.0e9 
+variable cunits index GPa
+# metal unit constants from LAMMPS
+# force->nktv2p = 1.6021765e6;
+# force->boltz = 8.617343e-5;
+variable boltz equal 8.617343e-5
+variable nktv2p equal 1.6021765e6
+variable vkt equal vol/(${boltz}*${temp})/${nktv2p}
+variable ffac equal ${pconv}*${vkt}
+
+variable F11 equal -(f_avssq[1]-f_avs[1]*f_avs[1])*${ffac}
+variable F22 equal -(f_avssq[2]-f_avs[2]*f_avs[2])*${ffac}
+variable F33 equal -(f_avssq[3]-f_avs[3]*f_avs[3])*${ffac}
+variable F44 equal -(f_avssq[4]-f_avs[4]*f_avs[4])*${ffac}
+variable F55 equal -(f_avssq[5]-f_avs[5]*f_avs[5])*${ffac}
+variable F66 equal -(f_avssq[6]-f_avs[6]*f_avs[6])*${ffac}
+
+variable F12 equal -(f_avssq[7]-f_avs[1]*f_avs[2])*${ffac}
+variable F13 equal -(f_avssq[8]-f_avs[1]*f_avs[3])*${ffac}
+variable F14 equal -(f_avssq[9]-f_avs[1]*f_avs[4])*${ffac}
+variable F15 equal -(f_avssq[10]-f_avs[1]*f_avs[5])*${ffac}
+variable F16 equal -(f_avssq[11]-f_avs[1]*f_avs[6])*${ffac}
+
+variable F23 equal -(f_avssq[12]-f_avs[2]*f_avs[3])*${ffac}
+variable F24 equal -(f_avssq[13]-f_avs[2]*f_avs[4])*${ffac}
+variable F25 equal -(f_avssq[14]-f_avs[2]*f_avs[5])*${ffac}
+variable F26 equal -(f_avssq[15]-f_avs[2]*f_avs[6])*${ffac}
+
+variable F34 equal -(f_avssq[16]-f_avs[3]*f_avs[4])*${ffac}
+variable F35 equal -(f_avssq[17]-f_avs[3]*f_avs[5])*${ffac}
+variable F36 equal -(f_avssq[18]-f_avs[3]*f_avs[6])*${ffac}
+
+variable F45 equal -(f_avssq[19]-f_avs[4]*f_avs[5])*${ffac}
+variable F46 equal -(f_avssq[20]-f_avs[4]*f_avs[6])*${ffac}
+
+variable F56 equal -(f_avssq[21]-f_avs[5]*f_avs[6])*${ffac}
+
+# Born term
+
+compute     virial all pressure NULL virial
+compute     born all born/matrix numdiff ${delta} virial
+fix 	    avborn all ave/time ${neveryborn} ${nrepeatborn} ${nfreq} c_born[*] ave running
+
+variable bfac equal ${pconv}*${nktv2p}/vol
+variable B vector f_avborn*${bfac}
+
+# Kinetic term
+
+variable kfac equal ${pconv}*${nktv2p}*atoms*${boltz}*${temp}/vol
+variable K11 equal 4.0*${kfac}
+variable K22 equal 4.0*${kfac}
+variable K33 equal 4.0*${kfac}
+variable K44 equal 2.0*${kfac}
+variable K55 equal 2.0*${kfac}
+variable K66 equal 2.0*${kfac}
+
+# Add F, K, and B together
+
+variable C11 equal v_F11+v_B[1]+v_K11
+variable C22 equal v_F22+v_B[2]+v_K22
+variable C33 equal v_F33+v_B[3]+v_K33
+variable C44 equal v_F44+v_B[4]+v_K44
+variable C55 equal v_F55+v_B[5]+v_K55
+variable C66 equal v_F66+v_B[6]+v_K66
+
+variable C12 equal v_F12+v_B[7] 
+variable C13 equal v_F13+v_B[8] 
+variable C14 equal v_F14+v_B[9] 
+variable C15 equal v_F15+v_B[10]
+variable C16 equal v_F16+v_B[11]
+
+variable C23 equal v_F23+v_B[12] 
+variable C24 equal v_F24+v_B[13] 
+variable C25 equal v_F25+v_B[14]
+variable C26 equal v_F26+v_B[15]
+
+variable C34 equal v_F34+v_B[16]
+variable C35 equal v_F35+v_B[17]
+variable C36 equal v_F36+v_B[18]
+
+variable C45 equal v_F45+v_B[19]
+variable C46 equal v_F46+v_B[20]
+
+variable C56 equal v_F56+v_B[21]
+
+thermo		${nthermo}
+thermo_style custom step temp pe press density f_avt f_avp f_avpe v_F11 v_F22 v_F33 v_F44 v_F55 v_F66 v_F12 v_F13 v_F23 v_B[1] v_B[2] v_B[3] v_B[4] v_B[5] v_B[6] v_B[7] v_B[8] v_B[12]
+thermo_modify norm no

examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in

@@ -0,0 +1,21 @@
+# 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 MD
+
+timestep ${timestep}
+fix 4 all nve
+if "${thermostat} == 1" then &
+   "fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed}"
+
+

examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in

@@ -0,0 +1,26 @@
+# this generates a 2-atom triclinic cell
+# due to rotation on to x-axis,
+# elastic constant analysis is not working yet
+
+# unit lattice vectors are
+# a1 = (1 0 0)
+# a2 = (1/2 sqrt3/2 0)
+# a3 = (1/2 1/(2sqrt3) sqrt2/sqrt3)
+
+variable        a1x equal 1
+variable 	a2x equal 1/2
+variable 	a2y equal sqrt(3)/2
+variable 	a3x equal 1/2
+variable 	a3y equal 1/(2*sqrt(3))
+variable	a3z equal sqrt(2/3)
+variable 	l equal $a/sqrt(2)
+
+lattice         custom ${l}             &
+                a1 ${a1x} 0 0           &
+                a2 ${a2x} ${a2y} 0.0    &
+                a3 ${a3x} ${a3y} ${a3z} &
+                basis 0 0 0             &
+                basis 0.25 0.25 0.25    &
+		spacing 1 1 1
+
+region		box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y}