diff --git a/doc/src/Commands_bond.txt b/doc/src/Commands_bond.txt index fbf292aab2..3b6d708863 100644 --- a/doc/src/Commands_bond.txt +++ b/doc/src/Commands_bond.txt @@ -37,6 +37,7 @@ OPT. "harmonic (iko)"_bond_harmonic.html, "harmonic/shift (o)"_bond_harmonic_shift.html, "harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html, +"mm3"_bond_mm3.html, "morse (o)"_bond_morse.html, "nonlinear (o)"_bond_nonlinear.html, "oxdna/fene"_bond_oxdna.html, @@ -67,10 +68,12 @@ OPT. "cosine/shift (o)"_angle_cosine_shift.html, "cosine/shift/exp (o)"_angle_cosine_shift_exp.html, "cosine/squared (o)"_angle_cosine_squared.html, +"cross"_angle_cross.html, "dipole (o)"_angle_dipole.html, "fourier (o)"_angle_fourier.html, "fourier/simple (o)"_angle_fourier_simple.html, "harmonic (iko)"_angle_harmonic.html, +"mm3"_angle_mm3.html, "quartic (o)"_angle_quartic.html, "sdk (o)"_angle_sdk.html, "table (o)"_angle_table.html :tb(c=4,ea=c) @@ -120,8 +123,10 @@ OPT. "cossq (o)"_improper_cossq.html, "cvff (io)"_improper_cvff.html, "distance"_improper_distance.html, +"distharm"_improper_distharm.html, "fourier (o)"_improper_fourier.html, "harmonic (iko)"_improper_harmonic.html, "inversion/harmonic"_improper_inversion_harmonic.html, "ring (o)"_improper_ring.html, -"umbrella (o)"_improper_umbrella.html :tb(c=4,ea=c) +"umbrella (o)"_improper_umbrella.html, +"sqdistharm"_improper_sqdistharm.html :tb(c=4,ea=c) diff --git a/doc/src/Commands_pair.txt b/doc/src/Commands_pair.txt index 46f4950ad7..b498e6bfda 100644 --- a/doc/src/Commands_pair.txt +++ b/doc/src/Commands_pair.txt @@ -154,6 +154,7 @@ OPT. "lj/sf/dipole/sf (go)"_pair_dipole.html, "lj/smooth (o)"_pair_lj_smooth.html, "lj/smooth/linear (o)"_pair_lj_smooth_linear.html, +"lj/switch3/coulgauss/long"_pair_lj_switch3_coulgauss.html, "lj96/cut (go)"_pair_lj96.html, "lubricate (o)"_pair_lubricate.html, "lubricate/poly (o)"_pair_lubricate.html, diff --git a/doc/src/Eqs/angle_cross.jpg b/doc/src/Eqs/angle_cross.jpg new file mode 100644 index 0000000000..3a5522382b Binary files /dev/null and b/doc/src/Eqs/angle_cross.jpg differ diff --git a/doc/src/Eqs/angle_cross.tex b/doc/src/Eqs/angle_cross.tex new file mode 100644 index 0000000000..9d1fdcb7f8 --- /dev/null +++ b/doc/src/Eqs/angle_cross.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} +\thispagestyle{empty} +$$ + E = K_{SS} \left(r_{12}-r_{12,0}\right)\left(r_{32}-r_{32,0}\right) + K_{BS0}\left(r_{12}-r_{12,0}\right)\left(\theta-\theta_0\right) + K_{BS1}\left(r_{32}-r_{32,0}\right)\left(\theta-\theta_0\right) +$$ + +\end{document} diff --git a/doc/src/Eqs/angle_mm3.jpg b/doc/src/Eqs/angle_mm3.jpg new file mode 100644 index 0000000000..ef234c3ed0 Binary files /dev/null and b/doc/src/Eqs/angle_mm3.jpg differ diff --git a/doc/src/Eqs/angle_mm3.tex b/doc/src/Eqs/angle_mm3.tex new file mode 100644 index 0000000000..e2d96f2d69 --- /dev/null +++ b/doc/src/Eqs/angle_mm3.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} +\thispagestyle{empty} +$$ + E = K (\theta - \theta_0)^2 \left[ 1 - 0.014(\theta - \theta_0) + 5.6(10)^{-5} (\theta - \theta_0)^2 - 7.0(10)^{-7} (\theta - \theta_0)^3 + 9(10)^{-10} (\theta - \theta_0)^4 \right] +$$ + +\end{document} diff --git a/doc/src/Eqs/bond_mm3.jpg b/doc/src/Eqs/bond_mm3.jpg new file mode 100644 index 0000000000..b8bbbbf7b0 Binary files /dev/null and b/doc/src/Eqs/bond_mm3.jpg differ diff --git a/doc/src/Eqs/bond_mm3.tex b/doc/src/Eqs/bond_mm3.tex new file mode 100644 index 0000000000..549500ebac --- /dev/null +++ b/doc/src/Eqs/bond_mm3.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} +\thispagestyle{empty} +$$ + E = K (r - r_0)^2 \left[ 1 - 2.55(r-r_0) + (7/12) 2.55^2(r-r_0)^2 \right] +$$ + +\end{document} diff --git a/doc/src/Eqs/improper_distharm.jpg b/doc/src/Eqs/improper_distharm.jpg new file mode 100644 index 0000000000..cef334351d Binary files /dev/null and b/doc/src/Eqs/improper_distharm.jpg differ diff --git a/doc/src/Eqs/improper_distharm.tex b/doc/src/Eqs/improper_distharm.tex new file mode 100644 index 0000000000..218de6a0dd --- /dev/null +++ b/doc/src/Eqs/improper_distharm.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} +\thispagestyle{empty} +$$ + E = K (d - d_0)^2 +$$ + +\end{document} diff --git a/doc/src/Eqs/improper_sqdistharm.jpg b/doc/src/Eqs/improper_sqdistharm.jpg new file mode 100644 index 0000000000..6f5cc97f43 Binary files /dev/null and b/doc/src/Eqs/improper_sqdistharm.jpg differ diff --git a/doc/src/Eqs/improper_sqdistharm.tex b/doc/src/Eqs/improper_sqdistharm.tex new file mode 100644 index 0000000000..1b50a309a0 --- /dev/null +++ b/doc/src/Eqs/improper_sqdistharm.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} +\thispagestyle{empty} +$$ + E = K (d^2 - d_0^2)^2 +$$ + +\end{document} diff --git a/doc/src/Eqs/pair_coulgauss.jpg b/doc/src/Eqs/pair_coulgauss.jpg new file mode 100644 index 0000000000..c4eb346043 Binary files /dev/null and b/doc/src/Eqs/pair_coulgauss.jpg differ diff --git a/doc/src/Eqs/pair_coulgauss.tex b/doc/src/Eqs/pair_coulgauss.tex new file mode 100644 index 0000000000..216d3b3360 --- /dev/null +++ b/doc/src/Eqs/pair_coulgauss.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} + \thispagestyle{empty} +\begin{eqnarray*} + E &=& \frac{q_i q_j \mathrm{erf}\left( r/\sqrt{\gamma_1^2+\gamma_2^2} \right) }{\epsilon r_{ij}} +\end{eqnarray*} + +\end{document} diff --git a/doc/src/Eqs/pair_lj_switch3.jpg b/doc/src/Eqs/pair_lj_switch3.jpg new file mode 100644 index 0000000000..1cc04de6a8 Binary files /dev/null and b/doc/src/Eqs/pair_lj_switch3.jpg differ diff --git a/doc/src/Eqs/pair_lj_switch3.tex b/doc/src/Eqs/pair_lj_switch3.tex new file mode 100644 index 0000000000..29161fb2ef --- /dev/null +++ b/doc/src/Eqs/pair_lj_switch3.tex @@ -0,0 +1,11 @@ +\documentclass[12pt]{article} + +\begin{document} + \thispagestyle{empty} + +\begin{eqnarray*} + E = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6} \right] +% \qquad r < r_c \\ +\end{eqnarray*} + +\end{document} diff --git a/doc/src/Eqs/pair_mm3_switch3.jpg b/doc/src/Eqs/pair_mm3_switch3.jpg new file mode 100644 index 0000000000..caad9bbe99 Binary files /dev/null and b/doc/src/Eqs/pair_mm3_switch3.jpg differ diff --git a/doc/src/Eqs/pair_mm3_switch3.tex b/doc/src/Eqs/pair_mm3_switch3.tex new file mode 100644 index 0000000000..a80278cce2 --- /dev/null +++ b/doc/src/Eqs/pair_mm3_switch3.tex @@ -0,0 +1,11 @@ +\documentclass[12pt]{article} + +\begin{document} + \thispagestyle{empty} +\begin{eqnarray*} + E &=& \epsilon_{ij} \left[ -2.25 \left(\frac{r_{v,ij}}{r_{ij}}\right)^6 + 1.84(10)^5 \exp\left[-12.0 r_{ij}/r_{v,ij}\right] \right] S_3(r_{ij}) \\ + r_{v,ij} &=& r_{v,i} + r_{v,j} \\ + \epsilon_{ij} &=& \sqrt{\epsilon_i \epsilon_j} +\end{eqnarray*} + +\end{document} diff --git a/doc/src/Eqs/pair_switch3.jpg b/doc/src/Eqs/pair_switch3.jpg new file mode 100644 index 0000000000..1b0528c080 Binary files /dev/null and b/doc/src/Eqs/pair_switch3.jpg differ diff --git a/doc/src/Eqs/pair_switch3.tex b/doc/src/Eqs/pair_switch3.tex new file mode 100644 index 0000000000..7ae67a5e65 --- /dev/null +++ b/doc/src/Eqs/pair_switch3.tex @@ -0,0 +1,14 @@ +\documentclass[12pt]{article} + +\begin{document} +\thispagestyle{empty} + +\begin{eqnarray*} + S_3(r) = \left\lbrace \begin{array}{ll} + 1 & \quad\mathrm{if}\quad r < r_\mathrm{c} - w \\ + 3x^2 - 2x^3 & \quad\mathrm{if}\quad r < r_\mathrm{c} \quad\mathrm{with\quad} x=\frac{r_\mathrm{c} - r}{w} \\ + 0 & \quad\mathrm{if}\quad r >= r_\mathrm{c} + \end{array} \right. +\end{eqnarray*} + +\end{document} diff --git a/doc/src/Packages_details.txt b/doc/src/Packages_details.txt index 4768c50ca4..ebfe01d8a7 100644 --- a/doc/src/Packages_details.txt +++ b/doc/src/Packages_details.txt @@ -100,7 +100,8 @@ as contained in the file name. "USER-SPH"_#PKG-USER-SPH, "USER-TALLY"_#PKG-USER-TALLY, "USER-UEF"_#PKG-USER-UEF, -"USER-VTK"_#PKG-USER-VTK :tb(c=6,ea=c) +"USER-VTK"_#PKG-USER-VTK, +"USER-YAFF"_#PKG-USER-YAFF, :tb(c=6,ea=c) :line @@ -2067,3 +2068,39 @@ lib/vtk/README "dump vtk"_dump_vtk.html :ul +:line + +USER-YAFF package :link(PKG-USER-YAFF),h4 + +[Contents:] + +Some potentials that are also implemented in the Yet Another Force Field ("YAFF"_yaff) code. +The expressions and their use are discussed in the following papers + +Vanduyfhuys et al., J. Comput. Chem., 36 (13), 1015-1027 (2015) "link"_vanduyfhuys2015 +Vanduyfhuys et al., J. Comput. Chem., 39 (16), 999-1011 (2018) "link"_vanduyfhuys2018 :ul + +which discuss the "QuickFF"_quickff methodology. + + +:link(vanduyfhuys2015,http://dx.doi.org/10.1002/jcc.23877) +:link(vanduyfhuys2018,http://dx.doi.org/10.1002/jcc.25173) +:link(quickff,http://molmod.github.io/QuickFF) +:link(yaff,https://github.com/molmod/yaff) + + +[Author:] Steven Vandenbrande. + +[Supporting info:] + +src/USER-YAFF/README +"angle_style cross"_angle_cross.html +"angle_style mm3"_angle_mm3.html +"bond_style mm3"_bond_mm3.html +"improper_style distharm"_improper_distharm.html +"improper_style sqdistharm"_improper_sqdistharm.html +"pair_style mm3/switch3/coulgauss/long"_pair_mm3_switch3_coulgauss.html +"pair_style lj/switch3/coulgauss/long"_pair_lj_switch3_coulgauss.html +examples/USER/yaff :ul + + diff --git a/doc/src/Packages_user.txt b/doc/src/Packages_user.txt index cea3db4377..4a702e971f 100644 --- a/doc/src/Packages_user.txt +++ b/doc/src/Packages_user.txt @@ -75,7 +75,8 @@ Package, Description, Doc page, Example, Library "USER-SPH"_Packages_details.html#PKG-USER-SPH, smoothed particle hydrodynamics,"SPH User Guide"_PDF/SPH_LAMMPS_userguide.pdf, USER/sph, no "USER-TALLY"_Packages_details.html#PKG-USER-TALLY, pairwise tally computes,"compute XXX/tally"_compute_tally.html, USER/tally, no "USER-UEF"_Packages_details.html#PKG-USER-UEF, extensional flow,"fix nvt/uef"_fix_nh_uef.html, USER/uef, no -"USER-VTK"_Packages_details.html#PKG-USER-VTK, dump output via VTK, "compute vtk"_dump_vtk.html, n/a, ext :tb(ea=c,ca1=l) +"USER-VTK"_Packages_details.html#PKG-USER-VTK, dump output via VTK, "compute vtk"_dump_vtk.html, n/a, ext +"USER-YAFF"_Packages_details.html#PKG-USER-YAFF, additional styles implemented in YAFF, "angle_style cross"_angle_cross.html, USER/yaff, no :tb(ea=c,ca1=l) :link(MOFplus,https://www.mofplus.org/content/show/MOF-FF) :link(PLUMED,http://www.plumed.org) diff --git a/doc/src/angle_cross.txt b/doc/src/angle_cross.txt new file mode 100644 index 0000000000..d9d83ed4b6 --- /dev/null +++ b/doc/src/angle_cross.txt @@ -0,0 +1,62 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +angle_style cross command :h3 + +[Syntax:] + +angle_style cross :pre + +[Examples:] + +angle_style cross +angle_coeff 1 200.0 100.0 100.0 1.25 1.25 107.0 :pre + +[Description:] + +The {cross} angle style uses a potential that couples the bond stretches of +a bend with the angle stretch of that bend: + +:c,image(Eqs/angle_cross.jpg) + +where r12,0 is the rest value of the bond length between atom 1 and 2, +r32,0 is the rest value of the bond length between atom 2 and 2, +and theta0 is the rest value of the angle. KSS is the force constant of +the bond stretch-bond stretch term and KBS0 and KBS1 are the force constants +of the bond stretch-angle stretch terms. + +The following coefficients must be defined for each angle type via the +"angle_coeff"_angle_coeff.html command as in the example above, or in +the data file or restart files read by the "read_data"_read_data.html +or "read_restart"_read_restart.html commands: + +KSS (energy/distance^2) +KBS0 (energy/distance/rad) +KBS1 (energy/distance/rad) +r12,0 (distance) +r32,0 (distance) +theta0 (degrees) :ul + +Theta0 is specified in degrees, but LAMMPS converts it to radians +internally; hence the units of KBS0 and KBS1 are in energy/distance/radian. + +[Restrictions:] + +This angle style can only be used if LAMMPS was built with the +USER_YAFF package. See the "Build package"_Build_package.html doc +page for more info. + +[Related commands:] + +"angle_coeff"_angle_coeff.html + +[Default:] none + +:line + + diff --git a/doc/src/angle_mm3.txt b/doc/src/angle_mm3.txt new file mode 100644 index 0000000000..9ae032c4ff --- /dev/null +++ b/doc/src/angle_mm3.txt @@ -0,0 +1,55 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +angle_style mm3 command :h3 + +[Syntax:] + +angle_style mm3 :pre + +[Examples:] + +angle_style mm3 +angle_coeff 1 100.0 107.0 :pre + +[Description:] + +The {mm3} angle style uses the potential that is anharmonic in the angle +as defined in "(Allinger)"_#mm3-allinger1989 + +:c,image(Eqs/angle_mm3.jpg) + +where theta0 is the equilibrium value of the angle, and K is a +prefactor. The anharmonic prefactors have units deg^(-n), for example +-0.014 deg^(-1), 5.6(10)^(-5) deg^(-2), ... + +The following coefficients must be defined for each angle type via the +"angle_coeff"_angle_coeff.html command as in the example above, or in +the data file or restart files read by the "read_data"_read_data.html +or "read_restart"_read_restart.html commands: + +K (energy/radian^2) +theta0 (degrees) :ul + +Theta0 is specified in degrees, but LAMMPS converts it to radians +internally; hence the units of K are in energy/radian^2. + +[Restrictions:] + +This angle style can only be used if LAMMPS was built with the +USER_YAFF package. See the "Build package"_Build_package.html doc +page for more info. + +[Related commands:] + +"angle_coeff"_angle_coeff.html + +[Default:] none + +:line + diff --git a/doc/src/angle_style.txt b/doc/src/angle_style.txt index 2c3c1f2bd1..2f2da678d8 100644 --- a/doc/src/angle_style.txt +++ b/doc/src/angle_style.txt @@ -81,10 +81,12 @@ of (g,i,k,o,t) to indicate which accelerated styles exist. "cosine/shift"_angle_cosine_shift.html - angle cosine with a shift "cosine/shift/exp"_angle_cosine_shift_exp.html - cosine with shift and exponential term in spring constant "cosine/squared"_angle_cosine_squared.html - angle with cosine squared term +"cross"_angle_cross.html - cross term coupling angle and bond lengths "dipole"_angle_dipole.html - angle that controls orientation of a point dipole "fourier"_angle_fourier.html - angle with multiple cosine terms "fourier/simple"_angle_fourier_simple.html - angle with a single cosine term "harmonic"_angle_harmonic.html - harmonic angle +"mm3"_angle_mm3.html - anharmonic angle "quartic"_angle_quartic.html - angle with cubic and quartic terms "sdk"_angle_sdk.html - harmonic angle with repulsive SDK pair style between 1-3 atoms "table"_angle_table.html - tabulated by angle :ul diff --git a/doc/src/angles.txt b/doc/src/angles.txt index 3e91ba9917..3d8a47b2eb 100644 --- a/doc/src/angles.txt +++ b/doc/src/angles.txt @@ -14,11 +14,13 @@ Angle Styles :h1 angle_cosine_shift angle_cosine_shift_exp angle_cosine_squared + angle_cross angle_dipole angle_fourier angle_fourier_simple angle_harmonic angle_hybrid + angle_mm3 angle_none angle_quartic angle_sdk diff --git a/doc/src/bond_mm3.txt b/doc/src/bond_mm3.txt new file mode 100644 index 0000000000..bfc802036e --- /dev/null +++ b/doc/src/bond_mm3.txt @@ -0,0 +1,58 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +bond_style mm3 command :h3 + +[Syntax:] + +bond_style mm3 :pre + +[Examples:] + +bond_style mm3 +bond_coeff 1 100.0 107.0 :pre + +[Description:] + +The {mm3} bond style uses the potential that is anharmonic in the bond +as defined in "(Allinger)"_#mm3-allinger1989 + +:c,image(Eqs/bond_mm3.jpg) + +where r0 is the equilibrium value of the bond, and K is a +prefactor. The anharmonic prefactors have units angstrom^(-n): +-2.55 angstrom^(-1) and (7/12)2.55$2 angstrom^(-2). The code takes +care of the necessary unit conversion for these factors internally. +Note that the MM3 papers contains an error in Eq (1): +(7/12)2.55 should be replaced with (7/12)2.55^2 + +The following coefficients must be defined for each bond type via the +"bond_coeff"_bond_coeff.html command as in the example above, or in +the data file or restart files read by the "read_data"_read_data.html +or "read_restart"_read_restart.html commands: + +K (energy/distance^2) +r0 (distance) :ul + +[Restrictions:] + +This bond style can only be used if LAMMPS was built with the +USER_YAFF package. See the "Build package"_Build_package.html doc +page for more info. + +[Related commands:] + +"bond_coeff"_bond_coeff.html + +[Default:] none + +:line + +:link(mm3-allinger1989) +[(Allinger)] Allinger, Yuh, Lii, JACS, 111,94, 8551-8566 +(1990), diff --git a/doc/src/bond_style.txt b/doc/src/bond_style.txt index ae19f2369d..aba6d3a778 100644 --- a/doc/src/bond_style.txt +++ b/doc/src/bond_style.txt @@ -86,6 +86,7 @@ accelerated styles exist. "harmonic"_bond_harmonic.html - harmonic bond "harmonic/shift"_bond_harmonic_shift.html - shifted harmonic bond "harmonic/shift/cut"_bond_harmonic_shift_cut.html - shifted harmonic bond with a cutoff +"mm3"_bond_mm3.html - MM3 anharmonic bond "morse"_bond_morse.html - Morse bond "nonlinear"_bond_nonlinear.html - nonlinear bond "oxdna/fene"_bond_oxdna.html - modified FENE bond suitable for DNA modeling diff --git a/doc/src/bonds.txt b/doc/src/bonds.txt index d33515eb88..48896e711c 100644 --- a/doc/src/bonds.txt +++ b/doc/src/bonds.txt @@ -13,6 +13,7 @@ Bond Styles :h1 bond_harmonic_shift bond_harmonic_shift_cut bond_hybrid + bond_mm3 bond_morse bond_none bond_nonlinear diff --git a/doc/src/improper_distharm.txt b/doc/src/improper_distharm.txt new file mode 100644 index 0000000000..52815e76aa --- /dev/null +++ b/doc/src/improper_distharm.txt @@ -0,0 +1,53 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +improper_style distharm command :h3 + +[Syntax:] + +improper_style distharm + +[Examples:] + +improper_style distharm +improper_coeff 1 25.0 0.5 :pre + +[Description:] + +The {distharm} improper style uses the potential + +:c,image(Eqs/improper_distharm.jpg) + +where d is the oriented distance between the central atom and the plane formed +by the other three atoms. If the 4 atoms in an improper quadruplet +(listed in the data file read by the "read_data"_read_data.html +command) are ordered I,J,K,L then the L-atom is assumed to be the +central atom. Note that this is different from the convention used +in the improper_style distance. The distance d is oriented and can take +on negative values. This may lead to unwanted behavior if d0 is not equal to zero. + +The following coefficients must be defined for each improper type via +the improper_coeff command as in the example above, or in the data +file or restart files read by the read_data or read_restart commands: + +K (energy/distance^2) +d0 (distance) :ul + +:line + +[Restrictions:] + +This improper style can only be used if LAMMPS was built with the +USER-YAFF package. See the "Build package"_Build_package.html doc +page for more info. + +[Related commands:] + +"improper_coeff"_improper_coeff.html + +[Default:] none diff --git a/doc/src/improper_sqdistharm.txt b/doc/src/improper_sqdistharm.txt new file mode 100644 index 0000000000..7473fd8c5d --- /dev/null +++ b/doc/src/improper_sqdistharm.txt @@ -0,0 +1,54 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +improper_style sqdistharm command :h3 + +[Syntax:] + +improper_style sqdistharm + +[Examples:] + +improper_style sqdistharm +improper_coeff 1 50.0 0.1 :pre + +[Description:] + +The {sqdistharm} improper style uses the potential + +:c,image(Eqs/improper_sqdistharm.jpg) + +where d is the distance between the central atom and the plane formed +by the other three atoms. If the 4 atoms in an improper quadruplet +(listed in the data file read by the "read_data"_read_data.html +command) are ordered I,J,K,L then the L-atom is assumed to be the +central atom. Note that this is different from the convention used +in the improper_style distance. + +The following coefficients must be defined for each improper type via +the improper_coeff command as in the example above, or in the data +file or restart files read by the read_data or read_restart commands: + +K (energy/distance^4) +d0^2 (distance^2) :ul + +Note that d0^2 (in units distance^2) has be provided and not d0. + +:line + +[Restrictions:] + +This improper style can only be used if LAMMPS was built with the +USER-MISC package. See the "Build package"_Build_package.html doc +page for more info. + +[Related commands:] + +"improper_coeff"_improper_coeff.html + +[Default:] none diff --git a/doc/src/improper_style.txt b/doc/src/improper_style.txt index c5e0be8a81..cd72da4d07 100644 --- a/doc/src/improper_style.txt +++ b/doc/src/improper_style.txt @@ -78,11 +78,13 @@ more of (g,i,k,o,t) to indicate which accelerated styles exist. "cossq"_improper_cossq.html - improper with a cosine squared term "cvff"_improper_cvff.html - CVFF improper "distance"_improper_distance.html - improper based on distance between atom planes +"distharm"_improper_distharm.html - improper that is harmonic in the out-of-plane distance "fourier"_improper_fourier.html - improper with multiple cosine terms "harmonic"_improper_harmonic.html - harmonic improper "inversion/harmonic"_improper_inversion_harmonic.html - harmonic improper with Wilson-Decius out-of-plane definition "ring"_improper_ring.html - improper which prevents planar conformations "umbrella"_improper_umbrella.html - DREIDING improper :ul +"sqdistharm"_improper_sqdistharm.html - improper that is harmonic in the square of the out-of-plane distance :line diff --git a/doc/src/impropers.txt b/doc/src/impropers.txt index ca6c839c95..ce829197fe 100644 --- a/doc/src/impropers.txt +++ b/doc/src/impropers.txt @@ -9,6 +9,7 @@ Improper Styles :h1 improper_cossq improper_cvff improper_distance + improper_distharm improper_fourier improper_harmonic improper_hybrid @@ -16,6 +17,7 @@ Improper Styles :h1 improper_none improper_ring improper_umbrella + improper_sqdistharm improper_zero END_RST --> diff --git a/doc/src/lammps.book b/doc/src/lammps.book index a1bbaf896e..73b2723dc3 100644 --- a/doc/src/lammps.book +++ b/doc/src/lammps.book @@ -687,11 +687,13 @@ angle_cosine_periodic.html angle_cosine_shift.html angle_cosine_shift_exp.html angle_cosine_squared.html +angle_cross.html angle_dipole.html angle_fourier.html angle_fourier_simple.html angle_harmonic.html angle_hybrid.html +angle_mm3.html angle_none.html angle_quartic.html angle_sdk.html @@ -725,6 +727,7 @@ improper_class2.html improper_cossq.html improper_cvff.html improper_distance.html +improper_distharm.html improper_fourier.html improper_harmonic.html improper_hybrid.html @@ -732,6 +735,7 @@ improper_inversion_harmonic.html improper_none.html improper_ring.html improper_umbrella.html +improper_sqdistharm.html improper_zero.html lammps_commands_kspace.html diff --git a/doc/src/pair_lj_switch3_coulgauss.txt b/doc/src/pair_lj_switch3_coulgauss.txt new file mode 100644 index 0000000000..b13b78a744 --- /dev/null +++ b/doc/src/pair_lj_switch3_coulgauss.txt @@ -0,0 +1,86 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +pair_style lj/switch3/coulgauss/long command :h3 + +[Syntax:] + +pair_style style args :pre + +style = {lj/switch3/coulgauss/long} +args = list of arguments for a particular style :ul + {lj/switch3/coulgauss/long} args = cutoff (cutoff2) width + cutoff = global cutoff for LJ (and Coulombic if only 1 arg) (distance units) + cutoff2 = global cutoff for Coulombic (optional) (distance units) + width = width parameter of the smoothing function (distance units) :pre + +[Examples:] + +pair_style lj/switch3/coulgauss/long 12.0 3.0 +pair_coeff 1 0.2 2.5 1.2 :pre + +pair_style lj/switch3/coulgauss/long 12.0 10.0 3.0 +pair_coeff 1 0.2 2.5 1.2 :pre + +[Description:] + +The {lj/switch3/coulgauss} style evaluates the LJ +vdW potentia + +:c,image(Eqs/pair_lj_switch3.jpg) + +, which goes smoothly to zero at the cutoff r_c as defined +by the switching function + +:c,image(Eqs/pair_switch3.jpg) + +where w is the width defined in the arguments. This potential +is combined with Coulomb interaction between Gaussian charge densities: + +:c,image(Eqs/pair_coulgauss.jpg) + +where qi and qj are the +charges on the 2 atoms, epsilon is the dielectric constant which +can be set by the "dielectric"_dielectric.html command, gamma_i and gamma_j +are the widths of the Gaussian charge distribution and erf() is the error-function. +This style has to be used in conjunction with the "kspace_style"_kspace_style.html command + +If one cutoff is specified it is used for both the vdW and Coulomb +terms. If two cutoffs are specified, the first is used as the cutoff +for the vdW terms, and the second is the cutoff for the Coulombic term. + +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: + +epsilon (energy) +sigma (distance) +gamma (distance) :ul + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +Shifting the potential energy is not necessary because the switching +function ensures that the potential is zero at the cut-off. + + +[Restrictions:] + +These styles are part of the USER-YAFF package. They are only +enabled if LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] none + diff --git a/doc/src/pair_mm3_switch3_coulgauss.txt b/doc/src/pair_mm3_switch3_coulgauss.txt new file mode 100644 index 0000000000..3e0e24150e --- /dev/null +++ b/doc/src/pair_mm3_switch3_coulgauss.txt @@ -0,0 +1,88 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +pair_style mm3/switch3/coulgauss/long command :h3 + +[Syntax:] + +pair_style style args :pre + +style = {mm3/switch3/coulgauss/long} +args = list of arguments for a particular style :ul + {mm3/switch3/coulgauss/long} args = cutoff (cutoff2) width + cutoff = global cutoff for MM3 (and Coulombic if only 1 arg) (distance units) + cutoff2 = global cutoff for Coulombic (optional) (distance units) + width = width parameter of the smoothing function (distance units) :pre + +[Examples:] + +pair_style mm3/switch3/coulgauss/long 12.0 3.0 +pair_coeff 1 0.2 2.5 1.2 :pre + +pair_style mm3/switch3/coulgauss/long 12.0 10.0 3.0 +pair_coeff 1 0.2 2.5 1.2 :pre + +[Description:] + +The {mm3/switch3/coulgauss} style evaluates the MM3 +vdW potential "(Allinger)"_#mm3-allinger1989 + +:c,image(Eqs/pair_mm3_switch3.jpg) + +, which goes smoothly to zero at the cutoff r_c as defined +by the switching function + +:c,image(Eqs/pair_switch3.jpg) + +where w is the width defined in the arguments. This potential +is combined with Coulomb interaction between Gaussian charge densities: + +:c,image(Eqs/pair_coulgauss.jpg) + +where qi and qj are the +charges on the 2 atoms, epsilon is the dielectric constant which +can be set by the "dielectric"_dielectric.html command, gamma_i and gamma_j +are the widths of the Gaussian charge distribution and erf() is the error-function. +This style has to be used in conjunction with the "kspace_style"_kspace_style.html command + +If one cutoff is specified it is used for both the vdW and Coulomb +terms. If two cutoffs are specified, the first is used as the cutoff +for the vdW terms, and the second is the cutoff for the Coulombic term. + +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: + +epsilon (energy) +r_v (distance) +gamma (distance) :ul + +:line + +[Mixing, shift, table, tail correction, restart, rRESPA info]: + +Mixing rules are fixed for this style as defined above. + +Shifting the potential energy is not necessary because the switching +function ensures that the potential is zero at the cut-off. + + +[Restrictions:] + +These styles are part of the USER-YAFF package. They are only +enabled if LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +[Related commands:] + +"pair_coeff"_pair_coeff.html + +[Default:] none + diff --git a/doc/src/pair_style.txt b/doc/src/pair_style.txt index 714ec55d2e..9782772456 100644 --- a/doc/src/pair_style.txt +++ b/doc/src/pair_style.txt @@ -220,6 +220,7 @@ accelerated styles exist. "lj/sf/dipole/sf"_pair_dipole.html - LJ with dipole interaction with shifted forces "lj/smooth"_pair_lj_smooth.html - smoothed Lennard-Jones potential "lj/smooth/linear"_pair_lj_smooth_linear.html - linear smoothed LJ potential +"lj/switch3/coulgauss"_pair_lj_switch3_coulgauss - smoothed LJ vdW potential with Gaussian electrostatics "lj96/cut"_pair_lj96.html - Lennard-Jones 9/6 potential "lubricate"_pair_lubricate.html - hydrodynamic lubrication forces "lubricate/poly"_pair_lubricate.html - hydrodynamic lubrication forces with polydispersity @@ -232,6 +233,7 @@ accelerated styles exist. "meam/sw/spline"_pair_meam_sw_spline.html - splined version of MEAM with a Stillinger-Weber term "mgpt"_pair_mgpt.html - simplified model generalized pseudopotential theory (MGPT) potential "mie/cut"_pair_mie.html - Mie potential +"mm3/switch3/coulgauss"_pair_mm3_switch3_coulgauss - smoothed MM3 vdW potential with Gaussian electrostatics "momb"_pair_momb.html - Many-Body Metal-Organic (MOMB) force field "morse"_pair_morse.html - Morse potential "morse/smooth/linear"_pair_morse.html - linear smoothed Morse potential diff --git a/doc/src/pairs.txt b/doc/src/pairs.txt index 2810252402..62cbb7f56d 100644 --- a/doc/src/pairs.txt +++ b/doc/src/pairs.txt @@ -61,6 +61,7 @@ Pair Styles :h1 pair_lj_smooth pair_lj_smooth_linear pair_lj_soft + pair_lj_switch3_coulgauss pair_lubricate pair_lubricateU pair_mdf @@ -70,6 +71,7 @@ Pair Styles :h1 pair_meso pair_mgpt pair_mie + pair_mm3_switch3_coulgauss pair_momb pair_morse pair_multi_lucy diff --git a/src/USER-YAFF/README b/src/USER-YAFF/README new file mode 100644 index 0000000000..4fabe67390 --- /dev/null +++ b/src/USER-YAFF/README @@ -0,0 +1,8 @@ +This package implements the styles that are needed to use typical force fields +generated by QuickFF for the simulation of metal-organic frameworks. The +QuickFF methodology is detailed in following papers: + Vanduyfhuys et al., J. Comput. Chem., 36 (13), 1015-1027 (2015) http://dx.doi.org/10.1002/jcc.23877 + Vanduyfhuys et al., J. Comput. Chem., 39 (16), 999-1011 (2018) http://dx.doi.org/10.1002/jcc.25173 +The corresponding software package can be found on http://molmod.github.io/QuickFF + + diff --git a/src/USER-YAFF/angle_cross.cpp b/src/USER-YAFF/angle_cross.cpp new file mode 100644 index 0000000000..a31931a23b --- /dev/null +++ b/src/USER-YAFF/angle_cross.cpp @@ -0,0 +1,344 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "angle_cross.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleCross::AngleCross(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleCross::~AngleCross() +{ + if (copymode) return; + + if (allocated) { + memory->destroy(setflag); + memory->destroy(kss); + memory->destroy(kbs0); + memory->destroy(kbs1); + memory->destroy(r00); + memory->destroy(r01); + memory->destroy(theta0); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleCross::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double dtheta,dtheta2,dtheta3,dtheta4,de_angle; + double dr1,dr2,tk1,tk2,aa1,aa2,aa11,aa12,aa21,aa22; + double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22,b1,b2; + double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22; + + eangle = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + // force & energy for bond-bond term + dr1 = r1 - r00[type]; + dr2 = r2 - r01[type]; + tk1 = kss[type] * dr1; + tk2 = kss[type] * dr2; + + f1[0] = -delx1*tk2/r1; + f1[1] = -dely1*tk2/r1; + f1[2] = -delz1*tk2/r1; + + f3[0] = -delx2*tk1/r2; + f3[1] = -dely2*tk1/r2; + f3[2] = -delz2*tk1/r2; + + if (eflag) eangle = kss[type]*dr1*dr2; + + // force & energy for bond-angle term + dtheta = acos(c) - theta0[type]; + + aa1 = s * dr1 * kbs0[type]; + aa2 = s * dr2 * kbs1[type]; + + aa11 = aa1 * c / rsq1; + aa12 = -aa1 / (r1 * r2); + aa21 = aa2 * c / rsq1; + aa22 = -aa2 / (r1 * r2); + + vx11 = (aa11 * delx1) + (aa12 * delx2); + vx12 = (aa21 * delx1) + (aa22 * delx2); + vy11 = (aa11 * dely1) + (aa12 * dely2); + vy12 = (aa21 * dely1) + (aa22 * dely2); + vz11 = (aa11 * delz1) + (aa12 * delz2); + vz12 = (aa21 * delz1) + (aa22 * delz2); + + aa11 = aa1 * c / rsq2; + aa21 = aa2 * c / rsq2; + + vx21 = (aa11 * delx2) + (aa12 * delx1); + vx22 = (aa21 * delx2) + (aa22 * delx1); + vy21 = (aa11 * dely2) + (aa12 * dely1); + vy22 = (aa21 * dely2) + (aa22 * dely1); + vz21 = (aa11 * delz2) + (aa12 * delz1); + vz22 = (aa21 * delz2) + (aa22 * delz1); + + b1 = kbs0[type] * dtheta / r1; + b2 = kbs1[type] * dtheta / r2; + + f1[0] -= vx11 + b1*delx1 + vx12; + f1[1] -= vy11 + b1*dely1 + vy12; + f1[2] -= vz11 + b1*delz1 + vz12; + + f3[0] -= vx21 + b2*delx2 + vx22; + f3[1] -= vy21 + b2*dely2 + vy22; + f3[2] -= vz21 + b2*delz2 + vz22; + + if (eflag) eangle += kbs0[type]*dr1*dtheta + kbs1[type]*dr2*dtheta; + + // apply force to each of 3 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleCross::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(kss,n+1,"angle:kss"); + memory->create(kbs0,n+1,"angle:kbs0"); + memory->create(kbs1,n+1,"angle:kbs1"); + memory->create(r00,n+1,"angle:r00"); + memory->create(r01,n+1,"angle:r01"); + memory->create(theta0,n+1,"angle:theta0"); + memory->create(setflag,n+1,"angle:setflag"); + + for (int i = 1; i <= n; i++) + setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs +------------------------------------------------------------------------- */ + +void AngleCross::coeff(int narg, char **arg) +{ + if (narg != 7) error->all(FLERR,"Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(FLERR,arg[0],atom->nangletypes,ilo,ihi); + + int count = 0; + + double kss_one = force->numeric(FLERR,arg[1]); + double kbs0_one = force->numeric(FLERR,arg[2]); + double kbs1_one = force->numeric(FLERR,arg[3]); + double r0_one = force->numeric(FLERR,arg[4]); + double r1_one = force->numeric(FLERR,arg[5]); + double theta0_one = force->numeric(FLERR,arg[6]); + + for (int i = ilo; i <= ihi; i++) { + kss[i] = kss_one; + kbs0[i] = kbs0_one; + kbs1[i] = kbs1_one; + r00[i] = r0_one; + r01[i] = r1_one; + // Convert theta0 from degrees to radians + theta0[i] = theta0_one*MY_PI/180.0; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); +} + +/* ---------------------------------------------------------------------- */ + +double AngleCross::equilibrium_angle(int i) +{ + return theta0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleCross::write_restart(FILE *fp) +{ + fwrite(&kss[1],sizeof(double),atom->nangletypes,fp); + fwrite(&kbs0[1],sizeof(double),atom->nangletypes,fp); + fwrite(&kbs1[1],sizeof(double),atom->nangletypes,fp); + fwrite(&r00[1],sizeof(double),atom->nangletypes,fp); + fwrite(&r01[1],sizeof(double),atom->nangletypes,fp); + fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleCross::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&kss[1],sizeof(double),atom->nangletypes,fp); + fread(&kbs0[1],sizeof(double),atom->nangletypes,fp); + fread(&kbs1[1],sizeof(double),atom->nangletypes,fp); + fread(&r00[1],sizeof(double),atom->nangletypes,fp); + fread(&r01[1],sizeof(double),atom->nangletypes,fp); + fread(&theta0[1],sizeof(double),atom->nangletypes,fp); + } + + MPI_Bcast(&kss[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&kbs0[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&kbs1[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&r00[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&r01[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void AngleCross::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nangletypes; i++) + fprintf(fp,"%d %g %g %g %g\n", + i,kss[i],kbs0[i],kbs1[i],r00[i],r01[i],theta0[i]/MY_PI*180.0); +} + +/* ---------------------------------------------------------------------- */ + +double AngleCross::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + double dtheta = acos(c) - theta0[type]; + double dr1 = r1 - r00[type]; + double dr2 = r2 - r01[type]; + double energy = kss[type]*dr1*dr2+kbs0[type]*dr1*dtheta + kbs1[type]*dr2*dtheta; + return energy; +} diff --git a/src/USER-YAFF/angle_cross.h b/src/USER-YAFF/angle_cross.h new file mode 100644 index 0000000000..7709c10414 --- /dev/null +++ b/src/USER-YAFF/angle_cross.h @@ -0,0 +1,57 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS + +AngleStyle(cross,AngleCross) + +#else + +#ifndef LMP_ANGLE_CROSS_H +#define LMP_ANGLE_CROSS_H + +#include +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleCross : public Angle { + public: + AngleCross(class LAMMPS *); + virtual ~AngleCross(); + virtual void compute(int, int); + void coeff(int, char **); + double equilibrium_angle(int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_data(FILE *); + double single(int, int, int, int); + + protected: + double *kss,*kbs0,*kbs1,*r00,*r01,*theta0; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args for angle coefficients + +Self-explanatory. Check the input script or data file. + +*/ diff --git a/src/USER-YAFF/angle_mm3.cpp b/src/USER-YAFF/angle_mm3.cpp new file mode 100644 index 0000000000..ccd9b2e20e --- /dev/null +++ b/src/USER-YAFF/angle_mm3.cpp @@ -0,0 +1,288 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "angle_mm3.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +AngleMM3::AngleMM3(LAMMPS *lmp) : Angle(lmp) {} + +/* ---------------------------------------------------------------------- */ + +AngleMM3::~AngleMM3() +{ + if (copymode) return; + + if (allocated) { + memory->destroy(setflag); + memory->destroy(theta0); + memory->destroy(k2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleMM3::compute(int eflag, int vflag) +{ + int i1,i2,i3,n,type; + double delx1,dely1,delz1,delx2,dely2,delz2; + double eangle,f1[3],f3[3]; + double dtheta,dtheta2,dtheta3,dtheta4,de_angle; + double dr1,dr2,tk1,tk2,aa1,aa2,aa11,aa12,aa21,aa22; + double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22,b1,b2; + double vx11,vx12,vy11,vy12,vz11,vz12,vx21,vx22,vy21,vy22,vz21,vz22; + + eangle = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **anglelist = neighbor->anglelist; + int nanglelist = neighbor->nanglelist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nanglelist; n++) { + i1 = anglelist[n][0]; + i2 = anglelist[n][1]; + i3 = anglelist[n][2]; + type = anglelist[n][3]; + + // 1st bond + + delx1 = x[i1][0] - x[i2][0]; + dely1 = x[i1][1] - x[i2][1]; + delz1 = x[i1][2] - x[i2][2]; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + r1 = sqrt(rsq1); + + // 2nd bond + + delx2 = x[i3][0] - x[i2][0]; + dely2 = x[i3][1] - x[i2][1]; + delz2 = x[i3][2] - x[i2][2]; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + r2 = sqrt(rsq2); + + // angle (cos and sin) + + c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + // force & energy for angle term + + dtheta = acos(c) - theta0[type]; + dtheta2 = dtheta*dtheta; + dtheta3 = dtheta2*dtheta; + dtheta4 = dtheta3*dtheta; + // MM3 angle term, taking into account that dtheta is expressed in rad + de_angle = 2.0*k2[type]*dtheta*(1.0-1.203211*dtheta+0.367674*dtheta2-0.3239159*dtheta3+0.711270*dtheta4); + + a = -de_angle*s; + a11 = a*c / rsq1; + a12 = -a / (r1*r2); + a22 = a*c / rsq2; + + f1[0] = a11*delx1 + a12*delx2; + f1[1] = a11*dely1 + a12*dely2; + f1[2] = a11*delz1 + a12*delz2; + + f3[0] = a22*delx2 + a12*delx1; + f3[1] = a22*dely2 + a12*dely1; + f3[2] = a22*delz2 + a12*delz1; + // MM3 angle term, taking into account that dtheta is expressed in rad + if (eflag) eangle = k2[type]*dtheta2*(1.0-0.802141*dtheta+0.183837*dtheta2-0.131664*dtheta3+0.237090*dtheta4); + + // apply force to each of 3 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= f1[0] + f3[0]; + f[i2][1] -= f1[1] + f3[1]; + f[i2][2] -= f1[2] + f3[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, + delx1,dely1,delz1,delx2,dely2,delz2); + } +} + +/* ---------------------------------------------------------------------- */ + +void AngleMM3::allocate() +{ + allocated = 1; + int n = atom->nangletypes; + + memory->create(setflag,n+1,"angle:setflag"); + memory->create(k2,n+1,"angle:k2"); + memory->create(theta0,n+1,"angle:theta0"); + for (int i = 1; i <= n; i++) + setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs + else -> Angle coeffs +------------------------------------------------------------------------- */ + +void AngleMM3::coeff(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR,"Incorrect args for angle coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(FLERR,arg[0],atom->nangletypes,ilo,ihi); + + int count = 0; + + double k2_one = force->numeric(FLERR,arg[1]); + double theta0_one = force->numeric(FLERR,arg[2]); + + // convert theta0 from degrees to radians + + for (int i = ilo; i <= ihi; i++) { + k2[i] = k2_one; + theta0[i] = theta0_one/180.0 * MY_PI; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); + +} + +/* ---------------------------------------------------------------------- */ + +double AngleMM3::equilibrium_angle(int i) +{ + return theta0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void AngleMM3::write_restart(FILE *fp) +{ + fwrite(&k2[1],sizeof(double),atom->nangletypes,fp); + fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void AngleMM3::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k2[1],sizeof(double),atom->nangletypes,fp); + fread(&theta0[1],sizeof(double),atom->nangletypes,fp); + } + + MPI_Bcast(&k2[1],atom->nangletypes,MPI_DOUBLE,0,world); + MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void AngleMM3::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nangletypes; i++) + fprintf(fp,"%d %g %g\n", + i,k2[i],theta0[i]/MY_PI*180.0); +} + +/* ---------------------------------------------------------------------- */ + +double AngleMM3::single(int type, int i1, int i2, int i3) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + double dtheta = acos(c) - theta0[type]; + double dtheta2 = dtheta*dtheta; + double dtheta3 = dtheta2*dtheta; + double dtheta4 = dtheta3*dtheta; + + double energy = k2[type]*dtheta2*(1.0-0.802141*dtheta+0.183837*dtheta2-0.131664*dtheta3+0.237090*dtheta4); + + return energy; +} diff --git a/src/USER-YAFF/angle_mm3.h b/src/USER-YAFF/angle_mm3.h new file mode 100644 index 0000000000..2d19b4d1b4 --- /dev/null +++ b/src/USER-YAFF/angle_mm3.h @@ -0,0 +1,57 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef ANGLE_CLASS + +AngleStyle(mm3,AngleMM3) + +#else + +#ifndef LMP_ANGLE_MM3_H +#define LMP_ANGLE_MM3_H + +#include +#include "angle.h" + +namespace LAMMPS_NS { + +class AngleMM3 : public Angle { + public: + AngleMM3(class LAMMPS *); + virtual ~AngleMM3(); + virtual void compute(int, int); + void coeff(int, char **); + double equilibrium_angle(int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_data(FILE *); + double single(int, int, int, int); + + protected: + double *theta0,*k2; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args for angle coefficients + +Self-explanatory. Check the input script or data file. + +*/ diff --git a/src/USER-YAFF/bond_mm3.cpp b/src/USER-YAFF/bond_mm3.cpp new file mode 100644 index 0000000000..f0e4197c6b --- /dev/null +++ b/src/USER-YAFF/bond_mm3.cpp @@ -0,0 +1,220 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande +------------------------------------------------------------------------- */ + +#include +#include +#include "bond_mm3.h" +#include "atom.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondMM3::BondMM3(LAMMPS *lmp) : Bond(lmp) {} + +/* ---------------------------------------------------------------------- */ + +BondMM3::~BondMM3() +{ + if (copymode) return; + + if (allocated) { + memory->destroy(setflag); + memory->destroy(r0); + memory->destroy(k2); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondMM3::compute(int eflag, int vflag) +{ + int i1,i2,n,type; + double delx,dely,delz,ebond,fbond; + double rsq,r,dr,dr2,de_bond,K3,K4; + + ebond = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + /* + E = K(r-r0)^2 [1-2.55*(r-r0)+(7/12)*2.55^(2)*(r-r0)^2] + with -2.55 in angstrom^(-1) and (7/12)*2.55^(2) in angstrom^(-2) + These prefactors are converted here to the correct units + */ + K3 = -2.55/force->angstrom; + K4 = 7.0/12.0*2.55*2.55/force->angstrom/force->angstrom; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + + delx = x[i1][0] - x[i2][0]; + dely = x[i1][1] - x[i2][1]; + delz = x[i1][2] - x[i2][2]; + + rsq = delx*delx + dely*dely + delz*delz; + r = sqrt(rsq); + dr = r - r0[type]; + dr2 = dr*dr; + + // force & energy + + de_bond = 2.0*k2[type]*dr*(1.0 + 1.5*K3*dr + 2.0*K4*dr2); + if (r > 0.0) fbond = -de_bond/r; + else fbond = 0.0; + + if (eflag) ebond = k2[type]*dr2*(1.0+K3*dr+K4*dr2); + + // apply force to each of 2 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += delx*fbond; + f[i1][1] += dely*fbond; + f[i1][2] += delz*fbond; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] -= delx*fbond; + f[i2][1] -= dely*fbond; + f[i2][2] -= delz*fbond; + } + + if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondMM3::allocate() +{ + allocated = 1; + int n = atom->nbondtypes; + + memory->create(r0,n+1,"bond:r0"); + memory->create(k2,n+1,"bond:k2"); + + memory->create(setflag,n+1,"bond:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs from one line in input script or data file +------------------------------------------------------------------------- */ + +void BondMM3::coeff(int narg, char **arg) +{ + if (narg != 3) error->all(FLERR,"Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi); + + double k2_one = force->numeric(FLERR,arg[1]); + double r0_one = force->numeric(FLERR,arg[2]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k2[i] = k2_one; + r0[i] = r0_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + return an equilbrium bond length +------------------------------------------------------------------------- */ + +double BondMM3::equilibrium_distance(int i) +{ + return r0[i]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondMM3::write_restart(FILE *fp) +{ + fwrite(&k2[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondMM3::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k2[1],sizeof(double),atom->nbondtypes,fp); + fread(&r0[1],sizeof(double),atom->nbondtypes,fp); + } + MPI_Bcast(&k2[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondMM3::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) + fprintf(fp,"%d %g %g\n",i,k2[i],r0[i]); +} + +/* ---------------------------------------------------------------------- */ + +double BondMM3::single(int type, double rsq, int i, int j, double &fforce) +{ + /* + E = K(r-r0)^2 [1-2.55*(r-r0)+(7/12)*2.55^(2)*(r-r0)^2] + with -2.55 in angstrom^(-1) and (7/12)*2.55^(2) in angstrom^(-2) + These prefactors are converted here to the correct units + */ + double K3 = -2.55/force->angstrom; + double K4 = 7.0/12.0*2.55*2.55/force->angstrom/force->angstrom; + double r = sqrt(rsq); + double dr = r - r0[type]; + double dr2 = dr*dr; + double de_bond = 2.0*k2[type]*dr*(1.0 + 1.5*K3*dr + 2.0*K4*dr2); + if (r > 0.0) fforce = -de_bond/r; + else fforce = 0.0; + return k2[type]*dr2*(1.0+K3*dr+K4*dr2); +} diff --git a/src/USER-YAFF/bond_mm3.h b/src/USER-YAFF/bond_mm3.h new file mode 100644 index 0000000000..9711d89529 --- /dev/null +++ b/src/USER-YAFF/bond_mm3.h @@ -0,0 +1,57 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef BOND_CLASS + +BondStyle(mm3,BondMM3) + +#else + +#ifndef LMP_BOND_MM3_H +#define LMP_BOND_MM3_H + +#include +#include "bond.h" + +namespace LAMMPS_NS { + +class BondMM3 : public Bond { + public: + BondMM3(class LAMMPS *); + virtual ~BondMM3(); + virtual void compute(int, int); + void coeff(int, char **); + double equilibrium_distance(int); + void write_restart(FILE *); + void read_restart(FILE *); + void write_data(FILE *); + double single(int, double, int, int, double &); + + protected: + double *r0,*k2; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Incorrect args for bond coefficients + +Self-explanatory. Check the input script or data file. + +*/ diff --git a/src/USER-YAFF/improper_distharm.cpp b/src/USER-YAFF/improper_distharm.cpp new file mode 100644 index 0000000000..9a54afed9a --- /dev/null +++ b/src/USER-YAFF/improper_distharm.cpp @@ -0,0 +1,269 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande, heavily based on the + improper_distance code by Paolo Raiteri (Curtin University) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "improper_distharm.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define TOLERANCE 0.05 +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +ImproperDistHarm::ImproperDistHarm(LAMMPS *lmp) : Improper(lmp) {} + +/* ---------------------------------------------------------------------- */ + +ImproperDistHarm::~ImproperDistHarm() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(chi); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperDistHarm::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,n,type; + double xab, yab, zab; // bond 1-2 + double xac, yac, zac; // bond 1-3 + double xad, yad, zad; // bond 1-4 + double xbc, ybc, zbc; // bond 2-3 + double xbd, ybd, zbd; // bond 2-4 + double xcd, ycd, zcd; // bond 3-4 + double xna, yna, zna, rna; // normal + double da; + + double eimproper,f1[3],f2[3],f3[3],f4[3]; + double domega,a; + + eimproper = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **improperlist = neighbor->improperlist; + int nimproperlist = neighbor->nimproperlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nimproperlist; n++) { + i1 = improperlist[n][0]; + i2 = improperlist[n][1]; + i3 = improperlist[n][2]; + i4 = improperlist[n][3]; + type = improperlist[n][4]; + + // geometry of 4-body + // 4 is the central atom + // 1-2-3 are ment to be equivalent + // I need the bonds between 2-3 and 3-4 to get the plane normal + // Then I need the bond 1-4 to project it onto the normal to the plane + + // bond 1->2 + xab = x[i2][0] - x[i1][0]; + yab = x[i2][1] - x[i1][1]; + zab = x[i2][2] - x[i1][2]; + domain->minimum_image(xab,yab,zab); + + // bond 1->3 + xac = x[i3][0] - x[i1][0]; + yac = x[i3][1] - x[i1][1]; + zac = x[i3][2] - x[i1][2]; + domain->minimum_image(xac,yac,zac); + + // bond 1->4 + xad = x[i4][0] - x[i1][0]; + yad = x[i4][1] - x[i1][1]; + zad = x[i4][2] - x[i1][2]; + domain->minimum_image(xad,yad,zad); + + // bond 2-3 + xbc = x[i3][0] - x[i2][0]; + ybc = x[i3][1] - x[i2][1]; + zbc = x[i3][2] - x[i2][2]; + domain->minimum_image(xbc,ybc,zbc); + + // bond 2-4 + xbd = x[i4][0] - x[i2][0]; + ybd = x[i4][1] - x[i2][1]; + zbd = x[i4][2] - x[i2][2]; + domain->minimum_image(xbd,ybd,zbd); + + // bond 3-4 + xcd = x[i4][0] - x[i3][0]; + ycd = x[i4][1] - x[i3][1]; + zcd = x[i4][2] - x[i3][2]; + domain->minimum_image(xcd,ycd,zcd); + + xna = ybc*zcd - zbc*ycd; + yna = -(xbc*zcd - zbc*xcd); + zna = xbc*ycd - ybc*xcd; + rna = 1.0 / sqrt(xna*xna+yna*yna+zna*zna); + xna *= rna; + yna *= rna; + zna *= rna; + + da = -(xna*xad + yna*yad + zna*zad); + + + domega = k[type]*(da - chi[type])*(da - chi[type]); + a = 2.0* k[type]*(da - chi[type]); + + if (eflag) eimproper = domega; + + f1[0] = a*( -xna); + f1[1] = a*( -yna); + f1[2] = a*( -zna); + f4[0] = a*( xna); + f4[1] = a*( yna); + f4[2] = a*( zna); + + f2[0] = a*( yad*zcd - zad*ycd )*rna + a*da*rna*( yna*zcd - zna*ycd); + f2[1] = a*( zad*xcd - xad*zcd )*rna + a*da*rna*( zna*xcd - xna*zcd); + f2[2] = a*( xad*ycd - yad*xcd )*rna + a*da*rna*( xna*ycd - yna*xcd); + + f3[0] = - a*( yad*zcd - zad*ycd )*rna - a*da*rna*( yna*zcd - zna*ycd); + f3[1] = - a*( zad*xcd - xad*zcd )*rna - a*da*rna*( zna*xcd - xna*zcd); + f3[2] = - a*( xad*ycd - yad*xcd )*rna - a*da*rna*( xna*ycd - yna*xcd); + + f3[0] += -a*( yad*zbc - zad*ybc )*rna - a*da*rna*( yna*zbc - zna*ybc); + f3[1] += -a*( zad*xbc - xad*zbc )*rna - a*da*rna*( zna*xbc - xna*zbc); + f3[2] += -a*( xad*ybc - yad*xbc )*rna - a*da*rna*( xna*ybc - yna*xbc); + f4[0] += a*( yad*zbc - zad*ybc )*rna + a*da*rna*( yna*zbc - zna*ybc); + f4[1] += a*( zad*xbc - xad*zbc )*rna + a*da*rna*( zna*xbc - xna*zbc); + f4[2] += a*( xad*ybc - yad*xbc )*rna + a*da*rna*( xna*ybc - yna*xbc); + + // apply force to each of 4 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += f2[0]; + f[i2][1] += f2[1]; + f[i2][2] += f2[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (newton_bond || i4 < nlocal) { + f[i4][0] += f4[0]; + f[i4][1] += f4[1]; + f[i4][2] += f4[2]; + } + + if (evflag) + ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f2,f3,f4, + xab,yab,zab,xac,yac,zac,xad-xac,yad-yac,zad-zac); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperDistHarm::allocate() +{ + allocated = 1; + int n = atom->nimpropertypes; + + memory->create(k,n+1,"improper:k"); + memory->create(chi,n+1,"improper:chi"); + + memory->create(setflag,n+1,"improper:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void ImproperDistHarm::coeff(int narg, char **arg) +{ +// if (which > 0) return; + if (narg != 3) error->all(FLERR,"Incorrect args for improper coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(FLERR,arg[0],atom->nimpropertypes,ilo,ihi); + + double k_one = force->numeric(FLERR,arg[1]); + double chi_one = force->numeric(FLERR,arg[2]); + + // convert chi from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + //chi[i] = chi_one/180.0 * PI; + chi[i] = chi_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void ImproperDistHarm::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&chi[1],sizeof(double),atom->nimpropertypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void ImproperDistHarm::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->nimpropertypes,fp); + fread(&chi[1],sizeof(double),atom->nimpropertypes,fp); + } + MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&chi[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1; +} diff --git a/src/USER-YAFF/improper_distharm.h b/src/USER-YAFF/improper_distharm.h new file mode 100644 index 0000000000..b8b9ae780e --- /dev/null +++ b/src/USER-YAFF/improper_distharm.h @@ -0,0 +1,47 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef IMPROPER_CLASS + +ImproperStyle(distharm,ImproperDistHarm) + +#else + +#ifndef LMP_IMPROPER_DISTHARM_H +#define LMP_IMPROPER_DISTHARM_H + +#include +#include "improper.h" + +namespace LAMMPS_NS { + +class ImproperDistHarm : public Improper { + public: + ImproperDistHarm(class LAMMPS *); + ~ImproperDistHarm(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + private: + double *k,*chi; + + void allocate(); +}; + +} + +#endif +#endif + diff --git a/src/USER-YAFF/improper_sqdistharm.cpp b/src/USER-YAFF/improper_sqdistharm.cpp new file mode 100644 index 0000000000..763d82f1c5 --- /dev/null +++ b/src/USER-YAFF/improper_sqdistharm.cpp @@ -0,0 +1,269 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande, heavily based on the + improper_distance code by Paolo Raiteri (Curtin University) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "improper_sqdistharm.h" +#include "atom.h" +#include "comm.h" +#include "neighbor.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define TOLERANCE 0.05 +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +ImproperSQDistHarm::ImproperSQDistHarm(LAMMPS *lmp) : Improper(lmp) {} + +/* ---------------------------------------------------------------------- */ + +ImproperSQDistHarm::~ImproperSQDistHarm() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(k); + memory->destroy(chi); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperSQDistHarm::compute(int eflag, int vflag) +{ + int i1,i2,i3,i4,n,type; + double xab, yab, zab; // bond 1-2 + double xac, yac, zac; // bond 1-3 + double xad, yad, zad; // bond 1-4 + double xbc, ybc, zbc; // bond 2-3 + double xbd, ybd, zbd; // bond 2-4 + double xcd, ycd, zcd; // bond 3-4 + double xna, yna, zna, rna; // normal + double da; + + double eimproper,f1[3],f2[3],f3[3],f4[3]; + double domega,a; + + eimproper = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + double **x = atom->x; + double **f = atom->f; + int **improperlist = neighbor->improperlist; + int nimproperlist = neighbor->nimproperlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + for (n = 0; n < nimproperlist; n++) { + i1 = improperlist[n][0]; + i2 = improperlist[n][1]; + i3 = improperlist[n][2]; + i4 = improperlist[n][3]; + type = improperlist[n][4]; + + // geometry of 4-body + // 4 is the central atom + // 1-2-3 are ment to be equivalent + // I need the bonds between 2-3 and 3-4 to get the plane normal + // Then I need the bond 1-4 to project it onto the normal to the plane + + // bond 1->2 + xab = x[i2][0] - x[i1][0]; + yab = x[i2][1] - x[i1][1]; + zab = x[i2][2] - x[i1][2]; + domain->minimum_image(xab,yab,zab); + + // bond 1->3 + xac = x[i3][0] - x[i1][0]; + yac = x[i3][1] - x[i1][1]; + zac = x[i3][2] - x[i1][2]; + domain->minimum_image(xac,yac,zac); + + // bond 1->4 + xad = x[i4][0] - x[i1][0]; + yad = x[i4][1] - x[i1][1]; + zad = x[i4][2] - x[i1][2]; + domain->minimum_image(xad,yad,zad); + + // bond 2-3 + xbc = x[i3][0] - x[i2][0]; + ybc = x[i3][1] - x[i2][1]; + zbc = x[i3][2] - x[i2][2]; + domain->minimum_image(xbc,ybc,zbc); + + // bond 2-4 + xbd = x[i4][0] - x[i2][0]; + ybd = x[i4][1] - x[i2][1]; + zbd = x[i4][2] - x[i2][2]; + domain->minimum_image(xbd,ybd,zbd); + + // bond 3-4 + xcd = x[i4][0] - x[i3][0]; + ycd = x[i4][1] - x[i3][1]; + zcd = x[i4][2] - x[i3][2]; + domain->minimum_image(xcd,ycd,zcd); + + xna = ybc*zcd - zbc*ycd; + yna = -(xbc*zcd - zbc*xcd); + zna = xbc*ycd - ybc*xcd; + rna = 1.0 / sqrt(xna*xna+yna*yna+zna*zna); + xna *= rna; + yna *= rna; + zna *= rna; + + da = -(xna*xad + yna*yad + zna*zad); + + domega = k[type]*(da*da - chi[type])*(da*da - chi[type]); + a = 4.0 * da* k[type]*(da*da - chi[type]); + + if (eflag) eimproper = domega; + + f1[0] = a*( -xna); + f1[1] = a*( -yna); + f1[2] = a*( -zna); + f4[0] = a*( xna); + f4[1] = a*( yna); + f4[2] = a*( zna); + + f2[0] = a*( yad*zcd - zad*ycd )*rna + a*da*rna*( yna*zcd - zna*ycd); + f2[1] = a*( zad*xcd - xad*zcd )*rna + a*da*rna*( zna*xcd - xna*zcd); + f2[2] = a*( xad*ycd - yad*xcd )*rna + a*da*rna*( xna*ycd - yna*xcd); + + f3[0] = - a*( yad*zcd - zad*ycd )*rna - a*da*rna*( yna*zcd - zna*ycd); + f3[1] = - a*( zad*xcd - xad*zcd )*rna - a*da*rna*( zna*xcd - xna*zcd); + f3[2] = - a*( xad*ycd - yad*xcd )*rna - a*da*rna*( xna*ycd - yna*xcd); + + f3[0] += -a*( yad*zbc - zad*ybc )*rna - a*da*rna*( yna*zbc - zna*ybc); + f3[1] += -a*( zad*xbc - xad*zbc )*rna - a*da*rna*( zna*xbc - xna*zbc); + f3[2] += -a*( xad*ybc - yad*xbc )*rna - a*da*rna*( xna*ybc - yna*xbc); + f4[0] += a*( yad*zbc - zad*ybc )*rna + a*da*rna*( yna*zbc - zna*ybc); + f4[1] += a*( zad*xbc - xad*zbc )*rna + a*da*rna*( zna*xbc - xna*zbc); + f4[2] += a*( xad*ybc - yad*xbc )*rna + a*da*rna*( xna*ybc - yna*xbc); + + + // apply force to each of 4 atoms + + if (newton_bond || i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += f2[0]; + f[i2][1] += f2[1]; + f[i2][2] += f2[2]; + } + + if (newton_bond || i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + + if (newton_bond || i4 < nlocal) { + f[i4][0] += f4[0]; + f[i4][1] += f4[1]; + f[i4][2] += f4[2]; + } + + if (evflag) + ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f2,f3,f4, + xab,yab,zab,xac,yac,zac,xad-xac,yad-yac,zad-zac); + } +} + +/* ---------------------------------------------------------------------- */ + +void ImproperSQDistHarm::allocate() +{ + allocated = 1; + int n = atom->nimpropertypes; + + memory->create(k,n+1,"improper:k"); + memory->create(chi,n+1,"improper:chi"); + + memory->create(setflag,n+1,"improper:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one type +------------------------------------------------------------------------- */ + +void ImproperSQDistHarm::coeff(int narg, char **arg) +{ +// if (which > 0) return; + if (narg != 3) error->all(FLERR,"Incorrect args for improper coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + force->bounds(FLERR,arg[0],atom->nimpropertypes,ilo,ihi); + + double k_one = force->numeric(FLERR,arg[1]); + double chi_one = force->numeric(FLERR,arg[2]); + + // convert chi from degrees to radians + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + k[i] = k_one; + //chi[i] = chi_one/180.0 * PI; + chi[i] = chi_one; + setflag[i] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void ImproperSQDistHarm::write_restart(FILE *fp) +{ + fwrite(&k[1],sizeof(double),atom->nimpropertypes,fp); + fwrite(&chi[1],sizeof(double),atom->nimpropertypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void ImproperSQDistHarm::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + fread(&k[1],sizeof(double),atom->nimpropertypes,fp); + fread(&chi[1],sizeof(double),atom->nimpropertypes,fp); + } + MPI_Bcast(&k[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + MPI_Bcast(&chi[1],atom->nimpropertypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1; +} diff --git a/src/USER-YAFF/improper_sqdistharm.h b/src/USER-YAFF/improper_sqdistharm.h new file mode 100644 index 0000000000..301b5066cb --- /dev/null +++ b/src/USER-YAFF/improper_sqdistharm.h @@ -0,0 +1,47 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef IMPROPER_CLASS + +ImproperStyle(sqdistharm,ImproperSQDistHarm) + +#else + +#ifndef LMP_IMPROPER_SQDISTHARM_H +#define LMP_IMPROPER_SQDISTHARM_H + +#include +#include "improper.h" + +namespace LAMMPS_NS { + +class ImproperSQDistHarm : public Improper { + public: + ImproperSQDistHarm(class LAMMPS *); + ~ImproperSQDistHarm(); + void compute(int, int); + void coeff(int, char **); + void write_restart(FILE *); + void read_restart(FILE *); + + private: + double *k,*chi; + + void allocate(); +}; + +} + +#endif +#endif + diff --git a/src/USER-YAFF/pair_lj_switch3_coulgauss_long.cpp b/src/USER-YAFF/pair_lj_switch3_coulgauss_long.cpp new file mode 100644 index 0000000000..b3a267807d --- /dev/null +++ b/src/USER-YAFF/pair_lj_switch3_coulgauss_long.cpp @@ -0,0 +1,714 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_lj_switch3_coulgauss_long.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairLJSwitch3CoulGaussLong::PairLJSwitch3CoulGaussLong(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; + respa_enable = 1; + writedata = 1; + ftable = NULL; + qdist = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +PairLJSwitch3CoulGaussLong::~PairLJSwitch3CoulGaussLong() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(gamma); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + memory->destroy(offset); + } + if (ftable) free_tables(); +} + +/* ---------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::compute(int eflag, int vflag) +{ + int i,ii,j,jj,inum,jnum,itype,jtype,itable,jtable,ktable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,ecoul2,fpair; + double fraction,fraction2,table; + double r,r2inv,r6inv,forcecoul,forcecoul2,forcelj,factor_coul,factor_lj,tr,ftr,trx; + double grij,expm2,prefactor,prefactor2,t,erfc1,erfc2,rrij,expn2,expb,g_ewald2i,g_ewaldi; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq, lookup_corr; + + 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; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + union_int_float_t rsq_lookup; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + g_ewaldi = 1.0/g_ewald; + g_ewald2i = g_ewaldi*g_ewaldi; + lookup_corr = 0.0; + + // 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]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc1 = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc1 + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + // Lennard-Jones potential + r = sqrt(rsq); + r6inv = r2inv*r2inv*r2inv; + forcelj = r6inv*(12.0*lj3[itype][jtype]*r6inv - 6.0*lj4[itype][jtype]); + // Correction for Gaussian radii + if (lj2[itype][jtype]==0.0) { + // This means a point charge is considerd, so the correction is zero + expn2 = 0.0; + erfc2 = 0.0; + forcecoul2 = 0.0; + } + else { + rrij = lj2[itype][jtype]*r; + expn2 = exp(-rrij*rrij); + erfc2 = erfc(rrij); + prefactor2 = -qqrd2e*qtmp*q[j]/r; + forcecoul2 = prefactor2*(erfc2+EWALD_F*rrij*expn2); + } + } else forcelj = 0.0; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*erfc1; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + ecoul += prefactor2*erfc2*factor_coul; + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + } else evdwl = 0.0; + + // Truncation, see Yaff Switch33 + if (truncw>0) { + if (rsq < cut_ljsq[itype][jtype]) { + if (r>cut_lj[itype][jtype]-truncw) { + trx = (cut_lj[itype][jtype]-r)*truncwi; + tr = trx*trx*(3.0-2.0*trx); + ftr = 6.0*trx*(1.0-trx)*r*truncwi; + forcelj = forcelj*tr + evdwl*ftr; + evdwl *= tr; + } + } + } + + fpair = (forcecoul + factor_coul*forcecoul2 + factor_lj*forcelj) * r2inv; + evdwl *= factor_lj; + + 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 (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); + memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); + memory->create(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(gamma,n+1,n+1,"pair:gamma"); + memory->create(lj1,n+1,n+1,"pair:lj1"); + memory->create(lj2,n+1,n+1,"pair:lj2"); + memory->create(lj3,n+1,n+1,"pair:lj3"); + memory->create(lj4,n+1,n+1,"pair:lj4"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::settings(int narg, char **arg) +{ + if (narg < 2 || narg > 3) error->all(FLERR,"Illegal pair_style command"); + + cut_lj_global = force->numeric(FLERR,arg[0]); + if (narg == 2) { + cut_coul = cut_lj_global; + truncw = force->numeric(FLERR,arg[1]); + } + else { + cut_coul = force->numeric(FLERR,arg[1]); + truncw = force->numeric(FLERR,arg[2]); + } + if (truncw>0.0) truncwi = 1.0/truncw; + else truncwi = 0.0; + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) + if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 6) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + double gamma_one = force->numeric(FLERR,arg[4]); + + double cut_lj_one = cut_lj_global; + if (narg == 6) cut_lj_one = force->numeric(FLERR,arg[5]); + + 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; + gamma[i][j] = gamma_one; + cut_lj[i][j] = cut_lj_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style lj/switch3/coulgauss/long requires atom attribute q"); + + // request regular or rRESPA neighbor list + + int irequest; + int respa = 0; + + if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } + + irequest = neighbor->request(this,instance_me); + + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; + + cut_coulsq = cut_coul * cut_coul; + + // set rRESPA cutoffs + + if (strstr(update->integrate_style,"respa") && + ((Respa *) update->integrate)->level_inner >= 0) + cut_respa = ((Respa *) update->integrate)->cutoff; + else cut_respa = NULL; + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); + g_ewald = force->kspace->g_ewald; + + // setup force tables + + if (ncoultablebits) init_tables(cut_coul,cut_respa); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJSwitch3CoulGaussLong::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]); + gamma[i][j] = 1.0/sqrt(gamma[i][i]*gamma[i][i]+gamma[j][j]*gamma[j][j]); + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + } + + double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + if (gamma[i][i]==0.0 && gamma[j][j]==0.0) lj2[i][j] = 0.0; + else lj2[i][j] = 1.0/sqrt(gamma[i][i]*gamma[i][i]+gamma[j][j]*gamma[j][j]); + 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); + + if (offset_flag && (cut_lj[i][j] > 0.0)) { + // Truncation is active, so offset is zero, except if truncw==0.0 + if (truncw==0.0) { + double r = cut_lj[i][j]; + double r2inv = 1.0/(r*r); + double r6inv = r2inv*r2inv*r2inv; + double r12inv = r6inv*r6inv; + offset[i][j] = lj3[i][j]*r12inv-lj4[i][j]*r6inv; + } + else {offset[i][j] = 0.0;} + } else offset[i][j] = 0.0; + + cut_ljsq[j][i] = cut_ljsq[i][j]; + cut_lj[j][i] = cut_lj[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]; + offset[j][i] = offset[i][j]; + + // check interior rRESPA cutoff + + if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) + error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + + // compute I,J contribution to long-range tail correction + // count total # of atoms of type I and J via Allreduce + + if (tail_flag) { + int *type = atom->type; + int nlocal = atom->nlocal; + + double count[2],all[2]; + count[0] = count[1] = 0.0; + for (int k = 0; k < nlocal; k++) { + if (type[k] == i) count[0] += 1.0; + if (type[k] == j) count[1] += 1.0; + } + MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + + double cg = epsilon[i][j]; + double cg1 = cut_lj[i][j]; + double cg3 = sigma[i][j]; + if (truncw > 0.0) { + double cg5 = truncw; + double t1 = pow(cg1, 0.2e1); + double t2 = t1 * cg1; + double t3 = t1 * t1; + double t4 = t3 * t2; + double t5 = -cg1 + cg5; + double t6 = t5 * t5; + double t8 = t6 * t6; + double t9 = t8 * t6 * t5; + double t10 = t4 * t9; + double t11 = log(-t5); + double t14 = log(cg1); + double t19 = pow(cg5, 0.2e1); + double t20 = t19 * t19; + double t21 = t20 * t19; + double t24 = pow(cg3, 0.2e1); + double t25 = t24 * t24; + double t26 = t25 * t24; + double t29 = t20 * cg5; + double t35 = t3 * t3; + double t41 = t19 * cg5; + double t58 = t21 * t3 * t1 - t21 * t26 / 0.84e2 - 0.6e1 * t29 * t4 + t29 * cg1 * t26 / 0.18e2 + 0.15e2 * t20 * t35 - t20 * t1 * t26 / 0.9e1 - 0.20e2 * t41 * t35 * cg1 + t41 * t2 * t26 / 0.9e1 + 0.15e2 * t19 * t35 * t1 - t19 * t3 * t26 / 0.18e2 - 0.6e1 * t35 * t2 * cg5 + t35 * t3; + double t71 = -0.4e1 * cg * (0.2e1 * t10 * t11 - 0.2e1 * t10 * t14 + (cg5 - 0.2e1 * cg1) * t58 * cg5) * t26 / t4 / t41 / t9; + etail_ij = 2.0*MY_PI*all[0]*all[1]*t71; + ptail_ij = 2.0*MY_PI*all[0]*all[1]*t71; + } + else { + double t1 = pow(cg3, 0.2e1); + double t2 = t1 * t1; + double t3 = t2 * t1; + double t5 = pow(cg1, 0.2e1); + double t6 = t5 * t5; + double t10 = t6 * t6; + double t16 = -0.4e1 / 0.9e1 * t3 * cg * (0.3e1 * t6 * t5 - t3) / t10 / cg1; + t1 = pow(cg3, 0.2e1); + t2 = t1 * t1; + t3 = t2 * t1; + t5 = pow(cg1, 0.2e1); + t6 = t5 * t5; + double t11 = t6 * t6; + double t17 = 0.8e1 / 0.3e1 * t3 * cg * (0.3e1 * t6 * t5 - 0.2e1 * t3) / t11 / cg1; + etail_ij = 2.0*MY_PI*all[0]*all[1]*t16; + ptail_ij = -2.0/3.0*MY_PI*all[0]*all[1]*t17; + } + } + + return cut; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::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(&gamma[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::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(&gamma[i][j],sizeof(double),1,fp); + fread(&cut_lj[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(&gamma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::write_restart_settings(FILE *fp) +{ + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,sizeof(double),1,fp); + fwrite(&truncw,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&ncoultablebits,sizeof(int),1,fp); + fwrite(&tabinner,sizeof(double),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,sizeof(double),1,fp); + fread(&truncw,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + fread(&tail_flag,sizeof(int),1,fp); + fread(&ncoultablebits,sizeof(int),1,fp); + fread(&tabinner,sizeof(double),1,fp); + } + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&truncw,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); + MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); +} + + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],gamma[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairLJSwitch3CoulGaussLong::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],gamma[i][j],cut_lj[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairLJSwitch3CoulGaussLong::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,grij,expm2,t,erfc1,prefactor,prefactor2,rrij,expn2,erfc2; + double fraction,table,forcecoul,forcecoul2,forcelj,phicoul,phicoul2,philj; + double expb, ecoul, evdwl, trx, tr, ftr; + + int itable; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc1 = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc1 + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup_single; + rsq_lookup_single.f = rsq; + itable = rsq_lookup_single.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r = sqrt(rsq); + r6inv = r2inv*r2inv*r2inv; + rrij = lj2[itype][jtype] * r; + if (rrij==0.0) { + expn2 = 0.0; + erfc2 = 0.0; + } + else { + expn2 = exp(-rrij*rrij); + erfc2 = erfc(rrij); + } + prefactor2 = -force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul2 = prefactor2 * (erfc2 + EWALD_F*rrij*expn2); + forcelj = expb*lj1[itype][jtype]*r-6.0*lj4[itype][jtype]*r6inv; + } else forcelj = 0.0; + + double eng = 0.0; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc1; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + ecoul += prefactor2*erfc2*factor_coul; + evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - + offset[itype][jtype]; + } else evdwl = 0.0; + + // Truncation, see Yaff Switch3 + if (truncw>0) { + if (rsq < cut_ljsq[itype][jtype]) { + if (r>cut_lj[itype][jtype]-truncw) { + trx = (cut_lj[itype][jtype]-r)*truncwi; + tr = trx*trx*(3.0-2.0*trx); + ftr = 6.0*trx*(1.0-trx)*r*truncwi; + forcelj = forcelj*tr + evdwl*ftr; + evdwl *= tr; + } + } + } + eng += evdwl*factor_lj; + fforce = (forcecoul + factor_coul*forcecoul2 + factor_lj*forcelj) * r2inv; + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairLJSwitch3CoulGaussLong::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"gamma") == 0) return (void *) gamma; + return NULL; +} diff --git a/src/USER-YAFF/pair_lj_switch3_coulgauss_long.h b/src/USER-YAFF/pair_lj_switch3_coulgauss_long.h new file mode 100644 index 0000000000..de6bc107f7 --- /dev/null +++ b/src/USER-YAFF/pair_lj_switch3_coulgauss_long.h @@ -0,0 +1,91 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(lj/switch3/coulgauss/long,PairLJSwitch3CoulGaussLong) + +#else + +#ifndef LMP_PAIR_LJ_SWITCH3_COULGAUSS_LONG_H +#define LMP_PAIR_LJ_SWITCH3_COULGAUSS_LONG_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairLJSwitch3CoulGaussLong : public Pair { + + public: + PairLJSwitch3CoulGaussLong(class LAMMPS *); + virtual ~PairLJSwitch3CoulGaussLong(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + virtual void init_style(); + virtual 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 *); + void write_data(FILE *); + void write_data_all(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + + virtual void *extract(const char *, int &); + + protected: + double cut_lj_global; + double truncw, truncwi; + double **cut_lj,**cut_ljsq; + double cut_coul,cut_coulsq; + double **epsilon,**sigma,**gamma; + double **lj1,**lj2,**lj3,**lj4,**offset; + double *cut_respa; + double qdist; // TIP4P distance from O site to negative charge + double g_ewald; + + virtual void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style lj/cut/coul/long requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair style requires a KSpace style + +No kspace style is defined. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +*/ diff --git a/src/USER-YAFF/pair_mm3_switch3_coulgauss_long.cpp b/src/USER-YAFF/pair_mm3_switch3_coulgauss_long.cpp new file mode 100644 index 0000000000..0d0e43a0ea --- /dev/null +++ b/src/USER-YAFF/pair_mm3_switch3_coulgauss_long.cpp @@ -0,0 +1,715 @@ +/* ---------------------------------------------------------------------- + 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: Steven Vandenbrande +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_mm3_switch3_coulgauss_long.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "kspace.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define EWALD_F 1.12837917 +#define EWALD_P 0.3275911 +#define A1 0.254829592 +#define A2 -0.284496736 +#define A3 1.421413741 +#define A4 -1.453152027 +#define A5 1.061405429 + +/* ---------------------------------------------------------------------- */ + +PairMM3Switch3CoulGaussLong::PairMM3Switch3CoulGaussLong(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; + respa_enable = 1; + writedata = 1; + ftable = NULL; + qdist = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +PairMM3Switch3CoulGaussLong::~PairMM3Switch3CoulGaussLong() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut_lj); + memory->destroy(cut_ljsq); + memory->destroy(epsilon); + memory->destroy(sigma); + memory->destroy(gamma); + memory->destroy(lj1); + memory->destroy(lj2); + memory->destroy(lj3); + memory->destroy(lj4); + memory->destroy(offset); + } + if (ftable) free_tables(); +} + +/* ---------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::compute(int eflag, int vflag) +{ + int i,ii,j,jj,inum,jnum,itype,jtype,itable,jtable,ktable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,ecoul2,fpair; + double fraction,fraction2,table; + double r,r2inv,r6inv,forcecoul,forcecoul2,forcelj,factor_coul,factor_lj,tr,ftr,trx; + double grij,expm2,prefactor,prefactor2,t,erfc1,erfc2,rrij,expn2,expb,g_ewald2i,g_ewaldi; + int *ilist,*jlist,*numneigh,**firstneigh; + double rsq, lookup_corr; + + 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; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + union_int_float_t rsq_lookup; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + g_ewaldi = 1.0/g_ewald; + g_ewald2i = g_ewaldi*g_ewaldi; + lookup_corr = 0.0; + + // 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]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + 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; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc1 = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc1 + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + // Repulsive exponential part + r = sqrt(rsq); + expb = lj3[itype][jtype]*exp(-lj1[itype][jtype]*r); + forcelj = expb*lj1[itype][jtype]*r; + // Attractive r^-6 part + r6inv = r2inv*r2inv*r2inv; + forcelj -= 6.0*lj4[itype][jtype]*r6inv; + // Correction for Gaussian radii + if (lj2[itype][jtype]==0.0) { + // This means a point charge is considered, so the correction is zero + expn2 = 0.0; + erfc2 = 0.0; + forcecoul2 = 0.0; + } + else { + rrij = lj2[itype][jtype]*r; + expn2 = exp(-rrij*rrij); + erfc2 = erfc(rrij); + prefactor2 = -qqrd2e*qtmp*q[j]/r; + forcecoul2 = prefactor2*(erfc2+EWALD_F*rrij*expn2); + } + } else forcelj = 0.0; + + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*erfc1; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } + if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; + } else ecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + ecoul += prefactor2*erfc2*factor_coul; + evdwl = expb-lj4[itype][jtype]*r6inv-offset[itype][jtype]; + } else evdwl = 0.0; + + // Truncation, see Yaff Switch3 + if (truncw>0) { + if (rsq < cut_ljsq[itype][jtype]) { + if (r>cut_lj[itype][jtype]-truncw) { + trx = (cut_lj[itype][jtype]-r)*truncwi; + tr = trx*trx*(3.0-2.0*trx); + ftr = 6.0*trx*(1.0-trx)*r*truncwi; + forcelj = forcelj*tr + evdwl*ftr; + evdwl *= tr; + } + } + } + + fpair = (forcecoul + factor_coul*forcecoul2 + factor_lj*forcelj) * r2inv; + evdwl *= factor_lj; + + 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 (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); + memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); + memory->create(epsilon,n+1,n+1,"pair:epsilon"); + memory->create(sigma,n+1,n+1,"pair:sigma"); + memory->create(gamma,n+1,n+1,"pair:gamma"); + memory->create(lj1,n+1,n+1,"pair:lj1"); + memory->create(lj2,n+1,n+1,"pair:lj2"); + memory->create(lj3,n+1,n+1,"pair:lj3"); + memory->create(lj4,n+1,n+1,"pair:lj4"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::settings(int narg, char **arg) +{ + if (narg < 2 || narg > 3) error->all(FLERR,"Illegal pair_style command"); + + cut_lj_global = force->numeric(FLERR,arg[0]); + if (narg == 2) { + cut_coul = cut_lj_global; + truncw = force->numeric(FLERR,arg[1]); + } + else { + cut_coul = force->numeric(FLERR,arg[1]); + truncw = force->numeric(FLERR,arg[2]); + } + if (truncw>0.0) truncwi = 1.0/truncw; + else truncwi = 0.0; + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) + if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 6) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + double epsilon_one = force->numeric(FLERR,arg[2]); + double sigma_one = force->numeric(FLERR,arg[3]); + double gamma_one = force->numeric(FLERR,arg[4]); + + double cut_lj_one = cut_lj_global; + if (narg == 6) cut_lj_one = force->numeric(FLERR,arg[5]); + + 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; + gamma[i][j] = gamma_one; + cut_lj[i][j] = cut_lj_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style mm3/switch3/coulgauss/long requires atom attribute q"); + + // request regular or rRESPA neighbor list + + int irequest; + int respa = 0; + + if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + } + + irequest = neighbor->request(this,instance_me); + + if (respa >= 1) { + neighbor->requests[irequest]->respaouter = 1; + neighbor->requests[irequest]->respainner = 1; + } + if (respa == 2) neighbor->requests[irequest]->respamiddle = 1; + + cut_coulsq = cut_coul * cut_coul; + + // set rRESPA cutoffs + + if (strstr(update->integrate_style,"respa") && + ((Respa *) update->integrate)->level_inner >= 0) + cut_respa = ((Respa *) update->integrate)->cutoff; + else cut_respa = NULL; + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); + g_ewald = force->kspace->g_ewald; + + // setup force tables + + if (ncoultablebits) init_tables(cut_coul,cut_respa); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairMM3Switch3CoulGaussLong::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + epsilon[i][j] = sqrt(epsilon[i][i]*epsilon[j][j]); + sigma[i][j] = sigma[i][i] + sigma[j][j]; + gamma[i][j] = 1.0/sqrt(gamma[i][i]*gamma[i][i]+gamma[j][j]*gamma[j][j]); + cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]); + } + + double cut = MAX(cut_lj[i][j],cut_coul+2.0*qdist); + cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; + lj1[i][j] = 12.0 / (sigma[i][j]); + if (gamma[i][i]==0.0 && gamma[j][j]==0.0) lj2[i][j] = 0.0; + else lj2[i][j] = 1.0/sqrt(gamma[i][i]*gamma[i][i]+gamma[j][j]*gamma[j][j]); + lj3[i][j] = 1.84e5 * epsilon[i][j]; + lj4[i][j] = 2.25 * epsilon[i][j] * pow(sigma[i][j],6.0); + + if (offset_flag && (cut_lj[i][j] > 0.0)) { + // Truncation is active, so offset is zero, except if truncw==0.0 + if (truncw==0.0) { + double r = cut_lj[i][j]; + double r2inv = 1.0/(r*r); + double r6inv = r2inv*r2inv*r2inv; + double expb = lj3[i][j]*exp(-lj1[i][j]*r); + offset[i][j] = expb-lj4[i][j]*r6inv; + } + else {offset[i][j] = 0.0;} + } else offset[i][j] = 0.0; + + cut_ljsq[j][i] = cut_ljsq[i][j]; + cut_lj[j][i] = cut_lj[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]; + offset[j][i] = offset[i][j]; + + // check interior rRESPA cutoff + + if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) + error->all(FLERR,"Pair cutoff < Respa interior cutoff"); + + // compute I,J contribution to long-range tail correction + // count total # of atoms of type I and J via Allreduce + + if (tail_flag) { + int *type = atom->type; + int nlocal = atom->nlocal; + + double count[2],all[2]; + count[0] = count[1] = 0.0; + for (int k = 0; k < nlocal; k++) { + if (type[k] == i) count[0] += 1.0; + if (type[k] == j) count[1] += 1.0; + } + MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + + double cg = epsilon[i][j]; + double cg1 = cut_lj[i][j]; + double cg3 = 2.0*sigma[i][j];//Mind the factor 2 here!!! + if (truncw > 0.0) { + double cg5 = truncw; + double t1 = pow(cg3, 0.2e1); + double t2 = t1 * cg3; + double t3 = 0.5e1 / 0.216e3 * t2; + double t5 = cg1 / 0.9e1; + double t8 = -cg1 + cg5; + double t14 = t8 * t8; + double t17 = 0.1e1 / cg3; + double t20 = exp(0.12e2 * t17 * cg5); + double t30 = pow(cg1, 0.2e1); + double t36 = exp(-0.12e2 * t17 * cg1); + double t37 = pow(cg5, 0.2e1); + double t39 = 0.1e1 / t37 / cg5; + double t43 = cg1 * t8; + double t44 = log(-t8); + double t47 = log(cg1); + double t54 = t1 * t1; + double t64 = cg * (0.6388888889e3 * ((-t3 + (0.7e1 / 0.36e2 * cg5 - t5) * t1 - 0.2e1 / 0.3e1 * t8 * (cg5 - cg1 / 0.4e1) * cg3 + cg5 * t14) * t20 + t3 + (cg5 / 0.12e2 + t5) * t1 + (cg5 + cg1 / 0.3e1) * cg1 * cg3 / 0.2e1 + t30 * cg5) * t2 * t36 * t39 - 0.225e1 * (0.2e1 * t43 * t44 - 0.2e1 * t43 * t47 + cg5 * (cg5 - 0.2e1 * cg1)) * t54 * t1 / cg1 / t8 * t39); + etail_ij = 2.0*MY_PI*all[0]*all[1]*t64; + ptail_ij = 2.0*MY_PI*all[0]*all[1]*t64; + } + else { + double t2 = pow(cg3, 0.2e1); + double t3 = t2 * t2; + double t7 = 0.12e2 / cg3 * cg1; + double t8 = exp(t7); + double t11 = pow(cg1, 0.2e1); + double t12 = t11 * t11; + double t17 = t11 * cg1; + double t21 = exp(-t7); + double t27 = -0.9259259259e-2 * cg3 * cg * (0.81e2 * t3 * cg3 * t8 - 0.1656000e7 * t12 * cg1 - 0.276000e6 * cg3 * t12 - 0.23000e5 * t2 * t17) * t21 / t17; + double t1 = pow(cg3, 0.2e1); + t2 = t1 * t1; + double t6 = 0.12e2 / cg3 * cg1; + t7 = exp(t6); + double t10 = pow(cg1, 0.2e1); + t11 = t10 * t10; + double t19 = t10 * cg1; + double t25 = exp(-t6); + double t29 = 0.5555555556e-1 * cg * (0.81e2 * t2 * t1 * t7 - 0.3312000e7 * t11 * t10 - 0.828000e6 * cg3 * t11 * cg1 - 0.138000e6 * t1 * t11 - 0.11500e5 * t19 * t1 * cg3) * t25 / t19; + etail_ij = 2.0*MY_PI*all[0]*all[1]*t27; + ptail_ij = -2.0/3.0*MY_PI*all[0]*all[1]*t29; + } + } + + return cut; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::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(&gamma[i][j],sizeof(double),1,fp); + fwrite(&cut_lj[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::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(&gamma[i][j],sizeof(double),1,fp); + fread(&cut_lj[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(&gamma[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::write_restart_settings(FILE *fp) +{ + fwrite(&cut_lj_global,sizeof(double),1,fp); + fwrite(&cut_coul,sizeof(double),1,fp); + fwrite(&truncw,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); + fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&ncoultablebits,sizeof(int),1,fp); + fwrite(&tabinner,sizeof(double),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_lj_global,sizeof(double),1,fp); + fread(&cut_coul,sizeof(double),1,fp); + fread(&truncw,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + fread(&tail_flag,sizeof(int),1,fp); + fread(&ncoultablebits,sizeof(int),1,fp); + fread(&tabinner,sizeof(double),1,fp); + } + printf("Reading from restart, trunc = %f\n",truncw); + MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&truncw,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); + MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); + MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); +} + + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],gamma[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairMM3Switch3CoulGaussLong::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],gamma[i][j],cut_lj[i][j]); +} + +/* ---------------------------------------------------------------------- */ + +double PairMM3Switch3CoulGaussLong::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r2inv,r6inv,r,grij,expm2,t,erfc1,prefactor,prefactor2,rrij,expn2,erfc2; + double fraction,table,forcecoul,forcecoul2,forcelj,phicoul,phicoul2,philj; + double expb, ecoul, evdwl, trx, tr, ftr; + + int itable; + + r2inv = 1.0/rsq; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc1 = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc1 + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup_single; + rsq_lookup_single.f = rsq; + itable = rsq_lookup_single.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = atom->q[i]*atom->q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + } else forcecoul = 0.0; + + if (rsq < cut_ljsq[itype][jtype]) { + r = sqrt(rsq); + r6inv = r2inv*r2inv*r2inv; + expb = lj3[itype][jtype]*exp(-lj1[itype][jtype]*r); + rrij = lj2[itype][jtype] * r; + if (rrij==0.0) { + expn2 = 0.0; + erfc2 = 0.0; + } + else { + expn2 = exp(-rrij*rrij); + erfc2 = erfc(rrij); + } + prefactor2 = -force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul2 = prefactor2 * (erfc2 + EWALD_F*rrij*expn2); + forcelj = expb*lj1[itype][jtype]*r-6.0*lj4[itype][jtype]*r6inv; + } else forcelj = 0.0; + + double eng = 0.0; + if (rsq < cut_coulsq) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc1; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + eng += phicoul; + } + + if (rsq < cut_ljsq[itype][jtype]) { + ecoul += prefactor2*erfc2*factor_coul; + evdwl = expb-lj4[itype][jtype]*r6inv-offset[itype][jtype]; + } else evdwl = 0.0; + + // Truncation, see Yaff Switch3 + if (truncw>0) { + if (rsq < cut_ljsq[itype][jtype]) { + if (r>cut_lj[itype][jtype]-truncw) { + trx = (cut_lj[itype][jtype]-r)*truncwi; + tr = trx*trx*(3.0-2.0*trx); + ftr = 6.0*trx*(1.0-trx)*r*truncwi; + forcelj = forcelj*tr + evdwl*ftr; + evdwl *= tr; + } + } + } + eng += evdwl*factor_lj; + fforce = (forcecoul + factor_coul*forcecoul2 + factor_lj*forcelj) * r2inv; + + return eng; +} + +/* ---------------------------------------------------------------------- */ + +void *PairMM3Switch3CoulGaussLong::extract(const char *str, int &dim) +{ + dim = 0; + if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul; + dim = 2; + if (strcmp(str,"epsilon") == 0) return (void *) epsilon; + if (strcmp(str,"sigma") == 0) return (void *) sigma; + if (strcmp(str,"gamma") == 0) return (void *) gamma; + return NULL; +} diff --git a/src/USER-YAFF/pair_mm3_switch3_coulgauss_long.h b/src/USER-YAFF/pair_mm3_switch3_coulgauss_long.h new file mode 100644 index 0000000000..76016c40de --- /dev/null +++ b/src/USER-YAFF/pair_mm3_switch3_coulgauss_long.h @@ -0,0 +1,91 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(mm3/switch3/coulgauss/long,PairMM3Switch3CoulGaussLong) + +#else + +#ifndef LMP_PAIR_MM3_SWITCH3_COULGAUSS_LONG_H +#define LMP_PAIR_MM3_SWITCH3_COULGAUSS_LONG_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairMM3Switch3CoulGaussLong : public Pair { + + public: + PairMM3Switch3CoulGaussLong(class LAMMPS *); + virtual ~PairMM3Switch3CoulGaussLong(); + virtual void compute(int, int); + virtual void settings(int, char **); + void coeff(int, char **); + virtual void init_style(); + virtual 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 *); + void write_data(FILE *); + void write_data_all(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); + + virtual void *extract(const char *, int &); + + protected: + double cut_lj_global; + double truncw, truncwi; + double **cut_lj,**cut_ljsq; + double cut_coul,cut_coulsq; + double **epsilon,**sigma,**gamma; + double **lj1,**lj2,**lj3,**lj4,**offset; + double *cut_respa; + double qdist; // TIP4P distance from O site to negative charge + double g_ewald; + + virtual void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style lj/cut/coul/long requires atom attribute q + +The atom style defined does not have this attribute. + +E: Pair style requires a KSpace style + +No kspace style is defined. + +E: Pair cutoff < Respa interior cutoff + +One or more pairwise cutoffs are too short to use with the specified +rRESPA cutoffs. + +*/