forked from lijiext/lammps
Implementation and documentation for USER-YAFF package
This commit is contained in:
Steven Vandenbrande
parent
bd3b712902
commit
794bf5084c
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 34 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 37 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 22 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 8.1 KiB |
|
|
@ -0,0 +1,9 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
\thispagestyle{empty}
|
||||
$$
|
||||
E = K (d - d_0)^2
|
||||
$$
|
||||
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 8.7 KiB |
|
|
@ -0,0 +1,9 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
\thispagestyle{empty}
|
||||
$$
|
||||
E = K (d^2 - d_0^2)^2
|
||||
$$
|
||||
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 17 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 20 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 40 KiB |
|
|
@ -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}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 25 KiB |
|
|
@ -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}
|
||||
|
|
@ -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
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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),
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -13,6 +13,7 @@ Bond Styles :h1
|
|||
bond_harmonic_shift
|
||||
bond_harmonic_shift_cut
|
||||
bond_hybrid
|
||||
bond_mm3
|
||||
bond_morse
|
||||
bond_none
|
||||
bond_nonlinear
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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 -->
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
||||
|
|
@ -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 <math.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <stdio.h>
|
||||
#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.
|
||||
|
||||
*/
|
||||
|
|
@ -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 <math.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <stdio.h>
|
||||
#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.
|
||||
|
||||
*/
|
||||
|
|
@ -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 <math.h>
|
||||
#include <stdlib.h>
|
||||
#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);
|
||||
}
|
||||
|
|
@ -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 <stdio.h>
|
||||
#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.
|
||||
|
||||
*/
|
||||
|
|
@ -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 <mpi.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <stdio.h>
|
||||
#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
|
||||
|
||||
|
|
@ -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 <mpi.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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 <stdio.h>
|
||||
#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
|
||||
|
||||
|
|
@ -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 <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
@ -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 <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#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;
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue