From fb702fa9d6587d92dc4783738e742028df939c03 Mon Sep 17 00:00:00 2001 From: Aidan Thompson <athomps@sandia.gov> Date: Sun, 25 Aug 2019 23:02:41 -0600 Subject: [PATCH] Added FCC, BCC, and ICOS examples --- doc/src/compute_orientorder_atom.txt | 4 +- examples/steinhardt/in.bcc | 74 +++++++++++++++---- examples/steinhardt/in.fcc | 69 ++++++++++++++--- examples/steinhardt/in.icos | 106 +++++++++++++++++++++++++++ 4 files changed, 227 insertions(+), 26 deletions(-) create mode 100644 examples/steinhardt/in.icos diff --git a/doc/src/compute_orientorder_atom.txt b/doc/src/compute_orientorder_atom.txt index 2e83f3ea07..f4e1cbb924 100644 --- a/doc/src/compute_orientorder_atom.txt +++ b/doc/src/compute_orientorder_atom.txt @@ -19,6 +19,8 @@ keyword = {cutoff} or {nnn} or {degrees} or {components} {cutoff} value = distance cutoff {nnn} value = number of nearest neighbors {degrees} values = nlvalues, l1, l2,... + {wl} value = yes or no + {wl/hat} value = yes or no {components} value = ldegree :pre :ule @@ -63,7 +65,7 @@ specified distance cutoff are used. The optional keyword {degrees} defines the list of order parameters to be computed. The first argument {nlvalues} is the number of order -parameters. This is followed by that number of integers giving the +parameters. This is followed by that number of non-negative integers giving the degree of each order parameter. Because {Q}2 and all odd-degree order parameters are zero for atoms in cubic crystals (see "Steinhardt"_#Steinhardt), the default order parameters are {Q}4, diff --git a/examples/steinhardt/in.bcc b/examples/steinhardt/in.bcc index acdfca3647..a0348fdcd9 100644 --- a/examples/steinhardt/in.bcc +++ b/examples/steinhardt/in.bcc @@ -22,26 +22,70 @@ mass 1 1.0 pair_style lj/cut ${rcut} pair_coeff * * 1.0 1.0 ${rcut} -# initial velocities +# 14 neighbors, perfect crystal -velocity all create 5.0 482748 -fix 1 all nve +compute qlwlhat all orientorder/atom degrees 6 2 4 6 8 10 12 nnn 14 wl/hat yes +compute avql all reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] c_qlwlhat[6] +compute avwlhat all reduce ave c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10] c_qlwlhat[11] c_qlwlhat[12] -# 8 neighbors, perfect crystal - -compute qn all orientorder/atom degrees 5 4 6 8 10 12 nnn 8 -compute avqn all reduce ave c_qn[*] - -thermo_style custom step temp epair etotal c_avqn[*] +thermo_style custom step temp epair etotal c_avql[*] c_avwlhat[*] run 0 -# 14 neighbors, perfect crystal dynamically melting +# check Q_l values -uncompute qn -compute qn all orientorder/atom degrees 5 4 6 8 10 12 nnn 14 +print " " +print "*******************************************************************" +print " " +print "Comparison with reference values of Q_l " +print " [Table I in W. Mickel, S. C. Kapfer," +print " G. E. Schroeder-Turkand, K. Mecke, " +print " J. Chem. Phys. 138, 044501 (2013).]" +print " " -timestep 0.003 -thermo 1 +variable q2ref equal 0.0 +variable q4ref equal 0.036 +variable q6ref equal 0.511 +variable q8ref equal 0.429 +variable q10ref equal 0.195 +variable q12ref equal 0.405 -run 20 +variable q2 equal c_avql[1] +variable q4 equal c_avql[2] +variable q6 equal c_avql[3] +variable q8 equal c_avql[4] +variable q10 equal c_avql[5] +variable q12 equal c_avql[6] + +print "q2 = $(v_q2:%10.6f) delta = $(v_q2-v_q2ref:%10.4f)" +print "q4 = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)" +print "q6 = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)" +print "q8 = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)" +print "q10 = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)" +print "q12 = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)" + +# check W_l_hat values + +print " " +print "Comparison with reference values of W_l_hat" +print " [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, " +print " Phys. Rev. B 28, 784 (1983).]" +print " " + +variable w4hatref equal 0.159317 +variable w6hatref equal 0.013161 +variable w8hatref equal -0.058455 +variable w10hatref equal -0.090130 + +variable w4hat equal c_avwlhat[2] +variable w6hat equal c_avwlhat[3] +variable w8hat equal c_avwlhat[4] +variable w10hat equal c_avwlhat[5] + +print "w4hat = $(v_w4hat:%10.6f) delta = $(v_w4hat-v_w4hatref:%10.6f)" +print "w6hat = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)" +print "w8hat = $(v_w8hat:%10.6f) delta = $(v_w8hat-v_w8hatref:%10.6f)" +print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)" +print " " +print "*******************************************************************" +print " " diff --git a/examples/steinhardt/in.fcc b/examples/steinhardt/in.fcc index 6c6a919d7b..6d2775d0bb 100644 --- a/examples/steinhardt/in.fcc +++ b/examples/steinhardt/in.fcc @@ -22,18 +22,67 @@ mass 1 1.0 pair_style lj/cut ${rcut} pair_coeff * * 1.0 1.0 ${rcut} -# initial velocities +# 12 neighbors, perfect crystal -velocity all create 5.0 482748 -fix 1 all nve +compute qlwlhat all orientorder/atom wl/hat yes +compute avql all reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] +compute avwlhat all reduce ave c_qlwlhat[6] c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10] -# 12 neighbors, perfect crystal dynamically melting +thermo_style custom step temp epair etotal c_avql[*] c_avwlhat[*] -compute qn all orientorder/atom # degrees 5 4 6 8 10 12 nnn 12 -compute avqn all reduce ave c_qn[*] +run 0 -timestep 0.003 -thermo_style custom step temp epair etotal c_avqn[*] -thermo 1 +# check Q_l values -run 20 +print " " +print "*******************************************************************" +print " " +print "Comparison with reference values of Q_l " +print " [Table I in W. Mickel, S. C. Kapfer," +print " G. E. Schroeder-Turkand, K. Mecke, " +print " J. Chem. Phys. 138, 044501 (2013).]" +print " " + +variable q4ref equal 0.190 +variable q6ref equal 0.575 +variable q8ref equal 0.404 +variable q10ref equal 0.013 +variable q12ref equal 0.600 + +variable q4 equal c_avql[1] +variable q6 equal c_avql[2] +variable q8 equal c_avql[3] +variable q10 equal c_avql[4] +variable q12 equal c_avql[5] + +print "q4 = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)" +print "q6 = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)" +print "q8 = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)" +print "q10 = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)" +print "q12 = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)" + +# check W_l_hat values + +print " " +print "Comparison with reference values of W_l_hat" +print " [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, " +print " Phys. Rev. B 28, 784 (1983).]" +print " " + +variable w4hatref equal -0.159316 +variable w6hatref equal -0.013161 +variable w8hatref equal 0.058454 +variable w10hatref equal -0.090128 + +variable w4hat equal c_avwlhat[1] +variable w6hat equal c_avwlhat[2] +variable w8hat equal c_avwlhat[3] +variable w10hat equal c_avwlhat[4] + +print "w4hat = $(v_w4hat:%10.6f) delta = $(v_w4hat-v_w4hatref:%10.6f)" +print "w6hat = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)" +print "w8hat = $(v_w8hat:%10.6f) delta = $(v_w8hat-v_w8hatref:%10.6f)" +print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)" +print " " +print "*******************************************************************" +print " " diff --git a/examples/steinhardt/in.icos b/examples/steinhardt/in.icos new file mode 100644 index 0000000000..d0d61a902d --- /dev/null +++ b/examples/steinhardt/in.icos @@ -0,0 +1,106 @@ +# Steinhardt-Nelson bond orientational order parameters for icosahedral cluster +# W_6_hat is sensitive to icosohedral order + +variable rcut equal 1.2 # a bit bigger than LJ Rmin +variable rcutred equal 0.75 # a bit bigger than 1/sqrt(2) + +# create a perfect fcc crystallite + +atom_style atomic +boundary s s s +lattice fcc 1.0 # neighbors at LJ Rmin +region box block 0 2 0 2 0 2 +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +region centralatom sphere 1 1 1 0.0 side in +group centralatom region centralatom + +region mysphere sphere 1 1 1 ${rcutred} side out +delete_atoms region mysphere + +# LJ potential + +pair_style lj/cut 100.0 +pair_coeff * * 1.0 1.0 100.0 + +# define output for central atom + +compute qlwlhat all orientorder/atom wl/hat yes cutoff ${rcut} nnn NULL +compute avql centralatom reduce ave c_qlwlhat[1] c_qlwlhat[2] c_qlwlhat[3] c_qlwlhat[4] c_qlwlhat[5] +compute avwlhat centralatom reduce ave c_qlwlhat[6] c_qlwlhat[7] c_qlwlhat[8] c_qlwlhat[9] c_qlwlhat[10] +variable q6 equal c_avql[2] +variable w6hat equal c_avwlhat[2] + +compute mype all pe/atom +compute centralatompe centralatom reduce ave c_mype + +# gently equilibrate the crystallite + +velocity all create 0.001 482748 +fix 1 all nve +neighbor 0.3 bin +neigh_modify every 1 check no delay 0 +timestep 0.003 +thermo_style custom step temp epair etotal c_centralatompe v_q6 v_w6hat +thermo 10 + +run 10 + +# quench to icosehedral cluster + +minimize 1.0e-10 1.0e-6 100 1000 + +# check Q_l values + +print " " +print "*******************************************************************" +print " " +print "Comparison with reference values of Q_l " +print " [Table I in W. Mickel, S. C. Kapfer," +print " G. E. Schroeder-Turkand, K. Mecke, " +print " J. Chem. Phys. 138, 044501 (2013).]" +print " " + +variable q4ref equal 0.0 +variable q6ref equal 0.663 +variable q8ref equal 0.0 +variable q10ref equal 0.363 +variable q12ref equal 0.585 + +variable q4 equal c_avql[1] +variable q6 equal c_avql[2] +variable q8 equal c_avql[3] +variable q10 equal c_avql[4] +variable q12 equal c_avql[5] + +print "q4 = $(v_q4:%10.6f) delta = $(v_q4-v_q4ref:%10.4f)" +print "q6 = $(v_q6:%10.6f) delta = $(v_q6-v_q6ref:%10.4f)" +print "q8 = $(v_q8:%10.6f) delta = $(v_q8-v_q8ref:%10.4f)" +print "q10 = $(v_q10:%10.6f) delta = $(v_q10-v_q10ref:%10.4f)" +print "q12 = $(v_q12:%10.6f) delta = $(v_q12-v_q12ref:%10.4f)" + +# check W_l_hat values + +print " " +print "Comparison with reference values of W_l_hat" +print " [Table I in P. Steinhardt, D. Nelson, and M. Ronchetti, " +print " Phys. Rev. B 28, 784 (1983).]" +print " " + +variable w6hatref equal -0.169754 +variable w10hatref equal -0.093967 + +variable w4hat equal c_avwlhat[1] +variable w6hat equal c_avwlhat[2] +variable w8hat equal c_avwlhat[3] +variable w10hat equal c_avwlhat[4] +variable w12hat equal c_avwlhat[5] + +print "w6hat = $(v_w6hat:%10.6f) delta = $(v_w6hat-v_w6hatref:%10.6f)" +print "w10hat = $(v_w10hat:%10.6f) delta = $(v_w10hat-v_w10hatref:%10.6f)" +print " " +print "*******************************************************************" +print " " +