From 111d350a223ffd4edf2b9c00e9ac35e72d962090 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Tue, 28 Mar 2017 12:34:46 -0600
Subject: [PATCH] fix gcmc units change for chemical potential
---
doc/src/Eqs/fix_gcmc1.jpg | Bin 0 -> 5617 bytes
doc/src/Eqs/fix_gcmc1.tex | 9 ++
doc/src/Eqs/fix_gcmc2.jpg | Bin 0 -> 10626 bytes
doc/src/Eqs/fix_gcmc2.tex | 10 ++
doc/src/Eqs/fix_gcmc3.jpg | Bin 0 -> 7437 bytes
doc/src/Eqs/fix_gcmc3.tex | 9 ++
doc/src/dihedral_charmm.txt | 26 +++--
doc/src/fix_gcmc.txt | 55 ++++++++--
doc/src/pair_charmm.txt | 6 +-
examples/gcmc/in.gcmc.lj | 41 +++++---
examples/gcmc/log.24Mar17.gcmc.lj.g++.1 | 124 +++++++++++++---------
examples/gcmc/log.24Mar17.gcmc.lj.g++.4 | 132 ++++++++++++++----------
src/MC/fix_gcmc.cpp | 16 ++-
src/library.cpp | 2 +-
src/update.cpp | 2 +-
15 files changed, 292 insertions(+), 140 deletions(-)
create mode 100644 doc/src/Eqs/fix_gcmc1.jpg
create mode 100644 doc/src/Eqs/fix_gcmc1.tex
create mode 100644 doc/src/Eqs/fix_gcmc2.jpg
create mode 100644 doc/src/Eqs/fix_gcmc2.tex
create mode 100644 doc/src/Eqs/fix_gcmc3.jpg
create mode 100644 doc/src/Eqs/fix_gcmc3.tex
diff --git a/doc/src/Eqs/fix_gcmc1.jpg b/doc/src/Eqs/fix_gcmc1.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..158cf8b61f98390acfc353032e882e26fdbcf8a3
GIT binary patch
literal 5617
zcmeHJc~leGo2?KqfFJ}>Sq+Au$dZJu5hO&2Y|;p#>;fWdfDj>IqC~=?$fgJ=ZnP+B
zqiktqRUqsjZQF_<ON2lIw1CnifDk2+33QJ=XU;6YGyjhFRlWME-uu4y?yWj?u5eH|
z2FM?>b+83MAONrzH9*J)JZ&iW5CCv;0(1ZXkOag)3IJFXfkX|shzG<sVgT4Jx&i?B
zF6f`pyJFu)Id{R|V$AwPAzNx!SZwSma}+9ygbWBGo(x6?5(y|uz$ugg5{&{Zu#{5)
zff2#68YhE8@lhD~TzwN<10RHedzd((olfC`!|?WJV}g&KJ$x+iY($`05FBf%VL>sc
z5Ka+-V*@lOgvcn8IRyjXC~hu_*B_(c8XHq$BQS7JCl?JIF(z2U1Zj*!!$sLKK_TX@
z2W|h(72RRre{YgZCL_s4NMcMV%D~Lb423pC85-(~X6Td7M8yVB^rJ|c|6bu>Fexww
ze<~JFjM7-I7;usp7mI<zQKA)4e=o5fe%q_1#&>*I;JX6f75J{ecLlyH@c&PNf80)R
zl*r<cMP3FF&I0>f!icd%QW)`+h7l4C?6Y=o5??o5B5cH38=-ukde%@j@Kk&{QR>Dg
zzIB4I8GuRx7=Qx?!2vNS2n+=Y+X0x!v`K(Au#qO}Kw@BViA@kmDe29k3Dxo<cLxTG
ziGw91)|UaD6O9AnPzi-y23DK49tnWJqZJL4a_&iLT0d#o=GwzW7@dqsmXhASLrGah
zYxf@Qy*kDwKbY<_GshjYv9+^zICS)wo4bdnmv>-La7btvKAaRAcbZH&b2jDth182_
z=^43sSFhz?zj5<c(ftPxe=dIXxa4VNRdr2m-LvPdZS5VMT`#(SecRjju76-~h{@)R
zPfY&tVQPAYH@~p>d1;yd<?A{xk>&i0mZ<-W*?-~%74Z@i7YB<&)_H-%$m_(R;u5<I
zHYr#gfdoWvg&QVGDq82<d(tAMY2?b?b~2_%dOO0Ht;JiXw!!RQBbNMsV)l*LKfK-p
zI6!OzU@_5K0)s_wN?Zg9h{QS|l8_BZ{tf93Y+i@#HwZ;VK%yF=<V~VORth5ZZR}qT
zgd-y7E+=FHGGLJCV**0~EU+pl%FqV>`!nLnu@E=p0fcs^gL<mHvgIF%>Ne3$a&fP>
z@o%?yTQ-hoZ!tY@^6`Cez7W7g3|4~!4E7kP-mvaMZrU;lrXL{rIgHV>SbKQ>v&ua1
zSi>#L+*{bB*z(<939h}a`JA~~mHy@Gn!1{J&)gm8Hofv1uLC%!jN>l{<nD(RWdPd$
z#VLIA*RIR(4L0nRbEn_8F01@h;Qea0Y9OCU8n&&?d++lUQE3MX(CkvR&-k;%%~2bn
zU{vW}0F%aYrfDs9mS};szG6O1TA*3Pm+r_R>s8J~F|V~<oYiV5HAaG7$i4OkR-tRX
zy|k*PCyW8eQc}vSf6&5?*Yl>43)F@`=a;GVPqKIg?`8!3d{Bue6}S3~fe{4tMJ)|d
z{8#p{{Zs$AW2H#M)R}r$7X+j&XgLnz63v(&l$A}s9z7tTcvLMh%c}<1;`1|vZg~~>
zsyiP*#jT(SM}&ZE?Km!RlXF&8R)Gvj2q^TtF5kbVfEpsNWk+<<_+EIJp!}{}zqf9?
z5D0xwhxiIis;)N|%G1}h@mTzGDUQCQOUlt_`OfI;>BTh}N-xg;R8kLeYesygN8T07
z(CxOw2y9`JYtobtoUpL-ER{DVvbIAZ`&}LDPLK#t>|#G54ddtnxfkXK=L&)JPlc}B
z{^HP5Z@Xwr>9FUp9Z6x>`$dlvlrnVHK{4aS$AVI^2V#=+Oj=ZFCMJIAOkIEM=I%Hl
zP<8wg3fE)Y&50+F+f}*=nm=A<zf}4;`PV-`&L%ig%P|d+_67hoo?YFtnnxSEhNfH@
zQN_5!z{3&Smv?%3X})W4d8ME6d`!LPxL*tO^}fRqQ-wpMv8JS0{~_JxQ87XAeLS=8
z2k$L~B#t`RBSY4DU!4BIORiym5+)Ti_L7Y(0tPOXc-%DxEvz!05O`P)vDJBKX1Sx5
z+0=%G@WDrUkn}jO54V^wWj{(d&%&bnoSF85*TYr^vd&FE)~v6$6(U|~5zBQ>7l08o
zg=!%nJ#4j~GS|WEY`U-iEI&{#W$}G&k8xbiDES&RUiBUxdO|$>qr&;L3oD`RuuTkk
ze&wN_Sd=0k)^_&y+w2$p`YNwo-V%`Ii-SE!w&}z3ciV4sHFnXI7Av`&7#Pt8lWUL|
z%WQOCX}P2N)>!qY(-W-lN_;-rtD6?%XmD%t-a#XW>w%3gfAB1r4+GQbYO8nGv}Ujn
zOKE-^Cp9m0%e^z~Mnj(|LX>juUSCPj@#0sP<2^{@zPL^mkJO5}28p5}e#-HzjwOP*
z5O5*ruui<y>1Iq5Dt?RC)Ei&cuY2TYRF^<W7+f=5Ev}%2$feb=5d5R;M(m;6wv;*K
zr-tz5+U;D6G`fz@%OYkO!#lTwI-zCt{@E{OpsyqwY8?x^T9JsLs&Uhsrz|y`^S<z7
zJ$l`p$yx;p{_R(MjqQ7F6jMFNuOvNlw|?a7_7tq>YPU~a8;m;6vrDC-V|g53zb*&!
z#zk|#!ljQlW*Gyt==$#Jl;46ELh~XbFB#6JB%X?}{^;6U7E)c?_&k~AGwsw^7W3lK
zz*cd*OmwP7{drqLX<oI%lob8<Ii@<69iA}aU%A}UaoOk6{ucM=^!COW`DlgJ!c|rK
zg0EE`JZQFqPGP67<6$0MplHM<Ui=(?+Ytj9>q<5e0%_(QBNL_SRW#;7{Pc(YS8`MX
z)jW0$Cq8m`p==ClJ03=ZGWH96_zN7lv}J};W(z~&-Sswelfq1TMMw69*~mzVSDjCp
zPL=PO!_^#L1ebN(4`;DU4UtoHK`ncBa|IRIpCIGkHAYM))8sAT)6Y?DO;UWl8I=p~
zt}+uIa<*sq$c88z?TkBTJ^EW??o>zsQ8Ifux{5)Z*^{0iEl}eXrV|9A93db^i8#)a
zxSM8GGUILYaAeVX)9Ch>=_k*4ybOGX9^E@}qPKc9SqMNHj+Bq`ueY-#xc%j&Gx>!l
z7=5U!Ovl@$r9bL?+3JLFLgWzKRF$qZgVm&gdE(F<o8CHikzv5U4r_CUMlc)7AK_JK
zX?oSV&yb5NvVL8M(<kD|xoE$odEK-M6XhQ(Eh_5_4eN9)A5_G*G8DM<q*z_zX38F$
zUIDgf=uBZqZ&kf#S^2I7-v>0LlAET}{QmtXWGWS694b;1rL~I^vBI8>pS;-R$B#Xy
zm;~07$G_F|pGG^7qKz*~si7Z}Ph*lFAyzylOh4`5jxL?&d9=xG6C`kp+x0bDp9n|;
z)in7L`~XYooa@}78?2G$D7~mH_X*Hks|Ag^nmh6eMw7oFVV_GocB=GPgtC9~Z(;4O
z+40=6K7jOd`1TH#fyu)Of=ticrE3G!C|kD_nF;wbB`g^(IhC=K?`+34O)2wmdb7A$
z+j!NoM|l`Wh)mM-oj0<d(4HK*IXHPEAMpm#i50`1^r0vF=_KPVHT)=rt<2DjYMScc
zQ26Ox)i(vW(&_OZBg&oKjGm69mHUl{cSS_KGG+zVN4n7)MIx*Vu@(F_OYV@jHZX0Y
z=mPRDXyP8-e0y{&btgBwS$5~7{%2TLg4r7j-9Ho!8ykz65`DGS=wr{V^w3{YbbhEs
z#L;_77G5h-_ws4uSQUyL2bR8gzKI7(LNIQ(Efe*A)qBua<cah=n5*Y{{n9Zx6IcK!
zFG6O0?}|-<a+x$WhMXr=iGQD$Yg9B<<DTOs1X{7|rcX_?x-Yj~(kS%yiqmbXIhdrp
z_9`-R+z5STqD$VnY)xKZ@@@_0P02jM^Q+yFCO9%K8|6#(?$a-8swBvahSol<N$&4^
z6~nV41mneZ5O+L6+;a+%2R_gI<wod**zi$telZ98=>rek%#ixU{7!v&Sg-e>SNh4s
zS4P)bzvMMUqYhWhU_MTnZ{&dIQEcbb1SB<NwT@wy@Mct|!Kau)bLK|!s0L3neS5U}
z%Te2Z-RJLm!Ce-2I3EUI62Ix?Tl-~rEXcGPx3YzjOoZ`EI};H6D^wh$z0TOafh@6>
zQ5qe;GaFm5<Q??Jw@W?sCRwG6U{dvlWlwl__4nUb20k@iY}(ua;a0r#=oNq-*gdeK
zN0>4+`jQ>=%1MuB^i()(X6WhR$zvzfOG>sKZ~^9~{-fan3*Y{4ChUJQXxd=m(7yqg
CKN|i3
literal 0
HcmV?d00001
diff --git a/doc/src/Eqs/fix_gcmc1.tex b/doc/src/Eqs/fix_gcmc1.tex
new file mode 100644
index 0000000000..c4b0d62527
--- /dev/null
+++ b/doc/src/Eqs/fix_gcmc1.tex
@@ -0,0 +1,9 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\mu &=&\mu^{id} + \mu^{ex}
+\end{eqnarray*}
+
+\end{document}
\ No newline at end of file
diff --git a/doc/src/Eqs/fix_gcmc2.jpg b/doc/src/Eqs/fix_gcmc2.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..d054f584305d763e5b7baa8881a1f1a0d96f0a62
GIT binary patch
literal 10626
zcmeHscUV)~w(p{MBV9lNDS}c$??pfaMWk$%CPcshAqJ#}VuA%hS_A|Hl%h1HNeL(*
zy@-l*kWdrYs1!+1aS0)L(R0o|^?mog_uc>Q9xHQYWY#st8X5C9#-zWeF9Q2dpFDLE
zU}6FQYsLl8iGZ^u8sP%~wzj}A001}vW+nlEg~2f~E+8HOu>Rr!;0WUz09Xo{{@z;1
z{I^!7!-Xt=<AK{hJJN|<hy6lBuNpw10U@fMmjk`vs@{PpDBANXR9#gK3K*H9uX=hT
z;i1xAa9>1#vCL*myNopAvayVlwym1&RSUQu!a6D#{#(@9bKX%%Z~e<MrY6!xXah9r
zDheLzDUC*52?#Mj8_WEf+<?LVtcJ=+|LPKoG?sC;wU@RC42DZ<t7@sL$uLF-U-mI@
zJaO``V;R4UW&Zk*5fKrp5t^!j!M;#+eSLkXng&!uLxs^pB_uK+)Dx`|5F-1x8BV}M
zyn_)}LlJ=i(m!YP^a>0MHI|WqGFE{8b%~$sAJ1wc{U`mYfj>3yrw0Diz@HlUQv?6|
zH1Kz{0}o(8oCpS%0q7fmp1ogSXkdt6;8kf&RW(4*{FE*0Pr${XU%dP;cK`UH<Y^pG
z$NDLj>;5`uj-ht|{G7mTro${uG5|9_6AM2Ry$6tDfHpR!U-WAj<Hp3y!pg?Z!O6wV
z!{~tD&p>x9EX=GdY-~T5VY<d>2Uz*p1P-g4u?sqQa>xV;X<Sb$;*>S7?G$z#A<Joc
z1>fM}77;xlCN6(OLGkD@Eo~iLJ$(a<6P72ftWTZ(?cDhbPR=f_-k0G%zJ3V*kkGL3
z2y|rB&0BG|;}h;Arf1yC%*wu>ll$m#aY<=e`ID!0^$qyOrspr3ySjUN`}$uEydE7J
zpO~C_KTRN#mX?2C`S9`6DrM`-*Kgb3!5=$6<zhh2zlg=S|3$KYk&B-p7c(m>3oFM@
zxtN$EehSXd%63?tUBJwN!!t-wM&mlCka=2BZ6}wkrXyL{D|m!kL{5t+Px&d@FUkIU
zg5CK4lI#z`{w~)HU;!}yA{J)GTVi2hyeU=&v2n2dBo0oFU&Q%W;{HWEKZ)-TqBAC8
zV$8u9%+C1madB|{(e~eF=yMF@E=VT;yev!%V`AY4Oo2V>qeMmEQDQem;f|IR=f@_G
z!xWt|@5i4Sub-J}&U_o5l_q-W&Wi0}&g%~|%@op)AH1u`5}f#DlhvQ5Tn7e~u3-;U
zhe{ZS^FSKbvJJFGRqqHYDsNe+ObaM8sYO_$KRgKWU|YndNJ1$K8*WvgL>FQHI>|0>
zQj+}{g&T%Y=kjS9$J;z+!?(=$hiN(SDC;!oat|E5B>RZixOoXuQ~oFsQ2Z~;(^;5{
zs||{6qQ`z9^5~6exv-PD_PC?gozg9wOA!Q(@CwP;yG;IxN?)>jP4UqFb~dzyIXFbB
z;d+lSIW%PyD0*wGBPul~0zulWpLKr&XB&L^>1A|u#ho)=zpKP0%67dK2<Mi+3K@B-
z_$Ze9Uy3Px_@^Sb0&%ANPDH6r*g>i--~POASN6x{Gs+u{`T>jA2U}0jfiagx>2Jzx
zKool%2?_!82fji1{KKR0#Swypm%S-H5$Y<vl_^3aqA58i(yS{}U%Y;H?6`Cyn72)7
z8-~&WUmKjn$M4}zP2p2EC4S<`-@6V8+K$JZ>_2|Z@$>Iosp_6=-`)*eyidBjoI*VP
zL$_Pv`mO@D^VlLqW9}O6TvEq(gSnSOnfS2_6n#Xc@c5F7(D!HNA@BA(O-R*dUk{}N
zaRv2O-J7?q--tYU^ijk~y!WPcO^|Q->gzC``z!XREI!PhYB(74Qm+t#GYYAs>J`w$
zklm8<<5K2(1);mwkb@b$P~ohc%BHxgM~&q1>?VN$As>l`<W-TouYyAJY`VhJq<nYz
z(Ni65N5BzEH;lZJa77y?fv#{ES8#I+tgA}J6{=*W<`N!OR62?1Hn%o8HQntuIwjX-
zcn5MUWCxCORUZUYgdC2u3HMPguN^Usr^#Q#dORmL-fy;~v`(QEY}`^*7Y48zpN}p(
zSr=aqkug+agR|aDc3SbJ%E9V3_avv1=>WNI)bxkjajD7%r{$@CtA`aAaQh7Vf<eSj
zF;G$xPW;o!zQafoPfPh`fVN`4WK7fk;-FXSR&QS?YQDGUC>*gAVp@VVoXF>C1cjFd
zZAeG5a;M_82AS!A<G9u*vAgb;WnEJ;YOh{(<~(9L{miiTChQ2Uv$2hBFOLp*nZ`HQ
z^t_C)B<G@Yb*HMs_*{4ROy)@_#5`$CW7fzIuH$ifN#9U!@jMPZ`mUu4GF^lb0UdSn
zdQ62|?TSF?=@a)(xJJ*UT|J+Zz^pIwzN?3FN1xQ{OZ*T^<D(6mh=B!U*mllX099#e
zlN%A`Ui}iOH90xyAxa)BO8jD%)cqY3Hp`D#JJjmv{;n6Ee<o+bdw=hKj`4hDP_?}Q
z+CTT33z)R`GH}rjKuo$`h%a-4&3)2o?yD-7Sy)*3kriyFGT5N@@?q-<Yj4exE;U0=
zQ*pG!;)k0ia&(~6oF+{*9^1etv`J9a$W^yQFND($Pw?U2wU__+%}TCnyZ*=t$05lM
zqJ`UylOY9DSjfkDNwUJtuM3cmomU#5`#wE$A;f1EJ1bWfIWO@kC%^5KG0XCiTawEd
z#Gn`psC=C%5C);7?2V%p&M)(eDA@ko{3jyc=lA6yWq)wCKHv-pP+gLIZRizi+)r#@
z-+LURvR46(F<azpwSW_Az~inwYtefa4&$Lu!i+7>Io;VPZp**g#BTwY&efaYI+%9o
z$qJlWM&q_3_jg9s+-w!3c&v6y9wXC%lo+TlyxdCmagtk#rA5c;zCPPgQ89_kSm|#Z
zUz^+Z(}5cabYOsc#B}`{8Ap3SDuh&PFWf2Mqe5%mqkrpZQyTu75_ldh^K2{-FM&k9
zd7@gZ*72-UE=ev$xv<#tTsAf{82z+O0<0(3^m!Z}H<Ben#6axDsn87R_A#wF$H~vz
z!(Fs<H_(SOmlKK$WCI=W&y6Us`~|Y(!fgyUrKG>Zp%M|kR*q~Jp{f%rsxEL&9y)`F
zXh#Ull#tjS%Faf_IRSdSCYLx{cD=!n?s))g**)U#!8wKs<@&mi5T1Xq-r({nzh?;1
z*J*Y>clMv#Zxrbsd+F$<=d-~~JJHS_@)WDvtrCxB!Hhofm#J~ZF6GF)z0)g(HqPhn
z77l2=O^8kBkPAkJx9eeqPP8{HY=lvuB&(#y2H7`g+^Y^~-OgOp>&;qR&xOZtC-6t&
z-W+Xf%D!TM2I$c^$`W<E5Hbwq0}Hw_8Y8F}8zP@;|HV=0xtEParL+AJ2!kib&V-l`
zNyK*CIG1mBsDe^Y<94SgkVH}_x}9yNWagwQmTmO-SBL=mkz#fJk*x1M5sp@fK$&A3
zU1paZ&E2lrnz4^o2YuCDx2ZCz*)^k0eBa@#K(ukh_lty2s`a(G<F@DSY0j)U@h`G|
z36NBj&$&J@n;dl`3dc!R#BHbUWttrESaaBD1$RwqNlF~!vW_ZMr6(H6J(UR|3ChY1
zKVE;3vF*QM3tXEhM&<@4gi1gKKuqYk*5287N^4@oOO7O-gq}15UGwTnPGdDB#cc`O
z*Um=u0-?kLf>P)9C>pu_aJ1E0*+_HIZsQ65<%Rxp`0>8HX!wyN#d8+&H?+1^`Z&x-
z(svV}-?FI!yQGJkJ=keYvk~ojy-AwX+#K4%-saX;v02$h^ttKS4ohDxHEPG9k;X{&
zc5{=Vzz=)HC-TGR+I&XZ14}<Y-!9}haap}6)BJ_Ki`gr_)Z<3J1^(&3;KEPH@JC^x
zkL*!v-kE|&Gbc9hnCt_s7Wd&o4mCce;Q9?X;$r%Klbs>tT@O(iPTueCnROr{`8xiG
zY!xQeDIfiqXxB{##K5sG)BRMtVtG1{V8kW=bzZ6^(#HXQq(q=IxJxe0>t3sWVPyF(
zCvTp!PppXs_;^<=<|4R6j_O}*>V+hy9TGiR*Km7`#!m<M=ztFv;r>1$i|boSwqwig
zMQd5;cC<HKEaz7I5xuT!aB@7cx{ACr)#F-_`0aVCB?)%`8H{Y6e2ntm%kvxE%jZ$@
zQ$d}Xf(nUDUbuVW(W&1f;0Up;Wr)W=?EGgSr(G49)Qpe%{iV_*2*mH2T2g$MIJ)fa
zM3vSrf5P(Ae6w|E(7|gSDT$6`4Vr%I^TpnL(y@P~Xlr@)loZtU3+^E1Rkn!|ShKVs
z@gCi12!T+Y#zJjn`&`d7d)r}+Q}0S1MOA7a_`E1`cU7dvELf5dL2!8->UG-`?^g(U
z&}K?K26j>+x|`dtQ-vu!XUI<RF=mf#i&vEt;-$`y4-szFz9YLoJUgj)!P?y(Qm<LW
zHhJ87`W@{x9jK=RJyO#xod!-_Gj!k}!bE|LiSrHFSAgP(mm9^iJ%8~+S4DXWMu>2*
zZYj=IR2jMdoZ)GY@_oA|1US8d;eFRXP^7kTV_F5gep38ZM*M|{rb&Zt^EOHGs^y)m
z0;M-R_?~QXk4mR{Q^t01@6y2C3Lz<vk>gUl)R4XVCOR=8cJa8URHu>9BwYN?VEi3Y
z`cUV)wXm8`b#EuXE2poiddHqvtqOe6r04rZr>_5?=|L>KAVo@|)q%X*GqoVx;4V$o
zA79trGnilK+|;nEBv^^&X>A;N=G-awYU3cMeZRyyLM}{TYVMvr2N$VE1b2IZ7mq!J
z9wjN&5>ebA&QFxpU#}wN?V3!NK1W}Q^zob6R<|a&aNiJW;eSr^rZKFOx!vzR)+ZXv
z?}rscwW=%ht~6788ler44coJmyXP{;VTh(XSj0iQh)Y>DjVtrZR%)uP>SjAROLiRS
z5R&LEJdPRc8r!w7n5rxt`^E#)v5r|o-a0B1dAH6w?_SzHxaa$uPE_MUzj>iytsx$C
z@Tei5!pZukv2_(0R;xfHQoNVfdf(Y`OK#&hASUN2=G6cZH(iDGshM8r^(D3ygdCjQ
zsmAN%=DsT{H;a@gS!-$7e%cgadHd5i<`erX0qNRE%nbQW2Sx{sCL>RSK@lV*9}G79
zHKXMH+D%K53pZa`9Lg9d?G9_GRFoC4I!NKtd048l0IhL2E>*Blk6RC;m~O`~%oLhH
zrk;N`bp<1{TDK;7WW=p&5Fym$XcMFsebnjvU4n3QJj+pydtG(4=iz%y{}<R!AN^a{
z@rMujM}P7kcrql}w87~nO#s;pi5=~Cp!nj9oSpqV1cx7|x}V<IPy8&8C{COo9W@e^
z@9jw|JRjWB+61tSLU5tpl*1U&PT?Lx-N$&?i!5!6K!pDkGhxohi#Tah=8_D^DycNg
z<eV<Y-0Tsr|L7`NRMqbp`BXcua!Sd&?9rQuU@!A7MgD2wk{HzSw}rK04A15N(X=5j
z*+c@pM$+wPZ2vS)+*z{a<cngE`-A&%*!7Cu!rr3x^r{fvtj0j?@(FF5o?H^&mX$^W
zJ3D`I9W4o5=-si4jS;K?zq@%g`^Iav40=%RUf#>?6U{g02&vtt$g;e^hMAM(265$U
zJEpNT$r#JYnpUMKLI9psTuyOPwESHPYajCY_~|MM_Tx6Z!5m`Wq~hoR9?}g-^w4Sp
z;a76|Z^s0d&g2bHwMo}Bt|>U`TyNw_$g=m!mx<)##q|$<w^sP_oyOHf>0Wf|u9<HB
z(rYS&IR-k}QQ}i$Y^uNt_l61WJSEggr8_>mI+An0yVM_Z;~GZEx6YPJC$7^Z(cZm`
z6*)E*N-$I60l0E>AQ9L7klCH3dy4S#L8C6_XznV{vne&Vi*mp9`>;P+uX|3IcV&|x
zoWI5LHH1{wV=G8h=nfR5o~7)M7;A-=xUbgtik;7b<5P<F_g#GYBmEiI0XgP%c7Lf`
z1p*CH$ynZ2b&5^5zPoPE4!p61H{kWwai8NJT-*X}S0K&$fE<UTq0@EDM4RMTo7w`J
zW&84kSMHJ(WxL9xf`lgWmeKgzd1dC$7%jv7&NM&By&`dMb4Cf(#{(T)AfH|kSS4`8
zVh+>Z(3Bq7R8tR7!nSj^PCCkZMjO3NPL90#{q9^_v}~Cd7qC4<ILT#XxHLeMkMRdH
zgNV=@Q4zO%8t|iEjggMS-<}v<@QoyVlCn{=3$O4LRKBV;TDs`AzQi_!-4#T~lgg6X
z;kbB=5;znUCMAqC@b(vUZ_ik<gNGT_ecCde?!ZMmH?Td2Wb6*j2b4@|4-r^8O~qoI
zo-X+17U0?>R+}$8qywDunMAMwU9Hsup3nNBy4>FSEQdBa<Ep_s5^L?fZ)pB9rbmjs
z6-72&v~9mdjY)I%*WGC=t17*TGMFWRQ45au2PJcOUol%hyIfxfH=4Q@m(AP3@X(e*
zfWu!36MgKz__MYmw-Yit!xxxaBcH^a7#(vZujCZMR$XsDl-1xVZvb7FPl|~h(9}A9
z^a4k?8ghBzz0Ooedp&mDiv0aNm|gv`P0&v8<j%KuMHMi!vYA=iD#KY<AGi1v-qJS%
z7kw3c4{`&0qL12CDZZUydqo(N)!d;T*v1-`8O?IGZnVMJfI!<kg`d<$nTgP7bl`1$
zI#&RyC~kd&Bx51=P?i7sdoysXV}Zj+jG{o}7Nu|yKNp}XEy)YS__!<5&Gk)`edVa?
zusy5Sx%%>f_#6$HEc7-9xA?4p=7U(wI|rB`Rd=tTAuz5T58q5MViHZb(yGl+K9AY8
zKy10cL29OIS-F2A#ke7R?41Zn@k+6vm4v5kx-8!4ZTme&Xktw}07B1`BW_@n=A6BI
zJ)|f#aqH`N&gIY`WKelI`jFZCfoAh1+2Cial7+FDY=?^{TPT{u_PB4m5N>pKZ<J(O
z)oKmv1##5&V-zYfOn&_K0g(f8X_wPV?g=FyOW)k8wi9;h#)@Es!8bhy@wZ|QFPTEt
zazCeE7Txj5d+ItjPdEU3eb<gvqcu$^V;Qxo)5MXl$TE!lVngi16R`|WBh|2wOjDez
zx?NO?hIKpb2oZAB>(HksJ4+|zpMBw6@X5YgLpWpM^{zmDp}E}#`TG!IVR7;frG^*`
zUfAtwl$336ooGOrDCH&T1*kSunc9zUNKXW+$Y($Tm37z=EPWj}lDfSS0BAvKZDt@3
z$aF?$d#&%yN3>C}h2pN<YPF;%x?;y|Q3<O_(9YYlk%()xiMZ;OC<;%vyjmM@Es;i~
z18k;XSs7J*4@b2mx~ZHot46=9jM4TKFvwX>MF`lueDffI@cO7j-T}7pFi%(ns6ZTI
zWHXwtH0yRVKnhLdSw76f&0VC6&qPdKTJ^h4Xo?FNm^IU5$^D&bN9I_xOTn+<guo$U
z|1JEU|GLbqzd>o7cVwhSPm(xCzMscQ?!>1ryyJj<Kp~&hO%KjzAcCCT*}m`VOv0uL
zRlNwYV?Nv8bY%n1Y`#TpSG^pn{Qk3TlwtFU6uD&n{USczvB2u7ua!_9v`tqV`*2mO
z%3_s-9hsH{U%8fdc67t5uKwm=c$S?i=XqJ%Eki?;Lzb7FT<&CBZS90|fY#vId8HP^
z=lMKIkt3HA{Yo6GQ<J&tG}6VarElImE_H}1NXfowqCp$Q3S(r^(Qc#P?zf5|N<fZN
zZicl%azT$C+MR3WJP?<-`k5rbFoRcw2;Q0FunZ$umES^rJXP4qnsTV`F$%#@b*je!
zcg=HOK9A~cghgMfovDqpdp(LTd3th-=|y!3p`Taj>sTklboi*c<G)YR4uOzv6OnPB
zCY@mSkIlyU1y{G<zFrYYpI4mia~Dn&Y&y12gUtubV895DC6$0LJ@ZfyN#qN09_N;L
zY2i40MfLM%UO!no9YkTIl)s|G^9PE+RYs<aES;Qc_CX6@Q2fc@!u1lDkFLTG+a8vF
znxwe8OcPE#bDMGND3VE3ynD>*pkLAXb8P1snd2(%f65us$N$ZTv7~^BJ7=&E>ggTY
z%OE03Jv3=4vYK!<t4Z9wJ)H)h{q0qC+9~Fh6X6%QUvqIBsChdCTDQX&c8qJp=l~I3
zjNNK%cR{EbE3#|_(gdK>5i}gK4-2{=q^MOjqjqtcewtK0^4B&|ie-OF=wt8wcnaIx
zctn%<J#U<eTY{tBa=7dY+)`{^Y4~Rp@e}m0eCeifhg7dwZ2Q6D+JsU{DUjYk!H*c8
z(KLB}k^JrDi^4u$K&oY<TCI`SS#*sU7`bSA6Wai}jZvVA?v;klv~z>IUAZ3F$Qsu{
zZ9G0Zj;iVB)jCEZ?=3uxV0De^4;0D{sde$qI1=XEqi*@#w6|t@2byGhz(gR%W<t9_
z!3K24KJ`y2!avE;R+@PwoBGaMH}{4}U6!dfB>M=-`bX5%1N74Ni5C=DFIEzq@9wRA
zkH(VP@*W{1HAdlM;?Y0WLt`m6D9`f#2%sU}zEOtf?)F&Q=%B=7A@$g6p2R(g?s--+
zCXcvxKqV`VIH$v7Xus8v#Flwa>7ZvxTF!krYmK)(4=$%i)(p0T@@;t`Bt{VIxL2Z$
z6s;dOfr4Ltn~<5`?Hog0#vT^UxXw4&JL^4db;W6_K#C?Da|vuA&cQlp++bWOTB*BD
zdZctTOgN((_jJ{F!`%4;TJ>ESQb<WowVY!*XZXgNwA_^sSU7f?4v--0bBhm^h>Q>?
zRjWV{G&Qe(5pUOcvpXH(dSmwDY+hCtvS!1uyx4uz=AOvMw-<Hjz`fhRPTPOT@6@!c
zEkg<Q+aVxs*XmQQ-EA#M!el<@u-44lA<rThgzv)EH5<Fk-{S}*8_-K+@&3%fuD%y#
z`PKPE0P>)Hv80`sD!gQai+$`ox-n%xg;4>$yPaGs`6~_Xk8WZ{<(q7JcnqwbZ*o@y
zR}aK=@JDR-lRw0_GEpM7&y2V7ml8e-ern&pVWjU9`YL_o(kL`jCVx;sx?CE))dF-l
zb_PBFFhSF&y_J$M(J$-npTH=6Y_{z0-t+11p`7wv{1DbW^>x|@8J%@1EL*)pChS7f
z?92nrA9*g`sZyz?bwho)lzBz)8VR|XNKn0tIk*}wPJ0*_xcEIsVG<>qG?*Kn;jDZq
z{c?CJtD<+xwpO~q9Y2=OuPG=;7(v=U!rQWm>@r6Efq4~BHqN@j*Z<~=`3d_|^%6JI
zH6H|d;7j_SSd_Y2R$Vj;E0^Q1n%y4OZV-4U#3iZd;X;t}890Zye2kOreM;XRfVRrp
zwWR|?P1jAyvW|3MNh^|(T%zux1sRc>Hk796u?ahhA($S($iy&#Y}GL^Fo>~v@4BOd
z-4KC?DMGp-oyPdL6<2DzROI7&HjGA4$}^m=8)xI~_stCMmC}I^GaUsW{}U`F*Jyzd
zNW0d`MF*~nXEV|k)eKLNRY4kzVC>alLy`oy_vk!_k<i_jUZWF|dMe+zHby>srCg3Y
zC>LCR>~^@`!F1$@mv*EnM$8}@<3)kR(zwyy-MW%u<WqHH-%!||+_=H#87?mxZKhRJ
zWa@8+XS`~B(WoK1)m8s&Kvsz9V_V>cE*&tLXQcxjDK>N<d1Qq)hW)a!Xh-E7ZwCis
zD#%mY_xH-$j>Y(bY2=Wu`GX)<u7Fqa^}>aTg5Sp5`sL@gYON3RT)C{ABJ|`i+`dYT
zMPK+Uq+=*iN<8K~MS>*NjZ2=F7;Zq;uy+>Dx_E$5Tcqtb_O_m?bL8*pW-<*X<`qgO
zi$xSN72h(A|M+KJyZhdRE%ErL;}X<Be`CW}fr=kT=LUkdQ{FVFJx>vQGu^JI|89P|
z#^xU~DU7oeQ~JceIq3fo1^#oiSo1~-JZ9GeVgKYI1Uir_`lAX*sm4pNL(Uh$^pj40
zGu=_L%FVH!;w@Pa?G~A*Sw8ka&G;(yGbO$&X(*+=X;yS$f2*uc1p^2!bPrr=K>2jo
zJf24qe(yY$b=crQK>T~>%Gf~u)=O#}LQMEd58AAx&J>k}BH9V_Bp1Ztp+YWLe^6tG
zFto9qFJPF9e0X66YjAS>$aywYr7N^flwbeFrO~pfNpY(i7(vk2ZTw0R&U)1G!v)Hb
z3Qfc3Up9|q?d*`O(sEuzg!hF%icJs}K8FdR4dlXl3*tlx3+o$&V02HENe+q<(j$IP
z)b(~=C_yKUOcOe_;rXn+l@F;uu6m>E?#zJPM}_APGr7eeWtcfeg7<Vgo5^9SBnVr?
zii2f+;-7Q&5(=jCg(m_UI#EyMkYUo}0t@TusD>;{=A&}(JkE~&yT|olpnO2m6Ms6u
zZsbY%(i1*xRTTrHG}e=j`oveb6%TSI?^>Z%U*oL4*Ni-r>N1fW4y%gpaXnv?Q4Z#=
z@C-VJ$fb+@-q@Hps?6czJaIS<u<k0|7)trPzMzZg&30AIBnx&)@iH>{Amg;5jgL&r
zL_D>q2BRH@vuw@%a;oRq+nVB3buU)d9F8mIcDIpwcxq?zKODAA{=XZ`Kcp-F)v#FT
KQ~$&BP5&Rn=r-a2
literal 0
HcmV?d00001
diff --git a/doc/src/Eqs/fix_gcmc2.tex b/doc/src/Eqs/fix_gcmc2.tex
new file mode 100644
index 0000000000..fc4d90355d
--- /dev/null
+++ b/doc/src/Eqs/fix_gcmc2.tex
@@ -0,0 +1,10 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\mu^{id} &=& k T \ln{\rho \Lambda^3} \\
+&=& k T \ln{\frac{\phi P \Lambda^3}{k T}}
+\end{eqnarray*}
+
+\end{document}
\ No newline at end of file
diff --git a/doc/src/Eqs/fix_gcmc3.jpg b/doc/src/Eqs/fix_gcmc3.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..e87764afd9181dd6c4707c0a33b0a80abd2170d0
GIT binary patch
literal 7437
zcmeHLcT`i^x<4TxO+bnu9fQ(CFVd0`qzx)k1gROMgAh6aLJ?_#fQo>?PzFInL`0g@
zC{=VEq!@Zh00jg|lya1i+`!Cr)_d#SyY6~_y|tddlYR0fd+%?b-}ikv`<Fu-r_BRL
zel|8U1|SdsFr^;=jR0IS!uWUsz`1jP8~^|;00V>vfYL1p{Qy#Y0LB9w0Ho;e001q8
z{M}p1@V$4j6#CtE*&j$Fupajg4Gq*$QVIxBbah8x_E2;~`zv8w1C>-1m6ZTJBqq?+
z&CerL;<AUAPk=sbqYVd>@Nw6N*=d|pJ{M@{;q7A@73^Uhb>7A;%Fj*P9fm|m=wWm)
z{(=4;p{^1b|0@9@Iv9P}fpZ<Yz5lfmOyXcjsGmOU(zy!~hUj1q2@OScMP(R$cCfpr
zj@23Czvj|E>BIiI$cTst#RxS;bg-9_ing}4lCr9jsw$j50v-|>5bBD72ZTue)595$
z5Vv5Tz)&A_fW*E>*URXzP<<FoiN1o;UrX$p->(%R@dH01@FN00BJd*uKO*oW0{{Oa
z@VDCW2%t-x2)dR5XzReK3*P8Zbci=PP(n>n88~HNc8+mhaM5vKpFA*+jEdfv$O4`+
zuEew6T>}kbY0m*}7Qh8qghF5d12+W94Wab_;&jo*1UbOLH2Mj`0A*x4#LU9V#!es5
zc!aL*pil-zC==8EGLRU0KfuV%#B*Hb^dVkLS7ulcpX$|Hzp+Rf)V(-rHAFh0b~!kK
zm5pEEn4r)}DQOv5Idu(9ty9`MhG&e7O-#*xwzj!wYj?@s!Oh*n)63h(HzYJHJOUFL
zbuBR|IVJUa+U@K+Ik|av^9#!El~>%ad{Fi9X?;UuQ}eTy)~@cJ-oE~q1FwchM&IGb
z#wR8T#D&GB_aBy5R>_~gY;JAufM37u^P)@6KWWj=e=_@TytwJS7#JC$jLiGIAPf=v
z#JL%nj;kEvIc>@88pI1zy~@I8aO=0a7p#(MR-~hsgNNAoPpA`4lJ}_{F#FGlCH$|L
zeJA!euPML~U^oDjfu1EOl%6R@I+&Q5_QA};d;pfez;*!keQ<mSjqU_O_d%b0i2mYW
zWoG@}_a7Uy8M<=kqD=yap%8j8LAe1WutzCNlL7u^>kQg^pl4qTq6yTja*T<6%UcZ^
za17JkMHPI4tIre+aCJ7ZX%?vAu@D)`b=FsOgbOpVIf}UZ#dQ^uv-bFrVsUd~Vl-gs
zBoN(ube;x~?o-@ofKgU|$xeR??yH#zMHY9A1`N58KkShS8&s<c%<YygTr_|U-G|(~
zM2VyUSWGkxxO3<Knny#zTrxuH?N+t3!!8k>zV#A5k-qXRZMh5|X!H2Qu4WEe6QY|u
zV6(_+`V!Ic6BD21PvW$<|2E#xrIK_S;7=A<=Orn<2s|y)6y-m5y&aAH6l@uqp|vXo
z2#ikqJ(%DAJV@Bmz}_R|)RQgNix!H4PDWZ>)nhHiMA6fe8|3Uame-oGnmIMe?yP@T
z%pgZ!xTN~QL6&RpApON}Ai1dkg$duOL^*{J^~>VaZ*P7AW(v3bR?MDnJ0ThGV9w4u
z8F_EgxrsYRH>4xBw)Qm*0A+h-WyrW>>QU;OSx5r=95?nn`Pah==(7o3AI=T&t_z^9
zc*jLG4Vj&6e|nb&U?h#@a~a~11_(jlPKB{pVbIQryk60fa=lShz@iNghg&Wr%v~s+
z5qmF#H@VeiBCq<u%kAobRHreS2BhM=Yy4x?z<P7C{)cX1@#J7N%w_o*SnsB2_#gci
zX;k9J3$2|Aj`^;+)(RXaf6cABlE&$ogW?NAieY^~^IaRPAZYOzB)5lEQY;(EX9~7*
z@Y$}eoOiw%mCa1OsYo|nXi4V3dAG^M$)!0rQBcURJ+fronOx#TM)k<!62v4yBXYD8
zXm>x=(W#Z>l^mcLN04p#>b}b`lkZO{Qq;@RCrCL~K{6a<ZWYW!KV?m95RI$I*XHA{
zcFOpL%ido_NWWccl&|*A$>l4?ya+iG?scTzqip*3@S_9z7IkBVSkciiIEsKu3UjvW
zeoYXpRj;-m>HZkt&za4bD#R{*HQ{b~Q)0Q5)H}A-=W_ADW&t2mmiE2ThW4w8G&gD9
zH@OiX%Z4emC#8Cl^*5E5V?}JHo_%axB(Jpke0}L$9TCEdh>~!djMw2xED`J!!v>Mo
zchX{oR<A!?nGC3nD|7n!8YE!ztL7<@o0g;Hb<JcWn3mx4Ejdc)Ui_vR(t`#}rY|7Z
z?De<Sx3fevU6e>~dWmI`Pt-5;2lgbLJ~1fT+p1zT(TYC{xnZ~F4DeYR0ep6{M6@U<
z{=#zt>qZonKdn_}ctt~ZT;caFc&At12Oac8M}*s{X|+*1B^tm!IA5~IC(a%>ze+E*
ztP9LlLoSgOSZC0(*4OqMB4xpX%U5gXiq$NB<YHrS<Ryl2(_vfBQ;hS_=TaS<@oX|>
z@oe(1+wIBGFN!(9jCoX2RQqrVH`bugI?J^AI#emI+`eOJ^l6<0OzvT6pMB97dxNf$
z9}#`Tr<jK{TSKhzy?D)H=utozd&{NxDyB#N3OCH=fwXe~c4b$Eea3f}xYOeZUcLd|
zBic4{+%p-77UwWQ6UNiUcFTvCf|kWPD{BjpZgZT2gp!FALSxNF;r5=%>83FaaXoI?
z(Xi>&h@=;bDlM%yZ;g{tW^W-zd>0@n=p3llTa2Iq3F$tuD(^O(&m)<!ieQ}Tiq4Jg
z&ql7&sf`1Z+SE5{GG+r$f^Rc=j=zNdLhJ;~;*oNl$5%=O3%5tXlm%z`?$+v$O|vY+
z6k$^F*U<Fz3|4nMYl-E%jWEbzylwx_&%Wr5zjx0i@<?WC3*OXPEns-?@b^pc5Vv@<
zt%bBOKETrGzcdJKc#l$>gOD0PFVM4ZPDtsR89L`PzZ+vnG_GU6r_D1}HQeAmPP+E{
zv-Z}Ovx^&bjGU^7>e?c8t_kG4Rx(1_85=#Mryqh|c)jxJv0`I5KFrL)v)$@;`AxI%
z?REZ<#MJ&J<D+_wD7a=kR&6&r5#-%HPaVWHiFm*_Kf>+Xj}oJ+oYPDyQLGM{+?*Bg
z%Fz?>y~V7@hmcNw&k-82^g>~b^Nv(-T;NX%As%aZ*IQHfZcE<$j4O^OreJvWn9v2$
z_?(*q4LLK-X`34W$kmNHOp$k+uVpa`J9Mj3{$zt)-cHn=x_dlRu|v`7U!Fb=OdCFu
z#_)g*07NtPZ`oSdxVH^OzRpaI7FRD|y_Mnf8Sk@x%HEZZ%f5-5iqKq)sq`MqZri&s
zmz9UqjqBnjSbS5dMpBqb;xs^^deb)z07D44M2fEKLXpt+$hL74Nr5!&p=Tl@#%VWV
zW^ApjDe&1jlhu0hp5`%<vEL^hnNo%{nNqF++$0Tnu2sT;^v$|j(kQ>$x4mNprdME2
zm?*?r`AqWHcx(sxJ_~;8U$|?${OZ*Ssb!NL*x9}7yPnglIMN#zb<%d?b*YL4mPhTv
z-cLo83W;q~F~?MBK%k32f7#gexEEG>`kK|dqL{oMhq!{)%4=Kt0zE>ZPiVlS1S0!#
zsgKR}=lN$j_2Y-fV#s-@x1F+J`gY*k(2sey`+JMw&Ex_Fjuxqr1`gOqo22R5OZcyk
zIHUOGbmM<Q^fxp%*EJ@;8rDc-$KYk|StgzN>8K6YPlx1RQA|jKy*i&7LmrM}Ey&^5
z(1AAhiz-ciYkVYZBg<Xb0eg>C5zhSTwD{{Uzi>I#rc&i7Zev@6lZD;LrzOc)&sgz+
z7(!LMWki+3rS6#O)r|LVL>yv^(Z96RHFEc(HF|qSy46|}O>@rg0AlnCVT<T1*|cw`
zD~VXF77chMha{W_N2uvqS;&pSWE$|*2qlM6TR``r5|K?ADQhKUbdMjI*?{~$E2I`<
zfy(W}cVDdIUkQ=4pRKtQzCHY~;?7+ArHlFbT)^FhJy<Ch3k^upU9kA3Fsnxca0p?#
z-Z?5xs>J%xl};gD>3rypqm<P~l*Q{wQu#JAkjwR$FH{PMrpp?exUM2{<sLtpE|f|*
zXu#*Rj@dMTlm%LMs?&h7Hw%NiUkEs=rRD#?<EuQC&4>oDUMZ07%3vo9W>75VH~2B#
z{SFa{YlYq4-t^`3PWe{u!J;4dzsko?_ln_P-Q6h;udj~1={z&BV#|v{{ob^l8`v;%
zPL1G}c2ePwwHby7C$hInx%mDgN(G$)ATLoX6Cvg+-k{aHaRZ@AT2D#i`O=#!#Ge!5
z;pRrz;B=w^72{JsJN#jP;#tB=rQ)R=Lp3S9yBac++7-(KPVAnc4h~Tq!3uJAKi&=@
zF}lUF88#foek@m1=MRheLqUoqiYwF+zWMCWN7)mZwShfdl2=j~nr+g_(U+=#$6oqk
zDk_c7voJdjxVgtoY}=YO6rbPn95l?7udmIIbMGf(OW0~nZcr>^7TDWR0+?p==64!O
zO6ykk%aiZk=QH2EY-7+1*ZNgP@G|q!&36ZzR~&8RUurdIAI`hv&@|Q%WF`wg#|J%X
zj@n2HOBXf@n?}6TD6fY8;=be|*DbObgepNkL1uhHsLo==-$^DYJ;ADiij!wEBX4PN
zSlbCX9ClwvdXZi?SJ|id?)gdT*)2&^cyRL_SzXQo3!#CMU-2}6r_q9~d^p(iDQWD}
zZ~ikhATmgK@Ub^^E1-S1w+^d9avz%$((~}Ux02_;nXRDZa<Z-Ydcb;UgPU6Sc3HCP
zAsSHl%C=jZU@o*NM4{}}#vZ1QAq6OCf+nb4jp5{iV$hzILH96^=1z@nHTaydaxp~v
z=!G~sK>3=6@-8SsZX3ni@(12ip%hIrs-M0Y7+6eY6Fu{S2s6?aBjN_2Td50F^qY^B
z!ozPhs~nw9c|JkfY5a#TN-aLglokJjsu=rbgNPEtDi2ZBz(GemRhs_G$BsoMM%gw#
zu=g%7N#67`+t_GWm%~HZYRbKH@p0Zcw(sp9BZlqWgLBr!4<q=fQ>`ER#96@Z?E<pJ
zRWw#<zIjwbUKH779-jTW+IhmF(nw*>r#@X)v~lJ+#KSU|v08+CO$Yg6PTZL~CbZLE
zd(1an72{_?D#9gUMJpTz)Nl7RTbo2w&CK+}6oQ=8$qJSYid1)D!}(G^$QQe{wp~ez
z$Y@^;S#G7>BqFRFYeeQ*mPKCQvWL|-l0MbKbIOfX0<Wv0>(ws1%8M&kx5WeGl2jyX
zyW_*yN#tP)ha1s+hPS`-1WFiFPh?Kj(+xY-CdS#`!rzL>7k$klDv>sc*FCjewpYAZ
z_#oxivl%;cYZWAPKML9|N!(z@3|c-J^;e|$*ZLvE`~+?lz=f^z1CD1dn(>`h)Whc#
zNghsr<~mhSF^{{7P#D4a;Zku8sPx*k#|?XS-IK^9%5gOML`H~Lr2+qfrFpscZ4FhF
z#BErXlugrm2Mcq5u?(n9=+N3;BcWf2#}gL__?vTFid!YA!C?5(qsJO=B-iiH;sx(B
zRrrmB^yIbW?z}3>4}c~~wt*gm5>}9x%BE6HL4fl2Co$-*R;{jRA3clvE8KS#Jn$Z)
z*BtamEpswzI|pE`gAzUkQh{cS1>YnIUeB|5LDY63*h}8&-O-5|#!SN5&yuEXht1)=
z$@%vkuRJruz8rGfJ+pe7??&sA7H_tTAz@>}nwXK&C!YH4kJ_pgEiqv*V<ASmSEsl&
zuNgH}?>)(t7gPyiFBQ7_^QlZr=fGdKg}{rK+jG1)jxfK$H8c0p0N!O>Q%QFafr}HD
zI!2IA+_K<x`MGIK3m<i6RsvqhVQAg_G}E!NVdUcR_FKnM8rezi13{a(F66{LLURTs
zjMTQ>N<#G%S@#L#M#Wa?Q(NTo9eX~n)*U<cx+7&AzApZ|&B*19<4vwhkAidtaHG^C
z6e!s)k*bFc1m7(<zR`*7Jn8IXTcB<$W#%$yIt=G%kM%PXiS{UzbI;I{P)QK^kfpwN
zn+A9y)8-B}#4&(cMEG_hRSa}?d)!@SQa0l}{kqyWYJDwRvm(5ie|x~&-`_oVR_4fb
zvX+|mErYU~9JQb<S5FpbTT<8DqdOH@mbrs^n$;a%_bTW^+!fQNF}dQ$P5!!BJ8F|1
z8{2`PAAM(uo|leC2o5mX?zo7P_L6Eo8&Yc;xny0){d=J;J|e|S)eBjZ)$;=DUc2Tl
z7^G?8dC<?5pu~)MIM-8C>HIOOX0-LSJ-_MYXC6Zfj+j~2lV(R5+XuQuaazkt2?5xM
zh|+$o9)Dfa;jLPgp4@Lzzac6-C7BmJiiyL020bxL%JEzV-q4!9fBwCpjsDA<FHp`V
WPzsgSiz;yVHxDlUe5e7XP5c7`dG?3^
literal 0
HcmV?d00001
diff --git a/doc/src/Eqs/fix_gcmc3.tex b/doc/src/Eqs/fix_gcmc3.tex
new file mode 100644
index 0000000000..6aac2b9e1f
--- /dev/null
+++ b/doc/src/Eqs/fix_gcmc3.tex
@@ -0,0 +1,9 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\Lambda &=& \sqrt{ \frac{h^2}{2 \pi m k T}}
+\end{eqnarray*}
+
+\end{document}
\ No newline at end of file
diff --git a/doc/src/dihedral_charmm.txt b/doc/src/dihedral_charmm.txt
index 802cf2c4ad..f4637e776d 100644
--- a/doc/src/dihedral_charmm.txt
+++ b/doc/src/dihedral_charmm.txt
@@ -40,8 +40,11 @@ field.
NOTE: The newer {charmmfsh} style was released in March 2017. We
recommend it be used instead of the older {charmm} style when running
-a simulation with the CHARMM force field. See the discussion below
-and more details on the "pair_style charmm"_pair_charmm.html doc page.
+a simulation with the CHARMM force field and Coulomb cutoffs, via the
+"pair_style lj/charmmfsw/coul/charmmfsh"_pair_charmm.html command.
+Otherwise the older {charmm} style is fine to use. See the discussion
+below and more details on the "pair_style charmm"_pair_charmm.html doc
+page.
The following coefficients must be defined for each dihedral type via the
"dihedral_coeff"_dihedral_coeff.html command as in the example above, or in
@@ -82,13 +85,18 @@ special_bonds 1-4 scaling factor to 0.0 (which is the
default). Otherwise 1-4 non-bonded interactions in dihedrals will be
computed twice.
-For simulations using the CHARMM force field, the difference between
-the {charmm} and {charmmfsh} styles is in the computation of the 1-4
-non-bond interactions, if the distance between the two atoms is within
-the switching distance of the pairwise potential defined by the
-corresponding CHARMM pair style, i.e. between the inner and outer
-cutoffs specified for the pair style. See the discussion on the
-"CHARMM pair_style"_pair_charmm.html doc page for details.
+For simulations using the CHARMM force field with a Coulomb cutoff,
+the difference between the {charmm} and {charmmfsh} styles is in the
+computation of the 1-4 non-bond interactions, though only if the
+distance between the two atoms is within the switching region of the
+pairwise potential defined by the corresponding CHARMM pair style,
+i.e. between the inner and outer cutoffs specified for the pair style.
+The {charmmfsh} style should only be used when using the "pair_style
+lj/charmmfsw/coul/charmmfsh"_pair_charmm.html to make the Coulombic
+pairwise calculations consistent. Use the {charmm} style with
+long-range Coulombics or the older "pair_style
+lj/charmm/coul/charmm"_pair_charmm.html command. See the discussion
+on the "CHARMM pair_style"_pair_charmm.html doc page for details.
Note that for AMBER force fields, which use pair styles with "lj/cut",
the special_bonds 1-4 scaling factor should be set to the AMBER
diff --git a/doc/src/fix_gcmc.txt b/doc/src/fix_gcmc.txt
index 1a419628e9..723e0ec6d9 100644
--- a/doc/src/fix_gcmc.txt
+++ b/doc/src/fix_gcmc.txt
@@ -122,11 +122,11 @@ If used with "fix nvt"_fix_nh.html, the temperature of the imaginary
reservoir, T, should be set to be equivalent to the target temperature
used in fix nvt. Otherwise, the imaginary reservoir
will not be in thermal equilibrium with the simulation cell. Also,
-it is important that the temperature used by fix nvt be dynamic,
+it is important that the temperature used by fix nvt be dynamic/dof,
which can be achieved as follows:
compute mdtemp mdatoms temp
-compute_modify mdtemp dynamic yes
+compute_modify mdtemp dynamic/dof yes
fix mdnvt mdatoms nvt temp 300.0 300.0 10.0
fix_modify mdnvt temp mdtemp :pre
@@ -204,10 +204,43 @@ atoms/molecules are assigned to two groups: the default group "all"
and the group specified in the fix gcmc command (which can also be
"all").
-The gas reservoir pressure can be specified using the {pressure}
+The chemical potential is a user-specified input parameter defined
+as:
+
+:c,image(Eqs/fix_gcmc1.jpg)
+
+The second term mu_ex is the excess chemical potential due to
+energetic interactions and is formally zero for the fictitious gas
+reservoir but is non-zero for interacting systems. So, while the chemical
+potential of the reservoir and the simulation cell are equal, mu_ex is not,
+and as a result, the densities of the two are generally quite different.
+The first term mu_id is the ideal gas contribution to the chemical potential.
+mu_id can be related to the density or pressure of the fictitious gas reservoir by:
+
+:c,image(Eqs/fix_gcmc2.jpg)
+
+where k is Boltzman's constant,
+T is the user-specified temperature, rho is the number density,
+P is the pressure, and phi is the fugacity coefficient.
+The constant Lambda is required for dimensional consistency.
+For all unit styles except {lj} it is defined as the thermal
+de Broglie wavelength
+
+:c,image(Eqs/fix_gcmc3.jpg)
+
+where h is Planck's constant, and m is the mass of the exchanged atom or molecule.
+For unit style {lj}, Lambda is simply set to the unity. Note that prior to March 2017
+Lambda for unit style {lj} was calculated using the above formula with h set to
+the rather specific value of 0.18292026. Chemical potential
+under the old definition can be converted to an equivalent value under the new
+definition by subtracting 3kTln(Lambda_old).
+
+As an alternative to specifying mu directly, the ideal gas reservoir
+can be defined by its pressure P using the {pressure}
keyword, in which case the user-specified chemical potential is
-ignored. For non-ideal gas reservoirs, the user may also specify the
-fugacity coefficient using the {fugacity_coeff} keyword.
+ignored. The user may also specify the
+fugacity coefficient phi using the {fugacity_coeff} keyword, which
+defaults to unity.
The {full_energy} option means that fix GCMC will compute the total
potential energy of the entire simulated system. The total system
@@ -224,7 +257,8 @@ potential energy calculations, including the following:
many-body pair styles
hybrid pair styles
eam pair styles
- triclinic systems
+ tail corrections
+
need to include potential energy contributions from other fixes :ul
In these cases, LAMMPS will automatically apply the {full_energy}
@@ -276,6 +310,13 @@ therefore, you will want to use the
current number of atoms is used as a normalizing factor each time
temperature is computed. Here is the necessary command:
+NOTE: If the density of the cell is initially very small or zero,
+and increases to a much larger density after a period of equilibration,
+then certain quantities that are only calculated once at the start
+(kspace parameters, tail corrections) may no longer be accurate.
+The solution is to start a new simulation after the equilibrium
+density has been reached.
+
With some pair_styles, such as "Buckingham"_pair_buck.html,
"Born-Mayer-Huggins"_pair_born.html and "ReaxFF"_pair_reax_c.html,
two atoms placed close to each other may have an arbitrary large,
@@ -366,7 +407,7 @@ referenced by the user for each subsequent fix gcmc command.
[Default:]
The option defaults are mol = no, maxangle = 10, overlap_cutoff = 0.0,
-and full_energy = no,
+fugacity_coeff = 1, and full_energy = no,
except for the situations where full_energy is required, as
listed above.
diff --git a/doc/src/pair_charmm.txt b/doc/src/pair_charmm.txt
index a3c93c7589..65b9f381c5 100644
--- a/doc/src/pair_charmm.txt
+++ b/doc/src/pair_charmm.txt
@@ -84,9 +84,9 @@ CHARMM force field.
The styles with {charmm} (not {charmmfsw} or {charmmfsh}) in their
name are the older, original LAMMPS implementations. They compute the
LJ and Coulombic interactions with an energy switching function (esw,
-a cubic polynomial, shown in the formula below), which ramps the
-energy smoothly to zero between the inner and outer cutoff. This can
-cause irregularities in pair-wise forces (due to the discontinuous 2nd
+shown in the formula below as S(r)), which ramps the energy smoothly
+to zero between the inner and outer cutoff. This can cause
+irregularities in pair-wise forces (due to the discontinuous 2nd
derivative of energy at the boundaries of the switching region), which
in some cases can result in detectable artifacts in an MD simulation.
diff --git a/examples/gcmc/in.gcmc.lj b/examples/gcmc/in.gcmc.lj
index 9b1f465f04..fc9afdb7f8 100644
--- a/examples/gcmc/in.gcmc.lj
+++ b/examples/gcmc/in.gcmc.lj
@@ -1,18 +1,23 @@
# GCMC for LJ simple fluid, no dynamics
+# T = 2.0
+# rho ~ 0.5
+# p ~ 1.5
+# mu_ex ~ 0.0
+# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8
-# variables available on command line
+# variables modifiable using -var command line switch
-variable mu index -21.0
-variable disp index 1.0
+variable mu index -1.25
variable temp index 2.0
-variable lbox index 10.0
+variable disp index 1.0
+variable lbox index 5.0
# global model settings
units lj
atom_style atomic
-pair_style lj/cut 3.0
-pair_modify tail yes
+pair_style lj/cut 3.0
+pair_modify tail no # turn of to avoid triggering full_energy
# box
@@ -28,15 +33,27 @@ mass * 1.0
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
+# averaging
+
+variable rho equal density
+variable p equal press
+variable nugget equal 1.0e-8
+variable lambda equal 1.0
+variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
+fix ave all ave/time 10 100 1000 v_rho v_p v_muex ave one file rho_vs_p.dat
+variable rhoav equal f_ave[1]
+variable pav equal f_ave[2]
+variable muexav equal f_ave[3]
+
# output
-variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
-variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
-variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
compute_modify thermo_temp dynamic yes
-thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
-thermo 100
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
+thermo 1000
# run
-run 1000
+run 10000
diff --git a/examples/gcmc/log.24Mar17.gcmc.lj.g++.1 b/examples/gcmc/log.24Mar17.gcmc.lj.g++.1
index 5f967192e4..36a9fe885d 100644
--- a/examples/gcmc/log.24Mar17.gcmc.lj.g++.1
+++ b/examples/gcmc/log.24Mar17.gcmc.lj.g++.1
@@ -1,28 +1,35 @@
LAMMPS (17 Mar 2017)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
+ using 1 OpenMP thread(s) per MPI task
# GCMC for LJ simple fluid, no dynamics
+# T = 2.0
+# rho ~ 0.5
+# p ~ 1.5
+# mu_ex ~ 0.0
+# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8
-# variables available on command line
+# variables modifiable using -var command line switch
-variable mu index -21.0
-variable disp index 1.0
+variable mu index -1.25
variable temp index 2.0
-variable lbox index 10.0
+variable disp index 1.0
+variable lbox index 5.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
-pair_modify tail yes
+pair_modify tail no # turn of to avoid triggering full_energy
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
-region box block 0 10.0 0 ${lbox} 0 ${lbox}
-region box block 0 10.0 0 10.0 0 ${lbox}
-region box block 0 10.0 0 10.0 0 10.0
+region box block 0 5.0 0 ${lbox} 0 ${lbox}
+region box block 0 5.0 0 5.0 0 ${lbox}
+region box block 0 5.0 0 5.0 0 5.0
create_box 1 box
-Created orthogonal box = (0 0 0) to (10 10 10)
+Created orthogonal box = (0 0 0) to (5 5 5)
1 by 1 by 1 MPI processor grid
# lj parameters
@@ -34,70 +41,89 @@ mass * 1.0
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
-fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
-fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 1.0
+
+# averaging
+
+variable rho equal density
+variable p equal press
+variable nugget equal 1.0e-8
+variable lambda equal 1.0
+variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
+variable muex equal -1.25-${temp}*ln(density*${lambda}+${nugget})
+variable muex equal -1.25-2.0*ln(density*${lambda}+${nugget})
+variable muex equal -1.25-2.0*ln(density*1+${nugget})
+variable muex equal -1.25-2.0*ln(density*1+1e-08)
+fix ave all ave/time 10 100 1000 v_rho v_p v_muex ave one file rho_vs_p.dat
+variable rhoav equal f_ave[1]
+variable pav equal f_ave[2]
+variable muexav equal f_ave[3]
# output
-variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
-variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
-variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+1e-08)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+1e-08)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+1e-08)
compute_modify thermo_temp dynamic yes
-thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
-thermo 100
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
+thermo 1000
# run
-run 1000
+run 10000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.3
ghost atom cutoff = 3.3
- binsize = 1.65, bins = 7 7 7
+ binsize = 1.65, bins = 4 4 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
-Per MPI rank memory allocation (min/avg/max) = 0.4369 | 0.4369 | 0.4369 Mbytes
-Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
- 0 0 0 0 -0 0 0 0 0 0
- 100 1.9042848 0.39026453 -1.7692765 2.8466449 0.292 292 0.3619855 0.30247792 0.40278761
- 200 1.8651924 0.47815517 -1.8494955 2.7886155 0.305 305 0.34021109 0.30357196 0.37759189
- 300 2.0626994 0.52068504 -1.8197295 3.0834166 0.291 291 0.32055605 0.3003043 0.36103862
- 400 2.0394818 0.53751435 -1.7636699 3.0482184 0.278 278 0.31698808 0.29995864 0.35441275
- 500 1.9628066 0.54594742 -1.7145336 2.9339513 0.287 287 0.31211861 0.29724228 0.35161407
- 600 1.9845913 0.40846162 -1.8199325 2.9669308 0.299 299 0.30976643 0.29612711 0.34933559
- 700 1.8582606 0.53445462 -1.7869306 2.777974 0.296 296 0.30642103 0.29446478 0.34633665
- 800 2.0340641 0.66057698 -1.7075279 3.0403148 0.283 283 0.30730979 0.29746793 0.34768045
- 900 2.0830765 0.63731971 -1.894775 3.114911 0.322 322 0.30636338 0.29737705 0.34737644
- 1000 1.9688933 0.50024802 -1.7013944 2.9428299 0.281 281 0.3053174 0.29772245 0.34788254
-Loop time of 3.98286 on 1 procs for 1000 steps with 281 atoms
-
-Performance: 108464.750 tau/day, 251.076 timesteps/s
-99.9% CPU use with 1 MPI tasks x no OpenMP threads
+Per MPI rank memory allocation (min/avg/max) = 0.433 | 0.433 | 0.433 Mbytes
+Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
+ 0 0 0 0 -0 0 0 0 0 0 0 0 0
+ 1000 2.4038954 2.1735585 -2.7041368 3.5476844 0.496 62 0.064790036 0.06313096 0.1081294 0.54304 1.4513524 -0.025479219
+ 2000 2.0461168 1.1913842 -2.9880181 3.0212194 0.512 64 0.067416408 0.066335853 0.11306166 0.52736 1.3274665 0.034690004
+ 3000 1.7930436 1.3788681 -3.2212667 2.6505861 0.552 69 0.067733191 0.066877836 0.1133516 0.5344 1.3834744 0.0070582537
+ 4000 1.981449 1.2541054 -2.8222868 2.9217977 0.472 59 0.068546991 0.067856412 0.11442807 0.52504 1.3815629 0.043309657
+ 5000 2.0946818 1.0701629 -3.5213291 3.0977688 0.568 71 0.06813743 0.067567891 0.11342906 0.53824 1.4049567 -0.0054539777
+ 6000 1.9793484 0.68224187 -3.410211 2.9247088 0.536 67 0.067797628 0.067420108 0.11295333 0.5384 1.401683 -0.0066894359
+ 7000 2.1885798 1.6745012 -3.185499 3.2345922 0.544 68 0.068630201 0.068261832 0.11403705 0.5244 1.449239 0.045987399
+ 8000 2.2175324 1.5897263 -3.078898 3.2759002 0.528 66 0.068180395 0.067899629 0.11332691 0.53928 1.5488388 -0.01075766
+ 9000 1.8610779 1.0396231 -2.923262 2.7465908 0.496 62 0.068346453 0.068028117 0.1134132 0.52912 1.4352871 0.027082544
+ 10000 2.1079271 1.1746643 -2.9112062 3.1091925 0.48 60 0.068352878 0.068054948 0.11335434 0.5316 1.4462327 0.018503094
+Loop time of 13.05 on 1 procs for 10000 steps with 60 atoms
+
+Performance: 331035.016 tau/day, 766.285 timesteps/s
+100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
-Pair | 0.10563 | 0.10563 | 0.10563 | 0.0 | 2.65
-Neigh | 0.33428 | 0.33428 | 0.33428 | 0.0 | 8.39
-Comm | 0.027969 | 0.027969 | 0.027969 | 0.0 | 0.70
-Output | 0.00017285 | 0.00017285 | 0.00017285 | 0.0 | 0.00
-Modify | 3.5096 | 3.5096 | 3.5096 | 0.0 | 88.12
-Other | | 0.005197 | | | 0.13
-
-Nlocal: 281 ave 281 max 281 min
+Pair | 0.37239 | 0.37239 | 0.37239 | 0.0 | 2.85
+Neigh | 0.94764 | 0.94764 | 0.94764 | 0.0 | 7.26
+Comm | 0.092473 | 0.092473 | 0.092473 | 0.0 | 0.71
+Output | 0.00023365 | 0.00023365 | 0.00023365 | 0.0 | 0.00
+Modify | 11.627 | 11.627 | 11.627 | 0.0 | 89.09
+Other | | 0.01054 | | | 0.08
+
+Nlocal: 60 ave 60 max 60 min
Histogram: 1 0 0 0 0 0 0 0 0 0
-Nghost: 977 ave 977 max 977 min
+Nghost: 663 ave 663 max 663 min
Histogram: 1 0 0 0 0 0 0 0 0 0
-Neighs: 5902 ave 5902 max 5902 min
+Neighs: 2133 ave 2133 max 2133 min
Histogram: 1 0 0 0 0 0 0 0 0 0
-Total # of neighbors = 5902
-Ave neighs/atom = 21.0036
-Neighbor list builds = 1000
+Total # of neighbors = 2133
+Ave neighs/atom = 35.55
+Neighbor list builds = 10000
Dangerous builds = 0
-Total wall time: 0:00:03
+Total wall time: 0:00:13
diff --git a/examples/gcmc/log.24Mar17.gcmc.lj.g++.4 b/examples/gcmc/log.24Mar17.gcmc.lj.g++.4
index db4cc7fd9a..8694d8b95e 100644
--- a/examples/gcmc/log.24Mar17.gcmc.lj.g++.4
+++ b/examples/gcmc/log.24Mar17.gcmc.lj.g++.4
@@ -1,28 +1,35 @@
LAMMPS (17 Mar 2017)
+OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:90)
+ using 1 OpenMP thread(s) per MPI task
# GCMC for LJ simple fluid, no dynamics
+# T = 2.0
+# rho ~ 0.5
+# p ~ 1.5
+# mu_ex ~ 0.0
+# comparable to Frenkel and Smit GCMC Case Study, Figure 5.8
-# variables available on command line
+# variables modifiable using -var command line switch
-variable mu index -21.0
-variable disp index 1.0
+variable mu index -1.25
variable temp index 2.0
-variable lbox index 10.0
+variable disp index 1.0
+variable lbox index 5.0
# global model settings
units lj
atom_style atomic
pair_style lj/cut 3.0
-pair_modify tail yes
+pair_modify tail no # turn of to avoid triggering full_energy
# box
region box block 0 ${lbox} 0 ${lbox} 0 ${lbox}
-region box block 0 10.0 0 ${lbox} 0 ${lbox}
-region box block 0 10.0 0 10.0 0 ${lbox}
-region box block 0 10.0 0 10.0 0 10.0
+region box block 0 5.0 0 ${lbox} 0 ${lbox}
+region box block 0 5.0 0 5.0 0 ${lbox}
+region box block 0 5.0 0 5.0 0 5.0
create_box 1 box
-Created orthogonal box = (0 0 0) to (10 10 10)
+Created orthogonal box = (0 0 0) to (5 5 5)
1 by 2 by 2 MPI processor grid
# lj parameters
@@ -34,70 +41,89 @@ mass * 1.0
fix mygcmc all gcmc 1 100 100 1 29494 ${temp} ${mu} ${disp}
fix mygcmc all gcmc 1 100 100 1 29494 2.0 ${mu} ${disp}
-fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 ${disp}
-fix mygcmc all gcmc 1 100 100 1 29494 2.0 -21.0 1.0
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 ${disp}
+fix mygcmc all gcmc 1 100 100 1 29494 2.0 -1.25 1.0
+
+# averaging
+
+variable rho equal density
+variable p equal press
+variable nugget equal 1.0e-8
+variable lambda equal 1.0
+variable muex equal ${mu}-${temp}*ln(density*${lambda}+${nugget})
+variable muex equal -1.25-${temp}*ln(density*${lambda}+${nugget})
+variable muex equal -1.25-2.0*ln(density*${lambda}+${nugget})
+variable muex equal -1.25-2.0*ln(density*1+${nugget})
+variable muex equal -1.25-2.0*ln(density*1+1e-08)
+fix ave all ave/time 10 100 1000 v_rho v_p v_muex ave one file rho_vs_p.dat
+variable rhoav equal f_ave[1]
+variable pav equal f_ave[2]
+variable muexav equal f_ave[3]
# output
-variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+0.1)
-variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+0.1)
-variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+0.1)
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+${nugget})
+variable tacc equal f_mygcmc[2]/(f_mygcmc[1]+1e-08)
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+${nugget})
+variable iacc equal f_mygcmc[4]/(f_mygcmc[3]+1e-08)
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+${nugget})
+variable dacc equal f_mygcmc[6]/(f_mygcmc[5]+1e-08)
compute_modify thermo_temp dynamic yes
-thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc
-thermo 100
+thermo_style custom step temp press pe ke density atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
+thermo 1000
# run
-run 1000
+run 10000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.3
ghost atom cutoff = 3.3
- binsize = 1.65, bins = 7 7 7
+ binsize = 1.65, bins = 4 4 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
-Per MPI rank memory allocation (min/avg/max) = 0.434 | 0.434 | 0.434 Mbytes
-Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc
- 0 0 0 0 -0 0 0 0 0 0
- 100 2.0328045 0.58661762 -1.6812724 3.0385824 0.287 287 0.35917318 0.30067507 0.38663622
- 200 1.9594279 0.50682399 -1.7308396 2.9287927 0.284 284 0.33788365 0.30337335 0.37300293
- 300 2.0602937 0.7028247 -1.9278541 3.0806296 0.315 315 0.31882007 0.29697498 0.36167185
- 400 1.995183 0.4328246 -1.8715454 2.983026 0.307 307 0.31527654 0.29681901 0.35673374
- 500 2.1390101 0.48232215 -1.554138 3.1960306 0.257 257 0.31372975 0.30003067 0.35558858
- 600 2.0584244 0.4929049 -1.6995569 3.0767263 0.283 283 0.31114213 0.29801665 0.35160109
- 700 1.9155066 0.49654243 -1.5770611 2.8624174 0.265 265 0.31056419 0.29944173 0.35157337
- 800 2.0883562 0.52731947 -1.8261112 3.1220925 0.3 300 0.30730979 0.29704354 0.34898892
- 900 2.0470677 0.5605993 -2.0130053 3.0610656 0.322 322 0.30484441 0.29586719 0.34678883
- 1000 2.004135 0.50642204 -1.6956257 2.9955798 0.283 283 0.30396929 0.29634309 0.34770304
-Loop time of 3.688 on 4 procs for 1000 steps with 283 atoms
-
-Performance: 117136.751 tau/day, 271.150 timesteps/s
-99.2% CPU use with 4 MPI tasks x no OpenMP threads
+Per MPI rank memory allocation (min/avg/max) = 0.4477 | 0.4477 | 0.4477 Mbytes
+Step Temp Press PotEng KinEng Density Atoms v_iacc v_dacc v_tacc v_rhoav v_pav v_muexav
+ 0 0 0 0 -0 0 0 0 0 0 0 0 0
+ 1000 1.956397 1.7699101 -2.7889468 2.8864874 0.488 61 0.068894746 0.067229075 0.1141726 0.53288 1.3832798 0.013392866
+ 2000 2.040943 0.56060899 -2.8001647 3.0077055 0.456 57 0.069858594 0.068831934 0.11629114 0.5232 1.3587174 0.049995794
+ 3000 2.0004866 1.5736515 -3.3098044 2.9572411 0.552 69 0.069594029 0.068727791 0.11592543 0.53096 1.4129434 0.020022578
+ 4000 2.1127942 2.642809 -2.8865084 3.1211733 0.528 66 0.070268697 0.069533235 0.11693806 0.52424 1.3444615 0.046884078
+ 5000 2.3663648 1.354269 -3.1917346 3.4957662 0.528 66 0.070519633 0.069960064 0.11710321 0.52688 1.3595814 0.036270867
+ 6000 1.9224136 0.82756699 -3.1965 2.839257 0.52 65 0.06985018 0.069474645 0.11628632 0.536 1.47062 0.00141549
+ 7000 2.0266192 1.5593811 -2.9972341 2.9931606 0.52 65 0.070244693 0.069880791 0.11666541 0.52528 1.3246332 0.040754793
+ 8000 1.7790467 1.8680568 -2.8028819 2.6275151 0.52 65 0.070454494 0.070172368 0.11736806 0.524 1.4213649 0.047985191
+ 9000 1.7968847 1.3195587 -3.261001 2.6550983 0.536 67 0.069952011 0.069618327 0.11650087 0.53904 1.4624201 -0.01069837
+ 10000 2.1566109 1.1015729 -3.4999837 3.1880335 0.552 69 0.069603309 0.069284134 0.11625548 0.53128 1.3587249 0.02075238
+Loop time of 13.0611 on 4 procs for 10000 steps with 69 atoms
+
+Performance: 330753.007 tau/day, 765.632 timesteps/s
+99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
-Pair | 0.024644 | 0.026027 | 0.027483 | 0.6 | 0.71
-Neigh | 0.085449 | 0.088998 | 0.092893 | 0.9 | 2.41
-Comm | 0.045756 | 0.051296 | 0.056578 | 1.7 | 1.39
-Output | 0.00028491 | 0.00030857 | 0.00035262 | 0.0 | 0.01
-Modify | 3.5189 | 3.5191 | 3.5194 | 0.0 | 95.42
-Other | | 0.002221 | | | 0.06
-
-Nlocal: 70.75 ave 77 max 68 min
-Histogram: 1 2 0 0 0 0 0 0 0 1
-Nghost: 514.25 ave 520 max 507 min
-Histogram: 1 0 0 0 1 0 0 1 0 1
-Neighs: 1483.5 ave 1715 max 1359 min
-Histogram: 2 0 0 1 0 0 0 0 0 1
-
-Total # of neighbors = 5934
-Ave neighs/atom = 20.9682
-Neighbor list builds = 1000
+Pair | 0.08888 | 0.09443 | 0.099889 | 1.4 | 0.72
+Neigh | 0.27721 | 0.28532 | 0.29177 | 1.1 | 2.18
+Comm | 0.27648 | 0.28875 | 0.30268 | 1.9 | 2.21
+Output | 0.00029635 | 0.00043058 | 0.00048113 | 0.0 | 0.00
+Modify | 12.384 | 12.384 | 12.384 | 0.0 | 94.82
+Other | | 0.008055 | | | 0.06
+
+Nlocal: 17.25 ave 23 max 10 min
+Histogram: 1 0 0 0 0 0 2 0 0 1
+Nghost: 506.5 ave 519 max 490 min
+Histogram: 1 0 1 0 0 0 0 0 0 2
+Neighs: 705.75 ave 998 max 369 min
+Histogram: 1 0 0 0 0 1 1 0 0 1
+
+Total # of neighbors = 2823
+Ave neighs/atom = 40.913
+Neighbor list builds = 10000
Dangerous builds = 0
-Total wall time: 0:00:03
+Total wall time: 0:00:13
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 780df25bdf..eae97deab5 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -432,7 +432,8 @@ void FixGCMC::init()
(force->pair == NULL) ||
(force->pair->single_enable == 0) ||
(force->pair_match("hybrid",0)) ||
- (force->pair_match("eam",0))
+ (force->pair_match("eam",0)) ||
+ (force->pair->tail_flag)
) {
full_flag = true;
if (comm->me == 0)
@@ -618,13 +619,18 @@ void FixGCMC::init()
}
// compute beta, lambda, sigma, and the zz factor
-
+ // For LJ units, lambda=1
beta = 1.0/(force->boltz*reservoir_temperature);
- double lambda = sqrt(force->hplanck*force->hplanck/
- (2.0*MY_PI*gas_mass*force->mvv2e*
+ if (strcmp(update->unit_style,"lj") == 0)
+ zz = exp(beta*chemical_potential);
+ else {
+ double lambda = sqrt(force->hplanck*force->hplanck/
+ (2.0*MY_PI*gas_mass*force->mvv2e*
force->boltz*reservoir_temperature));
+ zz = exp(beta*chemical_potential)/(pow(lambda,3.0));
+ }
+
sigma = sqrt(force->boltz*reservoir_temperature*tfac_insert/gas_mass/force->mvv2e);
- zz = exp(beta*chemical_potential)/(pow(lambda,3.0));
if (pressure_flag) zz = pressure*fugacity_coeff*beta/force->nktv2p;
imagezero = ((imageint) IMGMAX << IMG2BITS) |
diff --git a/src/library.cpp b/src/library.cpp
index ca3276ee8e..f9c5ad4892 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -790,7 +790,7 @@ void lammps_gather_atoms(void *ptr, char *name,
for (i = 0; i < nlocal; i++) {
offset = count*(tag[i]-1);
for (j = 0; j < count; j++)
- copy[offset++] = array[i][0];
+ copy[offset++] = array[i][j];
}
MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
diff --git a/src/update.cpp b/src/update.cpp
index 35492006cc..5599dc6c88 100644
--- a/src/update.cpp
+++ b/src/update.cpp
@@ -128,7 +128,7 @@ void Update::set_units(const char *style)
if (strcmp(style,"lj") == 0) {
force->boltz = 1.0;
- force->hplanck = 0.18292026; // using LJ parameters for argon
+ force->hplanck = 1.0;
force->mvv2e = 1.0;
force->ftm2v = 1.0;
force->mv2d = 1.0;
--
GitLab