From 4f71701e4ef23bc3fb74abfac1e88497fe586e22 Mon Sep 17 00:00:00 2001 From: athomps Date: Thu, 5 Nov 2015 01:33:46 +0000 Subject: [PATCH] Added hexatic bond orientational order parameter git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14234 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Eqs/hexorder.jpg | Bin 8635 -> 6808 bytes doc/Eqs/hexorder.tex | 2 +- doc/compute_hexorder_atom.txt | 17 ++-- src/compute_hexorder_atom.cpp | 149 ++++++++++++++++++++++++++++------ src/compute_hexorder_atom.h | 11 +-- 5 files changed, 135 insertions(+), 44 deletions(-) diff --git a/doc/Eqs/hexorder.jpg b/doc/Eqs/hexorder.jpg index cebb8d9d9e6e37800088e2e37d320a795fa46d03..466d6b80902cbd6323391d801b4d6389300d0d23 100644 GIT binary patch delta 5459 zcmZu#X*kqf*q@?oPl@b8SyEJz?MWU}X)t*-A^S3wWNS=hZ_K|)G4^$eGDVUsV;duD z)+s{P$ufo^J2NAuso7r7dtL9h_rvd;^W{GGbuZ_4pZjjIcy3hHxc4Ba>$MO_7z6_C z*@a!!=e^lNcA&0WxX{!dpaGxCH({nd;>r=>La2kS=H$D$FSh2k0V?XEm#oxXE47?; z6?`D(9wC2I(^Fz~dw|&tpuPmtHv3wwItAN>Pi|CAnrjK+z{f+EMhuBjM{6)6*9#{9 z{+im=JEwkq_`rvaT+v7I;^*wMFce*GC3Z0MJfuelc6Z}vEszQLj)dAl^i?%Dmj$5q z%dLw1ekFQ>AvHniDLyJ{c0;TgM)Mgq%uj)!42n}o8xeHzoJVJ&it5?}Ikz*4Qabci z2a8wl2fFzX7O_Y3#y7R{6!+g>{v@km{yGkxCGVxojCkxT0L5Zamz&tw%1B>0!zw|W zY}l`J;uO056mUHL!n;t-Jwk59iAyKpZl;Bb%R%6CPD)g zR=h37xK?u%utV#}6Qn{_L7?A4=AlZVGd|KBcY0T>k(42gqTnNg1TEkC)+`K)0Q((GQut-q^99C0llR0?-up^Pe@lQN(tV*yfJ^ z9Dz`uQk>`iQZraRxz<1x1d@iXcvSP=nC%5=lT;72uLwYg^Me59pGhY9@pf7w&?+`= zX2KML5yrOl4!PQv1Q?r~Y*p^_PbZzSKe<|-+-vm)zhu%W0M$@C0rpXOt}G38In3f|pZZ zqY7MLMQcRWyN~{+N$x?-O_cy^l{iTThnZ_d8nA-cee+FV8I)TM+lU62`fl~t^#C^< zD`sq1(;$jb@_eD%Ws$dfQd5K5_jIvElCLEHNp0JV4HkezjHz7DD4#a2G?sOvnFvJy zcc~LBm*VToD-4dc&Z^}lVfPSGx)*|DwHwe@YK_Q--t5Qi6abig7H4+SAok{H^D%F= zXLQ7o(mi7pveN(h9@ByzN*8l@d<8Fir}dJuEv4lTAFi>UWYd7(xLwfWvJ@LCZKZin zQPcjMqV;U-!4RFR41ldI@(HZK?^UZ=#4^fgn}vDX;DcwhzL}X+AmJu1-)NKfR3pi} z8$*}Z?yY?zR0jiF4fx2iJPdlikrEmnT~?XaxV)`?fI!y;Z*j@KV*%bM3wFw!0Q3Yc zI?A8}JbQ|j$V+WjVXXeNsvc>! zV44S(O%MXDT5DX_x68`%?`96&NgIs(8q$Br_j_%pP)&n}a{;@IMPh>2DS;es*SZ4f zAtQR+B=vVIP1#U-nW)b}^V8P;%L{D}+=MsOkt;Q=6?}ak1 z$(j{_I(ZH>WXguB@oosB^`xdmLtkjXYFjmOdlr|jBXB0!pW7cH%A4`)B((&z4!T!O zXosF$#f$A^z2+4cGbxKW<%U?>^WJS^q2Gt@II1hET^YG6{z=}s@2@8DWfp+uBt@It z`9-8V8pYtlKS15X4w#^-A6gYm>3+GlUfNsaWvgWySX)qEgWS0=x0Ec|rM7oPX>G-{ zC!UA*SD78zagmc)JG-i2Wp?$eM*3cp<4}k8c44P}G-@}~Z=k{4aJH`iv`?0o%Rhqd z#mO9F+>K)n*!o(HpWfu22fRy*pM9?wfAH|$7HKi2sILL>e)Aih_^3wqTzLVLV;V`N zEv+*-23&K-uZH zuFB>350!q)_+I=Bi%5I5?S!t^Rsv!JJQUM^|& z6=mL1Q~UD3ZiVXYL!o#n2w^ zE#!1g%4X}D7aQdteqK3wQq~*mKanV^!%;*p;Q!12An^Cf`FVIZH6E^*Orv7eu%C#QCzG)wRa;9XGFO zN)+-ACF6b;CM&6~e4TT_kLL`#;MM1_FnGwX^$0HyJX=>F$H2v+kn{MzxLWMf>Ll?flytgFwR5Sj8$wG2GzYxDWtQ$?}#*TL#t%ij22y;wjDb&KPN z?gcX(29K=pF9<;KADL?1*k9%n!-`QMU>>?z1&2lk z6-rmRE=U;O6Sqj%#4Ikv0X|oGYNlm&p)d;OiF2U4W01#MrSnOvrZm0;3kZ%OH&5g8 z{V68+<;^1PoRR^W-v*|KG;m+Mn>2G?<1zXD3Ouij0??Ug zo?fWy_3}sOs@jF`NO@#(_b`;+OU&X#AN>vTV{w7`9J%oya2qvT#jZ5OTJ3+OeH$2+ zIV8Tv&L3FOzd?pT5>#nbcbHM6obX6$+Q(|(wKOQGNSG0wx0pXG`G{hT6lPzEYZtL1 z?rJ%9il&PXQJWg$(Z?beXWV^W z;kOm*LhIv`&OLX52N-W5_gEY`=3zRQI4v)6iC1FgTpkGII$_)PnCq>iQ?8|di3)o# z*SW8-7UYu-DZj9OJwhd+|Iz)ovx;B@Yzo{i09BQ=TQPCn)D7i90`*dnU!yS=)G%gd zRXT0l7btN8s(NRu>ge8!AOn5zw`rX}?&ILrA*Ie#df;U`q;~K%c;QIrmc;jsW}kQapIJ6j#02#^|HzLP}OcB#GR{cXFPFnJQ&?!}k+^#} z_*nGs@bo+D_6BBeyde)7{#xyv#NW-Cn(%w|%LxfVY))P!BQz%}PX2B4>K-AcVw^{| z37hP1Qquz3YjE|F-F78h2rm=3N0nS1Vv|_+{I)rsH@k-2D{;-;2VRbN2Q7#z0hFJd zDHK8>p%z_5PC<97RfC_l-GarQoP!H{c&@n(lti4VrKL%)il9qqEp87Qm##ObXSzX^ zUu!r^99V*~V_ag0^nd{clsKVj3r0|3!xs zJ@3=`Tw~c~7F70gnp)-C1&(X0Tk2Lzz*ReUwTVUAOD5O4uBvJbaRKMq749Wu;QOhI zr@Mc|O$9i?)LuwevWw5ft}0%Qb*lbRH(E^84MeyTbt5Pt_2m5IT(Vo=T(;9MDE;B_ zO|0)jF9FC08Z;c%b{GELa;*obH7=SPoYJ)&t?+JSR*;go+Jz)qSyKyaZqz+{(tRpW z|LSq*yn4iOMYD!(U^I7mQ1mxX70K2{cI@VBZ0lKu9C2=K!^V=yFKi?wRN{?xKSKz& z-< zgWrf<^c(l5wfK)mnV7wqcwSH1c*QY+@0NeDgEC|1l0&9Vdp4}$P>p$(Er_=+GiU45jn^}MOia{W1^ma5CIsP`Yr zzdkwvpOR<6nXsO=WI&W_!)W45g(kJT0du2az=rF)RH);zoYSsyC#V}@!bE8g-t#(j z=$m?PorEgtWjkAQrkHD0n)7~u-SUk1 zZnuF<0ElwcOlR8S%f>S)vn@)1NTf0>@JSwRQTN3goc);Z)SXL7AlVo5;C+5nDXuQ? z24yj9*#Gd#*e&BM>EPkBY>PcmG!Ts*u2p5Nnlj45=!!hNz!x{F26P+|%7of5{qs2?+LgxWY zbm6MfcnjXVEf`hRlcds$6S4|azzpnj_|YM)1Q;ihI@W_2MrcmdWyMnYYiR+bbjR@0>0Vw;Mf-<<9099y}&ABPnDuv zd2#5YkxMK~W>4+x8b!!r47ovoHXD8k(IV8n=t%Oswok{bI~kzsr(Zl8Hf_*#&2Z)F z+wlx%bi#M5mxYjlyM`a!6VOmwWr($W!9J<6xKqFxrf?0LTE+6rl zjun9Pr9Ju#VyEYx22PeG@$z>#R9XI@>fSLk=mV|IqSD9*u*SxQEajzK8@~xk73j-t zD5)YmXDTsSt;j7nccrDCT7vaXI&fkjSxY4k=DFjE;>63?xhThHp^qHNNx#5n&^O8nIb<``ne~ zz7XZW{a-jHb|VPFLlfdy#cU%$%u+6VAyjm8!3QRiuH$nZv-c5bcF$zCE&Be*~Z+0yS% zWFkxE2HOV@Lzp(*UZ<;DhklUFH3Pn_==mZekf$u~E&qXG&nGk@o&leeB2PEoyf`={ zb7M5NOw1v--}R5|)wk~lj5Hb}yjs*=5y#gW%}A_&800jR#O_=S+1Xa-^2ZNJKzX+NTJ76p4y}h%`l7D54@IhA71Z32X#WI!aX#P^3tW z0wD^4NEcCoP(lw0f^-r}2$FL0&77P6T%4J6bN0-0vG>h-XYF@AWxdjCH2c`F_s|}r zw*nwR5D2t~kKpugdou)r_ZZcB35@SaN6Mc_N64S#0}NQ)iqXEHR1zhfaPsZm|L_HP zUZ6mbVLV=QjG;mk**diNO{?bK4_!rzA3d9%=H}_!9)EPEe!kdU|Ii)pAzP#92Zq>g zMtR#VnA=YTLlMS36$t28NS%6|9{3>JrpZ7|^A-t-gm5K1ZeL$kd2$T@?7OhfzxMcR zkJGFG=`W=8n#W$dy;rbS=XC~9#jOEX-;_!{Y6uJX&9(9zZ+iLn4*d@^nl96M)RPCY zV_D4r2hb%vkahBjrh1ZIZumlF|7gFwWBf>At>y9PuX`P$=HXbZ5Dyd`f)48A{t{;m z^B;}gMEPPKD5WIF*c_kbr>9V+t(j#H0@P?V;r|9r&?IuU6UxlnCbO0R&*sPH^=j8) z0{(bJ^~yp&dmC?qn>bn7tAH-i5kMO_4X+{e&R;sJ^c_B`6*6w9qZ8>(*kp(w_6p$373S;f0HjSMrPZj@K6;_Qi~K8U4CO3qg>F-(87)@ftJ% zmYUB9Okd+ZV)xZQ5;<8I_5Td^Fuy@NY~Vz(T#pJ~v$X z@vsYh7o{B6TE4e^6Lgp1GK9FAp28&+%v?Xk-T!hp?=BDYXVQ*y+;Uz^YJWjF01IHY zB2LppKUjT57ab~PD5A=i?)+1ja8e+hMWuIc%&51nz}no@lT2<*TjBiL)!pc8Bb~(o zuNB(0!VnJ+XM22eX?Ha{9P6@QDVUZzBGf{WfZhr(7$wtwy{~3rTbwC3J z)^|ZLBL%o$6$uvc;?Zq#DHpd#FB#x_Nl-<>EA?aYO9WlTS<&exQ{O|k+ZAF z<6+r*j(%f&gWe-1uk6ZVT0y znp0G+ll>57;xF#n)7t5Dq4M>=*NgrV*#?jmVIIih2E&X85~Umf)`#;9OCkPJA>)lG z|7O20>Jdtxu3D1=eF1ctV3k$Hx%cqE%%Vav;>0pLd%Ysew+lLg@)`!biP9`zs!Ut2 zLe;O@8SA4LL;OP{6zrO<4Z>l@($kog{qNSP$llLNi)H8 zP5s&$PM)U)@OFy_Vx3g)>OHuy?!E~9q|>pW5^pU5p8<<{e6y^3uz2l9MAV5-#OWL! zh7SM;2r*Tz z_YjnT0dhwvzq(aPTgj^{p(u?}m!Lq2;PK*J{Z#_Krv9h@@SM;JmIo5#GMu}>M=;Jq zn^+O`xbaqskTPApwT>{Tj55j`aBqJ1SlZQ5cxFB1j$(QZ#l22cICukVD2|n@(3)GSr&3v8OA9hgJXvB`S=kiKbRBOMI;djF(Hr_s=C7s+WR)I5NN4^~jM9^JtTIHJn5nDg6npBnc zei0Oq|gP=kf6bIWl?BwN)qL8qy#%W|GZg;7r3ZyDCDQl zG3$pk8v*-`=B^Niko`~%T zAiQ!gmppK?Fah)*G`dXg@Id=SSPJxA+M!Y^QW4`qS4zkFRWkyd@NVsAcp%@a_3Q3+ z&WB)k&+dr#=wvJ_7Ty1P>xesUV}IbHer}L8n$LD~rPA`>;JcP!beh zU&S~=^IFk)hgRw0{>y?d=i#bm5;mj*}xTe4m{EMH0N1z#+)Z3rh z5GLd^nLd7e&%jmzI2F|M+Vfsxi-mKfQ<1><{fp5XK>~A~jxRsee!ev$c+paeRSrS2RX0v2Kmey(K&)5$bqYa@8Yf3gl4=hmK6ja7|3{$t4` zCz-MCvQMNra*xH1Oxad%99Mn_<3Sb81@A+rY_~p6>JWRJS6tQf<0I(nM~MjwM&4H` z!y8dS|2%+Ob!CnZtaFLe_kE>--?$An4}@hI=hiv%jYyj@#*PPiAI(VR9&wZ5fs`Hi z1hK+^a#KV*p)4CeFLvy8hn1VNL`kA|Z*+ko0d5c&*}C|49_a7CI$Xi zA;Fm|H%8j?EOofN}_3!aOWQB zNOR?KPT=$R=Jo&Tgr0EnAnTgXjTM*mxVJkRrktDVRJxiMYQ>g1T99sqQ=7;BL&N>==n6g*OEXtNh}=!>0H!kvgqjVx zeUu|azqQn6DOM)^j~T-ac<5G5XCvae5v?S4q`(qo=-;%!jl9Z+tRpzr-sA zHL2el!XgE4y*_$9A&Am@m~Qm#=X@gdl1x=V!>{CXCvpA3Qr$y3&69@2*ZJbdC&C92P0MdDTjbaWww7)4ZSnlpE`L7Vb(HrV94#q;`Z=gLGg3UpD3orytmo? zny(F&Pj}!Zl@S8WSk|c-upE~JNU`l`*y#Mz*)!md-dL^-eY-6KAHBY0MNjVOdfbpC z)3Ehc9Hy4p{&C_eiD2QApuG3$?nMr)Fk3Yut#F7`(?l^i4MfLcCr9g7n8k2+@>tL7 z!Cl+6@2i;Y+U8u1P%Do{IGy~I--V+FOi4~5@pkI|5f3@GAvIM9olc#^u~%?dJy*oJ zrS{w4J#(p72Je17e{nAX%BLZ4xRwcp*9o})Mn;Lhvc zd~V;bvr_C@g+-|P_*lKmlpqH+=!+N`>8d=+dVy=Ibm-OiHUGrTSw#MlkmjYRvybVZ z&AsJzV?Q0sDr{{h<y*kUm`%&4wDNYUR7OFQh$d*)a z>fh279tf$BD!60>i9e?RfYDCdH}Rt7ux{TE7f-yCDsq;OEy$d^sM3@-;F|Het&N$L z3ET;aYa+JWjXrN}HepqGGxj44p;Y6(&AE@HwYe`}^y36fUT{i&8K!BxldL%xKX>Zx zZC5!3F+*T9oOz#N*+z`_)Mdr*UPeE^%usI!E4uj0oBzns&(liD^Y$q<)5#W-=s8!7 zn)D1XP>Yw{z1p=Ac5$*b@cNApTkd4;HehhpDSXPvAi@Y9UPOw{xK+1(eI56uZmf68 z+G)+VuZ^SOS6;)`ztD;<$5*uhZ{Z?4gzD?0DUJ-a` zrs!8*z7K^H-}di3irxNhYquKYALbIgU2gk+3E6PQ-s|a&NEAQczKZ_eD+;mHI3-ZO zsRmlXh2Gnqr`_dF=YirhX9u>lOxKmQnpOxOu3srKYtzsWw`_{E<%^mUI{Yhz?p83$ z139xC=&8(DuV<9~%Qz2ZfM8W!(nM1*`_Rjc@qr^xynm?MQZin$GmOr>rOoOqKr4SFDJ5C_WdtR07&$T>y2YQo)2y{-j;j7*pLPBIx%=BTbWDX=VK)YF4^Ry{*GAVm zOf~75+%HsFg5_H)z86Mgf0e1iv^L59g*p0@VaH)c2*XX?DJ+nv9gdUp!aCS#tp9U8nS{dz$pt&WvTlnRc#sn zRBwUERuKT|!tX;2otTwjn;;M=6SOj@6;gDIuZ*N8FP&b zyl*Po588fP{|y9^zVUA+rrdWaphkcCu43UT|m z;+(geKb^c{joKSBfaaqF`Fj<};5()T;iVHh4nI^rPJC&_o$M!lT^ zG1YXpDrPBMXlf82w^&RXSVXP7>GFZg6RItL%wEz_ty$pv&lV`u2O`34< z3{9r#Bq3oQjNT-iXk6@eLa3`yuLkI$Olz&0YO4VoTf^|t0nJr(hUN`2sd--m?nK$qQTCzrX;Sz(u}&hEU_bRxZ-SapLy`Gg81 zdVBpPJK7ucx}DLo=cZ$ z{ed|&ns1EC%Dy|_|7ec0zr^VRMCz_dkQ8QOQ zKzp9JIodZ3eXr+%w3?~iVqZ-C1Nix8{Xfo3Z;_Qu|0#RmKLGn~vn}YE0MeBbbJ@6c+>IJ}ncDmP)=5VI|1C{&Cb8UOl+ zNwQ0g_X#Qh#L2A(r?Wo0T8{o6_&GLtuz#gH^X_$lKdYavgJrj_p3ZpV?lxaZ$b5s? z)^8Hrnj?;(A_Et_a)`O6*9QTggD8Gbh%TGJTu*GxF;SAFBinXk>09G<#}=E_2Un$y zm!6AU`@k4ka6=GtzUwN^I)1z^^P9h9t;vTA57;Otjamvno}0Y|7Z&|7cB z7bd_2oBpV?$;4kK1rgP2ll5zlB7-7HJkY1!B-7K}&U`~svjXNI!y^45jUwUC#Ml2V z1(K?8lsNg}SiN=$`2B72Sj_l1Vc0>z&zsx#;NxLvkX>3P;NKf(s?B)=_<3l|tR-Pi zQN0uw0RgViqEo|K491w&WHcHH2|xu|q7*k0ciy8{)MM<@`U7t}XR_1l_d6*=t>8Z} zz^)WqnN`J*V8a*&)Pg7Mi!|5dTEP0lqWn1--Cu}TxKA(hz^}Y-U-PN4U;kYBMoqye zu-7BPzac!d22os`i_rTN;6TnXaO1CykY&jJMzsi}Pe^Y9@h)=YPUqk4R-KVxJ67`B zL?W2wWOIYIow|zk6z`#h_k^j}$G*KcMG|caVBY#3dq*FhR(dZYmn%F(2k^ADD5xJa z)eYa(D`HONK)aL-2DgR>wY68OVXhaM@T02ydiTCON*UAAwTmy_z4q|?z{JYUoAHU9 zaPyGw@Ubr~uBJYLS?Q^?jN9o>{VGbXS4DJDX|+TTL*o9~6WgfpSV;%hEDS$*?^9#7 zFwCeri6$ab0z~=>@fkp{=5HFsST}6?Pkc(?BBiFh0$c{AH*MtBwW{joee!h$gyrpj zESE`l2Sg5e(e0OjCxDV3MRbze57eNd#%w^A@-d7DSgaI})k+FQ};TPuiIxs{e_f7`7|f=^3Au{A^dz!O6spgyY6$DRkx6oM;?cFO13<0O)6~5CRiKH zYy1(2fw99`OGT{P1I*`~4DK<#&EPNnCQ=VImYkK9QIiJ{7T&SrUuxxcZZ?D#bFv70 z5rc6>1%O>^SJv&;9-BHl0km9}sEf30%UO6)UOA{6`>~>ptH>7)<7OH5Db{uv!?@)E ziwRW-xE#^IU86zbknaZ_oXk@+a^Z>F(dgSptY^(s5sGTv$@UvgISoYLzkegW%b${+ k;g}qtl=k0mu-{`Y`9A>spCA4|zia-F)BN9x3-Lz(7sJc7hX4Qo diff --git a/doc/Eqs/hexorder.tex b/doc/Eqs/hexorder.tex index 09fd4e2706..fd7dac717b 100644 --- a/doc/Eqs/hexorder.tex +++ b/doc/Eqs/hexorder.tex @@ -2,7 +2,7 @@ \begin{document} $$ - q_6 = \frac{1}{N_{neigh}}\sum_{j = 1}^{N_{neigh}} e^{6 i \theta({\bf r}_{ij})} + q_6 = \frac{1}{6}\sum_{j = 1}^{6} e^{6 i \theta({\bf r}_{ij})} $$ \end{document} diff --git a/doc/compute_hexorder_atom.txt b/doc/compute_hexorder_atom.txt index 4b7e6c4109..90df71d749 100644 --- a/doc/compute_hexorder_atom.txt +++ b/doc/compute_hexorder_atom.txt @@ -10,28 +10,27 @@ compute hexorder/atom command :h3 [Syntax:] -compute ID group-ID hexorder/atom cutoff type1 type2 ... :pre +compute ID group-ID hexorder/atom :pre ID, group-ID are documented in "compute"_compute.html command hexorder/atom = style name of this compute command -cutoff = distance within which to count neighbors (distance units) [Examples:] -compute 1 all hexorder/atom 2.0 :pre +compute 1 all hexorder/atom :pre [Description:] Define a computation that calculates {q}6 the hexatic bond-orientational order parameter for each atom in a group. This order parameter was introduced by "Nelson and Halperin"_#Nelson as a way to detect -hexagonal symmetry in two-dimensional systems. For a each atoms, {q}6 +hexagonal symmetry in two-dimensional systems. For each atom, {q}6 is a complex number (stored as two real numbers) defined as follows: :c,image(Eqs/hexorder.jpg) -where the sum is over all atoms that are within -the specified cutoff distance from the central atom. The angle theta +where the sum is over the six 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 central atom is calculated using all three @@ -46,9 +45,9 @@ 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. An order parameter for atoms that have no -neighbors of the specified atom type within the cutoff distance will -be zero. +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 neighbor list needed to compute this quantity is constructed each time the calculation is performed (i.e. each time a snapshot of atoms diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp index 9555bfc4ed..54b018512c 100644 --- a/src/compute_hexorder_atom.cpp +++ b/src/compute_hexorder_atom.cpp @@ -38,10 +38,7 @@ using namespace LAMMPS_NS; ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg != 4) error->all(FLERR,"Illegal compute hexorder/atom command"); - - double cutoff = force->numeric(FLERR,arg[3]); - cutsq = cutoff*cutoff; + if (narg != 3) error->all(FLERR,"Illegal compute hexorder/atom command"); ncol = 2; @@ -50,6 +47,10 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) : nmax = 0; q6array = NULL; + maxneigh = 0; + distsq = NULL; + nearest = NULL; + nnn = 6; } /* ---------------------------------------------------------------------- */ @@ -57,6 +58,8 @@ ComputeHexOrderAtom::ComputeHexOrderAtom(LAMMPS *lmp, int narg, char **arg) : ComputeHexOrderAtom::~ComputeHexOrderAtom() { memory->destroy(q6array); + memory->destroy(distsq); + memory->destroy(nearest); } /* ---------------------------------------------------------------------- */ @@ -65,9 +68,6 @@ void ComputeHexOrderAtom::init() { if (force->pair == NULL) error->all(FLERR,"Compute hexorder/atom requires a pair style be defined"); - if (sqrt(cutsq) > force->pair->cutforce) - error->all(FLERR, - "Compute hexorder/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list @@ -102,7 +102,7 @@ void ComputeHexOrderAtom::compute_peratom() invoked_peratom = update->ntimestep; - // grow coordination array if necessary + // grow order parameter array if necessary if (atom->nlocal > nmax) { memory->destroy(q6array); @@ -120,16 +120,16 @@ void ComputeHexOrderAtom::compute_peratom() numneigh = list->numneigh; firstneigh = list->firstneigh; - // compute order parameter(s) for each atom in group + // compute order parameter for each atom in group // use full neighbor list to count atoms less than cutoff double **x = atom->x; int *mask = atom->mask; + double cutsq = force->pair->cutforce * force->pair->cutforce; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; double* q6 = q6array[i]; - q6[0] = q6[1] = 0.0; if (mask[i] & groupbit) { xtmp = x[i][0]; ytmp = x[i][1]; @@ -137,31 +137,62 @@ void ComputeHexOrderAtom::compute_peratom() jlist = firstneigh[i]; jnum = numneigh[i]; + // insure distsq and nearest arrays are long enough + + if (jnum > maxneigh) { + memory->destroy(distsq); + memory->destroy(nearest); + maxneigh = jnum; + memory->create(distsq,maxneigh,"hexcoord/atom:distsq"); + memory->create(nearest,maxneigh,"hexcoord/atom:nearest"); + } + + // loop over list of all neighbors within force cutoff + // distsq[] = distance sq to each + // nearest[] = atom indices of neighbors + + int ncount = 0; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + if (rsq < cutsq) { + distsq[ncount] = rsq; + nearest[ncount++] = j; + } + } + + // if not nnn neighbors, order parameter = 0; + + if (ncount < nnn) { + q6[0] = q6[1] = 0.0; + continue; + } + + // store nnn nearest neighs in 1st nnn locations of distsq and nearest + + select2(nnn,ncount,distsq,nearest); + double usum = 0.0; double vsum = 0.0; - int ncount = 0; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; + for (jj = 0; jj < nnn; jj++) { + j = nearest[jj]; j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - if (rsq < cutsq) { - double u, v; - calc_q6(delx, dely, u, v); - usum += u; - vsum += v; - ncount++; - } - } - if (ncount > 0) { - double ninv = 1.0/ncount ; - q6[0] = usum*ninv; - q6[1] = vsum*ninv; + double u, v; + calc_q6(delx, dely, u, v); + usum += u; + vsum += v; } + q6[0] = usum/nnn; + q6[1] = vsum/nnn; } } } @@ -178,6 +209,70 @@ inline void ComputeHexOrderAtom::calc_q6(double delx, double dely, double &u, do v = ((6*a - 20*b1)*a + 6*b2)*x*y; } +/* ---------------------------------------------------------------------- + select2 routine from Numerical Recipes (slightly modified) + find k smallest values in array of length n + sort auxiliary array at same time +------------------------------------------------------------------------- */ + +#define SWAP(a,b) tmp = a; a = b; b = tmp; +#define ISWAP(a,b) itmp = a; a = b; b = itmp; + +/* ---------------------------------------------------------------------- */ + +void ComputeHexOrderAtom::select2(int k, int n, double *arr, int *iarr) +{ + int i,ir,j,l,mid,ia,itmp; + double a,tmp; + + arr--; + iarr--; + l = 1; + ir = n; + for (;;) { + if (ir <= l+1) { + if (ir == l+1 && arr[ir] < arr[l]) { + SWAP(arr[l],arr[ir]) + ISWAP(iarr[l],iarr[ir]) + } + return; + } else { + mid=(l+ir) >> 1; + SWAP(arr[mid],arr[l+1]) + ISWAP(iarr[mid],iarr[l+1]) + if (arr[l] > arr[ir]) { + SWAP(arr[l],arr[ir]) + ISWAP(iarr[l],iarr[ir]) + } + if (arr[l+1] > arr[ir]) { + SWAP(arr[l+1],arr[ir]) + ISWAP(iarr[l+1],iarr[ir]) + } + if (arr[l] > arr[l+1]) { + SWAP(arr[l],arr[l+1]) + ISWAP(iarr[l],iarr[l+1]) + } + i = l+1; + j = ir; + a = arr[l+1]; + ia = iarr[l+1]; + for (;;) { + do i++; while (arr[i] < a); + do j--; while (arr[j] > a); + if (j < i) break; + SWAP(arr[i],arr[j]) + ISWAP(iarr[i],iarr[j]) + } + arr[l+1] = arr[j]; + arr[j] = a; + iarr[l+1] = iarr[j]; + iarr[j] = ia; + if (j >= k) ir = j-1; + if (j <= k) l = i; + } + } +} + /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ diff --git a/src/compute_hexorder_atom.h b/src/compute_hexorder_atom.h index 99cd987161..38f34cab89 100644 --- a/src/compute_hexorder_atom.h +++ b/src/compute_hexorder_atom.h @@ -34,13 +34,15 @@ class ComputeHexOrderAtom : public Compute { double memory_usage(); private: - int nmax,ncol; - double cutsq; + int nmax,maxneigh,ncol,nnn; class NeighList *list; + double *distsq; + int *nearest; double **q6array; void calc_q6(double, double, double&, double&); + void select2(int, int, double *, int *); }; } @@ -60,11 +62,6 @@ E: Compute hexorder/atom requires a pair style be defined Self-explantory. -E: Compute hexorder/atom cutoff is longer than pairwise cutoff - -Cannot compute order parameter at distances longer than the pair cutoff, -since those atoms are not in the neighbor list. - W: More than one compute hexorder/atom It is not efficient to use compute hexorder/atom more than once.