diff -Naur lammps-15May08/doc/Eqs/pair_gromacs.jpg lammps-16May08/doc/Eqs/pair_gromacs.jpg --- lammps-15May08/doc/Eqs/pair_gromacs.jpg 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/doc/Eqs/pair_gromacs.jpg 2008-05-15 14:44:18.000000000 -0600 @@ -0,0 +1,49 @@ +JFIFHHC  !"$"$ ~E !1"AQa#2q$BRu3457bDWdrs?˥)JR)JR)JR)JR)JR)JR)JR)J^.5VIJTu[HT )J֩1HIqhrck.:zybŪ7"ԕWB0rGj mO 7;HLn|ď2O'9A}o2[!+qCwy=R!.$jWZZhnudUJ {-q3cV5Rʊ-J$J#%J=Vk]e!=kK冟Bw*3)I[b ' 85Q^aYTqЖ,djYR ?x zV?vl T a²P+q;Tv4bB,"N\uq)tH5.vf^ >x\֍j-`qwmoM)JR)JR)JRqNtwV7mK|28.:T1{puf"\iW7L\Qܢ ڬ)#v9 +*,˃] ܵ(vZ$g^#\H},ԩYrOl$mDJu uQa.~s>O?Ajg"xs:ݘ+Bv3qiJOB9=VK>Ra<7 $,C8>U6a(f6pԓsap։*![δ^ëJ_mE(q VsV~jpǻrp7ΔcsFkg$ +}z5#k*R;a:TRJTpyh#AQ$sI$y$I'I'Y?WӶ8v|HABߙe@[:{wkEUV;wخ=1)*IùW<ΨqOu.12܃uڋh< +k]W26Gzbnu/ \yn6Pl:IBB| VѩX;W\6[QYknS;i;2@V* 8 Z[0u\IP.BfK2o(xqhFG "wYm3]cdC%zU{U)TQ۴ͷ.a$o-$(@*H'$ 8KWS1-ěp 8,%?6Gb+TFŠġ<$oVTDJWٲ_,'~ѵ* RQuͪpȗ)/Cզ,~[M8R)))\'qRdoK6-̓ -i,ҝZX' HυWw 3q[HZ%R eҢQCSk<FTYiN9%#vvvEpjm +#[ G^JcSK{,#q[2+\TcIh ~4)->M ҺIՖ۾0Zݍw=c3th}3lvu);V2AH^2f +xEG|ڝǮVγ*ڼ,ƃS#3!.wOi1.D$7&䇢b,t1|_ Z_ؼ:<%%6q%OcH P CT'W^)۞[t(mնZNrXh6Iwt;sT ̽Erf#Aǘmե;CIF7v!Fn"Rj%giu !Nd2}ǥWn>3>.lE8I}n!i$۸q*}|vK)OBx; +.Ѡ0 &"S9IV\GsT;V-}7t.bV3iRL~+BdzC#LFn%<u$)Yc.|' U~HE+ r)PKVZ/zS;"ÈYi?U-=sb>\2ViiFLrB[/vԒrx$$)JR)JRKWoJ $MS Tk'Oiڐ~2e? Ę&}!Qa8.,甓O=KXnqIKAn-~R;% M_rw^bvg?eATE}ֶduNۮ]bKv_ mO[P_/SLf Z_yi)e*9OQH<{z P?j_6#HyהhC|.*@H쑊kIW Hm +'$vHWnnzve| +w D2ެdhPgx.,z%{N=Ѓw.) [Tv M)JR)JR)JR)JR)TTU|dm. ےXn-% Fx睧ڦٮn[+2I)X~6* DA U&iHJ炠ke +Cl=m$cas +S@VMAydXv#Pw O ג8 Mm33szBU, ʖx4-Rv9FsU#޴.z>ez9!ܸ<ҲN+ST2ytWD@2Pe3+*=lhn + )G H($QnN[X}bkMVۉl%e*$(ؔ 8'FTGS%h%+H>*'؂;KPL4\?u%Gjs +&)JR)JR)JR|öWkّ֬MJ")ĂVo̦9V%#$٬1*Jn6.*G.)+)N9 QL7 QJe]BB~yF%@q9ujӔyh?/Csܾ"km|/v~;xq\ګxsPsdoLiwya[)ų#I9Pd閣}gR:2w[oؒ.6ƺ%*IR'`AWXC~C!TLۈ8'VIq]m/!i֥tqN!NSa9q̀WS&;?=A`Êฅm .%8899#>|T!Sk>Q(urʹa\~)֤nTZgunT[wHDS!nmv:%n +Nj-QZ#T,(qZYP[* +󝅷~:)WirzYf)ˉ %$F}c]4I, Ŷ0t%J"Bl0Bꖵ F9#S-ShSnQ=q^C +meQ$zU6.t\XoaxVԥc@4 J)JR)JR6SK+ P(YI ʳh}6Ɵ6#LjPȸ)Q%@B99]GPLeљ!# o);NP#ptK}G#$eJ$ 8)JjPLYmG2XId'$spqXsjr" Q.O9Df]SQ"B҆I(pT@HjiP/69Kq NgTa+II Fpx56eXrnʻ9hRym)AIeJ8ZܵH +MSc- 8n +݂Ix,>ZJ2Y;3$O=,=p܈ KlBx%;| .XeC0CpIHae HJ ` +JTI62nPOCgNrKP1< - B -DHn3:PnY'h(Z +R)JR)JR)JR)JRjk֡іYɍt,9*IS +R'#?Vj;Dڃd-%(y-JqgԸ{mXq8t`[:Xr9ȨZ:oc)̦֒$SI9J^e6%;v)TZP <5vzaCBׄ[)Km,w;S#֙7vz-(D,ZP|*e)AC*M)\J ++Hf4vVBG'^ċ%Y)kD6QFs8?zi.x[e hdO4)}B\FI$+8rao%x;Aǥ{&<W}k +IpAu)JR)JR)JR)Uڒtm˃oz_e%@+(SJ e(*W٬L=d5싱^}r+3ZS+J[[d%+pH+$nsXz#^x3! +H(P;K3cx!R +F7$)ȟj)jYC1LK +}+aUW%.|;)xX>\~~m<&w2x#Opj^#Z;4Nv[lV)m[ҽ}hS]Bne 5 +;#ֶTfOwMDCa&pRD~(_Cg&J6p}ՠ0;%DPgoZJEz +. +qt,[lǙ~aGްj1g&4[mle oHZgܞy})JR)JR)JR)JRlo 7(#R@6 RBiR}ԙM)βaN :7 sAlGd3;im'j@6TQiR㰷Iq ҽœ[ ~+>ʉ q$8 ʩ]( \`"P)>^qT[ҵq1" JB} [l崀2Wۓk8 rKl$ 8dY)n-2ȭPTZlZf7MOuj;+iƒՔ$)9D%R,z7ki5M:!EqJ<'́sd9.yLyʴ8m!Em%ƟJprv1PSrܟdn\SɌڣ0.w +K8k)J2ThQ\2C12]XB=< +d9.yL\~~m<&w2x#Opj橵uM$:K1[Rsh;[Fw(N2p+ڷO˲ e|ǶeL} ?ԭ`T+?jYlV;'WkuPSg8Ϲ?"H_%J=%29AW?h}NeJ`%ԶX?JWs^t[pyt-Ҭ$01횑kI )8 @L.ċ橹m(~-NTӀrd!A'ex\ȷe(AqNYkI:Oyd]}R8 w 5a6YOzL +xnܐ=s+m4ȩ7e-'%Is[t%; +L%ZK/6*0n#8< )P7΅۔X) !|se>jLgYC0RoHFGqGҹ Ķ[` ` Oy5")P[cIø{rKn J9JPHBOp*m)JY-nrqr(2HlA;B98Jyܕ%hdyFv8 RI*<$ )PcELGmT6-iHHp})ϨJG*-P&YpRa.׀)ۻn H=[lkL6 D2x 6*oݶ]Sq=IMiJRTǹ<5ꔥ)QnVwncӁJTe%C'2]bj+ a$TI$II$)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR)JR \ No newline at end of file diff -Naur lammps-15May08/doc/Eqs/pair_gromacs.tex lammps-16May08/doc/Eqs/pair_gromacs.tex --- lammps-15May08/doc/Eqs/pair_gromacs.tex 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/doc/Eqs/pair_gromacs.tex 2008-05-15 14:44:18.000000000 -0600 @@ -0,0 +1,14 @@ +\documentstyle[12pt]{article} + +\begin{document} + +\begin{eqnarray*} +E_{LJ} & = & 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - + \left(\frac{\sigma}{r}\right)^6 \right] + S_{LJ}(r) + \qquad r < r_c \\ +E_C & = & \frac{C q_i q_j}{\epsilon r} + S_C(r) \qquad r < r_c \\ +S(r) & = & 0 \qquad r < r_1 \\ +S(r) & = & A (r - r_1)^2 + B (r - r_1)^3 \qquad r_1 < r < r_c +\end{eqnarray*} + +\end{document} diff -Naur lammps-15May08/doc/Section_commands.html lammps-16May08/doc/Section_commands.html --- lammps-15May08/doc/Section_commands.html 2008-05-14 15:40:02.000000000 -0600 +++ lammps-16May08/doc/Section_commands.html 2008-05-15 14:48:06.000000000 -0600 @@ -363,9 +363,10 @@ lj/charmm/coul/charmm/implicitlj/charmm/coul/longlj/charmm/coul/long/optlj/class2 lj/class2/coul/cutlj/class2/coul/longlj/cutlj/cut/opt lj/cut/coul/cutlj/cut/coul/debyelj/cut/coul/longlj/cut/coul/long/tip4p -lj/expandlj/smoothlubricatemeam -morsemorse/optresquaredsoft -swtabletersoffyukawa +lj/expandlj/gromacslj/gromacs/coul/gromacslj/smooth +lubricatemeammorsemorse/opt +resquaredsoftswtable +tersoffyukawa

These are pair styles contributed by users, which can be used if diff -Naur lammps-15May08/doc/Section_commands.txt lammps-16May08/doc/Section_commands.txt --- lammps-15May08/doc/Section_commands.txt 2008-05-14 15:40:02.000000000 -0600 +++ lammps-16May08/doc/Section_commands.txt 2008-05-15 14:48:06.000000000 -0600 @@ -512,6 +512,8 @@ "lj/cut/coul/long"_pair_lj.html, "lj/cut/coul/long/tip4p"_pair_lj.html, "lj/expand"_pair_lj_expand.html, +"lj/gromacs"_pair_gromacs.html, +"lj/gromacs/coul/gromacs"_pair_gromacs.html, "lj/smooth"_pair_lj_smooth.html, "lubricate"_pair_lubricate.html, "meam"_pair_meam.html, diff -Naur lammps-15May08/doc/pair_coeff.html lammps-16May08/doc/pair_coeff.html --- lammps-15May08/doc/pair_coeff.html 2008-02-29 18:13:20.000000000 -0700 +++ lammps-16May08/doc/pair_coeff.html 2008-05-15 14:48:06.000000000 -0600 @@ -122,6 +122,8 @@

  • pair_style lj/cut/coul/long - LJ with long-range Coulomb
  • pair_style lj/cut/coul/long/tip4p - LJ with long-range Coulomb for TIP4P water
  • pair_style lj/expand - Lennard-Jones for variable size particles +
  • pair_style lj/gromacs - GROMACS-style Lennard-Jones potential +
  • pair_style lj/gromacs/coul/gromacs - GROMACS-style LJ and Coulombic potential
  • pair_style lj/smooth - smoothed Lennard-Jones potential
  • pair_style lubricate - hydrodynamic lubrication forces
  • pair_style meam - modified embedded atom method (MEAM) diff -Naur lammps-15May08/doc/pair_coeff.txt lammps-16May08/doc/pair_coeff.txt --- lammps-15May08/doc/pair_coeff.txt 2008-02-29 18:13:20.000000000 -0700 +++ lammps-16May08/doc/pair_coeff.txt 2008-05-15 14:48:06.000000000 -0600 @@ -118,6 +118,8 @@ "pair_style lj/cut/coul/long"_pair_lj.html - LJ with long-range Coulomb "pair_style lj/cut/coul/long/tip4p"_pair_lj.html - LJ with long-range Coulomb for TIP4P water "pair_style lj/expand"_pair_lj_expand.html - Lennard-Jones for variable size particles +"pair_style lj/gromacs"_pair_gromacs.html - GROMACS-style Lennard-Jones potential +"pair_style lj/gromacs/coul/gromacs"_pair_gromacs.html - GROMACS-style LJ and Coulombic potential "pair_style lj/smooth"_pair_lj_smooth.html - smoothed Lennard-Jones potential "pair_style lubricate"_pair_lubricate.html - hydrodynamic lubrication forces "pair_style meam"_pair_meam.html - modified embedded atom method (MEAM) diff -Naur lammps-15May08/doc/pair_gromacs.html lammps-16May08/doc/pair_gromacs.html --- lammps-15May08/doc/pair_gromacs.html 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/doc/pair_gromacs.html 2008-05-15 14:48:06.000000000 -0600 @@ -0,0 +1,134 @@ + +
    LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Commands +
    + + + + + + +
    + +

    pair_style lj/gromacs command +

    +

    pair_style lj/gromacs/coul/gromacs command +

    +

    Syntax: +

    +
    pair_style style args 
+
    +
    • style = lj/gromacs or lj/gromacs/coul/gromacs +
    • args = list of arguments for a particular style +
    +
      lj/gromacs args = inner outer
+    inner, outer = global switching cutoffs for Lennard Jones
+  lj/gromacs/coul/gromacs args = inner outer (inner2) (outer2)
+    inner, outer = global switching cutoffs for Lennard Jones (and Coulombic if only 2 args)
+    inner2, outer2 = global switching cutoffs for Coulombic (optional) 
+
    +

    Examples: +

    +
    pair_style lj/gromacs 9.0 12.0
+pair_coeff * * 100.0 2.0
+pair_coeff 2 2 100.0 2.0 8.0 10.0 
+
    +
    pair_style lj/gromacs/coul/gromacs 9.0 12.0
+pair_style lj/gromacs/coul/gromacs 8.0 10.0 7.0 9.0
+pair_coeff * * 100.0 2.0 
+
    +

    Description: +

    +

    The lj/gromacs styles compute LJ and Coulombic interactions with an +additional switching function S(r) that ramps the energy and force +smoothly to zero between an inner and outer cutoff. It is a commonly +used potential in the GROMACS MD code and for +the coarse-grained models of (Marrink). +

    +
    +
    +

    R1 is the inner cutoff; Rc is the outer cutoff. The coefficients A +and B are computed by LAMMPS to perform the smoothing. The function +S(r) is actually applied once to each term of the LJ formula and once +to the Coulombic formula, so there are 2 or 3 sets of A,B coefficients +depending on which pair_style is used. The boundary conditions +applied to the smoothing function are as follows: S(r1) = S'(r1) = 0, +S(rc) = -F(rc), S'(rc) = -F'(rc), where F(r) is the correpsonding term +in the LJ or Coulombic function and a single quote represents a +derivative with respect to r. +

    +

    The inner and outer cutoff for the LJ and Coulombic terms can be the +same or different depending on whether 2 or 4 arguments are used in +the pair_style command. The inner LJ cutoff must be > 0, but the +inner Coulombic cutoff can be >= 0. +

    +

    The following coefficients must be defined for each pair of atoms +types via the pair_coeff command as in the examples +above, or in the data file or restart files read by the +read_data or read_restart +commands, or by mixing as described below: +

    +
    • epsilon (energy units) +
    • sigma (distance units) +
    • inner (distance units) +
    • outer (distance units) +
    +

    Note that sigma is defined in the LJ formula as the zero-crossing +distance for the potential, not as the energy minimum at 2^(1/6) +sigma. +

    +

    The last 2 coefficients are optional inner and outer cutoffs for style +lj/gromacs. If not specified, the global inner and outer values +are used. +

    +

    The last 2 coefficients cannot be used with style +lj/gromacs/coul/gromacs because this force field does not allow +varying cutoffs for individual atom pairs; all pairs use the global +cutoff(s) specified in the pair_style command. +

    +
    + +

    Mixing, shift, table, tail correction, restart, rRESPA info: +

    +

    For atom type pairs I,J and I != J, the epsilon and sigma coefficients +and cutoff distance for all of the lj/cut pair styles can be mixed. +The default mix value is geometric. See the "pair_modify" command +for details. +

    +

    None of the GROMACS pair styles support the +pair_modify shift option, since the Lennard-Jones +portion of the pair interaction is already smoothed to 0.0 at the +cutoff. +

    +

    The pair_modify table option is not relevant +for this pair style. +

    +

    None of the GROMACS pair styles support the +pair_modify tail option for adding long-range tail +corrections to energy and pressure, since there are no corrections for +a potential that goes to 0.0 at the cutoff. +

    +

    All of the GROMACS pair styles write their information to binary +restart files, so pair_style and pair_coeff commands do +not need to be specified in an input script that reads a restart file. +

    +

    All of the GROMACS pair styles can only be used via the pair +keyword of the run_style respa command. They do not +support the inner, middle, outer keywords. +

    +
    + +

    Restrictions: none +

    +

    Related commands: +

    +

    pair_coeff +

    +

    Default: none +

    +
    + + + +

    (Marrink) Marrink, de Vries, Mark, J Phys Chem B, 108, 750-760 (2004). +

    + diff -Naur lammps-15May08/doc/pair_gromacs.txt lammps-16May08/doc/pair_gromacs.txt --- lammps-15May08/doc/pair_gromacs.txt 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/doc/pair_gromacs.txt 2008-05-15 14:48:06.000000000 -0600 @@ -0,0 +1,126 @@ +"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 lj/gromacs command :h3 +pair_style lj/gromacs/coul/gromacs command :h3 + +[Syntax:] + +pair_style style args :pre + +style = {lj/gromacs} or {lj/gromacs/coul/gromacs} +args = list of arguments for a particular style :ul + {lj/gromacs} args = inner outer + inner, outer = global switching cutoffs for Lennard Jones + {lj/gromacs/coul/gromacs} args = inner outer (inner2) (outer2) + inner, outer = global switching cutoffs for Lennard Jones (and Coulombic if only 2 args) + inner2, outer2 = global switching cutoffs for Coulombic (optional) :pre + +[Examples:] + +pair_style lj/gromacs 9.0 12.0 +pair_coeff * * 100.0 2.0 +pair_coeff 2 2 100.0 2.0 8.0 10.0 :pre + +pair_style lj/gromacs/coul/gromacs 9.0 12.0 +pair_style lj/gromacs/coul/gromacs 8.0 10.0 7.0 9.0 +pair_coeff * * 100.0 2.0 :pre + +[Description:] + +The {lj/gromacs} styles compute LJ and Coulombic interactions with an +additional switching function S(r) that ramps the energy and force +smoothly to zero between an inner and outer cutoff. It is a commonly +used potential in the "GROMACS"_http://www.gromacs.org MD code and for +the coarse-grained models of "(Marrink)"_#Marrink. + +:c,image(Eqs/pair_gromacs.jpg) + +R1 is the inner cutoff; Rc is the outer cutoff. The coefficients A +and B are computed by LAMMPS to perform the smoothing. The function +S(r) is actually applied once to each term of the LJ formula and once +to the Coulombic formula, so there are 2 or 3 sets of A,B coefficients +depending on which pair_style is used. The boundary conditions +applied to the smoothing function are as follows: S(r1) = S'(r1) = 0, +S(rc) = -F(rc), S'(rc) = -F'(rc), where F(r) is the correpsonding term +in the LJ or Coulombic function and a single quote represents a +derivative with respect to r. + +The inner and outer cutoff for the LJ and Coulombic terms can be the +same or different depending on whether 2 or 4 arguments are used in +the pair_style command. The inner LJ cutoff must be > 0, but the +inner Coulombic cutoff can be >= 0. + +The following coefficients must be defined for each pair of atoms +types via the "pair_coeff"_pair_coeff.html command as in the examples +above, or in the data file or restart files read by the +"read_data"_read_data.html or "read_restart"_read_restart.html +commands, or by mixing as described below: + +epsilon (energy units) +sigma (distance units) +inner (distance units) +outer (distance units) :ul + +Note that sigma is defined in the LJ formula as the zero-crossing +distance for the potential, not as the energy minimum at 2^(1/6) +sigma. + +The last 2 coefficients are optional inner and outer cutoffs for style +{lj/gromacs}. If not specified, the global {inner} and {outer} values +are used. + +The last 2 coefficients cannot be used with style +{lj/gromacs/coul/gromacs} because this force field does not allow +varying cutoffs for individual atom pairs; all pairs use the global +cutoff(s) specified in the pair_style command. + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +For atom type pairs I,J and I != J, the epsilon and sigma coefficients +and cutoff distance for all of the lj/cut pair styles can be mixed. +The default mix value is {geometric}. See the "pair_modify" command +for details. + +None of the GROMACS pair styles support the +"pair_modify"_pair_modify.html shift option, since the Lennard-Jones +portion of the pair interaction is already smoothed to 0.0 at the +cutoff. + +The "pair_modify"_pair_modify.html table option is not relevant +for this pair style. + +None of the GROMACS pair styles support the +"pair_modify"_pair_modify.html tail option for adding long-range tail +corrections to energy and pressure, since there are no corrections for +a potential that goes to 0.0 at the cutoff. + +All of the GROMACS pair styles write their information to "binary +restart files"_restart.html, so pair_style and pair_coeff commands do +not need to be specified in an input script that reads a restart file. + +All of the GROMACS pair styles can only be used via the {pair} +keyword of the "run_style respa"_run_style.html command. They do not +support the {inner}, {middle}, {outer} keywords. + +:line + +[Restrictions:] none + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] none + +:line + +:link(Marrink) +[(Marrink)] Marrink, de Vries, Mark, J Phys Chem B, 108, 750-760 (2004). diff -Naur lammps-15May08/doc/pair_style.html lammps-16May08/doc/pair_style.html --- lammps-15May08/doc/pair_style.html 2008-02-29 18:13:20.000000000 -0700 +++ lammps-16May08/doc/pair_style.html 2008-05-15 14:48:06.000000000 -0600 @@ -124,6 +124,8 @@
  • pair_style lj/cut/coul/long - LJ with long-range Coulomb
  • pair_style lj/cut/coul/long/tip4p - LJ with long-range Coulomb for TIP4P water
  • pair_style lj/expand - Lennard-Jones for variable size particles +
  • pair_style lj/gromacs - GROMACS-style Lennard-Jones potential +
  • pair_style lj/gromacs/coul/gromacs - GROMACS-style LJ and Coulombic potential
  • pair_style lj/smooth - smoothed Lennard-Jones potential
  • pair_style lubricate - hydrodynamic lubrication forces
  • pair_style meam - modified embedded atom method (MEAM) diff -Naur lammps-15May08/doc/pair_style.txt lammps-16May08/doc/pair_style.txt --- lammps-15May08/doc/pair_style.txt 2008-02-29 18:13:20.000000000 -0700 +++ lammps-16May08/doc/pair_style.txt 2008-05-15 14:48:06.000000000 -0600 @@ -120,6 +120,8 @@ "pair_style lj/cut/coul/long"_pair_lj.html - LJ with long-range Coulomb "pair_style lj/cut/coul/long/tip4p"_pair_lj.html - LJ with long-range Coulomb for TIP4P water "pair_style lj/expand"_pair_lj_expand.html - Lennard-Jones for variable size particles +"pair_style lj/gromacs"_pair_gromacs.html - GROMACS-style Lennard-Jones potential +"pair_style lj/gromacs/coul/gromacs"_pair_gromacs.html - GROMACS-style LJ and Coulombic potential "pair_style lj/smooth"_pair_lj_smooth.html - smoothed Lennard-Jones potential "pair_style lubricate"_pair_lubricate.html - hydrodynamic lubrication forces "pair_style meam"_pair_meam.html - modified embedded atom method (MEAM) diff -Naur lammps-15May08/src/pair_lj_gromacs.cpp lammps-16May08/src/pair_lj_gromacs.cpp --- lammps-15May08/src/pair_lj_gromacs.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/src/pair_lj_gromacs.cpp 2008-05-15 14:49:38.000000000 -0600 @@ -0,0 +1,425 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mark Stevens (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_gromacs.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairLJGromacs::PairLJGromacs(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairLJGromacs::~PairLJGromacs() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + + memory->destroy_2d_double_array(cut); + memory->destroy_2d_double_array(cut_inner); + memory->destroy_2d_double_array(cut_inner_sq); + memory->destroy_2d_double_array(epsilon); + memory->destroy_2d_double_array(sigma); + memory->destroy_2d_double_array(lj1); + memory->destroy_2d_double_array(lj2); + memory->destroy_2d_double_array(lj3); + memory->destroy_2d_double_array(lj4); + memory->destroy_2d_double_array(ljsw1); + memory->destroy_2d_double_array(ljsw2); + memory->destroy_2d_double_array(ljsw3); + memory->destroy_2d_double_array(ljsw4); + memory->destroy_2d_double_array(ljsw5); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJGromacs::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,r2inv,r6inv,forcelj,factor_lj; + double r,t,fswitch,eswitch; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + if (j < nall) factor_lj = 1.0; + else { + factor_lj = special_lj[j/nall]; + j %= nall; + } + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_inner_sq[itype][jtype]) { + r = sqrt(rsq); + t = r - cut_inner[itype][jtype]; + fswitch = r*t*t*(ljsw1[itype][jtype] + ljsw2[itype][jtype]*t); + forcelj += fswitch; + } + + fpair = factor_lj*forcelj*r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + + ljsw5[itype][jtype]; + if (rsq > cut_inner_sq[itype][jtype]) { + eswitch = t*t*t*(ljsw3[itype][jtype] + ljsw4[itype][jtype]*t); + evdwl += eswitch; + } + evdwl *= factor_lj; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJGromacs::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + cut = memory->create_2d_double_array(n+1,n+1,"pair:cut"); + cut_inner = memory->create_2d_double_array(n+1,n+1,"pair:cut_inner"); + cut_inner_sq = memory->create_2d_double_array(n+1,n+1,"pair:cut_inner_sq"); + epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon"); + sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); + lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1"); + lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2"); + lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3"); + lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4"); + ljsw1 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw1"); + ljsw2 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw2"); + ljsw3 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw3"); + ljsw4 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw4"); + ljsw5 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw5"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJGromacs::settings(int narg, char **arg) +{ + if (narg != 2) error->all("Illegal pair_style command"); + + cut_inner_global = atof(arg[0]); + cut_global = atof(arg[1]); + + if (cut_inner_global <= 0.0 || cut_inner_global > cut_global) + error->all("Illegal pair_style command"); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) { + cut_inner[i][j] = cut_inner_global; + cut[i][j] = cut_global; + } + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJGromacs::coeff(int narg, char **arg) +{ + if (narg != 4 && narg != 6) + error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = atof(arg[2]); + double sigma_one = atof(arg[3]); + + double cut_inner_one = cut_inner_global; + double cut_one = cut_global; + if (narg == 6) { + cut_inner_one = atof(arg[4]); + cut_one = atof(arg[5]); + } + + if (cut_inner_one <= 0.0 || cut_inner_one > cut_one) + error->all("Incorrect args for pair coefficients"); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + epsilon[i][j] = epsilon_one; + sigma[i][j] = sigma_one; + cut_inner[i][j] = cut_inner_one; + cut[i][j] = cut_one; + 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 PairLJGromacs::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j]; + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + + double r6inv = 1.0/pow(cut[i][j],6.0); + double r8inv = 1.0/pow(cut[i][j],8.0); + double t = cut[i][j] - cut_inner[i][j]; + double t2inv = 1.0/(t*t); + double t3inv = t2inv/t; + double t3 = 1.0/t3inv; + double a6 = ( 7.0*cut_inner[i][j] - 10.0*cut[i][j])*r8inv*t2inv; + double b6 = ( 9.0*cut[i][j] - 7.0*cut_inner[i][j])*r8inv*t3inv; + double a12 = (13.0*cut_inner[i][j] - 16.0*cut[i][j])*r6inv*r8inv*t2inv; + double b12 = (15.0*cut[i][j] - 13.0*cut_inner[i][j])*r6inv*r8inv*t3inv; + double c6 = r6inv - t3*(a6/3.0 + b6*t/4.0); + double c12 = r6inv*r6inv - t3*(a12/3.0 + b12*t/4.0); + ljsw1[i][j] = lj1[i][j]*a12 - lj2[i][j]*a6; + ljsw2[i][j] = lj1[i][j]*b12 - lj2[i][j]*b6; + ljsw3[i][j] =-lj3[i][j]*a12/3.0 + lj4[i][j]*a6/3.0; + ljsw4[i][j] =-lj3[i][j]*b12/4.0 + lj4[i][j]*b6/4.0; + ljsw5[i][j] =-lj3[i][j]*c12 + lj4[i][j]*c6; + + cut_inner[j][i] = cut_inner[i][j]; + cut_inner_sq[j][i] = cut_inner_sq[i][j]; + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + ljsw1[j][i] = ljsw1[i][j]; + ljsw2[j][i] = ljsw2[i][j]; + ljsw3[j][i] = ljsw3[i][j]; + ljsw4[j][i] = ljsw4[i][j]; + ljsw5[j][i] = ljsw5[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJGromacs::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&epsilon[i][j],sizeof(double),1,fp); + fwrite(&sigma[i][j],sizeof(double),1,fp); + fwrite(&cut_inner[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJGromacs::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&epsilon[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + fread(&cut_inner[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJGromacs::write_restart_settings(FILE *fp) +{ + fwrite(&cut_inner_global,sizeof(double),1,fp); + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJGromacs::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&cut_inner_global,sizeof(double),1,fp); + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJGromacs::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,forcelj,philj; + double r,t,fswitch,phiswitch; + + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_inner_sq[itype][jtype]) { + r = sqrt(rsq); + t = r - cut_inner[itype][jtype]; + fswitch = r*t*t*(ljsw1[itype][jtype] + ljsw2[itype][jtype]*t); + forcelj += fswitch; + } + fforce = factor_lj*forcelj*r2inv; + + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + + ljsw5[itype][jtype]; + if (rsq > cut_inner_sq[itype][jtype]) { + phiswitch = t*t*t*(ljsw3[itype][jtype] + ljsw4[itype][jtype]*t); + philj += phiswitch; + } + + return factor_lj*philj; +} diff -Naur lammps-15May08/src/pair_lj_gromacs.h lammps-16May08/src/pair_lj_gromacs.h --- lammps-15May08/src/pair_lj_gromacs.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/src/pair_lj_gromacs.h 2008-05-15 14:49:38.000000000 -0600 @@ -0,0 +1,47 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef PAIR_LJ_GROMACS_H +#define PAIR_LJ_GROMACS_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJGromacs : public Pair { + public: + PairLJGromacs(class LAMMPS *); + virtual ~PairLJGromacs(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double cut_inner_global,cut_global; + double **cut,**cut_inner,**cut_inner_sq; + double **epsilon,**sigma; + double **lj1,**lj2,**lj3,**lj4; + double **ljsw1,**ljsw2,**ljsw3,**ljsw4,**ljsw5; + + void allocate(); +}; + +} + +#endif diff -Naur lammps-15May08/src/pair_lj_gromacs_coul_gromacs.cpp lammps-16May08/src/pair_lj_gromacs_coul_gromacs.cpp --- lammps-15May08/src/pair_lj_gromacs_coul_gromacs.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/src/pair_lj_gromacs_coul_gromacs.cpp 2008-05-15 14:51:14.000000000 -0600 @@ -0,0 +1,491 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mark Stevens (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_lj_gromacs_coul_gromacs.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairLJGromacsCoulGromacs::PairLJGromacsCoulGromacs(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairLJGromacsCoulGromacs::~PairLJGromacsCoulGromacs() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + + memory->destroy_2d_double_array(epsilon); + memory->destroy_2d_double_array(sigma); + memory->destroy_2d_double_array(lj1); + memory->destroy_2d_double_array(lj2); + memory->destroy_2d_double_array(lj3); + memory->destroy_2d_double_array(lj4); + memory->destroy_2d_double_array(ljsw1); + memory->destroy_2d_double_array(ljsw2); + memory->destroy_2d_double_array(ljsw3); + memory->destroy_2d_double_array(ljsw4); + memory->destroy_2d_double_array(ljsw5); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; + double r,tlj,tc,fswitch,fswitchcoul,eswitch,ecoulswitch; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + if (j < nall) factor_coul = factor_lj = 1.0; + else { + factor_coul = special_coul[j/nall]; + factor_lj = special_lj[j/nall]; + j %= nall; + } + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq < cut_bothsq) { + r2inv = 1.0/rsq; + + // skip if qi or qj = 0.0 since this potential may be used as + // coarse-grain model with many uncharged atoms + + if (rsq < cut_coulsq && qtmp != 0.0 && q[j] != 0.0) { + forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); + if (rsq > cut_coul_innersq) { + r = sqrt(rsq); + tc = r - cut_coul_inner; + fswitchcoul = qqrd2e * qtmp*q[j]*r*tc*tc*(coulsw1 + coulsw2*tc); + forcecoul += fswitchcoul; + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq) { + r6inv = r2inv*r2inv*r2inv; + jtype = type[j]; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + r = sqrt(rsq); + tlj = r - cut_lj_inner; + fswitch = r*tlj*tlj*(ljsw1[itype][jtype] + + ljsw2[itype][jtype]*tlj); + forcelj += fswitch; + } + } else forcelj = 0.0; + + fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + if (rsq < cut_coulsq) { + ecoul = qqrd2e * qtmp*q[j] * (sqrt(r2inv) - coulsw5); + if (rsq > cut_coul_innersq) { + ecoulswitch = tc*tc*tc * (coulsw3 + coulsw4*tc); + ecoul += qqrd2e*qtmp*q[j]*ecoulswitch; + } + ecoul *= factor_coul; + } else ecoul = 0.0; + if (rsq < cut_ljsq) { + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + + ljsw5[itype][jtype]; + if (rsq > cut_lj_innersq) { + eswitch = tlj*tlj*tlj * (ljsw3[itype][jtype] + + ljsw4[itype][jtype]*tlj); + evdwl += eswitch; + } + evdwl *= factor_lj; + } else evdwl = 0.0; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); + + epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon"); + sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma"); + lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1"); + lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2"); + lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3"); + lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4"); + ljsw1 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw1"); + ljsw2 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw2"); + ljsw3 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw3"); + ljsw4 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw4"); + ljsw5 = memory->create_2d_double_array(n+1,n+1,"pair:ljsw5"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::settings(int narg, char **arg) +{ + if (narg != 2 && narg != 4) + error->all("Illegal pair_style command"); + + cut_lj_inner = atof(arg[0]); + cut_lj = atof(arg[1]); + if (narg == 2) { + cut_coul_inner = cut_lj_inner; + cut_coul = cut_lj; + } else { + cut_coul_inner = atof(arg[2]); + cut_coul = atof(arg[3]); + } + + if (cut_lj_inner <= 0.0 || cut_coul_inner < 0.0) + error->all("Illegal pair_style command"); + if (cut_lj_inner > cut_lj || cut_coul_inner > cut_coul) + error->all("Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::coeff(int narg, char **arg) +{ + if (narg != 4) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = atof(arg[2]); + double sigma_one = atof(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + epsilon[i][j] = epsilon_one; + sigma[i][j] = sigma_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::init_style() +{ + if (!atom->q_flag) + error->all("Pair style lj/gromacs/coul/gromacs requires atom attribute q"); + + int irequest = neighbor->request(this); + + cut_lj_innersq = cut_lj_inner * cut_lj_inner; + cut_ljsq = cut_lj * cut_lj; + cut_coul_innersq = cut_coul_inner * cut_coul_inner; + cut_coulsq = cut_coul * cut_coul; + cut_bothsq = MAX(cut_ljsq,cut_coulsq); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJGromacsCoulGromacs::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], + sigma[i][i],sigma[j][j]); + sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); + } + + double cut = MAX(cut_lj,cut_coul); + + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + + double r6inv = 1.0/pow(cut_lj,6.0); + double r8inv = 1.0/pow(cut_lj,8.0); + double t = cut_lj - cut_lj_inner; + double t2inv = 1.0/(t*t); + double t3inv = t2inv/t; + double t3 = 1.0/t3inv; + double a6 = (7.0*cut_lj_inner - 10.0*cut_lj)*r8inv*t2inv; + double b6 = (9.0*cut_lj - 7.0*cut_lj_inner)*r8inv*t3inv; + double a12 = (13.0*cut_lj_inner - 16.0*cut_lj)*r6inv*r8inv*t2inv; + double b12 = (15.0*cut_lj - 13.0*cut_lj_inner)*r6inv*r8inv*t3inv; + double c6 = r6inv - t3*(a6/3.0 + b6*t/4.0); + double c12 = r6inv*r6inv - t3*(a12/3.0 + b12*t/4.0); + ljsw1[i][j] = lj1[i][j]*a12 - lj2[i][j]*a6; + ljsw2[i][j] = lj1[i][j]*b12 - lj2[i][j]*b6; + ljsw3[i][j] = -lj3[i][j]*a12/3.0 + lj4[i][j]*a6/3.0; + ljsw4[i][j] = -lj3[i][j]*b12/4.0 + lj4[i][j]*b6/4.0; + ljsw5[i][j] = -lj3[i][j]*c12 + lj4[i][j]*c6; + + double r3inv = 1.0/pow(cut_coul,3.0); + t = cut_coul - cut_coul_inner; + t2inv = 1.0/(t*t); + t3inv = t2inv/t; + double a1 = (2.0*cut_coul_inner - 5.0*cut_coul) * r3inv*t2inv; + double b1 = (4.0*cut_coul - 2.0*cut_coul_inner) * r3inv*t3inv; + coulsw1 = a1; + coulsw2 = b1; + coulsw3 = -a1/3.0; + coulsw4 = -b1/4.0; + coulsw5 = 1.0/cut_coul - t*t*t*(a1/3.0 + b1*t/4.0); + + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + ljsw1[j][i] = ljsw1[i][j]; + ljsw2[j][i] = ljsw2[i][j]; + ljsw3[j][i] = ljsw3[i][j]; + ljsw4[j][i] = ljsw4[i][j]; + ljsw5[j][i] = ljsw5[i][j]; + + return cut; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&epsilon[i][j],sizeof(double),1,fp); + fwrite(&sigma[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&epsilon[i][j],sizeof(double),1,fp); + fread(&sigma[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::write_restart_settings(FILE *fp) +{ + fwrite(&cut_lj_inner,sizeof(double),1,fp); + fwrite(&cut_lj,sizeof(double),1,fp); + fwrite(&cut_coul_inner,sizeof(double),1,fp); + fwrite(&cut_coul,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJGromacsCoulGromacs::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_lj_inner,sizeof(double),1,fp); + fread(&cut_lj,sizeof(double),1,fp); + fread(&cut_coul_inner,sizeof(double),1,fp); + fread(&cut_coul,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_lj_inner,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul_inner,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJGromacsCoulGromacs::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,forcecoul,forcelj,phicoul,philj; + double r,tlj,tc,fswitch,phiswitch,fswitchcoul,phiswitchcoul; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv); + if (rsq > cut_coul_innersq) { + r = sqrt(rsq); + tc = r - cut_coul_inner; + fswitchcoul = force->qqrd2e * + atom->q[i]*atom->q[j] * r*tc*tc * (coulsw1 + coulsw2*tc); + forcecoul += fswitchcoul; + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq) { + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); + if (rsq > cut_lj_innersq) { + r = sqrt(rsq); + tlj = r - cut_lj_inner; + fswitch = r*tlj*tlj*(ljsw1[itype][jtype] + ljsw2[itype][jtype]*tlj); + forcelj += fswitch; + } + } else forcelj = 0.0; + + fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv; + + double eng = 0.0; + if (rsq < cut_coulsq) { + phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * (sqrt(r2inv)-coulsw5); + if (rsq > cut_coul_innersq) { + phiswitchcoul = force->qqrd2e * atom->q[i]*atom->q[j] * + tc*tc*tc * (coulsw3 + coulsw4*tc); + phicoul += phiswitchcoul; + } + eng += factor_coul*phicoul; + } + + if (rsq < cut_ljsq) { + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + + ljsw5[itype][jtype]; + if (rsq > cut_lj_innersq) { + phiswitch = tlj*tlj*tlj * (ljsw3[itype][jtype] + + ljsw4[itype][jtype]*tlj); + philj += phiswitch; + } + eng += factor_lj*philj; + } + + return eng; +} diff -Naur lammps-15May08/src/pair_lj_gromacs_coul_gromacs.h lammps-16May08/src/pair_lj_gromacs_coul_gromacs.h --- lammps-15May08/src/pair_lj_gromacs_coul_gromacs.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-16May08/src/pair_lj_gromacs_coul_gromacs.h 2008-05-15 14:49:38.000000000 -0600 @@ -0,0 +1,49 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef PAIR_LJ_GROMACS_COUL_GROMACS_H +#define PAIR_LJ_GROMACS_COUL_GROMACS_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJGromacsCoulGromacs : public Pair { + public: + PairLJGromacsCoulGromacs(class LAMMPS *); + virtual ~PairLJGromacsCoulGromacs(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double cut_lj_inner,cut_lj,cut_coul_inner,cut_coul; + double cut_lj_innersq,cut_ljsq,cut_coul_innersq,cut_coulsq,cut_bothsq; + double **epsilon,**sigma; + double **lj1,**lj2,**lj3,**lj4; + double **ljsw1,**ljsw2,**ljsw3,**ljsw4,**ljsw5; + double coulsw1,coulsw2,coulsw3,coulsw4,coulsw5; + + void allocate(); +}; + +} + +#endif diff -Naur lammps-15May08/src/style.h lammps-16May08/src/style.h --- lammps-15May08/src/style.h 2008-05-14 15:38:34.000000000 -0600 +++ lammps-16May08/src/style.h 2008-05-15 14:49:38.000000000 -0600 @@ -300,6 +300,8 @@ #include "pair_lj_cut_coul_cut.h" #include "pair_lj_cut_coul_debye.h" #include "pair_lj_expand.h" +#include "pair_lj_gromacs.h" +#include "pair_lj_gromacs_coul_gromacs.h" #include "pair_lj_smooth.h" #include "pair_morse.h" #include "pair_soft.h" @@ -318,6 +320,8 @@ PairStyle(lj/cut/coul/cut,PairLJCutCoulCut) PairStyle(lj/cut/coul/debye,PairLJCutCoulDebye) PairStyle(lj/expand,PairLJExpand) +PairStyle(lj/gromacs,PairLJGromacs) +PairStyle(lj/gromacs/coul/gromacs,PairLJGromacsCoulGromacs) PairStyle(lj/smooth,PairLJSmooth) PairStyle(morse,PairMorse) PairStyle(soft,PairSoft)