From 214d346bdc05c2f8bd0d606ad3cbb0f1b3d4fda4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 2 Jul 2009 16:38:31 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2952 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/Eqs/heat_flux_J.jpg | Bin 0 -> 5580 bytes doc/Eqs/heat_flux_J.tex | 14 +++ doc/Eqs/heat_flux_k.jpg | Bin 0 -> 5138 bytes doc/Eqs/heat_flux_k.tex | 11 ++ doc/Scripts/correlate.py | 89 +++++++++++++++ doc/Section_commands.html | 8 +- doc/Section_commands.txt | 1 + doc/compute.html | 1 + doc/compute.txt | 1 + doc/compute_heat_flux.html | 147 ++++++++++++++++++++++++ doc/compute_heat_flux.txt | 142 +++++++++++++++++++++++ src/compute_heat_flux.cpp | 225 +++++++++++++++++++++++++++++++++++++ src/compute_heat_flux.h | 39 +++++++ src/style.h | 2 + 14 files changed, 676 insertions(+), 4 deletions(-) create mode 100644 doc/Eqs/heat_flux_J.jpg create mode 100644 doc/Eqs/heat_flux_J.tex create mode 100644 doc/Eqs/heat_flux_k.jpg create mode 100644 doc/Eqs/heat_flux_k.tex create mode 100644 doc/Scripts/correlate.py create mode 100644 doc/compute_heat_flux.html create mode 100644 doc/compute_heat_flux.txt create mode 100644 src/compute_heat_flux.cpp create mode 100644 src/compute_heat_flux.h diff --git a/doc/Eqs/heat_flux_J.jpg b/doc/Eqs/heat_flux_J.jpg new file mode 100644 index 0000000000000000000000000000000000000000..c30ae8dcf888efc0dcdda2d9845f685c264f3999 GIT binary patch literal 5580 zcmb_gXHXMNyA2QsMG%nQ1Bn=_^eQDFCG>=%At)e74^5f~FIA){5E2k+QbLiUVh905 z0Vx3#5dun)8bJ^gq^XEp-h1D<^L=yg++X+X%=2gW*=P35b7ptX@rUCD05`(S(hR`B zzyPo~8Gz#%fC+$^2?zu-F`ooxW@Z-9DOS)4PIGXuo#Hvo%gb|`hv$rdnBW;c5q=&X zA+V5$xVWUG2w>)6;RS=vs7M?0@jG~^h9>5+%BaBwAYp=9CPb+_;8XvG5H} z>#=`Ll1FSuwp^|}>om8@KHHNoYSmn3Wg#b(ljCSK!k~R9@{9CQu^_rk(4W#lskz>0iV0cwT#j%WDUr3Y`iU$59VlCdo3S48y%yzLCTBT<3<+V|eq(nb zWNuvTe_EH4$rLXyfo* zX!Fc#g%EPFK;@F~tR}KJNT9Djg?$87Nmk8w3Cx^-hdbZ30llD*ai6$S8E@ab<4;oJBcC>AX z5WiO{H}t7zSjk$ns<%fiXNkMl+|kH9#(WEk`Smn0x!CdTn?1v$%45Kw)-m9d?^H|b zyG_ho*$dy{#vBX3fOMS!bi8_MC6LTFoetmsYR|z&g+E_V9av$L+pc&V0ItMq!_baZfZ!A_DWlS!g^MxKlb(XZVF{5aL0*#k=<+Z zLA+jaKYeG-euR$=1d`yib=#9@O5{Wet?dqHuF@L(@be`4{hx zd9o>DF@B2qa~{`HNlSLAL}N>tbL5w zhvK8zJT&l5@+pTWZD2YiP+==7)*zB!5L3{+gJ0AIvau`=Jra!`z>s;RlK3)iVce z=3QP?3Oo2682x%#R8P6?^-(xpR$?AqwfD+M-p%i}Qvl)nM?0U`f`vnKwi zGE${OY_?cOp(}qQ{vh?$&96;!d*gLvp^zVs?LY17m9y?8~c@Xz5IptUb`twB!*i? zw~)=l82kLYUINs0i&qq7)ou6N@q38>APY*(Q%g=4UtMBL=a+hz`|}*5A-ciebth03>NC~PWRqoDH z%vWl>nsdupFZpvKmT+F3)tmUds^;zoU>R+RqI~Io>HOAToc+~mA8`|pY`O=8kXZ#+ z3&jVV&p##Az2~benF8KZ%bgD7cK<-D?YwV;Q3~!*b57^G7u2QK-_0I1i?%u%A&9LZG!v1%@#j8%Uqmej-R$`}+ zea|?)G$dN%O}71>0_U}cUBXnRZ`+KNd6CDaDH;`KGPGH>X#4Xto@8bY5nEWe(SP+O zg!QFHL@QLGyoO&mH3ile+7niRU!5(z@Ww@k;<1|JfdIO%`b3Wbw5sR$ zyj?#->OR9+-~Kz0Prg$01PUO3JYmmA$U(^9Lbk-8 z0KNew*bUZDPyM+$#|?sJ_-+@RvON;8q*mZ$kb9Qz?DT3_Q9iLYbc5rX{?Qkfrt$p} zmujp?*ddSVW-ANd>)YAwb`*ml5Riu|#k17*)~IIP74l^O(c;rac~$qJD(BYbvNmLo0F~>fr@GU9 zF$e^*aB!{h@l3SvksQlGlL9%ma&j%CMeVU#YgqYOQC9UmB%fycFgu(gqwUa5)t1W0 z8;_#8{E#nX&MZ;{kG-=#^txl~1qkE*#a&TqUVF+%Mh7HT(-awypGgah=&%?-jWxc$ zt)bEREW6M(jqr0IqI}45UucAvHx>e8E#OSIDencZj4dIT*aaAbyDjMwF`zyd-FkzI zMb!^*i`bC<`aml-Vgt863cymUQRE;SIz^F?h|A#Y3R+6a#}uE%3;QiPj#j6ia-V+IrQTvo_s z;y(MZKZKKtJPl0lNiZ<~veFVLca@jxOV?l=IQoF@? zoUlcYZfy#NcEICxl?wD44qi2*QFnF+U2UM2h(xdgX*>oP4~U;&7{^xlD|&l0XWwq^ zrm?qMrNyOC>ZMt;EN+}$NzQSkf@CkyODaL)7t^gfrA;Gg74!(G5@Qko!hDMNe6+#T zqsh+BYG~cVo+j&drOAC7xx#?K^1+aE&ckBvDgqZ5AD7OXON&~JUQSD`ph0#CDvi=GIClGw%dWZUc*um@9yLXJf2 zMm;KvRa#b>Nm6BLR@+$4ny$bXlM_{^%TaQdPj34@$xbKb2Zxt0wcjb>4D+Ny=Wx@c z*R+991q+z77jeOD$Z=W+An9zT8mBQ{Np)p`z`l_;R$U$Buj}Sa4~0?^5M(ZMX7jy< z?Wz{3b`kjt1snmhsj@VHzE_yecU43GVMN6jgZ#Bkhgiy-7oyF$Fq~ARzSkSUv;-#b z!zw%3r@3qqn)MQU0psM@$g(HhPp4(hD!u^TFwYvY1hX0ePL*1AeUCA)5pSOQ7=i8Z zOK26SAi^9*5$WG82hKk0Nn{Z?1|+45X#Nh~AmOnCSqtIx_Xg{R_M~1nlRu!mMui7A zjs29}Aj54M;?#|ES?V9qiAt1`E}kRxPo_z(?6YYTT-no#*1c&oMS))JEF!t^`i|HJ z+>8;K(`*aal8>zsH-56m?4oniX<7Fz1HJFZRx0%L%hs^5gokXb%o2}*_1VIz`I{PC z3|0=_)>w8(6GDpUcZ%3Y!8ibAV?=4hAUT+Dmk@x{Mm)Hk{PAbUAy%6kbIfPCu-Qi2VtuXZuK( z73ycU^|NSR@XPcPerpkL?R&R*xmo zDbS4Fo{ljfo^RNvm2WuSeKbJc2+x|Z`wg5dU(1Rw^%suPlt$BDd(y*lu4{<dZco0Y;b31^Xu%{fT{^ zz-ar;+Lcmcp?3B&y~R+|H2n@1^(eIea>S>`{B#$Cd|JxF52+v9&i+mF>Q_oT14IQH zC?;P&2KvREFCGDKC2%bseVAqs|6yyfLQ)>Xw0$|`Cr^anIA2)dl&II?xD@sKbZ{%L zHZ&k-==R3e^wWjh4xd(PUXxP$Ju9sG&^4`~l++|Pgg{IC*V73tbuh6oVRm3D`xvcL z|JnN0lIkPXxb6nm8ACpF7y4NQSbr6r8V+i|Rz8*g11H8d&G!X%nQ0BD9~Q8rQ9U8N`%+|ztpHV0ACy8g zkrCrIAY5I?c0BWq4Rpd_ohf-(LOy-xhWGEC2Ourcl6nv*_7WTq}85=54U>{zy!}oY_NcL1Fcd(y&j) z_(AXaA(lVeS%de2$7zHXEkbtIk3KfHs$nW+Nawg(_k2b@}l_HSR>ly1mg^1RS` zV;(a>8v5wmD+QPtIoF=VjrUTYo6YMhGwE(~RS$YFADzE3jE`1mRXM;Y!3L68zB`85 z)$8N@uh-u^J1^beLCl?;RY^Hm>4^mrZ7d?}V89>jCYhmN!|y|WSZ)D{=HTZ#W z!ZoU2U^&?ra0Xf%RHl_cyK)ThX;oEdoYC@WJ`4`(H!q!2+brJ6a=hyj{6N|6O3U3$ zXAJlNOICk=-_{2JG6IGC*$86R`7H|xM+>egVfKzPX8b9ZEc4JoK`%A1C0xWhjC&?s zjJJjAZ4wzBbZN`tgnMl1!ELFNR6+6c?h; zj#TGuO#U_EwZ)LLmKnvj#^Uu+EWc)jnybP&6k)49tn|F4`JLSA*mO>@Bv4X9?#+oy z(+CsDv`)kCW@2RB9ES2KYsl;si(2Zo`p!X{W`9tm7QgWA*=Q;C8drH_61<2B2J-01 zENN{C;a>~RXe7QO0f}&A0JvFv_n(BETPmihqgTatuF|(IH){;k+pj{8+&VV^B@03G_ZZ#3yn6B>H?g)H;_U^@BgHbtJwl7O&$Ibsw0W+WM)bfj402a>q`o3O&U1;<%Wv|GUdpt|-87L_5eFs}W2cAmP=f**41A7DYSzzpoEu zqlQPKdq1AAML=biH?F6JWJbh=r_@(8DJse|lLfd$rmKoSk${?~nXdT>{dRqZnVC2r zbh`fz>8h7DyR^ZVAoXRA56x*3bg;XE)J_+`Yi&7JVYzS^V7!0%M dL9zfO5&QpYH~!N_|GRJb?=1h{X)lkb{sWlWNqhhR literal 0 HcmV?d00001 diff --git a/doc/Eqs/heat_flux_J.tex b/doc/Eqs/heat_flux_J.tex new file mode 100644 index 0000000000..2c2b232058 --- /dev/null +++ b/doc/Eqs/heat_flux_J.tex @@ -0,0 +1,14 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ +\mathbf{J} += \sum_i e_i \mathbf{v}_i ++ \sum_{itq$!}3(2J-bV4-(Gij)YU3MM4f zgD6F+^b$aX2ueGmRF%W^yEFH@bMK#fcV~8XW_R}4A3HnqJZE!fD*zYD1Ze^=FfagS zza2QE0fqn*BM1ayWcnSLn3&Epv$HV&HcmD+FgrIV1j5b9&CSCve36Iu5+66WfP}y$ z5fL#lF`kQ(P)Sj!u&9{mSvSDO0_ZWkU|`?|&fWsT00Y1XF#g8)w}F_>pJQZZVEJ9u z-~t#J85z(0h69~v{8tCVIYtn`#C;wj!K1`1X<)_6r|d3;2#({wsDiig@XRZlV!4Vv zTL4)8HSrp7nvcTf5a9n;C^KqzBrHOmR-HeZ6$qUH9@>@!wCazn3)VkpHD={PiOu4D zRHM5%55{e^Z@4vwomTNr*E&)EW3b&tXG075S=$`7gj_gE1I*Ye{T}D!KQkK~f0NdL zSvLa3|B9Gc^A!2@4vKf)$!EE!@wmXh_|k1}2xLD3yUKZpMe<$#U(4pMb5~KSI>Kha zVX(16`}z>6mJ+r`i7(uP%!^qAWA2%y)~<10lwa%-!jWpk=KPtbo6Hxq^o65$08o`q z0zrJE(~@B4_iWfg`$a(h3H?#aT)|wSeHp@KCC92KxQ-s>yfE)s7_};z$!wFS1|+>4T6ryltGzZ7!5|qL0NKo9uj% zH5;^^6wZiM{_pK`>Hn}J!%`u;X^%b{&Dm{{`2GlTNL!rc#KdsXQb8L!L@1rnm2e*F z!&Q*Tku;a`#=WIbzXy|{2#C%#Y+4z?U?_1w8M)u6J#iF^yf<-cY63?6V>Xnc01-DG zrta}pvA_jjF@>My?%U;RQ@%m%5XE(qm(A6gRpm0&;51O)io_K&yNVDTSkJ^QvT^npCG=0S)%yEucRulwryWzANHlvYb(mQT}3sixOxADp5Q!w zg7S1UTHz4iaR%fXAWv%^Cli+4Qku8%yXar#v>SfTEz$lNZn4vZ%7A$)kmPJ7I98xn zB$OHvndvo|hDu19)jVNGBncz2hGQK`7y`pP6mA83zdu(xbNTu5!L@e>Hwq_)RVAH{ zm^z3G^7*`#roQF4&v5iabZuJpXU2yVJ_TB6)-tu#ciuF{jb*xrCqX!eXl_`4dI*Qx z#uhUQvwwU4%}Zfor{z|Iw(_I|#ueWi;8X_iFIpX77p%>%g%y_XncJ%ePEZDHg(^}u zd0|}T8pC&M?2G>45n+`@JX^$(M4r$SrPGy>O_N%cOXS?c?#U|(s>M(1Wvw$85&1A* z>of@mGqX2-X=N!lt0Vj5jMeT4IO(a~_Mv@$Rkvx?F0z+BhIip`t=nx&Ds0zxO2(kD z1uQV5L*Qn%^%qp-q?YS+)r;jDP8loFd#cvnaKEvl`--JrU(vc#6`$FCOdi)tsaX}kTJC!jrI|Jgdn-4J8Bud*IQwuGhlJprsX?I5J8?Gs1n_Tam0ohtr)}IHouxnLU(hb#JOfQT{WQ#l$ z;?NSo7NPg~n01JZSw);xIy2vN`-A_uxLYpd0_-w|A+YV+q9Jln->lnvS_aT` zihamBT=f<#)^lrS)7lOGj|6{0A4k!?p66E)kH*@&XTaM&rB_3gh!V0ItupVRcL+VG zy>M*Oz*QtpC!icf0{?D8qF!x(7FfAy;BZagrt^-dfS>kG34I?t^URsJUEe~p$NMXi z32YZ&hVJ4w1$2UpZ(LPt`s>;X-vgV`JcWx}lS*6@#ol0l2uvd>lYU!TzRa9es_T8D zdi9j+|Ld9#-eEuaneS#qJB?3DE))d8NGKN{4O?2A^cnDCDyL2iU1OqCp|*ft-}zu^ zR1B+pB}1?n^u1;t*i!xwqucF>Q#q{bf8zMM{cdg-;^mU?O-z|gB6pwX0ru_2*QcMY zu)NEwgy-9KP`FN1mPJ~@^yQtp#+LI1QTV#;USDic+A43W9I6PKgVb zIo8-Zy^Z1RM6!@pzjG#=+CoUp@8qFF~ZNckd=`$XFkuLJmInQ2a~tgi40z?R z=JR!1Gazsy`nX}8z#{Iay?2~Sl!J+-t_vC1nhI=?!z{?3HNLt%eY`9_;HPqBAcLnS zkfT1Yz;@OaOV7zBZQN%~(>UKLjYcE}HXA1g-1Xc$r%s}`rlZ0PjA~)CIaPcIcxn+l z!+oGmv<=ZHs&gEB)3W#(?pEqI6j%Tx!Cj!YF(PqulF}eI33s45oP@Wwbz`&@MkLn< z$AViX79+9`!6?tnU1vLr;L*;dSJp23nNvZMM(H<74o305`Yqpy=Ba)DO`RwSU5ib} zmQNN{V4;`Ko&^x|gQf!E(11_4!tX-^dT%zBO$)> zt~eW*e*XfJS;#WrX`;GidyT5i9=FtJ*J9g%?>vK~&W!@|+b!hOdE1$pkKm<54gCQo z8#iJAg>B;@MpeFNYKdgf>bRR``e#E;VG947M)(KGa5zbR!g`G3L8SxNx1Zb4H{5KUIfcq;kjhVs1bSP#(7~;H{(TEE@P4pPABy000$JTQRWs&@% z)Su<#S3)maDRk|70Gd6{$JpoLJ1E8Cx%z7l4{JvnVhTsR7Tvt(kf>RD+ugUnu`~4Q^omye_A%>inE;f#JTx7{ ziCa1SdTZy$2gsJsUG@z^lVCSnnsMQ_DHeig7V$+o#)KrYr#yHq*<^bvR<_IIlE{hK zob@72xLQj&y2SrLL1EqOML9AL5>qNyRtitt$A2gs4_En9VgoiM2ZN=0qnr0jsOA{S;*y7MQ!uQTgzZ z)l$x0iGgTRtDL0|Sej*6PsB;b5MWm75v@BaxpIT@g(&^FxRB;pY_Oa?NtF4mnZHt?OP|C!mTS+Ifh>D9wAJ)2c@#-8?Du}u5F zZE3f|ccSOp+Z3!8&C=wD6O+bCb=0Z>H$7Zc+FEC{G4T~!FKyEWboq-<(l2kZQF}aH zC-aA43gq6jSeJVzI8;kw1fV8!B@s>v6&R*MhXV_elP1&yH$CNJHI~*#Lb|Z@GMR@~b3&v0?C~qR&p=0mRq-;f z_5r_*#iN*RS z*cUK3$Govq;Qoc6d_CT3k#=l~bhw=5$J3hcw@PDO235?|*+{JoUA);dtd35p6lOj9 z9ktv5V7FV4QQGY{ajH9dY*$6vz2CL$Q7P+euOKT4{-@%ThOgL?Ib5t0bL{bLb*}#K zcyd0KwU=SYtM+zI0(EE`mXUi=X@F(XwyagUI-XO|?ncf`Kd$>Ff##*3yf68)uD7JB z=+Q#qG8@%j+f+huKmr9Yz7s}RhwsgcLu>R3UYirJv0=*B(yAaWdPz=D)e=(1#I{*2jahnIZZPl& z9;6!3HiMyg2RcwHoHrjPLKm%c5%xxGTy5_ScJi%@&{>D!j$E!GaP6)kPPZUYA)Zpvmm3jZ>1N(b}*iSl} z&8R$b1hm8wVA%dmxdA=L=SF0jB<;Gg;A6g==ArA^Y9t+UP0z&v6N^Q(vu+wAQ6kyHdDGqInlhX-<9 z3}z>`-(kvc>JB;V*Dm><`{%U2?| z;`L+;O^)y_5X%K6h*GKuN2ugEK`)p%NIm1iR9sw!24md9%jhmY)5jOwCMIAe6PIjp zUo(RE_j98Z8|J^f4aLuCs2s4ij^__m9MK&DvUueC|FBZT1x34gxgOV{}T<*jFBMS?;W@2-i{A+pA%xvW?7D>{A}p!-6q2VaYq z8F594CO}xW{ftc$S+6P{y}8Qe7GYjmFypIRO2z2H!>{#+D}?4!v%F#aCX}*{q^~5x0(+}lYsv#+y74`qdGsvONgX&VNDR?`LOWbi;lGzG^}rIn*?2>2=C#3BqNQe z^NlZehK!qSzC3iugtN%O71V?rCyKk7c2H3VaHM8$xZz}?R+7l3*Neq)kg24hpNa@S z%`lNErkTJe_C}vkoIQL@(zwcEK|TT*p<#jiv$Y0QNo`HnoQ9*PvUtligfD2+MkhzM zbayz?7j`o*<+EC&d=FjFSQ2+Uhmj}RL2z>P%i+$>FEX8)?_edn83x8q5qMgEe<~af zjf8IKOtt&B!CJ4d)DvMx(T(x}@LTv8L@+@n3q4gC90_K8NyPnIxM};OI9xs2mo*{k zPUm0r_k=r|y#__c%jTJ6OimbmNRDsyD`pBmVatZi+Ex$!SYX9IM%iqUQZm4F*1UEo zuVwHC>(y`tS#NpRYP@v<6DU9*!|daDe>p{eRVCgrZr&9QgAS{}$o4o9XvwZfPz(ZcE3DV~Lt28EV~ zj!x(@Ie%X`7N^Y18Dv^r)2aTV&M|nQJ*;qH*sOuec5aa34%P=%XXzub_q9O;xB$Sq w(`>&Jj?1BW@|zmBoOi|!_`gJOu|DpxZnZaC3D8@w^ZmbMv46`sc4zbd0npoNy8r+H literal 0 HcmV?d00001 diff --git a/doc/Eqs/heat_flux_k.tex b/doc/Eqs/heat_flux_k.tex new file mode 100644 index 0000000000..a6aa2aeabd --- /dev/null +++ b/doc/Eqs/heat_flux_k.tex @@ -0,0 +1,11 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ +\kappa = \frac{V}{k_B T^2} \int_0^\infty \langle J_x(0) J_x(t) \rangle \, dt += \frac{V}{3 k_B T^2} \int_0^\infty \langle \mathbf{J}(0) \cdot \mathbf{J}(t) \rangle \, dt +$$ + +\end{document} + diff --git a/doc/Scripts/correlate.py b/doc/Scripts/correlate.py new file mode 100644 index 0000000000..7db1763fcd --- /dev/null +++ b/doc/Scripts/correlate.py @@ -0,0 +1,89 @@ +#!/usr/bin/python +""" + function: + parse the block of thermo data in a lammps logfile and perform auto- and + cross correlation of the specified column data. The total sum of the + correlation is also computed which can be converted to an integral by + multiplying by the timestep. + output: + standard output contains column data for the auto- & cross correlations + plus the total sum of each. Note, only the upper triangle of the + correlation matrix is computed. + usage: + correlate.py [-c col] <-c col2> <-s max_correlation_time> [logfile] +""" +import sys +import re +import array + +# parse command line + +maxCorrelationTime = 0 +cols = array.array("I") +nCols = 0 +args = sys.argv[1:] +index = 0 +while index < len(args): + arg = args[index] + index += 1 + if (arg == "-c"): + cols.append(int(args[index])-1) + nCols += 1 + index += 1 + elif (arg == "-s"): + maxCorrelationTime = int(args[index]) + index += 1 + else : + filename = arg +if (nCols < 1): raise RuntimeError, 'no data columns requested' +data = [array.array("d")] +for s in range(1,nCols) : data.append( array.array("d") ) + +# read data block from log file + +start = False +input = open(filename) +nSamples = 0 +pattern = re.compile('\d') +line = input.readline() +while line : + columns = line.split() + if (columns and pattern.match(columns[0])) : + for i in range(nCols): + data[i].append( float(columns[cols[i]]) ) + nSamples += 1 + start = True + else : + if (start) : break + line = input.readline() +print "# read :",nSamples," samples of ", nCols," data" +if( maxCorrelationTime < 1): maxCorrelationTime = int(nSamples/2); + +# correlate and integrate + +correlationPairs = [] +for i in range(0,nCols): + for j in range(i,nCols): # note only upper triangle of the correlation matrix + correlationPairs.append([i,j]) +header = "# " +for k in range(len(correlationPairs)): + i = str(correlationPairs[k][0]+1) + j = str(correlationPairs[k][1]+1) + header += " C"+i+j+" sum_C"+i+j +print header +nCorrelationPairs = len(correlationPairs) +sum = [0.0] * nCorrelationPairs +for s in range(maxCorrelationTime) : + correlation = [0.0] * nCorrelationPairs + nt = nSamples-s + for t in range(0,nt) : + for p in range(nCorrelationPairs): + i = correlationPairs[p][0] + j = correlationPairs[p][1] + correlation[p] += data[i][t]*data[j][s+t] + output = "" + for p in range(0,nCorrelationPairs): + correlation[p] /= nt + sum[p] += correlation[p] + output += str(correlation[p]) + " " + str(sum[p]) + " " + print output diff --git a/doc/Section_commands.html b/doc/Section_commands.html index b553dd271c..5a7602f8ff 100644 --- a/doc/Section_commands.html +++ b/doc/Section_commands.html @@ -343,10 +343,10 @@ description:

These are compute styles contributed by users, which can be used if diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt index ef0ad89c7d..8ba86af7d6 100644 --- a/doc/Section_commands.txt +++ b/doc/Section_commands.txt @@ -462,6 +462,7 @@ description: "erotate/asphere"_compute_erotate_asphere.html, "erotate/sphere"_compute_erotate_sphere.html, "group/group"_compute_group_group.html, +"heat/flux"_compute_heat_flux.html, "ke"_compute_ke.html, "ke/atom"_compute_ke_atom.html, "pe"_compute_pe.html, diff --git a/doc/compute.html b/doc/compute.html index ab5641d615..1273319a88 100644 --- a/doc/compute.html +++ b/doc/compute.html @@ -117,6 +117,7 @@ available in LAMMPS:

  • erotate/asphere - rotational energy of aspherical particles
  • erotate/sphere - rotational energy of spherical particles
  • group/group - energy/force between two groups of atoms +
  • heat/flux - heat flux through a group of atoms
  • ke - translational kinetic energy
  • ke/atom - kinetic energy for each atom
  • pe - potential energy diff --git a/doc/compute.txt b/doc/compute.txt index 741eb71536..c0045acb6c 100644 --- a/doc/compute.txt +++ b/doc/compute.txt @@ -114,6 +114,7 @@ available in LAMMPS: "erotate/asphere"_compute_erotate_asphere.html - rotational energy of aspherical particles "erotate/sphere"_compute_erotate_sphere.html - rotational energy of spherical particles "group/group"_compute_group_group.html - energy/force between two groups of atoms +"heat/flux"_compute_heat_flux.html - heat flux through a group of atoms "ke"_compute_ke.html - translational kinetic energy "ke/atom"_compute_ke_atom.html - kinetic energy for each atom "pe"_compute_pe.html - potential energy diff --git a/doc/compute_heat_flux.html b/doc/compute_heat_flux.html new file mode 100644 index 0000000000..a645707a8d --- /dev/null +++ b/doc/compute_heat_flux.html @@ -0,0 +1,147 @@ + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    compute heat/flux command +

    +

    Syntax: +

    +
    compute ID group-ID heat/flux pe-ID 
+
    +
    • ID, group-ID are documented in compute command +
    • heat/flux = style name of this compute command +
    • pe-ID = ID of a compute that calculates per-atom potential energy +
    +

    Examples: +

    +
    compute myFlux all heat/flux myPE 
+
    +

    Description: +

    +

    Define a computation that calculates the heat flux vector based on +interactions between atoms in the specified group. This can be used +by itself to measure the heat flux between a hot and cold reservoir of +particles or to calculate a thermal conductivity using the Green-Kubo +formalism. +

    +

    The compute takes a pe-ID argument which is the ID of a compute +pe/atom that calculates per-atom potential +energy. It should be defined for the same group used by compute +heat/flux, though LAMMPS does not check for this. +

    +

    The Green-Kubo formulas relate the ensemble average of the +auto-correlation of the heat flux J to the thermal conductivity kappa. +

    +
    +
    +
    +
    +

    Ei is the per-atom energy (potential and kinetic). The potential term +is calculated by the compute pe-ID specified as an argument to +the compute heat/flux command. +

    +

    IMPORTANT NOTE: The per-atom potential energy calculated by the +pe-ID compute should only include pairwise energy, to be consistent +with the full heat-flux calculation. Thus if any bonds, angles, etc +exist in the system, the compute should limit its calculation to only +the pair contribution. E.g. it could be defined as +

    +
    compute myPE all pe/atom pair 
+
    +

    Note that if pair is not listed as the last argument, it will be +included by default, but so will other contributions such as bond, +angle, etc. +

    +

    The heat flux J is calculated by this compute for pairwise +interactions for any I,J pair where one of the 2 atoms in is the +compute group. It can be output every so many timesteps (e.g. via the +thermo_style custom command). Then as post-processing steps, an +autocorrelation can be performed, its integral estimated, and the +Green-Kubo formula evaluated. +

    +

    Here is an example of this procedure. First a LAMMPS input script for +solid Ar is appended below. A Python script +correlate.py is also given, which calculates +the autocorrelation of the flux output in the logfile flux.log, +produced by the LAMMPS run. It is invoked as +

    +
    correlate.py flux.log -c 3 -s 200 
+
    +

    The resulting data lists the autocorrelation in column 1 and the +integral of the autocorrelation in column 2. The integral of the +correlation needs to be multiplied by V/(kB T^2) times the sample +interval and the appropriate unit conversion factors. For real +units in LAMMPS, this is 2917703220.0 in this case. The +final thermal conductivity value obtained is 0.25 W/mK. +

    +

    Output info: +

    +

    This compute calculates a vector of length 6. The 6 components are +the x, y, z components of the full heat flux, followed by the x, y, z +components of just the convective portion of the flux, which is the +energy per atom times the velocity of the atom. +

    +

    The vector values calculated by this compute are "extensive", meaning +they scale with the number of atoms in the simulation. They should be +divided by the appropriate volume to get a flux. +

    +

    Restrictions: +

    +

    Only pairwise interactions, as defined by the pair_style command, are +included in this calculation. +

    +

    To use this compute you must define an atom_style, such as dpd or +granular, that communicates the velocites of ghost atoms. +

    +

    Related commands: none +

    +

    Default: none +

    +
    + +

    Sample LAMMPS input script +

    +
    atom_style      dpd
+units 		real
+dimension	3
+boundary	p p p
+lattice 	fcc  5.376  orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
+region  	box block 0 4 0 4 0 4
+create_box 	1 box
+create_atoms 	1 box
+mass 		1 39.948
+pair_style	lj/cut 13.0
+pair_coeff	* * 0.2381 3.405
+group 		every region box
+velocity 	all create 70 102486 mom yes rot yes dist gaussian
+timestep 	4.0
+thermo	        10 
+
    +
    # ------------- Equilibration and thermalisation ---------------- 
+
    +
    fix 		NPT all npt 70 70 10 xyz 0.0 0.0 100.0 drag 0.2
+run 		8000
+unfix           NPT 
+
    +
    # --------------- Equilibration in nve ----------------- 
+
    +
    fix 		NVE all nve
+run 		8000 
+
    +
    # -------------- Flux calculation in nve --------------- 
+
    +
    reset_timestep 0
+compute 	flux all heat_flux
+log     	flux.log
+variable        J  equal c_flux1/vol
+thermo_style 	custom step temp v_J 
+run 	        100000 
+
    + diff --git a/doc/compute_heat_flux.txt b/doc/compute_heat_flux.txt new file mode 100644 index 0000000000..1cd1727ff1 --- /dev/null +++ b/doc/compute_heat_flux.txt @@ -0,0 +1,142 @@ +"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 heat/flux command :h3 + +[Syntax:] + +compute ID group-ID heat/flux pe-ID :pre + +ID, group-ID are documented in "compute"_compute.html command +heat/flux = style name of this compute command +pe-ID = ID of a compute that calculates per-atom potential energy :ul + +[Examples:] + +compute myFlux all heat/flux myPE :pre + +[Description:] + +Define a computation that calculates the heat flux vector based on +interactions between atoms in the specified group. This can be used +by itself to measure the heat flux between a hot and cold reservoir of +particles or to calculate a thermal conductivity using the Green-Kubo +formalism. + +The compute takes a {pe-ID} argument which is the ID of a "compute +pe/atom"_compute_pe_atom.html that calculates per-atom potential +energy. Normally, it should be defined for the same group used by +compute heat/flux, though LAMMPS does not check for this. + +The Green-Kubo formulas relate the ensemble average of the +auto-correlation of the heat flux J to the thermal conductivity kappa. + +:c,image(Eqs/heat_flux_k.jpg) + +:c,image(Eqs/heat_flux_J.jpg) + +Ei is the per-atom energy (potential and kinetic). The potential +portion is calculated by the compute {pe-ID} specified as an argument +to the compute heat/flux command. + +IMPORTANT NOTE: The per-atom potential energy calculated by the +{pe-ID} compute should only include pairwise energy, to be consistent +with the second virial-like term in the formula for J. Thus if any +bonds, angles, etc exist in the system, the compute should limit its +calculation to only the pair contribution. E.g. it could be defined +as follows. Note that if {pair} is not listed as the last argument, +it will be included by default, but so will other contributions such +as bond, angle, etc. + +compute myPE all pe/atom pair :pre + +The second term of the heat flux equation for J is calculated by +compute heat/flux for pairwise interactions for any I,J pair where one +of the 2 atoms in is the compute group. It can be output every so +many timesteps (e.g. via the thermo_style custom command). Then as +post-processing steps, an autocorrelation can be performed, its +integral estimated, and the Green-Kubo formula evaluated. + +Here is an example of this procedure. First a LAMMPS input script for +solid Ar is appended below. A Python script +"correlate.py"_Scripts/correlate.py is also given, which calculates +the autocorrelation of the flux output in the logfile flux.log, +produced by the LAMMPS run. It is invoked as + +correlate.py flux.log -c 3 -s 200 :pre + +The resulting data lists the autocorrelation in column 1 and the +integral of the autocorrelation in column 2. The integral of the +correlation needs to be multiplied by V/(kB T^2) times the sample +interval and the appropriate unit conversion factors. For real +"units"_units.html in LAMMPS, this is 2917703220.0 in this case. The +final thermal conductivity value obtained is 0.25 W/mK. + +[Output info:] + +This compute calculates a vector of length 6. The 6 components are +the x, y, z components of the full heat flux, followed by the x, y, z +components of just the convective portion of the flux, which is the +energy per atom times the velocity of the atom. + +The vector values calculated by this compute are "extensive", meaning +they scale with the number of atoms in the simulation. They should be +divided by the appropriate volume to get a flux. + +[Restrictions:] + +Only pairwise interactions, as defined by the pair_style command, are +included in this calculation. + +To use this compute you must define an atom_style, such as dpd or +granular, that communicates the velocites of ghost atoms. + +[Related commands:] none + +[Default:] none + +:line + +Sample LAMMPS input script :h4 + +atom_style dpd +units real +dimension 3 +boundary p p p +lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 +region box block 0 4 0 4 0 4 +create_box 1 box +create_atoms 1 box +mass 1 39.948 +pair_style lj/cut 13.0 +pair_coeff * * 0.2381 3.405 +group every region box +velocity all create 70 102486 mom yes rot yes dist gaussian +timestep 4.0 +thermo 10 :pre + +# ------------- Equilibration and thermalisation ---------------- :pre + +fix NPT all npt 70 70 10 xyz 0.0 0.0 100.0 drag 0.2 +run 8000 +unfix NPT :pre + +# --------------- Equilibration in nve ----------------- :pre + +fix NVE all nve +run 8000 :pre + +# -------------- Flux calculation in nve --------------- :pre + +reset_timestep 0 +compute flux all heat_flux +log flux.log +variable J equal c_flux[1]/vol +thermo_style custom step temp v_J +run 100000 :pre + diff --git a/src/compute_heat_flux.cpp b/src/compute_heat_flux.cpp new file mode 100644 index 0000000000..320cab75a2 --- /dev/null +++ b/src/compute_heat_flux.cpp @@ -0,0 +1,225 @@ +/* ---------------------------------------------------------------------- + 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-93AL85000 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 authors: Reese Jones (Sandia) + Philip Howell (Siemens) + Vikas Varsney (Air Force Research Laboratory) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "compute_heat_flux.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "force.h" +#include "pair.h" +#include "modify.h" +#include "group.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM}; + +/* ---------------------------------------------------------------------- */ + +ComputeHeatFlux::ComputeHeatFlux(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 4) error->all("Illegal compute heat/flux command"); + + vector_flag = 1; + size_vector = 6; + extvector = 1; + + // store pe/atom ID used by heat flux computation + // insure it is valid for pe/atom computation + + int n = strlen(arg[3]) + 1; + id_atomPE = new char[n]; + strcpy(id_atomPE,arg[3]); + + int icompute = modify->find_compute(id_atomPE); + if (icompute < 0) error->all("Could not find compute heat/flux pe/atom ID"); + if (modify->compute[icompute]->peatomflag == 0) + error->all("Compute heat/flux compute ID does not compute pe/atom"); + + vector = new double[6]; +} + +/* ---------------------------------------------------------------------- */ + +ComputeHeatFlux::~ComputeHeatFlux() +{ + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHeatFlux::init() +{ + // error checks + + if (atom->avec->ghost_velocity == 0) + error->all("Compute heat/flux requires ghost atoms store velocity"); + + if (force->pair == NULL || force->pair->single_enable == 0) + error->all("Pair style does not support compute heat/flux"); + + int icompute = modify->find_compute(id_atomPE); + if (icompute < 0) + error->all("Pe/atom ID for compute heat/flux does not exist"); + atomPE = modify->compute[icompute]; + + pair = force->pair; + cutsq = force->pair->cutsq; + + // need an occasional half neighbor list + + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHeatFlux::init_list(int id, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeHeatFlux::compute_vector() +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double rsq,eng,fpair,factor_coul,factor_lj,factor; + double fdotv,massone,ke,pe; + int *ilist,*jlist,*numneigh,**firstneigh; + + invoked_vector = update->ntimestep; + + double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + + neighbor->build_one(list->index); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // heat flux J = \sum_i e_i v_i + \sum_{isingle(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); + + if (newton_pair || j < nlocal) factor = 1.0; + else factor = 0.5; + + // symmetrize velocities + + double vx = 0.5*(v[i][0]+v[j][0]); + double vy = 0.5*(v[i][1]+v[j][1]); + double vz = 0.5*(v[i][2]+v[j][2]); + fdotv = factor * fpair * (delx*vx + dely*vy + delz*vz); + + Jv[0] += fdotv*delx; + Jv[1] += fdotv*dely; + Jv[2] += fdotv*delz; + } + } + } + + // energy convection contribution + // uses per-atom potential energy + + if (!(atomPE->invoked_flag & INVOKED_PERATOM)) { + atomPE->compute_peratom(); + atomPE->invoked_flag |= INVOKED_PERATOM; + } + + double *mass = atom->mass; + double *rmass = atom->rmass; + double mvv2e = force->mvv2e; + + double Jc[3] = {0.0,0.0,0.0}; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + massone = (rmass) ? rmass[i] : mass[type[i]]; + ke = mvv2e * 0.5 * massone * + (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + pe = atomPE->scalar_atom[i]; + eng = pe + ke; + Jc[0] += v[i][0]*eng; + Jc[1] += v[i][1]*eng; + Jc[2] += v[i][2]*eng; + } + } + + // total flux + + double data[6] = {Jv[0]+Jc[0],Jv[1]+Jc[1],Jv[2]+Jc[2], + Jc[0],Jc[1],Jc[2]}; + MPI_Allreduce(data,vector,6,MPI_DOUBLE,MPI_SUM,world); +} diff --git a/src/compute_heat_flux.h b/src/compute_heat_flux.h new file mode 100644 index 0000000000..54aac14d79 --- /dev/null +++ b/src/compute_heat_flux.h @@ -0,0 +1,39 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef COMPUTE_HEATFLUX_H +#define COMPUTE_HEATFLUX_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeHeatFlux : public Compute { + public: + ComputeHeatFlux(class LAMMPS *, int, char **); + ~ComputeHeatFlux(); + void init(); + void init_list(int, class NeighList *); + void compute_vector(); + + private: + double **cutsq; + class Pair *pair; + class NeighList *list; + class Compute *atomPE; + char *id_atomPE; +}; + +} + +#endif diff --git a/src/style.h b/src/style.h index c6c6560411..1a5afd7237 100644 --- a/src/style.h +++ b/src/style.h @@ -81,6 +81,7 @@ CommandStyle(write_restart,WriteRestart) #include "compute_coord_atom.h" #include "compute_displace_atom.h" #include "compute_group_group.h" +#include "compute_heat_flux.h" #include "compute_ke.h" #include "compute_ke_atom.h" #include "compute_pe.h" @@ -106,6 +107,7 @@ ComputeStyle(cna/atom,ComputeCNAAtom) ComputeStyle(coord/atom,ComputeCoordAtom) ComputeStyle(displace/atom,ComputeDisplaceAtom) ComputeStyle(group/group,ComputeGroupGroup) +ComputeStyle(heat/flux,ComputeHeatFlux) ComputeStyle(ke,ComputeKE) ComputeStyle(ke/atom,ComputeKEAtom) ComputeStyle(pe,ComputePE)