From f227080d703c85ee4c55cd06d07a670e8cd4cb86 Mon Sep 17 00:00:00 2001 From: athomps Date: Wed, 4 Nov 2015 05:18:21 +0000 Subject: [PATCH] Added hexatic bond orientational order parameter git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14231 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Eqs/hexorder.jpg | Bin 0 -> 8635 bytes doc/Eqs/hexorder.tex | 8 ++ doc/compute_hexorder_atom.txt | 117 +++++++++++++++ src/compute_hexorder_atom.cpp | 263 ++++++++++++++++++++++++++++++++++ src/compute_hexorder_atom.h | 73 ++++++++++ 5 files changed, 461 insertions(+) create mode 100644 doc/Eqs/hexorder.jpg create mode 100644 doc/Eqs/hexorder.tex create mode 100644 doc/compute_hexorder_atom.txt create mode 100644 src/compute_hexorder_atom.cpp create mode 100644 src/compute_hexorder_atom.h diff --git a/doc/Eqs/hexorder.jpg b/doc/Eqs/hexorder.jpg new file mode 100644 index 0000000000000000000000000000000000000000..cebb8d9d9e6e37800088e2e37d320a795fa46d03 GIT binary patch literal 8635 zcmdUTcT`hN*YBZM1(YI10s^AaR6wLi6a+-1iu58SQeucw3`mKBC>@oiAfQN*8l@@% zkuIXr2^|sy=_Hg8B;_8S=Y79tefN+1-S4iu?pk-}?3~}s*;DqOnLT@kF~XPu_^%oo z8v;yBOh7vL0~iG6i^lr;t`=72hQ`+nzz+bxeEQO5Jp*Qb06_Q!`db<4N!i&uNOAN6 zoB$7S7*GYY-8=&BT(YpZ27u&;^n>lc|1*(?0k9?j7*N`;>py?{FH_tecl-lEDJsZk zJz)VJ0Kjw_#BfyL9WbBi6^I4h;eX-x``F*wN*~N`1MxwxzwzKccK;jy+Q%DPunm zr2ieHzgr;86O{eCfBU5XfxiX!Eq3CJrsnyRf93vP{M#RgzpLA{{L8*f!r%YN`==p~ zApd*+N-;630>FQj{-=&tM*zUQ003C2|KzD90>EoU0N_jeC-3+f0N_^vfTmXrJP>FQ z>~jkMOih7PV5^(}E0Y+&2BQ5LpvxrwAJ~xT1iC1VD@_X(Vjf0g@^Eo* za&xk?a~lZyaPR={K@jrnSF&H0<*%}AJe=&DJpXNCd;tVFfdD2g7N%nW zvj7u|028APfPk8?fvaqPq5X@Pm|0laz#f7900nCJ!S1uLFoT_CV+9St6b^m|SOwS) zNuRmIe%R8DFnR##X z3kr+g6_-?2RoB$kef(74(%RPE(fPIOTmQh|&@gUfbPP|J{WUlLdtq^jw6VFhP2QpG z{@K@y31InKt^d&Mf6_|;)QcG$16GcGy_lF$AhQUtvPqv|KXl2G!|m?jV`m?73hF(3 z^P!nbM%9Wanji;&*Vz-Q&dfs*BFTHnZe3X@a z&g|H|Gqno^PUjAKxDAS~0?zJMD>4tI>K=dNf}OXhyLSJ%vGAwwg-K!7r z{Y5u+8E<^{$AkEr97z=EIO@Ux%qhmRXwKk524IqSys@@lJu7svtZ%eW%sO@?ugd6H zILbGjXDXy!5ZW6zwl7KEmq}$gPQ{g%zy1^r z6?i2uGPriG+-rKPh5@LV8iY*8F#w;V8ZFT3^(eZq%R-aX#~~^_Z|VMOP66|HZ7zZ> zYOGW1=QTng0yq5TO7zp0zyzPrLaJZNI{hBCxAq?QiM$8@ZKe)`=Gmnl&M@z%^HFAk zX9=>AJF~V?x6Pkm>m;Fk`KFb_RQ*q+?o2nA^4hPm>Li{0*>JWe z1Nc2@PCI6_pd`eeTjJe%Wa zpP+Mn#<<42P2PdHKGIR(`%=73DFku%P=@nYyEc3ML(z8Zk^zL|5sv0*0qD)p+)*6i z=i3SjrrFj-Qc?^)zEv4-Y2pbuYKgkVY4yF z`t|kjZwasDMayUF(x-nl&5i8FlY422Fv=_e`=H)7YHQL!M+0wu*oamW04uF|*tZp| zs4P3mZ!9n};WWNTM$WAvkA-9~9sNrB3cZC*T-`eoW=~Xk*l?zldOpR1lC(1xj(MV9 zt~A}7-nWPE_mv;CId45^(;BWiar0QQEE zGI+!ZVKgzx_rsC{;&Da=Ex$x1<1R|ahu^-trNic8*~{}*`G2rf4ah1N12DW!(q{lX z)12sy;cU$!h>uX<_-B+)llLe2aLEr>O>ll*!NsiQ#-(a+;eKiPc{s%J6>7#tX^2-R zbOhx(jCR8cQ@rpZtv%x9KdWX_z$>xfv5&UG<=cfb1gM&%^eVN|fKc`7 zD*Yh2{$^C6ljEd5Yl`Bzs!Ce6iy`{;HUprXknij{u(;v01pT1Wz96v0pVjOS}dKS{T#6Ujs+Y<;}FofU5! zUr6cv(;H)Ry7>VFn=IOM5ZVxyH9GyXOSXSIacwmSr7&ve?Jn<1N!1(nR@|>V`pSlvdh+Q`4LK!uHl&vm1dP5-F9C zvmWY)8{y)H@L!5n+A zI@o8PJzDXHx_+F#b9__>d$3+fWKw9?KqYPG_KCa(8FD+V+Gu5#zyO#i3SY)TIVnU0kInlH0M3w*!F(RJ4f}zl>ue^XYIrW zjr`6>H#d$uGGAjJL5o4*Dv}JKWs!m)n#|M`J=$Ub@l%SukoJfZjZJX~S3{yr2K(2! z)K-bw{>g?NBpDvGDa~=avQ8kx%=kPV66a4n^qKBbD{{7GYK(rsh~(FeHOV#Yb3H`i zp>wMerpGkMY81X8VuccqF{Oq?|ARO*ML6Lyd#$%UW9(pgse$05mu_(inhvb@U?F}= z0;6d!llvmh(kdhP!=wjm!XX4~`}z|x$#T4jcj54xVI5O9SWsaH{fe>QY{=oS2mKI4 z`y0kNI^o+tuIiK*_0=pRoa@JX5W|b3rJsEm0E?sj?z9Oy6I;}GqZ;Ac({}$t#T)$c zXZRk^f)!zGe(XlNc+3R5t&x?!0yC%7(@S4V)s@uW9`Y%Do!6j1wx>)J;L{1yhav>1 zp*a`y+g5=>bvHk&#FNk)%MV&-ZZt6xan@TxgfVpydJU#oV!yjRKq zT;N`#3)gM@E?_b3z6~!d4Fow)ZmO=>l6*JnQh;~A6++tYCXT}DtBx8 zYhkK2OsSDflS)Br1i_Xfnu5*S9Z*RPN`>$U{)Q(O{Xtn*Fl7~QG?V7JEOh)&u9cPL zhl=%Lx*0$Mw63T3^BvaDK2&|Y zIm>#(+1jO_YvZ$eMDb*KL(OhQor&z6@Ed5nH8rym-n2HYJd8AqqsWlS6!RHAk%)eW z?BB`ITe)MM1e?O)+V10xGf+_{9|xr`0uPEMxCq*EQ&HSMe(^rr$tn-VQBUwiL$+Yp zy5Dq(&sh5H5@~rk6fMFP7htb!>#^-0M?_NOZKDJEIH;^2t0)0P90A+%nQzr+vKkrV z(3V+&k&Qm~eOV_nk+fmQ#@!UbWVkC*yxkK+7aIz5#`9$Hu>~jXwA@c<=e?g@P~Q0c zJ@Dndz=RHmM8%L($HvwguykR_!oz2JkkL zluSSDAi@BoEkVbyN`lgpcsigIQ*T#lw3X%RjXA6!)~zQpSE2!~5gTzX3_@|7npnVBQ5e}QsPT))sE#yR-F9==5jpLDpwY-r3) zgNr{NrgxZWR^&C<(nZG!Xqtpb0-BiFi{0f;rg$M^8Nk*Z7jPA{!1WlwuTzEN6RFQf z$R!HxXnb~)JOk)!4(?I z=?_x~@$$JlT*5m{H)JIqT=u_N?W^z`OF6ek=kBba_v>y2Mvkf;^<$q zy`C!Gc>mkzciuH)#T}dPc_(9k5}1!}8?6N0p)1pJ!p_V>6X^n!#CK!H$XP!2 zOCd1WJ4a}mEE#273pYvNoBHDH_Ck8EvQ?b-X4J&NSA9k=wM*FqEiX!FMmsY4nrQp6j-W7M5FuJ%`WMoxX}HyBcK zt4l<&3%n4Q6Z*Bu#hv@w^}XxE1|DKJ*@3w81oWP+g302+@GCj=$(p_Zp{^m7rb$ig z%N+h=6S&9A{vFH;Fl#moCe!NTp)ngPm*>e0Kt{pUL^t7jRqEop_($xy+}-xYTbT}x zeY^1_74l774RU0*6J?WoAVS7x`e1}K)u*wN*pcCbpVUL>H&1a6IecH*8UDk{uSWiX zA}cJo>EEvWQeXCH7p_wl&O(l+oSfworT3$SsOAJrWX`FKS-$q3Xu1e-r!}=Ma%0(; znAqKUzdk{ve)~1Qqg-0s`-!Xl4Ti8dY38eYmuQZ88M5Ke^M?8>8>cl+MMh&LM{8Hf z1#lg0wuWdoufovnW*)pUe@Qf}gCR zuGV1G?Gb9rZSH(b^T}5RJ-=MI#M}TSSJTpTpnF%qnc60H;dduGN>Y{N)v||ugpy<* zR3Dg`G&Gg)#Ubz92|D7=`M_?nTV!mgrLCvxz7`ex7F!Mda)6JVPVs_14hv9z<0(8p zZ8=bL;oOZU7JVwK1@`3_ge(kZ{3cBR=U;!0TT4q-=3e^Kn#MBA9)+I^4;^f|#q>B7 z^&XtPPXxA@OU%c9SQnR?nN5m4M^Rx)O2{YSmI8j+8fJWf8JTm>pHnkte}#4B!?jda zcJ~aag0vhupZ9$R5FK_DpM7X)$Jp_J^EDA+8{t?Va_GmRIQv`GMuWE>Wb`E2)NdHd z-blq6;rWvP6eTeLq zsTMnb($n2uRGe26JsL{BO)_f5#=59dV)rfwKe^D&IQovpS)0%)Esbz3cxZlqWa`c89jlF+ zPt{{RQzka+UU|L0bdE@qjHtpKof%`?A+dP#Xb0tIB|nQd%6|+Oy>n`1x;r#{BvY>_ z&3^WRWYen(o*9mn)x3xvG37du#}(=q9BoKu<>x+QmL2j-mEV+4UlrmVs6X|ZOVFZK z!aG*tuD*nKNeLUOhJVMW<0xk5o0<8VzfXu=z)p$T+ht_^84K4(*CSBi3HC+gziuFS zlVc=-K3#=iaXVs9Tefni6OjSLD9#ORE9q`XD>be*ytBHJuivVmz;DzTVFvDFlEI;0 zril)@qYS{7Vo6LUN4q|rW?!jsmToYvC~PE%>%Z%vN5SC8<8SY(HW!W;?hXg1-4t*2 zteA&j)UFrOvh9YPpeB^;B&JkcnKc#0G7^@#qBTd---#X|%5jcSKJ* zop~@!ZqKqhBfvMaqdWi>f>-$5c<$)g-NZ7RN~IO|69(Vy(uI9Q^TJE_mj3L_n772c zf+r)^u%k2Khn7N@RoEBO9VA3lOUmj5#xFiC{b&K)5{SGwD<8MCP{23Y`dmA-bWU*4 zwo7AY*?*pvpPFleeO3OQ0j&NOoGwfUP4#f!@j_mpiLOY=upo+M@i;t}Xg(dQdm6M< zY?uQ6fH^@VM(v@j>cZxtmNnV3i~r#p#VQl>sO(SbS^IAvyT0x~$7IN*=3{6&Ctjmt zeRREDSCN=T|Aar%VD#D;jCzAHKZ|7@m9}s`d70-XLykFWAvAAZl-hJITA>vYS}t9{egS)! zevHydK*z|A-+c#nk$FDxee*@HmR>rs!DV{IoD=`;ck?h8D45ioLo#!BK5@g4bGG{q3qV1I;TI>WDvlf^Yx&oH#d-bAB!V zkz>2uz}a`#k-U}boC(VVlM#JdDqmD0b2yu5g$=Mgm*uhTn(_pc)`f%~Fx1h_>*YYN zH0naM2nf7qntROE?jy3lPO4uRUqpg83UicFr|^b#3b^Cu;)|E6ZChQ5)=}PeJFPL&Zs=I_nkPEVnfIvS5lCqS z2fdfhPkX)f!^SmQtF1n@>1c!4tx{xwhi-0Y(ZsIhcbShhpOQ4)Z^d~RIDci4jq6U6 z8i@NH|9oMUVlds`?>barGU20(ogZop484cajysX}6ruw5&NCcO(Rfu72XxmWFC|q*VT}REm6({;B0tyAkBlu9C_5^N%qsZOjYsd9G5}Rw z94VtSnT4v=FE^+wMd@kD6nez%;+Sn@TUYZ~{iumf5+L~t_XEm`ohhdupjBi2#$6vx zzb`C*8@$;%#W&y=go%+nRc9ncilYv0&QH7V?z%^D~U&brbU=27PT2NJ9NQm;#$p5cW$S}<_+ zHCNcz7uPt}^>OvV8qJx|Rt(o-Jx41|!K`=dg4me0_5p44S`)bA5gChO!81}kVy7N3 zzXO;aI#^pY4t=X-07^~xu1~r?zThcmQ-&y`I*_~=gOv5GX@RWFusl=#bxhI#- zt+F~aJtL#*d{)qNWe}go-l(K0PqU^HxNn^4Nrf{-UV*_HE82VPWAv6kt}!b*>LhJK zQF0&yy!t)tyXi>jAh&F;;0C^JIy<#1HKra~&#@Ps(5{c#5z^!+K#q)m`AzOOPmYE8 z=knuP4$P!~v^N_4HSl9>@<88eSDL34%kQ-hR(wacubxVM<>atX){yoJv2(7Gb$cE= zhKlf8a?QkM=~@k1^%6%TXWI}~=LKSC!=6UP&#OUdR)Zm`88G}`W^yU38# zSEk2AYEdheb0hoJuE&(qL8Ju$kPXVS=RTxfBG@NZnY>#PQxj4BiFi@-;aN7{3pe*`uu6-Kk7cjq#I^i-oWp(U zLxU<21qE3M^$)(5xJ(TP@TLe^jO_a?7mkDl_QYX5ks}@*e|B0_Mgq(!iLc|ad=wkg z>x7-;HH-^?HzBk;M7}oq^{uIXoG)VVyOx>~qI=P@&O+miFmk;E9_}j?zrh$E(Ea za8I#diL2Ph4JP2dA}{5-Y2Wb=adf26l&i1ZpzxNd=!Q~x^#UwMm4!>p;`>UmaF=hy zkSoz*8T}9~sXom!N$*2fq8^N!Zt#&I2^BC?zppw1UG%5y1)RHe>rRGPnIe}$a@w9I zHIViieVdOpxQN^ykI#8FgO1Xb#;Id@FX7{tKbK+zx&tI?@88BfozN@wB@f$bNZ#s# zZOT6t5-6FW1zo=J8Z%V4keQVHbIH6*Ev@_}a{8Y0q20o57c-qwv*HQLCixmT`J)`E zp_Ju(iu(Zh2`!Z_rM?yLsZU4fuEMgdv@~jx6Jh8UJ^ramboWMmPysEy0o=|U={ziZ zCiWEfmTuE(bCzJybRIk5Bh6>w`Pn5v7v_CwD_sKIL$sSjs3$2qAtdeQI}|ch-1l;L zJ$;=3i9x;{w6rluQpkeG??eW>A2ylOmqkd(btPJC+GN&az5e`(a4UI)vxSFcqKgRo QV_?B#Ao%|!b1+8#2g*Fs7XSbN literal 0 HcmV?d00001 diff --git a/doc/Eqs/hexorder.tex b/doc/Eqs/hexorder.tex new file mode 100644 index 0000000000..09fd4e2706 --- /dev/null +++ b/doc/Eqs/hexorder.tex @@ -0,0 +1,8 @@ +\documentclass[12pt]{article} +\begin{document} + +$$ + q_6 = \frac{1}{N_{neigh}}\sum_{j = 1}^{N_{neigh}} e^{6 i \theta({\bf r}_{ij})} +$$ + +\end{document} diff --git a/doc/compute_hexorder_atom.txt b/doc/compute_hexorder_atom.txt new file mode 100644 index 0000000000..5a5171dc50 --- /dev/null +++ b/doc/compute_hexorder_atom.txt @@ -0,0 +1,117 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +compute hexorder/atom command :h3 + +[Syntax:] + +compute ID group-ID hexorder/atom cutoff type1 type2 ... :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) +typeN = atom type for Nth order parameter (see asterisk form below) :ul + +[Examples:] + +compute 1 all hexorder/atom 2.0 +compute 1 all hexorder/atom 6.0 1 2 +compute 1 all hexorder/atom 6.0 2*4 5*8 * :pre + +[Description:] + +Define a computation that calculates one or more hexatic bond orientational +order parameters for each atom in a group. The hexatic bond orientational order +parameter {q}6 "(Nelson)"_#Nelson for an atom is a complex number (stored as two real numbers). +It is defined as follows: + +:c,image(Eqs/hexorder.jpg) + +where the sum is over atoms of the specified atom type(s) that are within +the specified cutoff distance from 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 +{x}, {y}, and {z} components of the bond vector. +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 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 symmettry, |{q}6| << 1. + +The {typeN} keywords allow you to specify which atom types contribute +to each order parameter. One order parameter is computed for +each of the {typeN} keywords listed. If no {typeN} keywords are +listed, a single order parameter is calculated, which includes +atoms of all types (same as the "*" format, see below). + +The {typeN} keywords can be specified in one of two ways. An explicit +numeric value can be used, as in the 2nd example above. Or a +wild-card asterisk can be used to specify a range of atom types. This +takes the form "*" or "*n" or "n*" or "m*n". If N = the number of +atom types, then an asterisk with no numeric values means all types +from 1 to N. A leading asterisk means all types from 1 to n +(inclusive). A trailing asterisk means all types from n to N +(inclusive). A middle asterisk means all types from m to n +(inclusive). + +The value of all order parameters will be 0.0+0.0i 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 0.0+0.0i. + +The neighbor list needed to compute this quantity is constructed each +time the calculation is performed (i.e. each time a snapshot of atoms +is dumped). Thus it can be inefficient to compute/dump this quantity +too frequently. + +IMPORTANT NOTE: If you have a bonded system, then the settings of +"special_bonds"_special_bonds.html command can remove pairwise +interactions between atoms in the same bond, angle, or dihedral. This +is the default setting for the "special_bonds"_special_bonds.html +command, and means those pairwise interactions do not appear in the +neighbor list. Because this fix uses the neighbor list, it also means +those pairs will not be included in the order parameter. One way +to get around this, is to write a dump file, and use the +"rerun"_rerun.html command to compute the order parameter for snapshots +in the dump file. The rerun script can use a +"special_bonds"_special_bonds.html command that includes all pairs in +the neighbor list. + +[Output info:] + +If single {type1} keyword is specified (or if none are specified), +this compute calculates a per-atom array with 2 columns, giving the +real and imaginary parts of the order parameter, respectively. +If multiple {typeN} keywords are specified, this compute calculates +a per-atom array with 2*N columns, with each consecutive pair of +columns giving the real and imaginary parts of one order parameter. + +These values can be accessed by any command that uses +per-atom values from a compute as input. See "Section_howto +15"_Section_howto.html#howto_15 for an overview of LAMMPS output +options. + +The per-atom array values will be floating point numbers, as +explained above. + +[Restrictions:] none + +[Related commands:] + +"compute coord/atom"_compute_coord_atom.html + +[Default:] none + +:line + +:link(Nelson) +[(Nelson)] Nelson, Halperin, Phys Rev B, 19, 2457 (1979). diff --git a/src/compute_hexorder_atom.cpp b/src/compute_hexorder_atom.cpp new file mode 100644 index 0000000000..1fedc33c23 --- /dev/null +++ b/src/compute_hexorder_atom.cpp @@ -0,0 +1,263 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "compute_hexorder_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "pair.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +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; + + ncol = narg-4 + 1; + int ntypes = atom->ntypes; + typelo = new int[ncol]; + typehi = new int[ncol]; + + if (narg == 4) { + ncol = 2; + typelo[0] = 1; + typehi[0] = ntypes; + } else { + ncol = 0; + int iarg = 4; + while (iarg < narg) { + force->bounds(arg[iarg],ntypes,typelo[ncol],typehi[ncol]); + if (typelo[ncol] > typehi[ncol]) + error->all(FLERR,"Illegal compute hexorder/atom command"); + typelo[ncol+1] = typelo[ncol]; + typehi[ncol+1] = typehi[ncol]; + ncol+=2; + iarg++; + } + } + + peratom_flag = 1; + size_peratom_cols = ncol; + + nmax = 0; + q6array = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeHexOrderAtom::~ComputeHexOrderAtom() +{ + delete [] typelo; + delete [] typehi; + memory->destroy(q6array); +} + +/* ---------------------------------------------------------------------- */ + +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 + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->occasional = 1; + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"hexorder/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one compute hexorder/atom"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHexOrderAtom::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHexOrderAtom::compute_peratom() +{ + int i,j,m,ii,jj,inum,jnum,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int *ilist,*jlist,*numneigh,**firstneigh; + double *count; + + invoked_peratom = update->ntimestep; + + // grow coordination array if necessary + + if (atom->nlocal > nmax) { + memory->destroy(q6array); + nmax = atom->nmax; + memory->create(q6array,nmax,ncol,"hexorder/atom:q6array"); + array_atom = q6array; + } + + // invoke full neighbor list (will copy or build if necessary) + + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // compute order parameter(s) for each atom in group + // use full neighbor list to count atoms less than cutoff + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + + if (ncol == 2) { + 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]; + ztmp = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + double usum = 0.0; + double vsum = 0.0; + int ncount = 0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + jtype = type[j]; + 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 && jtype >= typelo[0] && jtype <= typehi[0]) { + 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; + } + } + } + } else { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + double* q6 = q6array[i]; + for (m = 0; m < ncol; m++) q6[m] = 0.0; + + if (mask[i] & groupbit) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (m = 0; m < ncol; m+=2) { + + double usum = 0.0; + double vsum = 0.0; + int ncount = 0; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + jtype = type[j]; + 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) { + if (jtype >= typelo[m] && jtype <= typehi[m]) { + double u, v; + calc_q6(delx, dely, u, v); + usum += u; + vsum += v; + ncount++; + } + } + if (ncount > 0) { + double ninv = 1.0/ncount ; + q6[m] = usum*ninv; + q6[m+1] = vsum*ninv; + } + } + } + } + } + } +} + +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; + double y = dely*rinv; + double a = x*x; + double b1 = y*y; + double b2 = b1*b1; + double b3 = b2*b1; + u = (( a - 15*b1)*a + 15*b2)*a - b3; + v = ((6*a - 20*b1)*a + 6*b2)*x*y; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeHexOrderAtom::memory_usage() +{ + double bytes = ncol*nmax * sizeof(double); + return bytes; +} diff --git a/src/compute_hexorder_atom.h b/src/compute_hexorder_atom.h new file mode 100644 index 0000000000..b9430addc2 --- /dev/null +++ b/src/compute_hexorder_atom.h @@ -0,0 +1,73 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(hexorder/atom,ComputeHexOrderAtom) + +#else + +#ifndef LMP_COMPUTE_HEXORDER_ATOM_H +#define LMP_COMPUTE_HEXORDER_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeHexOrderAtom : public Compute { + public: + ComputeHexOrderAtom(class LAMMPS *, int, char **); + ~ComputeHexOrderAtom(); + void init(); + void init_list(int, class NeighList *); + void compute_peratom(); + double memory_usage(); + + private: + int nmax,ncol; + double cutsq; + class NeighList *list; + + int *typelo,*typehi; + double **q6array; + + void calc_q6(double, double, double&, double&); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +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. + +*/