From 3bc7350704d2ecb09b2f6a955fd9b6acd3ca1580 Mon Sep 17 00:00:00 2001 From: athomps Date: Thu, 5 Nov 2015 06:53:08 +0000 Subject: [PATCH] Added hexatic bond orientational order parameter git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14236 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Eqs/hexorder.jpg | Bin 6808 -> 6565 bytes doc/Eqs/hexorder.tex | 2 +- doc/compute_hexorder_atom.txt | 43 ++++++++++++++++++++-------------- src/compute_hexorder_atom.cpp | 37 +++++++++++++++++++++++++---- src/compute_hexorder_atom.h | 2 ++ 5 files changed, 61 insertions(+), 23 deletions(-) diff --git a/doc/Eqs/hexorder.jpg b/doc/Eqs/hexorder.jpg index 466d6b80902cbd6323391d801b4d6389300d0d23..ba4938b8cd83d8ad23a3704ec7274f58c56e6b30 100644 GIT binary patch delta 5305 zcmai%dpy%`{Qg%uIE0cMLph`zQw|Z9M6p=TF|4G}XrVEzz00W_!W=^=XF1GiCZ{E* zN?OjJO->_GVI#BH;k!S7kKgx??;qFS_v60q>wdnj*Gq%4^Y!g}MF(CBfcAkvpgnv+ zGCuCj6o7yRn&1MIJ%B(OYk1cJeUuqYPp0xf13Zu{`cAipmY~-Huf`kO(J@kGYOw z3eG8WZUsQMrMWcPZ=k4Ay7xB^bet7_c>>v-x*tKh)!J1LBrd^p&wcRXLfSGESCa}# z{rt%mp5g2x5UzvLt}#FXj{tvjHBC^)KOM!<8zih4^cvX`a|!K%COmz#Q6q}f zy0&pu8pvf;xWo_qZ$@UYq_9_7IoFek;^A(i2?2+tLgdyOQmyx03G3$+jk@+ZldE z50*;QIC|;bm=D1=a|$q;H2i+H!iYW(^ulcixVxQj^lSDH{l2{e_U@-5lp9)*8iv{O z{&^RDa~EvXZeB@wGyUYDTSWxYpOkM$*ckTbs;L(TTe%R8-Vqkt|BUBtM3Sy5zO`eH zIan!QPCWg3@BfnKB3yh{JMARjXYvag}TEt!}_Zl{+Bx! zrIPdqwb>JC4SBwCXWzEF6kcv1sn<+Y)|b~N8(1^Cx=2K&LOF*L8D7{Cvo?0t{{ zB%%a=WKUfm)$N@;s^#x(Equ_*MU|sD`!?@BFjV`=(SF4&6~qHYXR*Cg=ShFE7J@Wp zKg=HPfK{Uhd#L3v$|z{)hnyWY8KqsCAooRUJ`GlK{^`4H0<^66oEWYY z`z&MjY1bK+dOutYbMj6(!f@sc-MgX$CiKQZL?mkPrPRz%weUTXDPd{Xj!#k^IXguI zu|JKNH6*$s0(&m}U!aTiPZ>2Pk_XBNnWZNZgOHE46Bs%@@Nshm#Kf#HK}1x&cUunC z8)>z=)$oqc)UH<_@BZdc`L&75H6LXp_Q=ZO`vI1E4*M9GN#xXmA1Enezk&(UeW~$) zfVmFK+(xBg2o#DB3c_2m7xwX{agP|Hbyzc=t)F)O|uF05w@ z;jzs~Y7|pZddF>GHrmI`24<5iKNRv|nmvmP-W^Qxk-OF}gfsov@=mWZ$H2O{%IRhx zyWJ)d1+{LgmU{i$Vl+NuIt-XbHZAXfoqk#F8gl47I5(!rJ3FC-c49=6rE_S4X^``G>YJsaOE9mBeHr|MMB9c|vDV*M)J^G@M9SFnr+ z+7F}iK9|>#km`;ZDFCq46Cy41}$@A}g*UL0?yHc}O~a z$NHWBjNRfXvo8*Zl9u|7j`wK=3LO&K)~K#n&i1eyTBN!Z)HYL-X9L>N#m2KYNtK_-r}~V092J>Nk<1%Z@NBPEz^7 zBP}IWUY>);T_@&0SIS0sz1)|uX&cH+x2JX`1VnjJ5XDVl6T=>z@`8OIbkD%G}H%aynj(rVq+^bRK_KX_yDOe--$ zdEb+p5}HeRYF7}1D_pdzh%V~pfkaJ3fHHK;z(Dq_nCV@_<8VJ_s`yI(PTyA5i8h-S zh-4mOF(PQ`+#vB*nxq_7m6MM>i7w#^C*k1;ph!$V2r@TOLT3)X&)&`<8rGI{mz?d$ zcWU0182m|??)x+{?|hc?3TCMyD-xGLGfH>d)nQ`2$%W@OlO7l%P(wj31vaHcs5Tym z_|Id@267dl-)5$Rxy|>xWEsvU!aea6lEP}o-FdhI`e^rt`+sP-byIy2q=Pu-3$#2@ ztKI7B!O(t14Z4@)`donpq+jxNJib9{nqAGJtw-}f%|%IU<6V1(? z*U}1f?D>$x@cv#lUY10)P7d-=Q$D>W7xAR|`$+mCF@+_pcFS}B?WWozF<6z z?K6TsjgGOfgF+u$_HW<#?dS=3m(|2GR#6B7W zXB55XfvAIhGurf-j5%@wxVa)8ZwfABVW)#A811sY@}YnVcd}Z;*-s(sn%_P3?HSK9 zty7hI8}-LCD(qccf9cxsKtDzZSp@7z&^DF~yI2(@{4jID)^YgLeV~~K%1LOSMNaom z*O(vW%Q$!oG>=p-2D92N2-`{7vL8~W| z1&rXkH^)rF&nxX|C`<3W{(_6jx}9mYRqapo`(?`7|fAYri%c z@LM_eYqYa;FLmVe(uS_dh}9vVxO=%}@q)hLCocnOhP@t{)s%pP{-xWlqpn-F)Gp#W z$xLLT!iuq_8$$T=fG)s7vj7k*+x{cYQJdm)iBi*Y!shYYxBi20fnqJ${rkkK+ zAU~*mYJ#HW=OnSbrsIbze&I5+M>#{Z4F z_EB6#)K;WD_gGYIXI8ySGB!#y*&A0UY43sQxfu% z$M!$>uD8GIDtAH>PuKK^%pUid^kbb_*vInR!#Ka3Q(d~ZFBE+XP{)1bEPb`Y%JS3W zfiYvL^FRk?(eg}xmJn@PAYw!RytSnR+`!O zw5fAb5fUAL@~=!a#xCr?e|z7r5b{mCZU|d1V`;T5qpR~+d*&}&y`S-P$KcHjNRC}W z-Uw&bdbo)uKIO)h`AypfdXSXSIa2rfQFOUcJzNAGhPYT6s67G`HJJ#Hkb3y$(SO=? zL9scwT;b2x&m2nE8C?Y$@T=F02>2c|KUl`)jGjrZJj1YdUfXueD`_w^qQa}QW<4n6 zh$%~?f6BAiXx}kM!Ag`a#bt{Jn##D&8p0#`UA|$C1kWp|Tfu)WPl+mUL}wnP-jw<3cH(dZ{&6voSibIS>)LJ_bAIpPNt0Mn7LgWo)w-*ty~q|aIgC^rT}NFAdefG1 zy#;7ZIVCOKWjFo~hTgu|5v?k!b4LMquQWyEQ@0u1*BSoaT#3jCf7XGVU)7lOYWv29 z%l?<`Vn(7*P`m3M!} ze-;s-t5@g9|B>gw&&^dZkz@Ti^JH-gV9C}coE)yH&;WxCOt}9Dsl=tYzkIZqLh9!l zj9}%_=iM06sRR)(70YjdOJgYc-#J2}1s&IXv)yGlu(`HEV?iz9lzboB9F-81G)3 zs*dz;r>xCYu9{4p0ES%t`6B(Y&R4{d%zKeHo)=D{Pn&%Y?{O9AIZXcH)7A zMDOsoH^e7W0GMSO&h%m32IgTzPLV$OI2$_o#g2JFWz-&bYh=$x!kSKls=vMQHC6DR z@z1u0la)LFD@~vI6i6bsx1SoRGPZ}#WVAm;{@VAl7%r~mM)#l6>i_XjUdq- ziW1~zKiEa}wxuM6&kw%is@*r=Ph=c;(sHDjP{Yj0F5YrPf3X(;t_?T-`3Gs?fm(sN zOgUmN_&L0xFpeS3tn6&d>F0r@mNNQtB13I-DHt#^>@X#Onf4KilQaytUv@ym_8x9i zHD%~;AoIqwjr_EwWUA1#{JFh1cA_I)>@!N7X1q6C1O6hjNO|`sDb?HW^s40Qx4m=h z49ZrcybPlB*D`=6XjW$)R)q)J_ivmnTZK88u0$D_|M}cp2?K5%*HQjPY6(=Uf5%X@ zYq3GWl@~eqz|Bv364P;cLFGQeA8TfQT1nc>a>aNcF#Jlp>QglZiEpQi6LGIdJ*hE@ zaxAM#bo#=4e?BLRZ~wjOQ^cyCW4+~67{%pE(|~{;ka8gNbvpFS4I35HI=QRi4+_Ri zvK|Vm>{E=-PBgR;ZJ9!bMLbL7fqdbD{6{Y@;A2VV6U_|bG@_%8uX(hxHCR5p`e3!V zHD127e$y1EQMR~n@!`#qH17-3CP7!{Nb#uDWI1Fi>V+QC)TJb#DZm7^xM_1@v~~#< z5!V4+@4ry!hUyjbI3)91WrZ20mh(rEMk5PPC9FHp2)ZC7o>qingbt1V8_;ZpSf zx7&&j1#}6a;WRh)aoM>dL6(ulb2N+@Np~eOlvne@TU{971eCJZk4@QxwAb15d*&|) z9#-CC$+Y3T3B&~OI|uNBwaEbYngsFb%>0=chy|ggCS=CEaHTLipmy4%U9<%U+3in5 zai~wcvY^Xf6X?}FigFC@I~qwG9c9(PKcJWx<2L@fG3Da4{YBK>r#hzuFx4{Qf--4u zQuOkA=oNoD>67n^WVZ0n3x!y~T&K#wjLP%wx2$xsRSOW~GImEz$YyPLU*le%`42FPuenj+@?dJtPqg*M3aPC>Ih zCEhgjUjly0;^Tq4BuO+6WWjj@4(Jr=10*^Sv=|nKrHMmjD3?9`F6gYYMR~QcWZ4wH<$nPX73cmbedlUQBO7NDs91HU> zlm|Mytht=H);Hp;Iaa1nOx-O=i9C-w(+_=XVI}Ei{Mx}v67+v)&rg6qWqw9uPhc!q z;HYx$7e%9bR>Z=2R6wv_&8&Y_4NTc8_;`kec6HkWn_`=jS~YdWmp@yAt{U%3aX&j@ zKXXMGSJ`SaCwZXQ_I}{7DI|inE{TBLWu8Xr))}6FM1`ux%cF0_&j>CAC`4Ztx~H1Q zG&$3^F*E(>pa1F#8mI`W!?%JXCs*)f7r+ZXbL^mO1_2E z-b`n{2rBDFJ^cK3)JCCAxg*ZOB1GuqCGy_3ozvx0QE@erQVv!pgijYg`F{oGQT_*C CmD#@l delta 5550 zcmbVQXHe7YvyOt&jv`1$Kq;c2qSB7YAu5ItIT!)yB`PA#5J9>zKLt?;y+%Qaf`~{7 zEd;4jqaYwPQbHgg9g-l#Xxim}XYR*4cjkV0_S??RzWeO6?{1OdbK}a!-TOLUi-1Hy zAkeNI*kyj+oh4!i>a0PCjPC+K1@(j!fhjBX5l?}H5JB&6v7p?+f3dZ&4N%h*zX;WI zt&$D7s$^UzA*n$xix>-C95y0F9j?X> zU(28P`)f*P&#dOPp?x3LbHpFTNuIUO#8UNn6}W-WbC7O5_?`8iH2{#|I~-~UIis${ zy(9#+U20Jl^r`GwTvfaJtWsQC>0xZCI2{UXNf%W`Z1c_ddAtWD=HUA1aK_YA#$C3z}N4ZDNzY zkXYK-*2a?Fp`Tt&DSASXMc{f(^*Cw54kGvISOFGvtN{3r`QIUJi4fa)m3qWxlDjll z4ecld6*u;9MdonN%!BEKT{iyQe$gXUcxhP0#9K7OAyAE4s&GyTWdwri zbzpzD3M8HhL82=h{El5>W_d}CVvHjnkcA*XT9$il)b>fsZls~3qH{*=rhJR-73Vdu?ESh!*LDjSlfOA-pCr^j( zWJAJr7(4@HTkinL_CZWC*4%Tp^dC3WzGIo<-GR^XH%U=Ws@Z!ityms64rrNO{ONc- z|K%juxEvo?-V#yy?xVk1qI*yi&{)AHjH8xKcBC1S>Uf6*Vf|oK3!;*7N{tG(%3TRfQ2A26B^Gm%BPhli(_AJ zBEgWr9ojhCrRdtyGLvhqw_!pP!+scPlznd#bV|fsG;G}r4^|SOIw=zhzx!3CXeDf8sLq#C$MR$eCm(@Kc)9m@7Y_ zRl`k|EQ`R>abh5}rP_6EtF$!lPR8Kv)PcyaA$SvW$ zq_P_PvLeSP%2_|3+$&oHJJZb@^ZKwG*^64#byriYPv7MH)QK*W%|Bg&^!bPWD~& z!C5$|;$u{atA`7`5O#~9He%|8lX=VRQw6$vZbf!vlKjr&z><+9Pfi^7qF|)8+AD$J zg*K_qoDqUL_zrYb^18aoP6?t7WG2N!UueN=TC{Sy7Z$G}@updy+a4jynh0xTjd+Y6 zrbk0$n~_w>kLh5)<`)&Qs0(=2`WV}D-mRmd-v@6yYAS179=;>_Nzu9YuSQ8=iH+eV zMw{OLMPfJ_KO}^IfVqe5GeuWDfaXu?f4RF>(o^VVt798jlV4Yj+CD$Km?Yh)v3pr% zb=j;tj!*Dcn;G7AQIJ|ay`ls)zw%WpZMW%Bm_u8es8b&Xy~FhD7%(rK<0}O1k>}?K z4qrI&w}}DeVJ|`ioVRDm8}%ot@v*JrPn4`G^3qQBa5!qn+1+gImQd|96KM7FRpRzIu_0|3-nnn__l{)W*- zLs+;cp+8!O(RMR7&kIJ4`Yy@HEl7;!j&MEH5?LFZewq7O%Z?5~Pq>g^^gnM>BjcA*MI(MbVCeIfbqhU$e zlgl@khCLly2_LCZ<9>O*!?kvh-I5%XKuU^h!_AwsUlhA{+9c9Bm0h0wlGX4Nh>urFB zl;CirgAM2L>uFAE$Hco*(IeWmbvff=lin#-fyaGqKEHdaERps)ShI7%6OQQ?-%R*`hsHja8N8;_O2X9VmRu0eQ+`&%YG|iO#one%y zVzQ9lYf90m(MV}JhvbG~`X*4kgNxn21g1<(Ba`b?&%w`g3KkJ@M^f~?-ys4lLN_WA zu*jeS*-F=WDWkiRmhl_dh51;(=L%oLtkf@6@>s2KJ`Qb66#!D~;D?l^ zDSVzk)s(QbQK*|;+)wvg$99wZ?@4x%r_X6Uro3N9J_If`SS~nbElmc{9?FsMaV^&gIxP zF(_##;-QM_V>#lZM;tF-kax7&{rKB1$A>3uiU0l-f)c^1#*C-anla>dAxIK)BywTe z-RBiyOSv|*E-vxxa~DK_$tG%-<)I@UX0r)Xic%N(#pcdHSzwM6u637%!Eze)YTB2m zu=}$edkShmK53A$^J~{4)Z+Ud-FrKu3`W8y!EHiNWpNvnh3}%Rs}2xp7YqFwOmLw3 zQFCa?lu2)()G?U)?aj)=yDxwY&q%&aRp+?x2+t8#Ghw}W1{{MRlG~H@KUrvT@Uy)% znI-V|7C;ZzdGyQ4(>|n;FO{UH65)$`C4wA6(p>ghg=t_sN*`jj6^}Mwcbr%cS~~h$ z(Wq{C1CmYEU@9rem=299@YB|t@R;MX2vHBuRk!})h*LH6RM{0VOfkLL?S8}JwWhQT zH<;>c?I+jX|Lw+hTQ0R5`EfrY8)qo+55w8g9(pmHZi;uG@<@13-7IeKJh;Wp=~?#B z-=xs->&h&B_Rx+77NhN$7nx!}H&_WR$-h234#*-bX=?rT&312j;*}F=sp?bHGn?-$qw5Q;4R|h%im;ML9+0Ta=N6)-l(A zF(HM|d-Xn7TXmWTmHwQfRr+>RP}c&v|B?jB`+z5;?_Hx<@5{2 zcyM$B=lj4*2=aji4TZJdL43Db?FMR03TFo<^=(JWy&G8Ne7!eLKG2$s4ckPqL^&WujW ze#NLkq9>!nsY|ua{^}BQ2+mPx8$8}ZiP@h@2P8b3pT*s+sQs+453sAtbf{^3q1!|p z)xpjpuE#9+jrr4?{l}tA&EJeauOqL&;+i6M#=pe=Y$j&$>bEBEnRP32Xa4^BS0>|b z0Zy`IKppE}8JN4jqt}Ed*2U61&nno93^{H>Wt1}jqYeCr=szp}m;|hhkwQmvw*sN_ z=xvq=`v3>Y_2bmc(z*~6z$XjsA&eI1YL9Ygy>ud>Jjh)C49diM&dg`2?kq}2-Q`!* z`wwMbA00zXDzXtQcz0_OAkMR4HVR}y6M;6j*^w|{-Su4x%<)L}Nmqqqv~>wllB@^+ zd95bwO`W%1d?oFYovnq_%g~wyudBaW3HgEuIhdLD>$>vwkOAG-B7Q&8~~ySN5;Qz%W)(mzQWYNfB$vLR|Kbhae$;n!Jmbvy30KHpmvJuZ zuV{hrA7qA+y%$5LrN+Foi$>&WV(sUvlE??v&#C8xa3NtC?USXn(73D~NY^bb~1U zq(zJg#4#9e+ynDoX<5}(=hAg<^{O^m7iH~eEDV(Kfa2g=H0eKJ=+96V1@4S(Hyw%_3Jh zBYc{K@{Bk? zxVM)v2bsE01!DfrRUB}8gZf7kVg&Gk)~VUp+|ppev`QnZXhLMQ3=epYLpOTcTihXA zFJBVIPU-bU(I7|%EQ&hfk_Vr5-_&mH1d4P~whj?X^EnL~w>WfU0Ie;$3BPZBt8qX^ zC#$cht;7R{@-hBq`QyZ9cM7W$;dkeHxFk%H9Q*qpu9=_>!L%P3@1Wy50kx%R9exrI zHV0s%3szLdnh6%I!RX5FM70*Y2sBU$+rP)*N4wab?-^6x+cmCZC(X~l^$Y-1_@m9) zaY+Wj_k+!2{w|+`oY>fy(t9tO5eahsir522G9-9j912UiOQA)J^P*^}MUt74{owwr ztwYD-lgZihIERVV@G8rL)7HOa+F8F1>MpxjEy?6PPc^6w`ZKxV*4jPx+adv{Yv6mk zqM{|bh=X7RmRhq!-T-Y*)hhn=x^Ox4948^3XC zCFsj77`Z$=dlE=UDs&6ZS#GYQ72~`U_Z{m`(oxHWdv1H8xp8uKE~-1`!k+&FhZNV9 zxloQ68^P46Xyigy+3KpyjnWTpei9K6RgJ4YHf34<6;fRyQ)qbeE!rwEC)`n#(x<_^ z3YG#q1w$A${u{w*zlm<*Y`ow+rzez2h&THDTCwD{e_iEr?!zCi8=ksgt&tlQx6WQ( z>J3p1-1~)VYB!7|J}@PXmd`YhGF1y+h!ozK_koM0>G@p4?tTQC*)@^n_%d%L_qVCy zuj%vbaFX2@D&$l$_#@1H`1<;^zSpyUFtZ)-`bP*VB=QS|Af{>?jn_XS1ZjYMcNS+D zS1z)6cA#zGAcSSp<#n>EW$*{ZLObBwvVkuu0(HXb?$RGP?p%BW^4aIa$de5>E({FH zT_1@lm2k-EbNwTKe9 g^3t!(A*C0x{B7rW7ki-qdkpmdgw+2yiU`O53){JS!vFvP diff --git a/doc/Eqs/hexorder.tex b/doc/Eqs/hexorder.tex index fd7dac717b..feb4e313b8 100644 --- a/doc/Eqs/hexorder.tex +++ b/doc/Eqs/hexorder.tex @@ -2,7 +2,7 @@ \begin{document} $$ - q_6 = \frac{1}{6}\sum_{j = 1}^{6} e^{6 i \theta({\bf r}_{ij})} + q_n = \frac{1}{n}\sum_{j = 1}^{n} e^{n i \theta({\bf r}_{ij})} $$ \end{document} diff --git a/doc/compute_hexorder_atom.txt b/doc/compute_hexorder_atom.txt index dd6d1d98a1..bdd79dc7a0 100644 --- a/doc/compute_hexorder_atom.txt +++ b/doc/compute_hexorder_atom.txt @@ -10,26 +10,31 @@ compute hexorder/atom command :h3 [Syntax:] -compute ID group-ID hexorder/atom :pre +compute ID group-ID hexorder/atom keyword values ... :pre -ID, group-ID are documented in "compute"_compute.html command -hexorder/atom = style name of this compute command +ID, group-ID are documented in "compute"_compute.html command :ulb,l +hexorder/atom = style name of this compute command :l +zero or more keyword/value pairs may be appended :l +keyword = {degree} :l + {n} value = degree of order parameter :pre +:ule [Examples:] -compute 1 all hexorder/atom :pre +compute 1 all hexorder/atom +compute 1 all hexorder/atom n 4 :pre [Description:] -Define a computation that calculates {q}6 the hexatic bond-orientational -order parameter for each atom in a group. This order +Define a computation that calculates {qn} the bond-orientational +order parameter for each atom in a group. The hexatic ({n} = 6) order parameter was introduced by "Nelson and Halperin"_#Nelson as a way to detect -hexagonal symmetry in two-dimensional systems. For each atom, {q}6 +hexagonal symmetry in two-dimensional systems. For each atom, {qn} is a complex number (stored as two real numbers) defined as follows: :c,image(Eqs/hexorder.jpg) -where the sum is over the six nearest neighbors +where the sum is over the {n} nearest neighbors of the central atom. The angle theta is formed by the bond vector rij and the {x} axis. theta is calculated only using the {x} and {y} components, whereas the distance from the @@ -38,16 +43,16 @@ central atom is calculated using all three Neighbor atoms not in the group are included in the order parameter of atoms in the group. -If the neighbors of the central atom lie on a hexagonal lattice, -then |{q}6| = 1. +The optional keyword {n} sets the degree of the order parameter. +The default value is 6. If the neighbors of the central atom +lie on a hexagonal lattice, then |{q}6| = 1. The complex phase of {q}6 depends on the orientation of the lattice relative to the {x} axis. For a liquid in which the atomic neighborhood lacks orientational symmetry, |{q}6| << 1. -The value of all order parameters will be zero for atoms not in the -specified compute group. If the atom does not have 6 neighbors (within -the potential cutoff), then its centro-symmetry parameter is set to -zero. +The value of {qn} will be zero for atoms not in the +specified compute group. If the atom has less than {n} neighbors (within +the potential cutoff), then {qn} is set to zero. The neighbor list needed to compute this quantity is constructed each time the calculation is performed (i.e. each time a snapshot of atoms @@ -70,7 +75,7 @@ the neighbor list. [Output info:] This compute calculates a per-atom array with 2 columns, giving the -real and imaginary parts of {q}6, respectively. +real and imaginary parts of {qn}, respectively. These values can be accessed by any command that uses per-atom values from a compute as input. See "Section_howto @@ -78,8 +83,8 @@ per-atom values from a compute as input. See "Section_howto options. The per-atom array contain pairs of numbers representing the -real and imaginary parts of {q}6, a complex number subject to the -constraint |{q}6| <= 1. +real and imaginary parts of {qn}, a complex number subject to the +constraint |{qn}| <= 1. [Restrictions:] none @@ -87,7 +92,9 @@ constraint |{q}6| <= 1. "compute coord/atom"_compute_coord_atom.html, "compute centro/atom"_compute_centro_atom.html -[Default:] none +[Default:] + +The option default is {n} = 6. :line diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp index 54b018512c..dcecf50d7f 100644 --- a/src/compute_hexorder_atom.cpp +++ b/src/compute_hexorder_atom.cpp @@ -16,6 +16,7 @@ ------------------------------------------------------------------------- */ #include +#include #include #include #include "compute_hexorder_atom.h" @@ -38,10 +39,24 @@ using namespace LAMMPS_NS; ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg != 3) error->all(FLERR,"Illegal compute hexorder/atom command"); + if (narg < 3 ) error->all(FLERR,"Illegal compute hexorder/atom command"); + + nnn = 6; + + // process optional args + + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg],"degree") == 0) { + if (iarg+1 > narg) error->all(FLERR,"Illegal lattice command"); + nnn = force->numeric(FLERR,arg[iarg+1]); + if (nnn < 0) + error->all(FLERR,"Illegal lattice command"); + iarg += 2; + } + } ncol = 2; - peratom_flag = 1; size_peratom_cols = ncol; @@ -50,7 +65,6 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) : maxneigh = 0; distsq = NULL; nearest = NULL; - nnn = 6; } /* ---------------------------------------------------------------------- */ @@ -187,7 +201,7 @@ void ComputeHexOrderAtom::compute_peratom() delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; double u, v; - calc_q6(delx, dely, u, v); + calc_qn(delx, dely, u, v); usum += u; vsum += v; } @@ -197,6 +211,8 @@ void ComputeHexOrderAtom::compute_peratom() } } +// this might be faster than pow(std::complex) on some platforms + inline void ComputeHexOrderAtom::calc_q6(double delx, double dely, double &u, double &v) { double rinv = 1.0/sqrt(delx*delx+dely*dely); double x = delx*rinv; @@ -205,10 +221,23 @@ inline void ComputeHexOrderAtom::calc_q6(double delx, double dely, double &u, do double b1 = y*y; double b2 = b1*b1; double b3 = b2*b1; + + // (x + i y)^6 coeffs: 1, 6, -15, -20, 15, 6, -1 + u = (( a - 15*b1)*a + 15*b2)*a - b3; v = ((6*a - 20*b1)*a + 6*b2)*x*y; } +inline void ComputeHexOrderAtom::calc_qn(double delx, double dely, double &u, double &v) { + double rinv = 1.0/sqrt(delx*delx+dely*dely); + double x = delx*rinv; + double y = dely*rinv; + std::complex z = x + y*1i; + std::complex zn = pow(z,nnn); + u = real(zn); + v = imag(zn); +} + /* ---------------------------------------------------------------------- select2 routine from Numerical Recipes (slightly modified) find k smallest values in array of length n diff --git a/src/compute_hexorder_atom.h b/src/compute_hexorder_atom.h index 38f34cab89..9f64c48f88 100644 --- a/src/compute_hexorder_atom.h +++ b/src/compute_hexorder_atom.h @@ -42,6 +42,8 @@ class ComputeHexOrderAtom : public Compute { double **q6array; void calc_q6(double, double, double&, double&); + void calc_q4(double, double, double&, double&); + void calc_qn(double, double, double&, double&); void select2(int, int, double *, int *); };