diff -Naur lammps-3Nov06/doc/Eqs/Script.create lammps-4Nov06/doc/Eqs/Script.create --- lammps-3Nov06/doc/Eqs/Script.create 2006-09-21 10:22:34.000000000 -0600 +++ lammps-4Nov06/doc/Eqs/Script.create 2006-11-03 10:09:44.000000000 -0700 @@ -40,5 +40,7 @@ latex pair_lj_smooth.tex latex pair_morse.tex latex pair_soft.tex +latex pair_sw.tex +latex pair_tersoff.tex latex pair_yukawa.tex latex stress_tensor.tex diff -Naur lammps-3Nov06/doc/Eqs/pair_debye.tex lammps-4Nov06/doc/Eqs/pair_debye.tex --- lammps-3Nov06/doc/Eqs/pair_debye.tex 2006-09-21 10:22:34.000000000 -0600 +++ lammps-4Nov06/doc/Eqs/pair_debye.tex 2006-11-03 10:09:44.000000000 -0700 @@ -5,4 +5,4 @@ E = \frac{C q_i q_j}{\epsilon r} \exp(- \kappa r) \qquad r < r_c $$ -\end{document} \ No newline at end of file +\end{document} diff -Naur lammps-3Nov06/doc/Eqs/pair_sw.jpg lammps-4Nov06/doc/Eqs/pair_sw.jpg --- lammps-3Nov06/doc/Eqs/pair_sw.jpg 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/Eqs/pair_sw.jpg 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,47 @@ +JFIFHHC  !"$"$ J!"1AQa#2BqU37RWt$SVbCDrscu?˥)JR)JR)JR)JR)JR)JR)JR)JRcTHZmimQcȑ!% sc`h`' ?pvmO\7R)(+VAf"Ҕ.PUΰBԺ2w; L[ZGºJ~P0r>=5ΐ Q +kAa*S_Xp=ǨlZ$5i~T̫m"k u/pGX aKO3lCv!d֑o$$ X4jߦ@^Z~A{~uGitMC|['%*kUĜKP'pQ<$z֎%0H%! »;F \oɴ2QWHQ)G梕ڳ3>7$  1BW8gOfFʞu%CmHXʼCiW: Rˍ\/wm1nOj +(BX|֏:Bߪ,7=F)QN i~Jm`=rh`զkZrS2;'%ֺRҕa4k.?PͲZZV8wY !؇9[ZGU2FN`4[~k)z;jI\BA]׋5 lݮOMTGp:.:,!@D^iw d=uuj)Z +ӃxSH>HKrH +zT2HF5u kS\A!I?54ԝ+mט.?׌ɒ5W*l`$nP$W )DŤ,4T 979*~MBBJ?5&՜q)J]kL~KUzP:*[V%+%8U*"]A}) 9]JJgPYJR)JRi:҆В)G y}xX#Tkۦo;S;![9 Tdgwϳ[Yy2@ZQAλ3R}#Ԏϻv|?=3y\T񞴟}%۵;>+}qZZ}֘eo<i-kP J@$ +{SRJn\;b !iNoPGb=$Wm\볭w'f9>XsX4{7KL6} +aFHIVUo7d^yw"MQ\Ͷ;7&be\poA^lrR60Xt+Nۏ%.jeIy]:JJʓe*ƕеxY g1_uTs.Cu~P0ԅ6Q +d@uXb2KKm>#(H y$(nᑕBѝd޾1w9Ϯ7=yz|{o-ٗ޾>/rL+n0y4i*Z֠I>@W;#oUY ܸvŀCҜ[ ` VzI  'PG&)(/NPeM JV+rLl--J#/ *K}߂u13;>+J17ڣw#j[J{HK[J*9;@ƩgO)U[H#n232|$\#'%lbye) V7HCuB~)(i>oP'3'jĴ\yϭ9\-)Ï8~La@H5kPWydrPb.'K@N< ϟx]E{%bvCQⰧ]9NvHPGs+[zsm1iL%䂜گoW+Q&ʸǕeTG],'{ cyg_;H!3qc?=f)JR)JR+ڶ՚zegrm*`$%AP;mgZVF#-6tFaVF^;)A2N]FӺmS4NG۾1._|9''ky>N@Q݌T> +H,S8Ԯʹδ^.uc Pa5vb@$ye)G,iJ}2/(NvmN""t&*N +#ȜU+(uq:~i*^>N+98ƕn+!%qԵoJ]O,Asm7]Y[ݣѤ(/QJ3m89\ Í3juNJGZJ?d{VڳRƒO\"!.Q]}pP eC>,*sBњ_4kCZven2rԥ(#rFIZhw\ fwצ% $to;I +; +ہIQx2<Ŋgٶ֖l}A2J&BhV$O5Vik %J7++|kC9PîejjF7/w;qq!攡X?1D楱@ԩW򺋔~m2lHJSW+Q5r2\̥(-)OU?3N٦KrLAkux!շZڲP?d{ +j7vVPi u %HT|Vs; Ź%[g5Ջ-{w!i)P97Qb9+f TAi+K@,P9jшjK +xRG*>ި[uYK1Gxw-!. V{M>VɮgB# I#2JMghJbTmzHN #`+Gs&F?-nAfa%ŧhm! m yя#5"]>޶%*HƬO)JR)JR)J\PޥZI R0pVu pʶ'>I!IM#]^y 5ӏA7]F|*;d#*k7=?{ L/= wm*i9z(MXN,>ۑd/nmIZ>V,[R)JR)JR)JR)JR+٘fW%,[aX0q]ew7KPb=:gqQK}xRqdJQSϸ8I ܢH@$œH"ʜ#2%ʟ鎡@$ 'Y&R)JWbEv[3䰕%+l+U<+VL5.j\]GS*)t +N22{A>x[?5)j^ĥ)m>rZHH%̴%OSTG* }+ۄ+RE;”Q3]ڹC+j~q[qkScrQ*YTIҹi3Vvdq.H(JR<-Z)JR)JR)JR)JR)Xy&v -3s٘쇛WH5 *T^3)Humh8ird3@ɬƘ:iكrd1~vBo[1P((n;RԷHWlbu3%"\ YZ[ˁڐxTSzinW[T+_ģ.+g.W8V-W IX>/;iK3A-t)`7T}q[ZRSo!ƙgJy1ݷ$`BTbے&2@ '\J6%JÓdww˜+-mb! TbniM{WgYCDe$<:!“54S*nNɔ&"Sg)J#i@|'B}4Yu @5ڔ)JR)JR͒de xDԄ)DROwD&d9.tPZT -0^T6)ASHYVI>gg85rxú*1將0B3= L)J'vsivg$⧝ii S` +V3SK:7FszKikq~%,, +6dF9cԮ735Vdc +!禴YJЮ33,k 3lO Z\eG yM-&%Jx) H8#)yeJ'j) ʈTa*Se)q +#9\xyzV=be\F!ԽQ$ԀyvutMmh +~Jꏓl6R~dhHGźk}l*~\Gn40!oJ<8*9UC鲔 6@H^5s EfCm !8 c:r֝953)O/6nJa9 r3[Rlw&@zC&=CRSg6W\;ݖer,Tөc}A$hWμۃ%"AN4 Pm!!)H$e\dť)XWдݬ8ڛ%T&J @^i{$yN”)JR)JW˻J6G>W|TU\d훤-8|!,Q"btn阥'O3ڗ%@5dL[nwShBP\%N9Zmj+$tzfŧ6Gd=nWNԓ9$ry~ymyE]e!d8bZjKeD+۵K~#y-~ݼ=v/C-E(F :@V7̍ԝM,w?G9͛nrnW+by%ȐiMŵǝ~i=ygAkAv 7yS^tr:P[k W!a'ڝL8Ro-Zd;ɕon: _d 8 x ѧ1!ulƚu+\w[orPiPI;T*tTΚDd) +~!I7[A *[G#}8њK,q;PFO'<SUHwqϞO|`^KOvho|KO7 +qD3'O, TԂ\>Z]zfp=2'hZ˛vÇܦ@;KVH$Gq&5>!/ZgZ3݆F1h4Vƹǭv%-Hm) ۀpyD_&ݮ][nU̵أuWFc4Q%U smV Jܚƕ 2cdU׼Z^vҴ:8 +`' -m7]3/4]]r&@\R 9vsVABulNHTAnEKjCj +CeHZ'Q㭦=?I?Sٮ7~Yz-En ;`nߴs1K7Qk4gV} Ƴ+i >?cYCuCZ 7lj +'$` 0i0O)'*<g:jO-^mR?*Ҕ)JR)JW&Xmq>q|ԧ[ifvwE4.\hM;+'I皛pP.!0l=|'3` kq ([O9G4 [CM6! )H}Ҕ)XM,w?G*H܏P-l2U)*8NqZ/7pyo&:+BC6i PH'7([Ie*W`j}NIsF%Rl%I)@DG +q>`$1c{[orln-9QJRέ92:MPڹ8{NNҝ'y!F@f F2vg$O$II沚D_&ݭ*Q_Zԧ6H}kj +BA#8>JT ^Zv-YJRg8ښRHt"ʴϋ( SkhA-J*sŽ`iPW#˓UV*D,sVֹdYQw:r[رqX4k×hj`ԪD붥) ‚ly!#rR!j]adfdvYp=)#9!;h%)V`ݸghnm>{e!>ńJ(OXΩ{0".Lɬ.<۴s_qN% +8iJR)JR)JR)JR:KJmWk}Um/aB@P!điJdA"m4e[R<$~d}je)JT92..02N£K7dO$d`9JҦUSZ[ R-( |y>fԥ*5KEs!-}؂A - [ک(߬lCLR|>''p`7$ʱlЮ DAHZBH +m)w#nO}IZGܵ }IygJqPujsKEXTxaJuƣK8 +(@kaM'sSҵ:-:[OQ 2+=f鉓b5yq !r䴞%D)^.F5fO6dJe.)OURSM:ZSbe .#ru#rPU +8 M1z4ͷPBW%#r*j}4a0jAņ[㣂8*j|)4])-X +c#ֲYB8(i`IUE{/Ohۓߤžۉn,eyn) y +PN$tJR B(9-8!bmxD̼ZECԓӵIu^G5a +~ii#ý0jh-i'f]opڡ$ Th=sӒctxv3 +ʐ\=>z8( )D ytw+.3̶k܆%묕PՂ57@_ob$ܱ)(mkRQoZH"9b\Yi +!/!+UXI/zpz1ήMw\xFjv?zKM3JFmV^䣂~p+KT_$ +iOw$θNͧӜ!*Q%*?oMc~1cM^B!Tj%D$Ody 56t-y)DC#a99'brETvU>;\+iǻ7@eHզ)KoCN(rJ}Rjn8k\> 1ԩ[I=I.oYϹ6xVe0\&5-iCܞ]mp^tΥ`֕Z7Yj_N VVT4خ["Zvnn2ԂC:TB g4fηq k\o'BܠaIl $hUN&SKz wxFy. + d#%\dk;0tŠ뮣ҷzLk< :›KVR7 ȑRsU}IU6kMh\˃$݄<|NoYvrJVnP5Exӿ 3!u>e*)qAh->zcJ!(ZV򶒐Jp ;ZF};ov~zvΗϥ&~d]j}5n]oL}>y;o#q4-ggvM?fs3~M9YݿSO٦ߩN~woփMX o4KSr&[Zxqd8G6]JUL-pwx׫4w4ACApI-JU#.>F]'v:Te)vHx'yo8HMK{Eóc҂TܔQ8 +83|{/L7 }^s!R +r|JyKS`1KOPOT2?k%4ӟ4iZ~4-ggvM?fs3~M9YݿS]v{a9ժuR)9H\RAEp4*RYwU#j+y/RRzH$~{GWqv+KfZ:IBǘJ?)IǺGvZmli 8I'چRA#WMHJBBRtV&r$u;\AZ}R} ԥ)JW9L5&3NABӒ20G*fE[vX+qJQd !RJR)\'jt7"ۃ +,>V?'"-pb-mZEa[Pd9<ԊTŦwC(0^KqINRAT|/UXu^v=X$q藚i<vqۭVdf]"ChR WRYd +0oln1T +[+,%GPv0|Xϐݵ0Ztc#ʭoˑ .΅ܟ$垨sWo7Qnmf$[cv‚înIP+^#Q\+V&vMH.FnZ}թ\WP4q) )$M-?.:D.!%F7G $Wukiґqm-eID)RGZT|[ ̖:=@2q퓏3SmBVRd(0Gep4ĕ/O+vw8w=lQ# +zrJMVUs{Z[ڌ\hAԶsˣ!,8ӿX &WgƉԁ=X7ǾOVR-%+AP +_;Fy85__&-%2N$mGm %{PH9P$(z߫En%"Vmv2Ҥ]^ +qEyxژODNJu/YZVfWsdڙaNp5)JR)JR)JR)JR+زmAm֜HR0RA8 hMF˜ +ESX$Y9 R99'&ҕv[: S 8c:33>/>MubobcXNHZ9 kn=ثXEt+#Ԥ3p|"])JTEۢ*싨odĴY.$rI#>Dcr;2)1JCI Vy9 dҸMsAXPPCͅp}kV[Qݶ}%8rT3k8FDhqٌ3$kD"}SK)JR)JR)JR)JR)JR)JR)JR \ No newline at end of file diff -Naur lammps-3Nov06/doc/Eqs/pair_sw.tex lammps-4Nov06/doc/Eqs/pair_sw.tex --- lammps-3Nov06/doc/Eqs/pair_sw.tex 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/Eqs/pair_sw.tex 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,18 @@ +\documentclass[12pt]{article} + +\begin{document} + +\begin{eqnarray*} + E & = & \sum_i \sum_{j > i} \phi_2 (r_{ij}) + + \sum_i \sum_{j \neq i} \sum_{k > j} + \phi_3 (r_{ij}, r_{ik}, \theta_{ijk}) \\ + \phi_2(r) & = & A \epsilon \left[ B (\frac{\sigma}{r})^p - + (\frac{\sigma}{r})^q \right] + \exp \left( \frac{\sigma}{r - a \sigma} \right) \\ + \phi_3(r,s,\theta) & = & \lambda \epsilon \left[ \cos \theta - + \cos \theta_0 \right]^2 + \exp \left( \frac{\gamma \sigma}{r - a \sigma} \right) + \exp \left( \frac{\gamma \sigma}{s - a \sigma} \right) +\end{eqnarray*} + +\end{document} \ No newline at end of file diff -Naur lammps-3Nov06/doc/Eqs/pair_tersoff.jpg lammps-4Nov06/doc/Eqs/pair_tersoff.jpg --- lammps-3Nov06/doc/Eqs/pair_tersoff.jpg 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/Eqs/pair_tersoff.jpg 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,58 @@ +JFIFHHC  !"$"$ GV !1"AQa2q#BW3Rbr$%u45SU&'67CDEVcfs?˥)JR)JR)JR)JR)JR)JR)JR)JRMkR~=[ɹ~V;`qd[Ud .o. l^yC/3FN?*-"ۚ7 H4v' +<mf"%@9R.?\?Xo_.N% h X+AJR)JR)JR)JR)JRc*rx[,qC/̊شT%aV#61fb 1d;-JYU|x･)JƜ^R2q%ٵM]}>{~*g'/f#_\B)o G_-'>Iq_]s{)FO<*eXp~ȸ#[EV DA_qeCCq[YM2?> i.Gsb"}u\ck+hOw>RDOv} ݢ=߱f+~EZR)JR)JR)JR)JR|򙟊f~;5 iypeU(ףKN{@/yNVL66`AXUSr'dv~n :1;aM>3V_^ +HxGK|GU`/Ň|ˌt6~߷Y[aao1%ZyW*@K9I>X]:kui90^fzuF=Ccbd.# l<@'}g1YO4K(1?e8 2ܨD#fl@Mmd[*уh-l./.#gX R~$cT{N'|n,+ +͒/o=4_12SL8K49L$7I@> P@i#i'lhuAWOm[R)JR)JR)JR)JR)JRnQ佛 9 Yۻ;4*`f_:c|cxIvWs o_ cpH5> +ݥ|@3s TP<}]s+! K\%՛d7h!7Y+U+ZY:]{V ] ){~³i 讜+ a~A|*^]qeydWK nkL~INm+եBVJ&c@e߳hɨ\ )E詒VX0OUULkެJR)JR)JR)JR)JR)JG9e|8bnA~n ޑ=zUI.26rYqp4X(uy}†c}Qsxn#a$8ԃS#I։e+^rƻL9d|b5{Hm{դY]P$uPlF &jcf K8cF13d[K}R)JR)JR)JR)JR)JR)JTk}w66[ H-'U'_o]ن3pk'_+JJ]iJR)JR)JR)JR)JR]ZC A~avgHdSoS7C[irη9("k٬I.%I5HII GUV?X'saۮ{zPlVn;,t 5u HPXupHa e&Fy&XngsiccVDb Q}Roo}6RڽC9̲vpam~P$aH\s% @MD#zCF}:>Gno,f[)U=k kDZ~LFxcuݍ)w!UNz'^vٌ/t6KXѮ\/29dxl;Aٞ$q(&5,6{>تN1~`gԓ:4$ʁ> 'dRvJ\7!$kGVd'0Gkk UO%]:)H":s>yݭ*s Y$GGdlɓiĖdi!S/= + ?Yr2|R˓d)nP\Kiq+fnj55EBA ^vM/\Mc;جmRv~ŕFڀjVW=Kce}ta{X-F}nH"axŕ6)Wh2nBEW2&ؙ{y +I$߀(~1#DUv]&[``=eQ%wvhֺ)JR)JR+)s]N?/.Oo,J$ʒ+<(?[b-,R[*-̺iVUK.Kw_<^SZ++H_NHY+e-{p֩ ; ;7yl s$ͳ&7ǫ}5_ʈTm?5>wrU)&9G+uSF d{}S]:zw}lѿ!Y Xi?tHZ \e6ױ~ȣ׬d\pܯes3wY K2$,{;we]]1&WcxY伋}R8@vF%~Ͻט-~oQ7+d?t4RV8%Z󕿲sJ+[;\K<ҷU]OU-f\EʹZIZ h +۴J$mVkpA3vw3yefYZBщ߅:L1c&$FMHt $."u>c(8{KT٣"-=c}W3sxfB XoS&y)&m?gt=7sWLǺcgXHOD~ڴ8l+>:i@}FmgpZGq0br}[HƽHbDFy;+>yi̖%0m|PzccBk`21>Q_?xou_NUQ$I'7ZLmem[ʒI7DTO(>vAx:5ȥ#ղ]OwP}! Pu[)8D=ӡ'y# Z&~ߝ~KgqhY&9^ ؘ̲HW} ]vտn.~g1%ݴfZ%fˠ5:@kW>77%S7wmnblJR)JR)JR}2m16VXXtU{?ҔC13^OmdW"8Q}6n uN: (h(O$k1u }Ccmi%K;0>GEAh3:vMEnS7@$jLl zUF@_MzJu^5>Ŧk+kVFc JGf,Y<ē$1G KHiT{^Jd%Xd"Į6<cHf'$hў=b߾ڽрe*F?z$Q,QTE +}Wvw +:V GGkз$1c % m +/VLxD%¢gOI}l¥Ҕ)JR)JR)JR|yGWO]^\lI,m~UOھJR)JRgmsMŋ7BBZWMt5޵*4 > io-%V.Cu:߾MpyuCrM<,/d;;k**M[ޏ;*#;UQbt^?\xP7lN +{9u t.Clď)t rFXg%6x> %gM93aHXݏ:l(b%B$ j{(@WR)JR yl|)JRc1<]~]`J_N t >B&u^yU: }"ɐDI-U2"f+5"-m.@i!R4ݙ`+85-HRRJ#|޷ 2&U,BcAU>Qͥ*_ڋ6x^1q +{GqSz7+.]qÑJUCbAafޔ)JR)JR)JRys.sqm%i/={ |!؂Fυ:{?X[;,Wg?!͂_^Is !eU'L}a)l2|UջG{27I!S3nI[l_Ad-267-gIm*t|+g0[1nSj]ͽ"kY"t7 |g V?6[ }^p+k9E1 mI$qv^l~zv[)^3>E0-q+ABje\mZ>_[D +aX|; Cgta2ey/^~*f.kv}nԫ͙Lyxѕ!z ёIq˂r%yfqDYze44Ac:+9ĭ?xRfezoq#=5I{_;'u|O3v8W/wg%dDNT;zL)!|G("X/mϖrMCqB[Bbu,x#vhx<ƕ"ii cQpQ$t]S_ -ox,̈|"WV:?/q-¬+yVH:<1\srX*+wr?a4B +k P;՝味gL7'GI=H 0 aS8aNVF7WN;ɱŏЈUcolvbv_ Zπp.H<^﹀Ɨ.ޮ?R{[m1B=_#Ϊ/,:tNxl~35)y4>Bqf{BS*Aj`qV<5#!ICI}I$I5-DZ_>H(UNm`DHj녇1oz +;(v<}Kv~טsae.c5v~{k^=?sXi?M>;Ĩd"$ǿj3hFŐ6HqkY^cw^c[-CkRJb('+>1?ⶼ}+U`}^ViGv}<=s3eƶ .JFhPDI $y=U2YlvV(>+/E-nI/1X$DZXy\BaܻJ :*``ýњKF^AUQX@UQ6I,Yޔ(@#D~i@;JRS':52=S87Xel2,ś_cF oz'[ B,1)r*$ܒk*q<0凔M ~RƐ"՘l&}ɫRyE.3%F9%H{<ED`ݬki䳒],o񣪕Yqpl>c,T\q!,Ř3;I;bOȰAIKLwD_"rc1geBA* +HQTb-ܳ|: +g߾bM)JRK4HƃlMp]lS+?> JxXD4gv(j=Sw/i[_RJR)JR)JRcσ;߷jR)JR+69dBٽ[wioܢ +|} ݩP,rrfoqa +Dm.v,=HQݾکK3<Zs"=53D +@;$v_eccDvI l/*t<1qfK=ͣ]IgvR_ePZ3*ɨNImso5f<֒i']vI޶C~TTGv @D8FA;ly>Ӑe_qXme$gw\*QPK Glc,k tmtRGԔ[>^+ss f[rGbb `mJR)JW8Qw`?4X&D)2:0*F{Ie G +:׈P7 }:"Ҕr[g/VhLPo{VmoQ H8Px{+J+e1uiI68D]t>++31+$oҔ)JR)JR)JR)JGk[0Po>B[۔R;Bԍ ^@מk$4ad8elZ[7u#鷶Do(ZwoڴHcWrr̾K5H(leHcR,A>H +[λAuqiO#[QHd,F֦f/nlm[\MQ˅0ڼ*4Oc.kߑUgL Ĭk :'C+񸜟M7fZI)a`A'>}ͺ8k--nmm!x'ݙ?qJV -ynq,0\8$4~PFmtjF/\܃Wtw": +dNT|79738Xe3X=n'WGAAQ{!ddwe$uzD +O!'jC A'JK=N۪" Z|6[{oYx沎iU^N"|T'X^?c7]X~ib'4AI.qLkv(-'8rr=7@slx&%})JR)JR)JR)JR&)?XcZC3)v2\GV[{f_':*k3|"9~@E.[+gX[7O:/uӖspZ㮣@7}DO$X-;յykٯgyXi$m` +:Tlu +9zhG c S}eIob YN+3l|eFq"uʁ 0 މMVxd$Re`v!|Wa 6=ڪm6NFe=RF@uU oonaK0Q5]ix+dgfcwG) |vm$yWffk~MyGƈm;m/o$Cɨ)Y#ť!Ul1,rA-SI7<0-HgmK[ůT?}>Us`'$i h]ܒ{ȋRĝ*k渥.ӊpWx@-v*;ک* rm(~D^fk%Uf#;9*5XD5&|_ gnd.]zH?2oo[)JR)JR)JR)JR쨌UFى\nmYgxXF#ԲS# F^R(c +N^JRbK1Oڡ[7,6K)ecDVc56xXN8uG"1<Dљa)$#DCKn+:tE{)JR)JR)JR)JV#,rqv(d9|I +#2*hNޛ^\Gf[,UɷHIJZe}mj{m}xK^2M!7U@?IZf_%Zq$x6W2$1_eU젟$-]Ow8ߧĭ($2Y#kSP66-&(˜m^u']^53Y&H߄V` 5T\j~iQ%S= =2^AN9w|}O,H o 5("GS' ƛ& e. +8%pd, &URznH +k뼯ǯ𘅶!ⶌ ,qG*t’Ka*n>SMa!%Ry] w(2 }%{Fq9KxW&NfM(}Bu~y6gIrXݍüh#x: "7lP"5\,9^;KKRy}>uYy \:.[ުH*_m|[U9x q  =wbk2,Fӌ#dtY!bH(e )o:|SŦ<%oYG!D[?KZR)JR)JR)JRd5>kK?u=r~fWF~ D4ª~0a&,VXV摯of`Y<(UQ+׹6$}M'Ql}&f18ͲYċ *0#z$o5[ LCpKGLt! E_ҕIݿ/Y.3y[ĖL)ԫ(. \A.m-d?ѿX{$::KlXٚh/q4w0"QǐAVRA|wks}e,G=GR~H'cx)1v9 &ʒIOnƥaHb73DZ^usio ɡ ZLu]pIo#/ mSHa邏87Q 䘄m*[=1j[蛉"c{k5[rFҲݘw[L,wjH׷^,HuU*sktD>6>JR)JR)JR)JR)JR)JR)JR/wHbU1:\!o1e)JR)JR)JR)JRif3ȫ!qбo`kT3K'Q_=5E)Q.u*)Ul;kx s%NQdeAclmE>~H7٬# {V\Km":sj]C\Ggi鋉QYQTIN6>|V/q ZCCy]mg+";,e^"mS"?^5ĮR +ه5'5ŒId, !OoIrDż},?Ս$}_QD(H +rgLmd d[3*>| ^I~ߊ[࢞ώar\b%H}Y EO[r866܊ W +wG rL[*lerw{`h2cGԪT<{WjR)JR)JR)JR+)[s Oo=HԳ5GaΉ[V=}:c1̵X{H$z[ƻ^ 7]B+KpZ/T*"N?v&ʩ}C660b1XoԁYHV\s=I"N]| R |2xdyK{74 {}5/;[;;-K ݬ:6ς/VxG)웓&@ߟ?uYm^CVUީ{L +LDa#K UJWm5z#aqy{/,bHԺ) $lnǎW9$TI.l[1iNlG` yKmu4bn9Q`|o_ ʐ$6uHhUZ6VB0&߉f.&pEa;&{[]cP5ߑ32WABuцU] ܲ%ؕbOMO_Ui4 0NKJR)JR)JR)JR)PXR^Bd{unC+(oov;}'B;#yE>0{1 -.jޔ*-SIn +-̾?tHFzϹ$ʯ3F"tuYHpkr;hqf,MuWPwU4 *)JR)JR)JR)JRGh Ian")_399&7G{)n'5ISaxZK˯`5Yv,6ȱٝc_t]`>g\Eb$)??JkuL8gSF:#Tсj\\A{vU^ +^N ?5eEErcEފ1`#pmOV:5a^eN&7U'J6Nkͬ[Esm*M$D;WR6?pEt+3b7Wk7Awh[-YܸR񲶻"o*'ciܤW7cd^̪}u@ҝK"VXIj$r-ҴD2k_O߷bI}$JԾ ,""]Ԯ h쟿Hn_alKt~wj^AoY5,_lnԮWv^ZMiu sG,r.Ն}U߀8lvۯT+Įmwb627K1(eRAK`6V U98̋;+[ \~5N4=P:Mwܬ;؅Y:ו'~|kUfa, +4\Kl%z+{|̦3&)*}EP?+|pVykp2+Ű& 2:?\-$[Lѡ[ Y$ zAv, *E{\HJ9m2luݐ7߷[g1'wYߎwilb79kmnr +дodl_I-Z_ks-Q]E13\8X X=At#ܼnlBtk nȓ2;m/nz@Mk0s|%;Diegv:†: pyQKAk}aY)y([jA:ao}3${^A;\39?٩~^w$jP>F!azuE'~*119Gdߵuf* yR4DkVxۜvfJwq:L0fޅ~q$oagB&̑L =v<{0Y߾r+!^[@`h&!C%@ Y] ki X}kJR)JU#M%O#FFVU` YG̛{8,.Mi5[%u + n@${avw%15(њ8 ޙ+֚}NƎEJRXeܔҔ$--#ED"TbPSwXe7g b8bX,t=wORo+< #89y$*D6`g'dTxa⮤oq~Co& oI~b9%U ћK%#;%V09P4IIRͳǟ5S1wpi]J#pnሷ;5WDdž@џ>}ǃ$c\su/ym䃲 +lB2N{귶R"2A%ɋLlQNϷ` P0@cSH"9qm;6 ^@mưXO +H]bԂ>bIrOI=8)#[Wԕ +< +mָᎵ#1E7F R]\zZӯQhx9N%9r+Lj8Mg,Ogy>|>mq8, Y0cvvCUq,U`(>yoğY᱐_䯒 Ow%UDhml£8,c1w8zj8' #$րwعt+*W Ч{:y>2xb-U("#d4gϳq]s4ol L/dدo l5Z*T;e^-4,+y d|lK +\8JwAHcgu+'dJ##8Xќ/ .\^-^]\KwtP=IPOPBaGPAR)JR)JRcs'NS҉ʘv'e6o)782뱦u׊n\5V#۬<_LDatSx֍sSf2uJ'}>L)JR)JRsX)JR)JR+de P[Z8+q{hGodDbDVWj.kq<7K9m'\HFhBIn}*[,ťџx.[X-PR^`B.=1z45qF1-,Vt)JR)JR)JR)JR)JR)JR)Uq-SI=ǣJ[nb_7q"yRCtaaTR)JRVԶJ&ۨHZS+ ;vgAƂ)JR)JR)JR$@4 x)JRļ{[孅 e8_S coZEf ]:E --ױ+'͢魟5ۛ; {Cq=ɂ%;N@]JR)J9?9/8\f"Zs$*Ġ8эDגJWGr1\fGmɮy~VU#kXK95^RX_Tu)JR)JR)JR X4(boFR,ϲL6?q V8,rI'TZEZ|C}"ڨA~oǎr?#3aE@@,Uef0d uVM{X^-QJU5|<5ܑ 62@`ڰ߶UgKKKKKKKK稗~ ܥ:7IeZ3zߚ|,g|,g|,g|,g|,g|,g|,g|,g|,gne!9ҵ4uRǪ+t 5)JR)JR|e9y哔emms6ڲV:DFv3>Ԟ[{#g7Иoq폖Opz aA¢ c={lu:K[G'QVE:V}-*|E",z;~Bo\{-Bc-ou%h[٭z#Dy'^߁UqT3$w,2?B4j!R,KGfѴF80;G]X;kD~*LEl]7?og$K^=<}>դ_]ZN7c'{& +{3'?&mMinޔtm;+hmYNU+u1gəm 1%M +dsxluV/Y-F煯5Dt3ƭbBly^S +mk)PƱʑҗ6?_: +ᏼX[&fA> &=-!mr UWs)΁JUFw8;"'i>Z'A5IHEڂ@'ןxjב|O!Y)e ӄ#tueJR-q%V\iwq{ifnzqqR\|j4v}@5%dAwi(ŃN ;J=65N;l> 7^ zU@0D:T0oƶ>b[mwr9i9q]ClrmlHAUU/(37%\A>Gj8F3al]Vq۫,H kgk?kEe;y 3ĩߴ P߰ +r{052wwV0ek&kT.:PH钱lr\?0Q``HMf UH6TǃٔxF熴}K R83+\Z l̪@z*]Ԓ@t _o9x7QEąVTD,[>?} 䙜:]IŷYc3"$ *FF!V鬫EJR)JR)JR)JR)TMgJTFU!ףkُ`?)JᏳ[Y!a1술 +ƦA-9nmp +egC1֎JR)Jc{ +o2؏"Dž{mym&rܙYKT;9$EP7:jqmXa{"(Gy'195kX-؇WN#@(>M\R)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)_ \ No newline at end of file diff -Naur lammps-3Nov06/doc/Eqs/pair_tersoff.tex lammps-4Nov06/doc/Eqs/pair_tersoff.tex --- lammps-3Nov06/doc/Eqs/pair_tersoff.tex 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/Eqs/pair_tersoff.tex 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,24 @@ +\documentclass[12pt]{article} + +\begin{document} + +\begin{eqnarray*} + E & = & \frac{1}{2} \sum_i \sum_{j \neq i} V_{ij} \\ + V_{ij} & = & f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) \right] \\ + f_C(r) & = & \left\{ \begin{array} {r@{\quad:\quad}l} + 1 & r < R - D \\ + \frac{1}{2} - \frac{1}{2} \sin \left( \frac{\pi}{2} \frac{r-R}{D} \right) & + R-D < r < R + D \\ + 0 & r > R + D + \end{array} \right. \\ + f_R(r) & = & A \exp (-\lambda_1 r) \\ + f_A(r) & = & -B \exp (-\lambda_2 r) \\ + b_{ij} & = & \left( 1 + \beta^n {\zeta_{ij}}^n \right)^{-\frac{1}{2n}} \\ + \zeta_{ij} & = & \sum_{k \neq i,j} f_C(r_{ik}) g(\theta_{ijk}) + \exp \left[ {\lambda_3}^3 (r_{ij} - r_{ik})^3 \right] \\ + g(\theta) & = & 1 + \frac{c^2}{d^2} - + \frac{c^2}{\left[ d^2 + + (\cos \theta - \cos \theta_0)^2\right]} +\end{eqnarray*} + +\end{document} \ No newline at end of file diff -Naur lammps-3Nov06/doc/pair_style.html lammps-4Nov06/doc/pair_style.html --- lammps-3Nov06/doc/pair_style.html 2006-09-21 10:22:34.000000000 -0600 +++ lammps-4Nov06/doc/pair_style.html 2006-11-03 10:09:44.000000000 -0700 @@ -125,6 +125,8 @@
  • pair_style lj/smooth - smoothed Lennard-Jones potential
  • pair_style morse - Morse potential
  • pair_style soft - Soft (cosine) potential +
  • pair_style sw - Stillinger-Weber 3-body potential +
  • pair_style tersoff - Tersoff 3-body potential
  • pair_style table - tabulated pair potential
  • pair_style yukawa - Yukawa potential diff -Naur lammps-3Nov06/doc/pair_style.txt lammps-4Nov06/doc/pair_style.txt --- lammps-3Nov06/doc/pair_style.txt 2006-09-21 10:22:34.000000000 -0600 +++ lammps-4Nov06/doc/pair_style.txt 2006-11-03 10:09:44.000000000 -0700 @@ -122,6 +122,8 @@ "pair_style lj/smooth"_pair_style_lj_smooth.html - smoothed Lennard-Jones potential "pair_style morse"_pair_style_morse.html - Morse potential "pair_style soft"_pair_style_soft.html - Soft (cosine) potential +"pair_style sw"_pair_style_sw.html - Stillinger-Weber 3-body potential +"pair_style tersoff"_pair_style_tersoff.html - Tersoff 3-body potential "pair_style table"_pair_style_table.html - tabulated pair potential "pair_style yukawa"_pair_style_yukawa.html - Yukawa potential :ul diff -Naur lammps-3Nov06/doc/pair_style_eam.html lammps-4Nov06/doc/pair_style_eam.html --- lammps-3Nov06/doc/pair_style_eam.html 2006-11-01 16:01:32.000000000 -0700 +++ lammps-4Nov06/doc/pair_style_eam.html 2006-11-03 10:09:44.000000000 -0700 @@ -151,30 +151,34 @@ values in it for exponents. C only recognizes "e" or "E" for scientific notation.

    -

    Only one pair_coeff command can be used which specifies a single -setfl file. DYNAMO setfl files contain information for M -elements. These are mapped to LAMMPS atom types by specifying N -additional arguments after the filename, where N is the number of -LAMMPS atom types: +

    Only a single pair_coeff command is used with the eam/alloy style +which specifies a DYNAMO setfl file, which contains information for +M elements. These are mapped to LAMMPS atom types by specifying N +additional arguments after the filename in the pair_coeff command, +where N is the number of LAMMPS atom types:

    • filename
    • N element names = mapping of setfl elements to atom types
    -

    As an example, the nialhjea setfl file has tabulated EAM values for -3 elements and their alloy interactions: Ni, Al, and H. If your -LAMMPS simulation had 4 atoms types and you wanted the 1st 3 to be Ni, -and the 4th to be Al, you would use the following pair_coeff command: +

    As an example, the potentials/nialhjea setfl file has tabulated EAM +values for 3 elements and their alloy interactions: Ni, Al, and H. If +your LAMMPS simulation has 4 atoms types and you want the 1st 3 to be +Ni, and the 4th to be Al, you would use the following pair_coeff +command: 

    pair_coeff * * nialhjea.eam.alloy Ni Ni Ni Al

    The 1st 2 arguments must be * * so as to span all LAMMPS atom types. The first three Ni arguments map LAMMPS atom types 1,2,3 to the Ni element in the setfl file. The final Al argument maps LAMMPS atom -type 4 to the Al element in the setfl file. If a mapping value is -specified as NULL, the mapping is not performed. This can be used -when an eam/alloy potential is used as part of the hybrid pair -style. The NULL values are placeholders for atom types that will be -used with other potentials. +type 4 to the Al element in the setfl file. Note that there is no +requirement that your simulation use all the elements specified by the +setfl file. +

    +

    If a mapping value is specified as NULL, the mapping is not performed. +This can be used when an eam/alloy potential is used as part of the +hybrid pair style. The NULL values are placeholders for atom types +that will be used with other potentials.

    Setfl files in the potentials directory of the LAMMPS distribution have a ".eam.alloy" suffix. A DYNAMO multi-element setfl file is diff -Naur lammps-3Nov06/doc/pair_style_eam.txt lammps-4Nov06/doc/pair_style_eam.txt --- lammps-3Nov06/doc/pair_style_eam.txt 2006-11-01 16:01:32.000000000 -0700 +++ lammps-4Nov06/doc/pair_style_eam.txt 2006-11-03 10:09:44.000000000 -0700 @@ -146,30 +146,34 @@ values in it for exponents. C only recognizes "e" or "E" for scientific notation. -Only one pair_coeff command can be used which specifies a single -{setfl} file. DYNAMO {setfl} files contain information for M -elements. These are mapped to LAMMPS atom types by specifying N -additional arguments after the filename, where N is the number of -LAMMPS atom types: +Only a single pair_coeff command is used with the {eam/alloy} style +which specifies a DYNAMO {setfl} file, which contains information for +M elements. These are mapped to LAMMPS atom types by specifying N +additional arguments after the filename in the pair_coeff command, +where N is the number of LAMMPS atom types: filename N element names = mapping of {setfl} elements to atom types :ul -As an example, the nialhjea {setfl} file has tabulated EAM values for -3 elements and their alloy interactions: Ni, Al, and H. If your -LAMMPS simulation had 4 atoms types and you wanted the 1st 3 to be Ni, -and the 4th to be Al, you would use the following pair_coeff command: +As an example, the potentials/nialhjea {setfl} file has tabulated EAM +values for 3 elements and their alloy interactions: Ni, Al, and H. If +your LAMMPS simulation has 4 atoms types and you want the 1st 3 to be +Ni, and the 4th to be Al, you would use the following pair_coeff +command: pair_coeff * * nialhjea.eam.alloy Ni Ni Ni Al :pre The 1st 2 arguments must be * * so as to span all LAMMPS atom types. The first three Ni arguments map LAMMPS atom types 1,2,3 to the Ni element in the {setfl} file. The final Al argument maps LAMMPS atom -type 4 to the Al element in the {setfl} file. If a mapping value is -specified as NULL, the mapping is not performed. This can be used -when an {eam/alloy} potential is used as part of the {hybrid} pair -style. The NULL values are placeholders for atom types that will be -used with other potentials. +type 4 to the Al element in the {setfl} file. Note that there is no +requirement that your simulation use all the elements specified by the +{setfl} file. + +If a mapping value is specified as NULL, the mapping is not performed. +This can be used when an {eam/alloy} potential is used as part of the +{hybrid} pair style. The NULL values are placeholders for atom types +that will be used with other potentials. {Setfl} files in the {potentials} directory of the LAMMPS distribution have a ".eam.alloy" suffix. A DYNAMO multi-element {setfl} file is diff -Naur lammps-3Nov06/doc/pair_style_sw.html lammps-4Nov06/doc/pair_style_sw.html --- lammps-3Nov06/doc/pair_style_sw.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/pair_style_sw.html 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,112 @@ + +

    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    pair_style sw command +

    +

    Syntax: +

    +
    pair_style sw 
+
    +

    Examples: +

    +
    pair_style sw
+pair_coeff * * si.sw Si
+pair_coeff * * SiC.sw Si C Si 
+
    +

    Description: +

    +

    The sw style computes a 3-body Stillinger-Weber +potential for the energy E of a system of atoms as +

    +
    +
    +

    where phi2 is a two-body term and phi3 is a three-body term. The +summations in the formula are over all neighbors J and K of atom I +within a cutoff distance = a*sigma. +

    +

    Only a single pair_coeff command is used with the sw style which +specifies a Stillinger-Weber potential file with parameters for all +needed elements. These are mapped to LAMMPS atom types by specifying +N additional arguments after the filename in the pair_coeff command, +where N is the number of LAMMPS atom types: +

    +
    • filename +
    • N element names = mapping of SW elements to atom types +
    +

    As an example, imagine the SiC.sw file has Stillinger-Weber values for +Si and C. If your LAMMPS simulation has 4 atoms types and you want +the 1st 3 to be Si, and the 4th to be C, you would use the following +pair_coeff command: +

    +
    pair_coeff * * SiC.sw Si Si Si C 
+
    +

    The 1st 2 arguments must be * * so as to span all LAMMPS atom types. +The first three Si arguments map LAMMPS atom types 1,2,3 to the Si +element in the SW file. The final C argument maps LAMMPS atom type 4 +to the C element in the SW file. If a mapping value is specified as +NULL, the mapping is not performed. This can be used when a sw +potential is used as part of the hybrid pair style. The NULL values +are placeholders for atom types that will be used with other +potentials. +

    +

    Stillinger-Weber files in the potentials directory of the LAMMPS +distribution have a ".sw" suffix. Lines that are not blank or +comments (starting with #) define parameters for a triplet of +elements. The parameters in a single entry correspond to coefficients +in the formula above: +

    +
    • element 1 (the center atom in a 3-body interaction) +
    • element 2 +
    • element 3 +
    • epsilon (energy units) +
    • sigma (distnace units) +
    • a +
    • lambda +
    • gamma +
    • costheta0 +
    • A +
    • B +
    • p +
    • q +
    +

    The non-annotated parameters are unitless. +

    +

    The Stillinger-Weber potential file must contain entries for all the +elements listed in the pair_coeff command. It can also contain +entries for additional elements not being used in a particular +simulation; LAMMPS ignores those entries. +

    +

    For a single-element simulation, only a single entry is required +(e.g. SiSiSi). For a two-element simulation, the file must contain 8 +entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that +specify SW parameters for all permutations of the two elements +interacting in three-body configurations. Thus for 3 elements, 27 +entries would be required, etc. Note that due to symmetries, some +parameters will typically be the same in multiple entries. +

    +

    Restrictions: +

    +

    This pair potential requires the newton setting to be +"on" for pair interactions. +

    +

    Related commands: +

    +

    pair_coeff +

    +

    Default: none +

    +
    + + + +

    (Stillinger) Stillinger and Weber, Phys Rev B, 31, 5262 (1985). +

    + diff -Naur lammps-3Nov06/doc/pair_style_sw.txt lammps-4Nov06/doc/pair_style_sw.txt --- lammps-3Nov06/doc/pair_style_sw.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/pair_style_sw.txt 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,106 @@ +"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 + +pair_style sw command :h3 + +[Syntax:] + +pair_style sw :pre + +[Examples:] + +pair_style sw +pair_coeff * * si.sw Si +pair_coeff * * SiC.sw Si C Si :pre + +[Description:] + +The {sw} style computes a 3-body "Stillinger-Weber"_#Stillinger +potential for the energy E of a system of atoms as + +:c,image(Eqs/pair_sw.jpg) + +where phi2 is a two-body term and phi3 is a three-body term. The +summations in the formula are over all neighbors J and K of atom I +within a cutoff distance = a*sigma. + +Only a single pair_coeff command is used with the {sw} style which +specifies a Stillinger-Weber potential file with parameters for all +needed elements. These are mapped to LAMMPS atom types by specifying +N additional arguments after the filename in the pair_coeff command, +where N is the number of LAMMPS atom types: + +filename +N element names = mapping of SW elements to atom types :ul + +As an example, imagine the SiC.sw file has Stillinger-Weber values for +Si and C. If your LAMMPS simulation has 4 atoms types and you want +the 1st 3 to be Si, and the 4th to be C, you would use the following +pair_coeff command: + +pair_coeff * * SiC.sw Si Si Si C :pre + +The 1st 2 arguments must be * * so as to span all LAMMPS atom types. +The first three Si arguments map LAMMPS atom types 1,2,3 to the Si +element in the SW file. The final C argument maps LAMMPS atom type 4 +to the C element in the SW file. If a mapping value is specified as +NULL, the mapping is not performed. This can be used when a {sw} +potential is used as part of the {hybrid} pair style. The NULL values +are placeholders for atom types that will be used with other +potentials. + +Stillinger-Weber files in the {potentials} directory of the LAMMPS +distribution have a ".sw" suffix. Lines that are not blank or +comments (starting with #) define parameters for a triplet of +elements. The parameters in a single entry correspond to coefficients +in the formula above: + +element 1 (the center atom in a 3-body interaction) +element 2 +element 3 +epsilon (energy units) +sigma (distnace units) +a +lambda +gamma +costheta0 +A +B +p +q :ul + +The non-annotated parameters are unitless. + +The Stillinger-Weber potential file must contain entries for all the +elements listed in the pair_coeff command. It can also contain +entries for additional elements not being used in a particular +simulation; LAMMPS ignores those entries. + +For a single-element simulation, only a single entry is required +(e.g. SiSiSi). For a two-element simulation, the file must contain 8 +entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that +specify SW parameters for all permutations of the two elements +interacting in three-body configurations. Thus for 3 elements, 27 +entries would be required, etc. Note that due to symmetries, some +parameters will typically be the same in multiple entries. + +[Restrictions:] + +This pair potential requires the "newton"_newton.html setting to be +"on" for pair interactions. + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] none + +:line + +:link(Stillinger) +[(Stillinger)] Stillinger and Weber, Phys Rev B, 31, 5262 (1985). diff -Naur lammps-3Nov06/doc/pair_style_tersoff.html lammps-4Nov06/doc/pair_style_tersoff.html --- lammps-3Nov06/doc/pair_style_tersoff.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/pair_style_tersoff.html 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,114 @@ + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    pair_style tersoff command +

    +

    Syntax: +

    +
    pair_style tersoff 
+
    +

    Examples: +

    +
    pair_style tersoff
+pair_coeff * * si.tersoff Si
+pair_coeff * * SiC.tersoff Si C Si 
+
    +

    Description: +

    +

    The tersoff style computes a 3-body Tersoff potential +for the energy E of a system of atoms as +

    +
    +
    +

    where f_R is a two-body term and f_A includes three-body interactions. +The summations in the formula are over all neighbors J and K of atom I +within a cutoff distance = R + D. +

    +

    Only a single pair_coeff command is used with the tersoff style +which specifies a Tersoff potential file with parameters for all +needed elements. These are mapped to LAMMPS atom types by specifying +N additional arguments after the filename in the pair_coeff command, +where N is the number of LAMMPS atom types: +

    +
    • filename +
    • N element names = mapping of Tersoff elements to atom types +
    +

    As an example, imagine the SiC.tersoff file has Tersoff values for Si +and C. If your LAMMPS simulation has 4 atoms types and you want the +1st 3 to be Si, and the 4th to be C, you would use the following +pair_coeff command: +

    +
    pair_coeff * * SiC.tersoff Si Si Si C 
+
    +

    The 1st 2 arguments must be * * so as to span all LAMMPS atom types. +The first three Si arguments map LAMMPS atom types 1,2,3 to the Si +element in the Tersoff file. The final C argument maps LAMMPS atom +type 4 to the C element in the Tersoff file. If a mapping value is +specified as NULL, the mapping is not performed. This can be used +when a tersoff potential is used as part of the hybrid pair style. +The NULL values are placeholders for atom types that will be used with +other potentials. +

    +

    Tersoff files in the potentials directory of the LAMMPS distribution +have a ".tersoff" suffix. Lines that are not blank or comments +(starting with #) define parameters for a triplet of elements. The +parameters in a single entry correspond to coefficients in the formula +above: +

    +
    • element 1 (the center atom in a 3-body interaction) +
    • element 2 +
    • element 3 +
    • lambda3 (1/distance units) +
    • c +
    • d +
    • costheta0 +
    • n +
    • beta +
    • lambda2 (1/distance units) +
    • B (energy units) +
    • R (distance units) +
    • D (distance units) +
    • lambda1 (1/distance units) +
    • A (energy units) +
    +

    The non-annotated parameters are unitless. +

    +

    The Tersoff potential file must contain entries for all the elements +listed in the pair_coeff command. It can also contain entries for +additional elements not being used in a particular simulation; LAMMPS +ignores those entries. +

    +

    For a single-element simulation, only a single entry is required +(e.g. SiSiSi). For a two-element simulation, the file must contain 8 +entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that +specify Tersoff parameters for all permutations of the two elements +interacting in three-body configurations. Thus for 3 elements, 27 +entries would be required, etc. Note that due to symmetries, some +parameters will typically be the same in multiple entries. +

    +

    Restrictions: +

    +

    This pair potential requires the newton setting to be +"on" for pair interactions. +

    +

    Related commands: +

    +

    pair_coeff +

    +

    Default: none +

    +
    + + + +

    (Tersoff) Tersoff, Phys Rev B, 37, 6991 (1988). +

    + diff -Naur lammps-3Nov06/doc/pair_style_tersoff.txt lammps-4Nov06/doc/pair_style_tersoff.txt --- lammps-3Nov06/doc/pair_style_tersoff.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/doc/pair_style_tersoff.txt 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,108 @@ +"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 + +pair_style tersoff command :h3 + +[Syntax:] + +pair_style tersoff :pre + +[Examples:] + +pair_style tersoff +pair_coeff * * si.tersoff Si +pair_coeff * * SiC.tersoff Si C Si :pre + +[Description:] + +The {tersoff} style computes a 3-body "Tersoff"_#Tersoff potential +for the energy E of a system of atoms as + +:c,image(Eqs/pair_tersoff.jpg) + +where f_R is a two-body term and f_A includes three-body interactions. +The summations in the formula are over all neighbors J and K of atom I +within a cutoff distance = R + D. + +Only a single pair_coeff command is used with the {tersoff} style +which specifies a Tersoff potential file with parameters for all +needed elements. These are mapped to LAMMPS atom types by specifying +N additional arguments after the filename in the pair_coeff command, +where N is the number of LAMMPS atom types: + +filename +N element names = mapping of Tersoff elements to atom types :ul + +As an example, imagine the SiC.tersoff file has Tersoff values for Si +and C. If your LAMMPS simulation has 4 atoms types and you want the +1st 3 to be Si, and the 4th to be C, you would use the following +pair_coeff command: + +pair_coeff * * SiC.tersoff Si Si Si C :pre + +The 1st 2 arguments must be * * so as to span all LAMMPS atom types. +The first three Si arguments map LAMMPS atom types 1,2,3 to the Si +element in the Tersoff file. The final C argument maps LAMMPS atom +type 4 to the C element in the Tersoff file. If a mapping value is +specified as NULL, the mapping is not performed. This can be used +when a {tersoff} potential is used as part of the {hybrid} pair style. +The NULL values are placeholders for atom types that will be used with +other potentials. + +Tersoff files in the {potentials} directory of the LAMMPS distribution +have a ".tersoff" suffix. Lines that are not blank or comments +(starting with #) define parameters for a triplet of elements. The +parameters in a single entry correspond to coefficients in the formula +above: + +element 1 (the center atom in a 3-body interaction) +element 2 +element 3 +lambda3 (1/distance units) +c +d +costheta0 +n +beta +lambda2 (1/distance units) +B (energy units) +R (distance units) +D (distance units) +lambda1 (1/distance units) +A (energy units) :ul + +The non-annotated parameters are unitless. + +The Tersoff potential file must contain entries for all the elements +listed in the pair_coeff command. It can also contain entries for +additional elements not being used in a particular simulation; LAMMPS +ignores those entries. + +For a single-element simulation, only a single entry is required +(e.g. SiSiSi). For a two-element simulation, the file must contain 8 +entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that +specify Tersoff parameters for all permutations of the two elements +interacting in three-body configurations. Thus for 3 elements, 27 +entries would be required, etc. Note that due to symmetries, some +parameters will typically be the same in multiple entries. + +[Restrictions:] + +This pair potential requires the "newton"_newton.html setting to be +"on" for pair interactions. + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] none + +:line + +:link(Tersoff) +[(Tersoff)] Tersoff, Phys Rev B, 37, 6991 (1988). diff -Naur lammps-3Nov06/potentials/README lammps-4Nov06/potentials/README --- lammps-3Nov06/potentials/README 2006-09-27 13:51:23.000000000 -0600 +++ lammps-4Nov06/potentials/README 2006-11-03 10:09:44.000000000 -0700 @@ -4,7 +4,7 @@ of the file formats and the various styles in LAMMPS that read these files. -The suffix of each file tells what pair style it is used with: +The suffix of each file indicates the pair style it is used with: eam embedded atom method (EAM) single element, DYNAMO funcfl format eam.alloy EAM multi-element alloy, DYNAMO setfl format diff -Naur lammps-3Nov06/potentials/si.sw lammps-4Nov06/potentials/si.sw --- lammps-3Nov06/potentials/si.sw 2006-10-12 14:37:55.000000000 -0600 +++ lammps-4Nov06/potentials/si.sw 2006-11-03 10:09:44.000000000 -0700 @@ -1,10 +1,12 @@ # Stillinger-Weber parameters for various elements and mixtures # multiple entries can be added to this file, LAMMPS reads the ones it needs +# these entries are in LAMMPS "metal" units: +# epsilon = eV; sigma = Angstroms +# other quantities are unitless -# format of an entry (one or more lines): -# element 1, element 2, element 3 -# epsilon, sigma, littlea, lambda, gamma, costheta, biga, bigb, powerp, powerq - -Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333 - 7.049556277 0.6022245584 4.0 0.0 +# format of a single entry (one or more lines): +# element 1, element 2, element 3, +# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q +Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333 + 7.049556277 0.6022245584 4.0 0.0 diff -Naur lammps-3Nov06/potentials/si.tersoff lammps-4Nov06/potentials/si.tersoff --- lammps-3Nov06/potentials/si.tersoff 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/potentials/si.tersoff 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,12 @@ +# Tersoff parameters for various elements and mixtures +# multiple entries can be added to this file, LAMMPS reads the ones it needs +# these entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless + +# format of a single entry (one or more lines): +# element 1, element 2, element 3, +# lambda3, c, d, costheta0, n, beta, lambda2, B, R, D, lambda1, A + +Si Si Si 1.3258 4.8381 2.0417 0.0000 22.956 + 0.33675 1.3258 95.373 3.0 0.2 3.2394 3264.7 diff -Naur lammps-3Nov06/src/MANYBODY/Install.csh lammps-4Nov06/src/MANYBODY/Install.csh --- lammps-3Nov06/src/MANYBODY/Install.csh 2006-09-27 13:51:33.000000000 -0600 +++ lammps-4Nov06/src/MANYBODY/Install.csh 2006-11-03 10:09:44.000000000 -0700 @@ -9,10 +9,14 @@ cp pair_eam.cpp .. cp pair_eam_alloy.cpp .. cp pair_eam_fs.cpp .. + cp pair_sw.cpp .. + cp pair_tersoff.cpp .. # cp pair_eam.h .. cp pair_eam_alloy.h .. cp pair_eam_fs.h .. + cp pair_sw.h .. + cp pair_tersoff.h .. else if ($1 == 0) then @@ -22,9 +26,13 @@ rm ../pair_eam.cpp rm ../pair_eam_alloy.cpp rm ../pair_eam_fs.cpp + rm ../pair_sw.cpp + rm ../pair_tersoff.cpp # rm ../pair_eam.h rm ../pair_eam_alloy.h rm ../pair_eam_fs.h + rm ../pair_sw.h + rm ../pair_tersoff.h endif diff -Naur lammps-3Nov06/src/MANYBODY/pair_eam_alloy.cpp lammps-4Nov06/src/MANYBODY/pair_eam_alloy.cpp --- lammps-3Nov06/src/MANYBODY/pair_eam_alloy.cpp 2006-11-01 16:00:53.000000000 -0700 +++ lammps-4Nov06/src/MANYBODY/pair_eam_alloy.cpp 2006-11-03 10:09:44.000000000 -0700 @@ -54,7 +54,7 @@ memory->destroy_2d_double_array(setfl->frho); memory->destroy_2d_double_array(setfl->rhor); memory->destroy_3d_double_array(setfl->z2r); - memory->sfree(setfl); + delete setfl; } setfl = new Setfl(); read_file(arg[2]); @@ -139,7 +139,7 @@ if (nwords != file->nelements + 1) error->all("Incorrect element names in EAM potential file"); - char *words[file->nelements]; + char **words = new char*[file->nelements+1]; nwords = 0; char *first = strtok(line," \t\n\r\f"); while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; @@ -150,6 +150,7 @@ file->elements[i] = new char[n]; strcpy(file->elements[i],words[i]); } + delete [] words; if (me == 0) { fgets(line,MAXLINE,fp); diff -Naur lammps-3Nov06/src/MANYBODY/pair_eam_fs.cpp lammps-4Nov06/src/MANYBODY/pair_eam_fs.cpp --- lammps-3Nov06/src/MANYBODY/pair_eam_fs.cpp 2006-11-01 16:00:53.000000000 -0700 +++ lammps-4Nov06/src/MANYBODY/pair_eam_fs.cpp 2006-11-03 10:09:44.000000000 -0700 @@ -54,7 +54,7 @@ memory->destroy_2d_double_array(fs->frho); memory->destroy_3d_double_array(fs->rhor); memory->destroy_3d_double_array(fs->z2r); - memory->sfree(fs); + delete fs; } fs = new Fs(); read_file(arg[2]); @@ -139,7 +139,7 @@ if (nwords != file->nelements + 1) error->all("Incorrect element names in EAM potential file"); - char *words[file->nelements]; + char **words = new char*[file->nelements+1]; nwords = 0; char *first = strtok(line," \t\n\r\f"); while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; @@ -150,6 +150,7 @@ file->elements[i] = new char[n]; strcpy(file->elements[i],words[i]); } + delete [] words; if (me == 0) { fgets(line,MAXLINE,fp); diff -Naur lammps-3Nov06/src/MANYBODY/pair_sw.cpp lammps-4Nov06/src/MANYBODY/pair_sw.cpp --- lammps-3Nov06/src/MANYBODY/pair_sw.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/src/MANYBODY/pair_sw.cpp 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,557 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_sw.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "update.h" +#include "memory.h" +#include "neighbor.h" +#include "memory.h" +#include "error.h" + +#define MAXLINE 1024 +#define DELTA 4 + +/* ---------------------------------------------------------------------- */ + +PairSW::PairSW() +{ + neigh_half_every = 0; + neigh_full_every = 1; + single_enable = 0; + one_coeff = 1; + + nelements = 0; + elements = NULL; + nparams = 0; + maxparam = 0; + params = NULL; + elem2param = NULL; +} + +/* ---------------------------------------------------------------------- + free all arrays + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairSW::~PairSW() +{ + if (elements) + for (int i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + memory->sfree(params); + memory->destroy_3d_int_array(elem2param); + + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + delete [] map; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::compute(int eflag, int vflag) +{ + int i,j,k,m,n,itag,jtag,itype,jtype,ktype,iparam,numneigh; + double xtmp,ytmp,ztmp,delx,dely,delz; + double rsq,rsq1,rsq2,eng,fforce; + double delr1[3],delr2[3],fj[3],fk[3]; + int *neighs; + double **f; + + eng_vdwl = 0.0; + if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; + + if (vflag == 2) f = update->f_pair; + else f = atom->f; + double **x = atom->x; + int *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + + // loop over full neighbor list of my atoms + + for (i = 0; i < nlocal; i++) { + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // two-body interactions, skip half of them + + neighs = neighbor->firstneigh_full[i]; + numneigh = neighbor->numneigh_full[i]; + + for (m = 0; m < numneigh; m++) { + j = neighs[m]; + jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) + continue; + } + + jtype = map[type[j]]; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + iparam = elem2param[itype][jtype][jtype]; + if (rsq > params[iparam].cutsq) continue; + + twobody(¶ms[iparam],rsq,fforce,eflag,eng); + + if (eflag) eng_vdwl += eng; + f[i][0] += fforce*delx; + f[i][1] += fforce*dely; + f[i][2] += fforce*delz; + f[j][0] -= fforce*delx; + f[j][1] -= fforce*dely; + f[j][2] -= fforce*delz; + } + + // three-body interactions + // cannot test I-J distance against cutoff outside of 2nd loop + // b/c must use I-J-K cutoff for both rij and rik + + for (m = 0; m < numneigh-1; m++) { + j = neighs[m]; + jtype = map[type[j]]; + + for (n = m+1; n < numneigh; n++) { + k = neighs[n]; + ktype = map[type[k]]; + iparam = elem2param[itype][jtype][ktype]; + + delr1[0] = x[j][0] - xtmp; + delr1[1] = x[j][1] - ytmp; + delr1[2] = x[j][2] - ztmp; + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + if (rsq1 > params[iparam].cutsq) continue; + + delr2[0] = x[k][0] - xtmp; + delr2[1] = x[k][1] - ytmp; + delr2[2] = x[k][2] - ztmp; + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + if (rsq2 > params[iparam].cutsq) continue; + + threebody(¶ms[iparam],rsq1,rsq2,delr1,delr2,fj,fk,eflag,eng); + + if (eflag) eng_vdwl += eng; + f[i][0] -= fj[0] + fk[0]; + f[i][1] -= fj[1] + fk[1]; + f[i][2] -= fj[2] + fk[2]; + f[j][0] += fj[0]; + f[j][1] += fj[1]; + f[j][2] += fj[2]; + f[k][0] += fk[0]; + f[k][1] += fk[1]; + f[k][2] += fk[2]; + } + } + } + if (vflag == 2) virial_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + map = new int[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairSW::settings(int narg, char **arg) +{ + if (narg != 0) error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairSW::coeff(int narg, char **arg) +{ + int i,j,m,n; + + if (!allocated) allocate(); + + if (narg != 3 + atom->ntypes) + error->all("Incorrect args for pair coefficients"); + + // insure I,J args are * * + + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all("Incorrect args for pair coefficients"); + + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if NULL + // nelements = # of unique elements + // elements = list of element names + + if (elements) { + for (i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + } + elements = new char*[atom->ntypes]; + for (i = 0; i < atom->ntypes; i++) elements[i] = NULL; + + nelements = 0; + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < nelements; j++) + if (strcmp(arg[i],elements[j]) == 0) break; + map[i-2] = j; + if (j == nelements) { + n = strlen(arg[i]) + 1; + elements[j] = new char[n]; + strcpy(elements[j],arg[i]); + nelements++; + } + } + + // read potential file and initialize potential parameters + + read_file(arg[2]); + setup(); + + // clear setflag since coeff() called once with I,J = * * + + n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + + int count = 0; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + count++; + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairSW::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); + + return cutmax; +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::init_style() +{ + if (atom->tag_enable == 0) + error->all("Pair style Stillinger-Weber requires atom IDs"); + if (force->newton_pair == 0) + error->all("Pair style Stillinger-Weber requires newton pair on"); +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::read_file(char *file) +{ + int params_per_line = 13; + + if (params) delete [] params; + params = NULL; + nparams = 0; + + // open file on proc 0 + + FILE *fp; + if (comm->me == 0) { + fp = fopen(file,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open Stillinger-Weber potential file %s",file); + error->one(str); + } + } + + // read each set of params from potential file + // one set of params can span multiple lines + // store params if all 3 element tags are in element list + + int i,n,nwords,ielement,jelement,kelement; + char *words[params_per_line]; + char line[MAXLINE],*ptr; + int eof = 0; + + while (1) { + if (comm->me == 0) { + ptr = fgets(line,MAXLINE,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + + // strip comment, skip line if blank + + if (ptr = strchr(line,'#')) *ptr = '\0'; + nwords = atom->count_words(line); + if (nwords == 0) continue; + + // concatenate additional lines until have params_per_line words + + while (nwords < params_per_line) { + n = strlen(line); + if (comm->me == 0) { + ptr = fgets(&line[n],MAXLINE-n,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + if (ptr = strchr(line,'#')) *ptr = '\0'; + nwords = atom->count_words(line); + } + + if (nwords != params_per_line) + error->all("Incorrect format in Stillinger-Weber potential file"); + + // words = ptrs to all words in line + + nwords = 0; + words[nwords++] = strtok(line," \t\n\r\f"); + while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + + for (ielement = 0; ielement < nelements; ielement++) + if (strcmp(words[0],elements[ielement]) == 0) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (strcmp(words[1],elements[jelement]) == 0) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (strcmp(words[2],elements[kelement]) == 0) break; + if (kelement == nelements) continue; + + // load up parameter settings and error check their values + + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].epsilon = atof(words[3]); + params[nparams].sigma = atof(words[4]); + params[nparams].littlea = atof(words[5]); + params[nparams].lambda = atof(words[6]); + params[nparams].gamma = atof(words[7]); + params[nparams].costheta = atof(words[8]); + params[nparams].biga = atof(words[9]); + params[nparams].bigb = atof(words[10]); + params[nparams].powerp = atof(words[11]); + params[nparams].powerq = atof(words[12]); + + if (params[nparams].epsilon < 0.0 || params[nparams].sigma <= 0.0 || + params[nparams].littlea <= 0.0 || params[nparams].lambda < 0.0 || + params[nparams].gamma < 0.0 || params[nparams].costheta > 1.0 || + params[nparams].costheta < -1.0 || params[nparams].biga < 0.0 || + params[nparams].bigb < 0.0 || params[nparams].powerp < 0.0 || + params[nparams].powerq < 0.0) + error->all("Illegal Stillinger-Weber parameter"); + + nparams++; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::setup() +{ + int i,j,k,m,n; + + // set elem2param for all triplet combinations + // must be a single exact match to lines read from file + // do not allow for ACB in place of ABC + + if (elem2param) memory->destroy_3d_int_array(elem2param); + elem2param = memory->create_3d_int_array(nelements,nelements,nelements, + "pair:elem2param"); + + for (i = 0; i < nelements; i++) + for (j = 0; j < nelements; j++) + for (k = 0; k < nelements; k++) { + n = -1; + for (m = 0; m < nparams; m++) { + if (i == params[m].ielement && j == params[m].jelement && + k == params[m].kelement) { + if (n >= 0) error->all("Potential file has duplicate entry"); + n = m; + } + } + if (n < 0) error->all("Potential file is missing an entry"); + elem2param[i][j][k] = n; + } + + + // compute parameter values derived from inputs + + for (m = 0; m < nparams; m++) { + params[m].cut = params[m].sigma*params[m].littlea; + params[m].cutsq = params[m].cut*params[m].cut; + params[m].sigma_gamma = params[m].sigma*params[m].gamma; + params[m].lambda_epsilon = params[m].lambda*params[m].epsilon; + params[m].lambda_epsilon2 = 2.0*params[m].lambda*params[m].epsilon; + params[m].c1 = params[m].biga*params[m].epsilon * + params[m].powerp*params[m].bigb * + pow(params[m].sigma,params[m].powerp); + params[m].c2 = params[m].biga*params[m].epsilon*params[m].powerq * + pow(params[m].sigma,params[m].powerq); + params[m].c3 = params[m].biga*params[m].epsilon*params[m].bigb * + pow(params[m].sigma,params[m].powerp+1.0); + params[m].c4 = params[m].biga*params[m].epsilon * + pow(params[m].sigma,params[m].powerq+1.0); + params[m].c5 = params[m].biga*params[m].epsilon*params[m].bigb * + pow(params[m].sigma,params[m].powerp); + params[m].c6 = params[m].biga*params[m].epsilon * + pow(params[m].sigma,params[m].powerq); + } + + // set cutmax to max of all params + + cutmax = 0.0; + for (m = 0; m < nparams; m++) + if (params[m].cut > cutmax) cutmax = params[m].cut; +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::twobody(Param *param, double rsq, double &fforce, + int eflag, double &eng) +{ + double r,rinvsq,rp,rq,rainv,rainvsq,expsrainv; + double delr[3]; + + r = sqrt(rsq); + rinvsq = 1.0/rsq; + rp = pow(r,-param->powerp); + rq = pow(r,-param->powerq); + rainv = 1.0 / (r - param->cut); + rainvsq = rainv*rainv*r; + expsrainv = exp(param->sigma * rainv); + fforce = (param->c1*rp - param->c2*rq + + (param->c3*rp -param->c4*rq) * rainvsq) * expsrainv * rinvsq; + if (eflag) eng = (param->c5*rp - param->c6*rq) * expsrainv; +} + +/* ---------------------------------------------------------------------- */ + +void PairSW::threebody(Param *param, double rsq1, double rsq2, + double *delr1, double *delr2, + double *fj, double *fk, int eflag, double &eng) +{ + double r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1; + double r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2; + double rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2; + double facang,facang12,csfacang,csfac1,csfac2; + + r1 = sqrt(rsq1); + rinvsq1 = 1.0/rsq1; + rainv1 = 1.0/(r1 - param->cut); + gsrainv1 = param->sigma_gamma * rainv1; + gsrainvsq1 = gsrainv1*rainv1/r1; + expgsrainv1 = exp(gsrainv1); + + r2 = sqrt(rsq2); + rinvsq2 = 1.0/rsq2; + rainv2 = 1.0/(r2 - param->cut); + gsrainv2 = param->sigma_gamma * rainv2; + gsrainvsq2 = gsrainv2*rainv2/r2; + expgsrainv2 = exp(gsrainv2); + + rinv12 = 1.0/(r1*r2); + cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12; + delcs = cs - param->costheta; + delcssq = delcs*delcs; + + facexp = expgsrainv1*expgsrainv2; + facrad = param->lambda_epsilon * facexp*delcssq; + frad1 = facrad*gsrainvsq1; + frad2 = facrad*gsrainvsq2; + facang = param->lambda_epsilon2 * facexp*delcs; + facang12 = rinv12*facang; + csfacang = cs*facang; + csfac1 = rinvsq1*csfacang; + + fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12; + fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12; + fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12; + + csfac2 = rinvsq2*csfacang; + + fk[0] = delr2[0]*(frad2+csfac2)-delr1[0]*facang12; + fk[1] = delr2[1]*(frad2+csfac2)-delr1[1]*facang12; + fk[2] = delr2[2]*(frad2+csfac2)-delr1[2]*facang12; + + if (eflag) eng = facrad; +} diff -Naur lammps-3Nov06/src/MANYBODY/pair_sw.h lammps-4Nov06/src/MANYBODY/pair_sw.h --- lammps-3Nov06/src/MANYBODY/pair_sw.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/src/MANYBODY/pair_sw.h 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,58 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 PAIR_SW_H +#define PAIR_SW_H + +#include "pair.h" + +class PairSW : public Pair { + public: + PairSW(); + ~PairSW(); + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + + private: + struct Param { + double epsilon,sigma; + double littlea,lambda,gamma,costheta; + double biga,bigb; + double powerp,powerq; + double cut,cutsq; + double sigma_gamma,lambda_epsilon,lambda_epsilon2; + double c1,c2,c3,c4,c5,c6; + int ielement,jelement,kelement; + }; + + double cutmax; // max cutoff for all elements + int nelements; // # of unique elements + char **elements; // names of unique elements + int ***elem2param; // mapping from element triplets to parameters + int *map; // mapping from atom types to elements + int nparams; // # of stored parameter sets + int maxparam; // max # of parameter sets + Param *params; // parameter set for an I-J-K interaction + + void allocate(); + void read_file(char *); + void setup(); + void twobody(Param *, double, double &, int, double &); + void threebody(Param *, double, double, double *, double *, + double *, double *, int, double &); +}; + +#endif diff -Naur lammps-3Nov06/src/MANYBODY/pair_tersoff.cpp lammps-4Nov06/src/MANYBODY/pair_tersoff.cpp --- lammps-3Nov06/src/MANYBODY/pair_tersoff.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/src/MANYBODY/pair_tersoff.cpp 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,792 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_tersoff.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "update.h" +#include "memory.h" +#include "neighbor.h" +#include "memory.h" +#include "error.h" + +#define MAXLINE 1024 +#define DELTA 4 + +/* ---------------------------------------------------------------------- */ + +PairTersoff::PairTersoff() +{ + neigh_half_every = 0; + neigh_full_every = 1; + single_enable = 0; + one_coeff = 1; + + PI2 = 2.0*atan(1.0); + PI4 = atan(1.0); + + nelements = 0; + elements = NULL; + nparams = 0; + maxparam = 0; + params = NULL; + elem2param = NULL; +} + +/* ---------------------------------------------------------------------- + free all arrays + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairTersoff::~PairTersoff() +{ + if (elements) + for (int i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + delete [] params; + memory->destroy_3d_int_array(elem2param); + + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + delete [] map; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::compute(int eflag, int vflag) +{ + int i,j,k,m,n,itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk,numneigh; + double xtmp,ytmp,ztmp,delx,dely,delz; + double rsq,rsq1,rsq2,eng,fforce; + double delr1[3],delr2[3],fi[3],fj[3],fk[3]; + double zeta_ij,prefactor; + int *neighs; + double **f; + + eng_vdwl = 0.0; + if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; + + if (vflag == 2) f = update->f_pair; + else f = atom->f; + double **x = atom->x; + int *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + + // loop over full neighbor list of my atoms + + for (i = 0; i < nlocal; i++) { + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // two-body interactions, skip half of them + + neighs = neighbor->firstneigh_full[i]; + numneigh = neighbor->numneigh_full[i]; + + for (m = 0; m < numneigh; m++) { + j = neighs[m]; + jtag = tag[j]; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < x[i][2]) continue; + else if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) + continue; + } + + jtype = map[type[j]]; + + delx = x[j][0] - xtmp; + dely = x[j][1] - ytmp; + delz = x[j][2] - ztmp; + rsq = delx*delx + dely*dely + delz*delz; + + iparam_ij = elem2param[itype][jtype][jtype]; + if (rsq > params[iparam_ij].cutsq) continue; + + repulsive(¶ms[iparam_ij],rsq,fforce,eflag,eng); + + if (eflag) eng_vdwl += eng; + f[i][0] += fforce*delx; + f[i][1] += fforce*dely; + f[i][2] += fforce*delz; + f[j][0] -= fforce*delx; + f[j][1] -= fforce*dely; + f[j][2] -= fforce*delz; + } + + // three-body interactions + // skip immediately if I-J is not within cutoff + + for (m = 0; m < numneigh; m++) { + j = neighs[m]; + jtype = map[type[j]]; + iparam_ij = elem2param[itype][jtype][jtype]; + + delr1[0] = x[j][0] - xtmp; + delr1[1] = x[j][1] - ytmp; + delr1[2] = x[j][2] - ztmp; + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + if (rsq1 > params[iparam_ij].cutsq) continue; + + // accumulate bondorder zeta for each i-j interaction via loop over k + + zeta_ij = 0.0; + for (n = 0; n < numneigh; n++) { + if (n == m) continue; + k = neighs[n]; + ktype = map[type[k]]; + iparam_ijk = elem2param[itype][jtype][ktype]; + + delr2[0] = x[k][0] - xtmp; + delr2[1] = x[k][1] - ytmp; + delr2[2] = x[k][2] - ztmp; + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + if (rsq2 > params[iparam_ijk].cutsq) continue; + + zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2); + } + + // pairwise force due to zeta + + force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fforce,prefactor,eflag,eng); + + if (eflag) eng_vdwl += eng; + f[i][0] += fforce*delr1[0]; + f[i][1] += fforce*delr1[1]; + f[i][2] += fforce*delr1[2]; + f[j][0] -= fforce*delr1[0]; + f[j][1] -= fforce*delr1[1]; + f[j][2] -= fforce*delr1[2]; + + // attractive term via loop over k + + for (n = 0; n < numneigh; n++) { + if (n == m) continue; + k = neighs[n]; + ktype = map[type[k]]; + iparam_ijk = elem2param[itype][jtype][ktype]; + + delr2[0] = x[k][0] - xtmp; + delr2[1] = x[k][1] - ytmp; + delr2[2] = x[k][2] - ztmp; + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + if (rsq2 > params[iparam_ijk].cutsq) continue; + + attractive(¶ms[iparam_ijk],prefactor, + rsq1,rsq2,delr1,delr2,fi,fj,fk); + + f[i][0] += fi[0]; + f[i][1] += fi[1]; + f[i][2] += fi[2]; + f[j][0] += fj[0]; + f[j][1] += fj[1]; + f[j][2] += fj[2]; + f[k][0] += fk[0]; + f[k][1] += fk[1]; + f[k][2] += fk[2]; + } + } + } + if (vflag == 2) virial_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + map = new int[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairTersoff::settings(int narg, char **arg) +{ + if (narg != 0) error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairTersoff::coeff(int narg, char **arg) +{ + int i,j,m,n; + + if (!allocated) allocate(); + + if (narg != 3 + atom->ntypes) + error->all("Incorrect args for pair coefficients"); + + // insure I,J args are * * + + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all("Incorrect args for pair coefficients"); + + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if NULL + // nelements = # of unique elements + // elements = list of element names + + if (elements) { + for (i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + } + elements = new char*[atom->ntypes]; + for (i = 0; i < atom->ntypes; i++) elements[i] = NULL; + + nelements = 0; + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < nelements; j++) + if (strcmp(arg[i],elements[j]) == 0) break; + map[i-2] = j; + if (j == nelements) { + n = strlen(arg[i]) + 1; + elements[j] = new char[n]; + strcpy(elements[j],arg[i]); + nelements++; + } + } + + // read potential file and initialize potential parameters + + read_file(arg[2]); + setup(); + + // clear setflag since coeff() called once with I,J = * * + + n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + + int count = 0; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + count++; + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairTersoff::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all("All pair coeffs are not set"); + + return cutmax; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::init_style() +{ + if (atom->tag_enable == 0) + error->all("Pair style Tersoff requires atom IDs"); + if (force->newton_pair == 0) + error->all("Pair style Tersoff requires newton pair on"); +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::read_file(char *file) +{ + int params_per_line = 15; + + if (params) delete [] params; + params = NULL; + nparams = 0; + + // open file on proc 0 + + FILE *fp; + if (comm->me == 0) { + fp = fopen(file,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open Tersoff potential file %s",file); + error->one(str); + } + } + + // read each line out of file, skipping blank lines or leading '#' + // store line of params if all 3 element tags are in element list + + int i,n,nwords,ielement,jelement,kelement; + char *words[params_per_line]; + char line[MAXLINE],*ptr; + int eof = 0; + + while (1) { + if (comm->me == 0) { + ptr = fgets(line,MAXLINE,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + + // strip comment, skip line if blank + + if (ptr = strchr(line,'#')) *ptr = '\0'; + nwords = atom->count_words(line); + if (nwords == 0) continue; + + // concatenate additional lines until have params_per_line words + + while (nwords < params_per_line) { + n = strlen(line); + if (comm->me == 0) { + ptr = fgets(&line[n],MAXLINE-n,fp); + if (ptr == NULL) { + eof = 1; + fclose(fp); + } else n = strlen(line) + 1; + } + MPI_Bcast(&eof,1,MPI_INT,0,world); + if (eof) break; + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + if (ptr = strchr(line,'#')) *ptr = '\0'; + nwords = atom->count_words(line); + } + + if (nwords != params_per_line) + error->all("Incorrect format in Tersoff potential file"); + + // words = ptrs to all words in line + + nwords = 0; + words[nwords++] = strtok(line," \t\n\r\f"); + while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue; + + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next line + + for (ielement = 0; ielement < nelements; ielement++) + if (strcmp(words[0],elements[ielement]) == 0) break; + if (ielement == nelements) continue; + for (jelement = 0; jelement < nelements; jelement++) + if (strcmp(words[1],elements[jelement]) == 0) break; + if (jelement == nelements) continue; + for (kelement = 0; kelement < nelements; kelement++) + if (strcmp(words[2],elements[kelement]) == 0) break; + if (kelement == nelements) continue; + + // load up parameter settings and error check their values + + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].lam3 = atof(words[3]); + params[nparams].c = atof(words[4]); + params[nparams].d = atof(words[5]); + params[nparams].h = atof(words[6]); + params[nparams].powern = atof(words[7]); + params[nparams].beta = atof(words[8]); + params[nparams].lam2 = atof(words[9]); + params[nparams].bigb = atof(words[10]); + params[nparams].bigr = atof(words[11]); + params[nparams].bigd = atof(words[12]); + params[nparams].lam1 = atof(words[13]); + params[nparams].biga = atof(words[14]); + + if (params[nparams].lam3 < 0.0 || params[nparams].c <= 0.0 || + params[nparams].d <= 0.0 || params[nparams].h > 1.0 || + params[nparams].h < -1.0 || params[nparams].powern <= 0.0 || + params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 || + params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 || + params[nparams].bigd < 0.0 || + params[nparams].bigd >= params[nparams].bigr || + params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0) + error->all("Illegal Tersoff parameter"); + + nparams++; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::setup() +{ + int i,j,k,m,n; + + // set elem2param for all element triplet combinations + // must be a single exact match to lines read from file + // do not allow for ACB in place of ABC + + if (elem2param) memory->destroy_3d_int_array(elem2param); + elem2param = memory->create_3d_int_array(nelements,nelements,nelements, + "pair:elem2param"); + + for (i = 0; i < nelements; i++) + for (j = 0; j < nelements; j++) + for (k = 0; k < nelements; k++) { + n = -1; + for (m = 0; m < nparams; m++) { + if (i == params[m].ielement && j == params[m].jelement && + k == params[m].kelement) { + if (n >= 0) error->all("Potential file has duplicate entry"); + n = m; + } + } + if (n < 0) error->all("Potential file is missing an entry"); + elem2param[i][j][k] = n; + } + + + // compute parameter values derived from inputs + + for (m = 0; m < nparams; m++) { + params[m].cut = params[m].bigr + params[m].bigd; + params[m].cutsq = params[m].cut*params[m].cut; + + params[m].c1 = pow(2.0*params[m].powern*1.0e-16,-1.0/params[m].powern); + params[m].c2 = pow(2.0*params[m].powern*1.0e-8,-1.0/params[m].powern); + params[m].c3 = 1.0/params[m].c2; + params[m].c4 = 1.0/params[m].c1; + } + + // set cutmax to max of all params + + cutmax = 0.0; + for (m = 0; m < nparams; m++) + if (params[m].cut > cutmax) cutmax = params[m].cut; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::repulsive(Param *param, double rsq, double &fforce, + int eflag, double &eng) +{ + double r,tmp_fc,tmp_fc_d,tmp_exp; + + r = sqrt(rsq); + tmp_fc = ters_fc(r,param); + tmp_fc_d = ters_fc_d(r,param); + tmp_exp = exp(-param->lam1 * r); + fforce = param->biga * tmp_exp * (tmp_fc_d - tmp_fc*param->lam1) / r; + + if (eflag) eng = tmp_fc * param->biga * tmp_exp; +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::zeta(Param *param, double rsqij, double rsqik, + double *delrij, double *delrik) +{ + double rij,rik,costheta,arg,ex_delr; + + rij = sqrt(rsqij); + rik = sqrt(rsqik); + costheta = (delrij[0]*delrik[0] + delrij[1]*delrik[1] + + delrij[2]*delrik[2]) / (rij*rik); + + arg = pow(param->lam3 * (rij-rik),3); + if (arg > 69.0776) ex_delr = 1.e30; + else if (arg < -69.0776) ex_delr = 0.0; + else ex_delr = exp(arg); + + return ters_fc(rik,param) * ters_gijk(costheta,param) * ex_delr; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::force_zeta(Param *param, double rsq, double zeta_ij, + double &fforce, double &prefactor, + int eflag, double &eng) +{ + double r,fa,fa_d,bij; + + r = sqrt(rsq); + fa = ters_fa(r,param); + fa_d = ters_fa_d(r,param); + bij = ters_bij(zeta_ij,param); + fforce = 0.5*bij*fa_d / r; + prefactor = -0.5*fa * ters_bij_d(zeta_ij,param); + if (eflag) eng = 0.5*bij*fa; +} + +/* ---------------------------------------------------------------------- + attractive term + use param_ij cutoff for rij test + use param_ijk cutoff for rik test +------------------------------------------------------------------------- */ + +void PairTersoff::attractive(Param *param, double prefactor, + double rsqij, double rsqik, + double *delrij, double *delrik, + double *fi, double *fj, double *fk) +{ + double rij_hat[3],rik_hat[3]; + double rij,rijinv,rik,rikinv; + + rij = sqrt(rsqij); + rijinv = 1.0/rij; + vec3_scale(rijinv,delrij,rij_hat); + + rik = sqrt(rsqik); + rikinv = 1.0/rik; + vec3_scale(rikinv,delrik,rik_hat); + + ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,param); +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_fc(double r, Param *param) +{ + double ters_R = param->bigr; + double ters_D = param->bigd; + + if (r < ters_R-ters_D) return 1.0; + if (r > ters_R+ters_D) return 0.0; + return 0.5*(1.0 - sin(PI2*(r - ters_R)/ters_D)); +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_fc_d(double r, Param *param) +{ + double ters_R = param->bigr; + double ters_D = param->bigd; + + if (r < ters_R-ters_D) return 0.0; + if (r > ters_R+ters_D) return 0.0; + return -(PI4/ters_D) * cos(PI2*(r - ters_R)/ters_D); +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_fa(double r, Param *param) +{ + if (r > param->bigr + param->bigd) return 0.0; + return -param->bigb * exp(-param->lam2 * r) * ters_fc(r,param); +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_fa_d(double r, Param *param) +{ + if (r > param->bigr + param->bigd) return 0.0; + return param->bigb * exp(-param->lam2 * r) * + (param->lam2 * ters_fc(r,param) - ters_fc_d(r,param)); +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_bij(double zeta, Param *param) +{ + double tmp = param->beta * zeta; + if (tmp > param->c1) return 1.0/sqrt(tmp); + if (tmp > param->c2) + return (1.0 - pow(tmp,-param->powern) / (2.0*param->powern))/sqrt(tmp); + if (tmp < param->c4) return 1.0; + if (tmp < param->c3) + return 1.0 - pow(tmp,param->powern)/(2.0*param->powern); + return pow(1.0 + pow(tmp,param->powern), -1.0/(2.0*param->powern)); +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_bij_d(double zeta, Param *param) +{ + double tmp = param->beta * zeta; + if (tmp > param->c1) return param->beta * -0.5*pow(tmp,-1.5); + if (tmp > param->c2) + return param->beta * (-0.5*pow(tmp,-1.5) * + (1.0 - 0.5*(1.0 + 1.0/(2.0*param->powern)) * + pow(tmp,-param->powern))); + if (tmp < param->c4) return 0.0; + if (tmp < param->c3) + return -0.5*param->beta * pow(tmp,param->powern-1.0); + + double tmp_n = pow(tmp,param->powern); + return -0.5 * pow(1.0+tmp_n, -1.0-(1.0/(2.0*param->powern)))*tmp_n / zeta; +} + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_gijk(double costheta, Param *param) +{ + double ters_c = param->c; + double ters_d = param->d; + + return 1.0 + pow(ters_c/ters_d,2.0) - + pow(ters_c,2.0) / (pow(ters_d,2.0) + pow(param->h - costheta,2.0)); +}; + +/* ---------------------------------------------------------------------- */ + +double PairTersoff::ters_gijk_d(double costheta, Param *param) +{ + double numerator = -2.0 * pow(param->c,2) * (param->h - costheta); + double denominator = pow(pow(param->d,2.0) + + pow(param->h - costheta,2.0),2.0); + return numerator/denominator; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::ters_zetaterm_d(double prefactor, + double *rij_hat, double rij, + double *rik_hat, double rik, + double *dri, double *drj, double *drk, + Param *param) +{ + double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp; + double dcosdri[3],dcosdrj[3],dcosdrk[3]; + + fc = ters_fc(rik,param); + dfc = ters_fc_d(rik,param); + tmp = pow(param->lam3 * (rij-rik),3.0); + + if (tmp > 69.0776) ex_delr = 1.e30; + else if (tmp < -69.0776) ex_delr = 0.0; + else ex_delr = exp(tmp); + + ex_delr_d = (3.0*pow(param->lam3,3.0)) * pow(rij-rik,2.0)*ex_delr; + cos_theta = vec3_dot(rij_hat,rik_hat); + gijk = ters_gijk(cos_theta,param); + gijk_d = ters_gijk_d(cos_theta,param); + costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); + + // compute the derivative wrt Ri + // dri = -dfc*gijk*ex_delr*rik_hat; + // dri += fc*gijk_d*ex_delr*dcosdri; + // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat); + + vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri); + vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri); + vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri); + vec3_scale(prefactor,dri,dri); + + // compute the derivative wrt Rj + // drj = fc*gijk_d*ex_delr*dcosdrj; + // drj += fc*gijk*ex_delr_d*rij_hat; + + vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj); + vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj); + vec3_scale(prefactor,drj,drj); + + // compute the derivative wrt Rk + // drk = dfc*gijk*ex_delr*rik_hat; + // drk += fc*gijk_d*ex_delr*dcosdrk; + // drk += -fc*gijk*ex_delr_d*rik_hat; + + vec3_scale(dfc*gijk*ex_delr,rik_hat,drk); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk); + vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk); + vec3_scale(prefactor,drk,drk); +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoff::costheta_d(double *rij_hat, double rij, + double *rik_hat, double rik, + double *dri, double *drj, double *drk) +{ + // first element is devative wrt Ri, second wrt Rj, third wrt Rk + + double cos_theta = vec3_dot(rij_hat,rik_hat); + + vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj); + vec3_scale(1.0/rij,drj,drj); + vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk); + vec3_scale(1.0/rik,drk,drk); + vec3_add(drj,drk,dri); + vec3_scale(-1.0,dri,dri); +}; + +/* ---------------------------------------------------------------------- + vector routines +------------------------------------------------------------------------- */ + +/* +double PairTersoff::vec3_dot(double *x, double *y) +{ + return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; +} + +void PairTersoff::vec3_add(double *x, double *y, double *z) +{ + z[0] = x[0]+y[0]; + z[1] = x[1]+y[1]; + z[2] = x[2]+y[2]; +} + +void PairTersoff::vec3_scale(double k, double *x, double *y) +{ + y[0] = k*x[0]; + y[1] = k*x[1]; + y[2] = k*x[2]; +} + +void PairTersoff::vec3_scaleadd(double k, double *x, double *y, double *z) +{ + z[0] = k*x[0]+y[0]; + z[1] = k*x[1]+y[1]; + z[2] = k*x[2]+y[2]; +} + +*/ diff -Naur lammps-3Nov06/src/MANYBODY/pair_tersoff.h lammps-4Nov06/src/MANYBODY/pair_tersoff.h --- lammps-3Nov06/src/MANYBODY/pair_tersoff.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-4Nov06/src/MANYBODY/pair_tersoff.h 2006-11-03 10:09:44.000000000 -0700 @@ -0,0 +1,88 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories + + 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 PAIR_TERSOFF_H +#define PAIR_TERSOFF_H + +#include "pair.h" + +class PairTersoff : public Pair { + public: + PairTersoff(); + ~PairTersoff(); + void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + + private: + struct Param { + double lam1,lam2,lam3; + double c,d,h; + double powern,beta; + double biga,bigb,bigd,bigr; + double cut,cutsq; + double c1,c2,c3,c4; + int ielement,jelement,kelement; + }; + + double PI2,PI4; + double cutmax; // max cutoff for all elements + int nelements; // # of unique elements + char **elements; // names of unique elements + int ***elem2param; // mapping from element triplets to parameters + int *map; // mapping from atom types to elements + int nparams; // # of stored parameter sets + int maxparam; // max # of parameter sets + Param *params; // parameter set for an I-J-K interaction + + void allocate(); + void read_file(char *); + void setup(); + void repulsive(Param *, double, double &, int, double &); + double zeta(Param *, double, double, double *, double *); + void force_zeta(Param *, double, double, double &, double &, int, double &); + void attractive(Param *, double, double, double, double *, double *, + double *, double *, double *); + + double ters_fc(double, Param *); + double ters_fc_d(double, Param *); + double ters_fa(double, Param *); + double ters_fa_d(double, Param *); + double ters_bij(double, Param *); + double ters_bij_d(double, Param *); + double ters_gijk(double, Param *); + double ters_gijk_d(double, Param *); + void ters_zetaterm_d(double, double *, double, double *, double, + double *, double *, double *, Param *); + void costheta_d(double *, double, double *, double, + double *, double *, double *); + + // vector functions, inline for efficiency + + inline double vec3_dot(double *x, double *y) { + return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; + } + inline void vec3_add(double *x, double *y, double *z) { + z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2]; + } + inline void vec3_scale(double k, double *x, double *y) { + y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2]; + } + inline void vec3_scaleadd(double k, double *x, double *y, double *z) { + z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2]; + } +}; + +#endif diff -Naur lammps-3Nov06/src/MANYBODY/style_manybody.h lammps-4Nov06/src/MANYBODY/style_manybody.h --- lammps-3Nov06/src/MANYBODY/style_manybody.h 2006-09-27 13:51:33.000000000 -0600 +++ lammps-4Nov06/src/MANYBODY/style_manybody.h 2006-11-03 10:09:44.000000000 -0700 @@ -15,10 +15,14 @@ #include "pair_eam.h" #include "pair_eam_alloy.h" #include "pair_eam_fs.h" +#include "pair_sw.h" +#include "pair_tersoff.h" #endif #ifdef PairClass PairStyle(eam,PairEAM) PairStyle(eam/alloy,PairEAMAlloy) PairStyle(eam/fs,PairEAMFS) +PairStyle(sw,PairSW) +PairStyle(tersoff,PairTersoff) #endif diff -Naur lammps-3Nov06/src/library.cpp lammps-4Nov06/src/library.cpp --- lammps-3Nov06/src/library.cpp 2006-10-31 15:16:12.000000000 -0700 +++ lammps-4Nov06/src/library.cpp 2006-11-03 10:09:44.000000000 -0700 @@ -24,7 +24,9 @@ LAMMPS *lammps; -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + minimal library interface functions +------------------------------------------------------------------------- */ void lammps_open(int argc, char **argv, MPI_Comm communicator) { @@ -54,9 +56,9 @@ return lammps->input->one(str); } -/* ---------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + add application-specific library functions here +------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ diff -Naur lammps-3Nov06/src/library.h lammps-4Nov06/src/library.h --- lammps-3Nov06/src/library.h 2006-09-27 13:51:33.000000000 -0600 +++ lammps-4Nov06/src/library.h 2006-11-03 10:09:44.000000000 -0700 @@ -11,15 +11,20 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -// prototypes for calling LAMMPS as a library +// library interface to LAMMPS +// new application-specific functions can be added #include "mpi.h" +extern "C" { + void lammps_open(int, char **, MPI_Comm); // start LAMMPS w/ command-line args void lammps_close(); // shut-down LAMMPS void lammps_file(char *); // execute an input script char *lammps_command(char *); // execute a single LAMMPS command int lammps_get_natoms(); // return # of atoms -void lammps_get_coords(double *); // retrieve list of coords on all procs -void lammps_put_coords(double *); // overwrite list of coords on all procs +void lammps_get_coords(double *); // retrieve atom coords from all procs +void lammps_put_coords(double *); // overwrite atom coords on all procs + +}