forked from lijiext/lammps
Merge branch 'master' into collected-small-changes
This commit is contained in:
Axel Kohlmeyer
commit
ef5c327f74
|
|
@ -178,7 +178,7 @@ set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE GRANULAR
|
|||
USER-MEAMC USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF
|
||||
USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-REAXC USER-SCAFACOS
|
||||
USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK
|
||||
USER-QUIP USER-QMMM)
|
||||
USER-QUIP USER-QMMM USER-YAFF)
|
||||
set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU)
|
||||
set(OTHER_PACKAGES CORESHELL QEQ)
|
||||
foreach(PKG ${DEFAULT_PACKAGES})
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
"sqdistharm"_improper_sqdistharm.html,
|
||||
"umbrella (o)"_improper_umbrella.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: 7.8 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: 9.0 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: 5.7 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: 2.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: 2.2 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: 4.5 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: 3.6 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: 11 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: 6.6 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}
|
||||
Binary file not shown.
|
|
@ -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(23), 8551-8566
|
||||
(1989),
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -36,10 +36,12 @@ react = mandatory argument indicating new reaction specification :l
|
|||
template-ID(post-reacted) = ID of a molecule template containing post-reaction topology :l
|
||||
map_file = name of file specifying corresponding atom-IDs in the pre- and post-reacted templates :l
|
||||
zero or more individual keyword/value pairs may be appended to each react argument :l
|
||||
individual_keyword = {prob} or {stabilize_steps} or {update_edges} :l
|
||||
individual_keyword = {prob} or {max_rxn} or {stabilize_steps} or {update_edges} :l
|
||||
{prob} values = fraction seed
|
||||
fraction = initiate reaction with this probability if otherwise eligible
|
||||
seed = random number seed (positive integer)
|
||||
{max_rxn} value = N
|
||||
N = maximum number of reactions allowed to occur
|
||||
{stabilize_steps} value = timesteps
|
||||
timesteps = number of timesteps to apply the internally-created "nve/limit"_fix_nve_limit.html fix to reacting atoms
|
||||
{update_edges} value = {none} or {charges} or {custom}
|
||||
|
|
@ -142,7 +144,7 @@ modified to match the post-reaction template.
|
|||
|
||||
A bonding atom pair will be identified if several conditions are met.
|
||||
First, a pair of atoms I,J within the specified react-group-ID of type
|
||||
itype and jtype must separated by a distance between {Rmin} and
|
||||
itype and jtype must be separated by a distance between {Rmin} and
|
||||
{Rmax}. It is possible that multiple bonding atom pairs are
|
||||
identified: if the bonding atoms in the pre-reacted template are not
|
||||
1-2, 1-3, or 1-4 neighbors, the closest bonding atom partner is set as
|
||||
|
|
@ -211,9 +213,10 @@ mandatory keyword is 'equivalences' and the optional keywords are
|
|||
N {equivalences} = # of atoms N in the reaction molecule templates
|
||||
N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template
|
||||
N {deleteIDs} = # of atoms N that are specified for deletion
|
||||
N {customIDs} = # of atoms N that are specified for a custom update :pre
|
||||
N {customIDs} = # of atoms N that are specified for a custom update
|
||||
N {constraints} = # of specified reaction constraints N :pre
|
||||
|
||||
The body of the map file contains two mandatory sections and three
|
||||
The body of the map file contains two mandatory sections and four
|
||||
optional sections. The first mandatory section begins with the keyword
|
||||
'BondingIDs' and lists the atom IDs of the bonding atom pair in the
|
||||
pre-reacted molecule template. The second mandatory section begins
|
||||
|
|
@ -230,7 +233,10 @@ Edges' and allows for forcing the update of a specific atom's atomic
|
|||
charge. The first column is the ID of an atom near the edge of the
|
||||
pre-reacted molecule template, and the value of the second column is
|
||||
either 'none' or 'charges.' Further details are provided in the
|
||||
discussion of the 'update_edges' keyword.
|
||||
discussion of the 'update_edges' keyword. The fourth optional section
|
||||
begins with the keyword 'Constraints' and lists additional criteria
|
||||
that must be satisfied in order for the reaction to occur. Currently,
|
||||
there is one type of constraint available, as discussed below.
|
||||
|
||||
A sample map file is given below:
|
||||
|
||||
|
|
@ -263,6 +269,18 @@ Equivalences :pre
|
|||
|
||||
:line
|
||||
|
||||
Any number of additional constraints may be specified in the
|
||||
Constraints section of the map file. Currently there is one type of
|
||||
additional constraint, of type 'distance', whose syntax is as follows:
|
||||
|
||||
distance {ID1} {ID2} {rmin} {rmax} :pre
|
||||
|
||||
where 'distance' is the required keyword, {ID1} and {ID2} are
|
||||
pre-reaction atom IDs, and these two atoms must be separated by a
|
||||
distance between {rmin} and {rmax} for the reaction to occur. This
|
||||
constraint can be used to enforce a certain orientation between
|
||||
reacting molecules.
|
||||
|
||||
Once a reaction site has been successfully identified, data structures
|
||||
within LAMMPS that store bond topology are updated to reflect the
|
||||
post-reacted molecule template. All force fields with fixed bonds,
|
||||
|
|
@ -285,7 +303,8 @@ The {prob} keyword can affect whether an eligible reaction actually
|
|||
occurs. The fraction setting must be a value between 0.0 and 1.0. A
|
||||
uniform random number between 0.0 and 1.0 is generated and the
|
||||
eligible reaction only occurs if the random number is less than the
|
||||
fraction.
|
||||
fraction. Up to N reactions are permitted to occur, as optionally
|
||||
specified by the {max_rxn} keyword.
|
||||
|
||||
The {stabilize_steps} keyword allows for the specification of how many
|
||||
timesteps a reaction site is stabilized before being returned to the
|
||||
|
|
|
|||
|
|
@ -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 -->
|
||||
|
|
|
|||
|
|
@ -596,6 +596,7 @@ pair_lj_long.html
|
|||
pair_lj_smooth.html
|
||||
pair_lj_smooth_linear.html
|
||||
pair_lj_soft.html
|
||||
pair_lj_switch3_coulgauss.html
|
||||
pair_lubricate.html
|
||||
pair_lubricateU.html
|
||||
pair_mdf.html
|
||||
|
|
@ -605,6 +606,7 @@ pair_meam_sw_spline.html
|
|||
pair_meso.html
|
||||
pair_mgpt.html
|
||||
pair_mie.html
|
||||
pair_mm3_switch3_coulgauss.html
|
||||
pair_momb.html
|
||||
pair_morse.html
|
||||
pair_multi_lucy.html
|
||||
|
|
@ -668,6 +670,7 @@ bond_harmonic_shift.html
|
|||
bond_harmonic_shift_cut.html
|
||||
bond_hybrid.html
|
||||
bond_morse.html
|
||||
bond_mm3.html
|
||||
bond_none.html
|
||||
bond_nonlinear.html
|
||||
bond_oxdna.html
|
||||
|
|
@ -687,11 +690,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 +730,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 +738,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 potential
|
||||
|
||||
: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
|
||||
|
|
|
|||
|
|
@ -60,6 +60,7 @@ Alejandre
|
|||
alessandro
|
||||
Alessandro
|
||||
aliceblue
|
||||
Allinger
|
||||
allocaters
|
||||
allosws
|
||||
AlO
|
||||
|
|
@ -410,6 +411,7 @@ cossq
|
|||
costheta
|
||||
Couette
|
||||
coul
|
||||
coulgauss
|
||||
coulombic
|
||||
Coulombic
|
||||
Coulombics
|
||||
|
|
@ -566,6 +568,7 @@ discretized
|
|||
disp
|
||||
dissipative
|
||||
Dissipative
|
||||
distharm
|
||||
dl
|
||||
dlambda
|
||||
DLAMMPS
|
||||
|
|
@ -2493,6 +2496,7 @@ Springer
|
|||
springgreen
|
||||
spx
|
||||
spz
|
||||
sqdistharm
|
||||
sqrt
|
||||
src
|
||||
srd
|
||||
|
|
@ -2801,6 +2805,8 @@ Valone
|
|||
valuev
|
||||
Valuev
|
||||
Vanden
|
||||
Vandenbrande
|
||||
Vanduyfhuys
|
||||
varavg
|
||||
Varshalovich
|
||||
Varshney
|
||||
|
|
@ -2958,6 +2964,8 @@ xy
|
|||
xyz
|
||||
xz
|
||||
xzhou
|
||||
YAFF
|
||||
yaff
|
||||
Yamada
|
||||
Yazdani
|
||||
Ybar
|
||||
|
|
@ -2976,6 +2984,7 @@ ys
|
|||
ysu
|
||||
yu
|
||||
Yu
|
||||
Yuh
|
||||
yukawa
|
||||
Yukawa
|
||||
Yusof
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
LAMMPS (30 Mar 2018)
|
||||
LAMMPS (4 Jan 2019)
|
||||
|
||||
units metal
|
||||
atom_style full
|
||||
|
|
@ -45,10 +45,10 @@ fix 1 all nph x 1. 1. 10.
|
|||
fix 2 all temp/csvr 350. 350. 0.1 64582
|
||||
|
||||
run 1000
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:138)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:139)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:139)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:139)
|
||||
WARNING: More than one compute entropy/atom (../compute_entropy_atom.cpp:139)
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
|
|
@ -56,27 +56,27 @@ Neighbor list info ...
|
|||
ghost atom cutoff = 13.2
|
||||
binsize = 6.6, bins = 21 6 6
|
||||
5 neighbor lists, perpetual/occasional/extra = 5 0 0
|
||||
(1) pair eam/fs, perpetual
|
||||
(1) pair eam/fs, perpetual, half/full from (2)
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/newton/tri
|
||||
stencil: half/bin/3d/newton/tri
|
||||
bin: standard
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) compute entropy/atom, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
(3) compute entropy/atom, perpetual, copy from (2)
|
||||
attributes: full, newton on
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(4) compute entropy/atom, perpetual
|
||||
attributes: full, newton on, ghost
|
||||
pair build: full/bin/ghost
|
||||
stencil: full/ghost/bin/3d
|
||||
bin: standard
|
||||
(3) compute entropy/atom, perpetual, copy from (2)
|
||||
attributes: full, newton on, ghost
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(4) compute entropy/atom, perpetual, copy from (2)
|
||||
attributes: full, newton on, ghost
|
||||
pair build: copy
|
||||
stencil: none
|
||||
bin: none
|
||||
(5) compute entropy/atom, perpetual, copy from (2)
|
||||
(5) compute entropy/atom, perpetual, copy from (4)
|
||||
attributes: full, newton on, ghost
|
||||
pair build: copy
|
||||
stencil: none
|
||||
|
|
@ -85,34 +85,34 @@ Setting up Verlet run ...
|
|||
Unit style : metal
|
||||
Current step : 0
|
||||
Time step : 0.002
|
||||
Per MPI rank memory allocation (min/avg/max) = 25.68 | 25.69 | 25.69 Mbytes
|
||||
Per MPI rank memory allocation (min/avg/max) = 19.6 | 19.6 | 19.6 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press Volume
|
||||
0 346.29871 -4285.222 0 -4101.9191 594.65353 165399.75
|
||||
500 359.33758 -4285.247 0 -4095.0423 471.98587 165847.18
|
||||
1000 348.99659 -4276.2274 0 -4091.4964 149.27188 166966.18
|
||||
Loop time of 5.3437 on 4 procs for 1000 steps with 4096 atoms
|
||||
500 359.33769 -4285.2472 0 -4095.0424 472.02043 165847.09
|
||||
1000 348.99683 -4276.2282 0 -4091.4971 149.38771 166965.86
|
||||
Loop time of 4.4394 on 4 procs for 1000 steps with 4096 atoms
|
||||
|
||||
Performance: 32.337 ns/day, 0.742 hours/ns, 187.136 timesteps/s
|
||||
99.8% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
Performance: 38.924 ns/day, 0.617 hours/ns, 225.256 timesteps/s
|
||||
99.5% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 4.2832 | 4.3257 | 4.3839 | 1.8 | 80.95
|
||||
Bond | 0.00018309 | 0.00019825 | 0.00021418 | 0.0 | 0.00
|
||||
Neigh | 0.42195 | 0.42512 | 0.42739 | 0.3 | 7.96
|
||||
Comm | 0.051679 | 0.1101 | 0.14916 | 10.8 | 2.06
|
||||
Output | 0.40909 | 0.4091 | 0.40911 | 0.0 | 7.66
|
||||
Modify | 0.060869 | 0.061921 | 0.06327 | 0.4 | 1.16
|
||||
Other | | 0.01161 | | | 0.22
|
||||
Pair | 3.5298 | 3.5931 | 3.6669 | 2.6 | 80.94
|
||||
Bond | 8.2099e-05 | 0.00011444 | 0.00014673 | 0.0 | 0.00
|
||||
Neigh | 0.43384 | 0.4445 | 0.45558 | 1.4 | 10.01
|
||||
Comm | 0.048904 | 0.11122 | 0.16728 | 12.7 | 2.51
|
||||
Output | 0.22595 | 0.22595 | 0.22596 | 0.0 | 5.09
|
||||
Modify | 0.053103 | 0.053795 | 0.054549 | 0.3 | 1.21
|
||||
Other | | 0.01068 | | | 0.24
|
||||
|
||||
Nlocal: 1024 ave 1040 max 1001 min
|
||||
Histogram: 1 0 0 0 0 0 2 0 0 1
|
||||
Nghost: 4614.25 ave 4700 max 4540 min
|
||||
Histogram: 1 1 0 0 0 0 0 1 0 1
|
||||
Neighs: 121747 ave 126398 max 116931 min
|
||||
Neighs: 121747 ave 126398 max 116930 min
|
||||
Histogram: 1 0 0 1 0 0 1 0 0 1
|
||||
FullNghs: 243494 ave 252523 max 233842 min
|
||||
FullNghs: 243494 ave 252523 max 233841 min
|
||||
Histogram: 1 0 0 1 0 0 1 0 0 1
|
||||
|
||||
Total # of neighbors = 973974
|
||||
|
|
@ -121,4 +121,4 @@ Ave special neighs/atom = 0
|
|||
Neighbor list builds = 13
|
||||
Dangerous builds = 0
|
||||
|
||||
Total wall time: 0:00:05
|
||||
Total wall time: 0:00:04
|
||||
|
|
|
|||
|
|
@ -0,0 +1,7 @@
|
|||
mof5:
|
||||
NPT simulation of MOF5 using a QuickFF force field
|
||||
|
||||
mil53al:
|
||||
NPT simulation of MIL-53(Al) using a QuickFF force field.
|
||||
If the pressure is high enough (for instance 2000 atm),
|
||||
a phase transition from large pore to narrow pore is observed
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,42 @@
|
|||
#########################################
|
||||
#General settings
|
||||
#########################################
|
||||
units real
|
||||
atom_style full
|
||||
boundary p p p
|
||||
dielectric 1
|
||||
|
||||
#########################################
|
||||
#Force field and system specification
|
||||
#########################################
|
||||
special_bonds lj 0.000000 0.000000 1.000000 coul 1.000000 1.000000 1.000000
|
||||
pair_style mm3/switch3/coulgauss/long 12.0000 4.0000
|
||||
pair_modify table 16 # Accuracy of the table used for real space electrostatics
|
||||
pair_modify mix arithmetic
|
||||
pair_modify tail no
|
||||
bond_style harmonic
|
||||
angle_style hybrid cosine/periodic cross harmonic cosine/squared
|
||||
dihedral_style fourier
|
||||
improper_style distharm
|
||||
box tilt large
|
||||
|
||||
read_data lammps.data # Data file location
|
||||
kspace_style pppm 0.0000001 # Ewald accuracy
|
||||
neighbor 2.0 multi
|
||||
neigh_modify every 2 delay 4 check yes
|
||||
|
||||
#########################################
|
||||
#Output settings
|
||||
#########################################
|
||||
thermo 100 # Provide output every n steps
|
||||
thermo_style custom step time etotal ke temp pe emol evdwl ecoul elong etail vol press
|
||||
thermo_modify line multi format float %20.12f
|
||||
|
||||
#########################################
|
||||
#Sampling options
|
||||
#########################################
|
||||
timestep 0.5 # in femtosecond
|
||||
velocity all create 0.0 5 # initial temperature in Kelvin and random seed
|
||||
fix 1 all npt temp 300.0 300.0 100.0 tri 2000.0 2000.0 1000.0 tchain 3 mtk yes nreset 1000
|
||||
fix_modify 1 energy yes # Add thermo/barostat contributions to energy
|
||||
run 10000
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,42 @@
|
|||
#########################################
|
||||
#General settings
|
||||
#########################################
|
||||
units real
|
||||
atom_style full
|
||||
boundary p p p
|
||||
dielectric 1
|
||||
|
||||
#########################################
|
||||
#Force field and system specification
|
||||
#########################################
|
||||
special_bonds lj 0.000000 0.000000 1.000000 coul 1.000000 1.000000 1.000000
|
||||
pair_style mm3/switch3/coulgauss/long 12.0000 4.0000
|
||||
pair_modify table 16 # Accuracy of the table used for real space electrostatics
|
||||
pair_modify mix arithmetic
|
||||
pair_modify tail no
|
||||
bond_style mm3
|
||||
angle_style hybrid cross mm3
|
||||
dihedral_style fourier
|
||||
improper_style distharm
|
||||
box tilt large
|
||||
|
||||
read_data lammps.data # Data file location
|
||||
kspace_style pppm 0.0000001 # Ewald accuracy
|
||||
neighbor 2.0 multi
|
||||
neigh_modify every 2 delay 4 check yes
|
||||
|
||||
#########################################
|
||||
#Output settings
|
||||
#########################################
|
||||
thermo 10 # Provide output every n steps
|
||||
thermo_style custom step time etotal ke temp pe emol evdwl ecoul elong etail vol press
|
||||
thermo_modify line multi format float %20.12f
|
||||
|
||||
#########################################
|
||||
#Sampling options
|
||||
#########################################
|
||||
timestep 0.5 # in femtosecond
|
||||
velocity all create 0.0 5 # initial temperature in Kelvin and random seed
|
||||
fix 1 all npt temp 300.0 300.0 100.0 tri 1.0 1.0 1000.0 tchain 3 mtk yes nreset 1000
|
||||
fix_modify 1 energy yes # Add thermo/barostat contributions to energy
|
||||
run 100
|
||||
|
|
@ -256,7 +256,6 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
|
|||
|
||||
if (vflag_fdotr) pair_virial_fdotr_compute(this);
|
||||
|
||||
|
||||
if (eflag_atom) {
|
||||
k_eatom.template modify<DeviceType>();
|
||||
k_eatom.template sync<LMPHostType>();
|
||||
|
|
@ -275,8 +274,8 @@ void PairSNAPKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
|
|||
|
||||
// free duplicated memory
|
||||
if (need_dup) {
|
||||
dup_f = decltype(dup_f)();
|
||||
dup_vatom = decltype(dup_vatom)();
|
||||
dup_f = decltype(dup_f)();
|
||||
dup_vatom = decltype(dup_vatom)();
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -453,6 +452,13 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAP<NEIGHFLAG,EVFLAG>,const
|
|||
//t4 += timer.seconds(); timer.reset();
|
||||
team.team_barrier();
|
||||
|
||||
if (quadraticflag) {
|
||||
my_sna.compute_bi(team);
|
||||
team.team_barrier();
|
||||
my_sna.copy_bi2bvec(team);
|
||||
team.team_barrier();
|
||||
}
|
||||
|
||||
// for neighbors of I within cutoff:
|
||||
// compute dUi/drj and dBi/drj
|
||||
// Fij = dEi/dRj = -dEi/dRi => add to Fi, subtract from Fj
|
||||
|
|
@ -472,10 +478,6 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAP<NEIGHFLAG,EVFLAG>,const
|
|||
my_sna.compute_dbidrj(team);
|
||||
//t7 += timer2.seconds(); timer2.reset();
|
||||
my_sna.copy_dbi2dbvec(team);
|
||||
if (quadraticflag) {
|
||||
my_sna.compute_bi(team);
|
||||
my_sna.copy_bi2bvec(team);
|
||||
}
|
||||
|
||||
Kokkos::single(Kokkos::PerThread(team), [&] (){
|
||||
F_FLOAT fij[3];
|
||||
|
|
@ -536,10 +538,10 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAP<NEIGHFLAG,EVFLAG>,const
|
|||
a_f(j,1) -= fij[1];
|
||||
a_f(j,2) -= fij[2];
|
||||
|
||||
// tally per-atom virial contribution
|
||||
// tally global and per-atom virial contribution
|
||||
|
||||
if (EVFLAG) {
|
||||
if (vflag) {
|
||||
if (vflag_either) {
|
||||
v_tally_xyz<NEIGHFLAG>(ev,i,j,
|
||||
fij[0],fij[1],fij[2],
|
||||
-my_sna.rij(jj,0),-my_sna.rij(jj,1),
|
||||
|
|
@ -554,11 +556,13 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAP<NEIGHFLAG,EVFLAG>,const
|
|||
// tally energy contribution
|
||||
|
||||
if (EVFLAG) {
|
||||
if (eflag) {
|
||||
if (eflag_either) {
|
||||
|
||||
if (!quadraticflag) {
|
||||
my_sna.compute_bi(team);
|
||||
team.team_barrier();
|
||||
my_sna.copy_bi2bvec(team);
|
||||
team.team_barrier();
|
||||
}
|
||||
|
||||
// E = beta.B + 0.5*B^t.alpha.B
|
||||
|
|
@ -566,14 +570,15 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAP<NEIGHFLAG,EVFLAG>,const
|
|||
// coeff[k] = alpha_ii or
|
||||
// coeff[k] = alpha_ij = alpha_ji, j != i
|
||||
|
||||
if (team.team_rank() == 0)
|
||||
Kokkos::single(Kokkos::PerThread(team), [&] () {
|
||||
Kokkos::single(Kokkos::PerTeam(team), [&] () {
|
||||
|
||||
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
|
||||
// evdwl = energy of atom I, sum over coeffs_k * Bi_k
|
||||
|
||||
double evdwl = d_coeffi[0];
|
||||
double evdwl = d_coeffi[0];
|
||||
|
||||
// linear contributions
|
||||
// could use thread vector range on this loop
|
||||
|
||||
// linear contributions
|
||||
for (int k = 1; k <= ncoeff; k++)
|
||||
evdwl += d_coeffi[k]*my_sna.bvec[k-1];
|
||||
|
||||
|
|
@ -589,11 +594,10 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAP<NEIGHFLAG,EVFLAG>,const
|
|||
}
|
||||
}
|
||||
}
|
||||
// ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);
|
||||
if (eflag_either) {
|
||||
if (eflag_global) ev.evdwl += evdwl;
|
||||
if (eflag_atom) d_eatom[i] += evdwl;
|
||||
}
|
||||
|
||||
//ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);
|
||||
if (eflag_global) ev.evdwl += evdwl;
|
||||
if (eflag_atom) d_eatom[i] += evdwl;
|
||||
});
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -327,29 +327,40 @@ void SNAKokkos<DeviceType>::compute_bi(const typename Kokkos::TeamPolicy<DeviceT
|
|||
clock_gettime(CLOCK_REALTIME, &starttime);
|
||||
#endif
|
||||
|
||||
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,idxj_max),
|
||||
Kokkos::parallel_for(Kokkos::TeamThreadRange(team,idxj_full_max),
|
||||
[&] (const int& idx) {
|
||||
const int j1 = idxj(idx).j1;
|
||||
const int j2 = idxj(idx).j2;
|
||||
const int j = idxj(idx).j;
|
||||
double b_j1_j2_j = 0.0;
|
||||
const int j1 = idxj_full(idx).j1;
|
||||
const int j2 = idxj_full(idx).j2;
|
||||
const int j = idxj_full(idx).j;
|
||||
|
||||
for(int mb = 0; 2*mb < j; mb++)
|
||||
for(int ma = 0; ma <= j; ma++) {
|
||||
b_j1_j2_j +=
|
||||
const int bound = (j+2)/2;
|
||||
double b_j1_j2_j = 0.0;
|
||||
double b_j1_j2_j_temp = 0.0;
|
||||
Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(team,(j+1)*bound),
|
||||
[&] (const int mbma, double& sum) {
|
||||
//for(int mb = 0; 2*mb <= j; mb++)
|
||||
//for(int ma = 0; ma <= j; ma++) {
|
||||
const int ma = mbma%(j+1);
|
||||
const int mb = mbma/(j+1);
|
||||
if (2*mb == j) return;
|
||||
sum +=
|
||||
uarraytot_r(j,ma,mb) * zarray_r(j1,j2,j,mb,ma) +
|
||||
uarraytot_i(j,ma,mb) * zarray_i(j1,j2,j,mb,ma);
|
||||
} // end loop over ma, mb
|
||||
},b_j1_j2_j_temp); // end loop over ma, mb
|
||||
b_j1_j2_j += b_j1_j2_j_temp;
|
||||
|
||||
// For j even, special treatment for middle column
|
||||
|
||||
if (j%2 == 0) {
|
||||
const int mb = j/2;
|
||||
for(int ma = 0; ma < mb; ma++) {
|
||||
b_j1_j2_j +=
|
||||
Kokkos::parallel_reduce(Kokkos::ThreadVectorRange(team, mb),
|
||||
[&] (const int ma, double& sum) {
|
||||
//for(int ma = 0; ma < mb; ma++) {
|
||||
sum +=
|
||||
uarraytot_r(j,ma,mb) * zarray_r(j1,j2,j,mb,ma) +
|
||||
uarraytot_i(j,ma,mb) * zarray_i(j1,j2,j,mb,ma);
|
||||
}
|
||||
},b_j1_j2_j_temp); // end loop over ma
|
||||
b_j1_j2_j += b_j1_j2_j_temp;
|
||||
|
||||
const int ma = mb;
|
||||
b_j1_j2_j +=
|
||||
|
|
@ -357,11 +368,13 @@ void SNAKokkos<DeviceType>::compute_bi(const typename Kokkos::TeamPolicy<DeviceT
|
|||
uarraytot_i(j,ma,mb) * zarray_i(j1,j2,j,mb,ma))*0.5;
|
||||
}
|
||||
|
||||
b_j1_j2_j *= 2.0;
|
||||
if (bzero_flag)
|
||||
b_j1_j2_j -= bzero[j];
|
||||
Kokkos::single(Kokkos::PerThread(team), [&] () {
|
||||
b_j1_j2_j *= 2.0;
|
||||
if (bzero_flag)
|
||||
b_j1_j2_j -= bzero[j];
|
||||
|
||||
barray(j1,j2,j) = b_j1_j2_j;
|
||||
barray(j1,j2,j) = b_j1_j2_j;
|
||||
});
|
||||
});
|
||||
//} // end loop over j
|
||||
//} // end loop over j1, j2
|
||||
|
|
@ -1028,6 +1041,8 @@ void SNAKokkos<DeviceType>::create_team_scratch_arrays(const typename Kokkos::Te
|
|||
uarraytot_i_a = uarraytot_i = t_sna_3d(team.team_scratch(1),jdim,jdim,jdim);
|
||||
zarray_r = t_sna_5d(team.team_scratch(1),jdim,jdim,jdim,jdim,jdim);
|
||||
zarray_i = t_sna_5d(team.team_scratch(1),jdim,jdim,jdim,jdim,jdim);
|
||||
bvec = Kokkos::View<double*, Kokkos::LayoutRight, DeviceType>(team.team_scratch(1),ncoeff);
|
||||
barray = t_sna_3d(team.team_scratch(1),jdim,jdim,jdim);
|
||||
|
||||
rij = t_sna_2d(team.team_scratch(1),nmax,3);
|
||||
rcutij = t_sna_1d(team.team_scratch(1),nmax);
|
||||
|
|
@ -1046,6 +1061,8 @@ T_INT SNAKokkos<DeviceType>::size_team_scratch_arrays() {
|
|||
size += t_sna_3d::shmem_size(jdim,jdim,jdim); // uarraytot_i_a
|
||||
size += t_sna_5d::shmem_size(jdim,jdim,jdim,jdim,jdim); // zarray_r
|
||||
size += t_sna_5d::shmem_size(jdim,jdim,jdim,jdim,jdim); // zarray_i
|
||||
size += Kokkos::View<double*, Kokkos::LayoutRight, DeviceType>::shmem_size(ncoeff); // bvec
|
||||
size += t_sna_3d::shmem_size(jdim,jdim,jdim); // barray
|
||||
|
||||
size += t_sna_2d::shmem_size(nmax,3); // rij
|
||||
size += t_sna_1d::shmem_size(nmax); // rcutij
|
||||
|
|
@ -1062,8 +1079,6 @@ KOKKOS_INLINE_FUNCTION
|
|||
void SNAKokkos<DeviceType>::create_thread_scratch_arrays(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team)
|
||||
{
|
||||
int jdim = twojmax + 1;
|
||||
bvec = Kokkos::View<double*, Kokkos::LayoutRight, DeviceType>(team.thread_scratch(1),ncoeff);
|
||||
barray = t_sna_3d(team.thread_scratch(1),jdim,jdim,jdim);
|
||||
|
||||
dbvec = Kokkos::View<double*[3], Kokkos::LayoutRight, DeviceType>(team.thread_scratch(1),ncoeff);
|
||||
dbarray = t_sna_4d(team.thread_scratch(1),jdim,jdim,jdim);
|
||||
|
|
@ -1079,8 +1094,6 @@ inline
|
|||
T_INT SNAKokkos<DeviceType>::size_thread_scratch_arrays() {
|
||||
T_INT size = 0;
|
||||
int jdim = twojmax + 1;
|
||||
size += Kokkos::View<double*, Kokkos::LayoutRight, DeviceType>::shmem_size(ncoeff); // bvec
|
||||
size += t_sna_3d::shmem_size(jdim,jdim,jdim); // barray
|
||||
|
||||
size += Kokkos::View<double*[3], Kokkos::LayoutRight, DeviceType>::shmem_size(ncoeff); // dbvec
|
||||
size += t_sna_4d::shmem_size(jdim,jdim,jdim); // dbarray
|
||||
|
|
|
|||
|
|
@ -189,17 +189,17 @@ void ImproperUmbrella::compute(int eflag, int vflag)
|
|||
dahy = ary-c*hry;
|
||||
dahz = arz-c*hrz;
|
||||
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
|
||||
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
|
||||
|
||||
f4[0] = dahx*rhr;
|
||||
f4[1] = dahy*rhr;
|
||||
f4[2] = dahz*rhr;
|
||||
f4[0] = dahx*rhr*a;
|
||||
f4[1] = dahy*rhr*a;
|
||||
f4[2] = dahz*rhr*a;
|
||||
|
||||
f1[0] = -(f2[0] + f3[0] + f4[0]);
|
||||
f1[1] = -(f2[1] + f3[1] + f4[1]);
|
||||
|
|
@ -208,27 +208,27 @@ void ImproperUmbrella::compute(int eflag, int vflag)
|
|||
// apply force to each of 4 atoms
|
||||
|
||||
if (newton_bond || i1 < nlocal) {
|
||||
f[i1][0] += f1[0]*a;
|
||||
f[i1][1] += f1[1]*a;
|
||||
f[i1][2] += f1[2]*a;
|
||||
f[i1][0] += f1[0];
|
||||
f[i1][1] += f1[1];
|
||||
f[i1][2] += f1[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i2 < nlocal) {
|
||||
f[i2][0] += f3[0]*a;
|
||||
f[i2][1] += f3[1]*a;
|
||||
f[i2][2] += f3[2]*a;
|
||||
f[i2][0] += f3[0];
|
||||
f[i2][1] += f3[1];
|
||||
f[i2][2] += f3[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i3 < nlocal) {
|
||||
f[i3][0] += f2[0]*a;
|
||||
f[i3][1] += f2[1]*a;
|
||||
f[i3][2] += f2[2]*a;
|
||||
f[i3][0] += f2[0];
|
||||
f[i3][1] += f2[1];
|
||||
f[i3][2] += f2[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i4 < nlocal) {
|
||||
f[i4][0] += f4[0]*a;
|
||||
f[i4][1] += f4[1]*a;
|
||||
f[i4][2] += f4[2]*a;
|
||||
f[i4][0] += f4[0];
|
||||
f[i4][1] += f4[1];
|
||||
f[i4][2] += f4[2];
|
||||
}
|
||||
|
||||
if (evflag) {
|
||||
|
|
@ -247,7 +247,7 @@ void ImproperUmbrella::compute(int eflag, int vflag)
|
|||
vb3y = x[i4][1] - x[i3][1];
|
||||
vb3z = x[i4][2] - x[i3][2];
|
||||
|
||||
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
|
||||
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,
|
||||
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -64,7 +64,7 @@ PACKUSER = user-atc user-awpmd user-bocs user-cgdna user-cgsdk user-colvars \
|
|||
user-mgpt user-misc user-mofff user-molfile \
|
||||
user-netcdf user-omp user-phonon user-plumed user-ptm user-qmmm \
|
||||
user-qtb user-quip user-reaxc user-scafacos user-smd user-smtbq \
|
||||
user-sdpd user-sph user-tally user-uef user-vtk
|
||||
user-sdpd user-sph user-tally user-uef user-vtk user-yaff
|
||||
|
||||
PACKLIB = compress gpu kim kokkos latte message mpiio mscg poems \
|
||||
python voronoi \
|
||||
|
|
|
|||
|
|
@ -74,10 +74,10 @@ pair_oxdna_coaxstk.cpp: coaxial stacking interaction between nucleotides
|
|||
|
||||
|
||||
pair_oxdna2_excv.cpp, pair_oxdna2_stk.cpp, pair_oxdna2_coaxstk.cpp:
|
||||
corresponding pair styles in oxDNA2 (see [3])
|
||||
corresponding pair styles in oxDNA2 (see [3])
|
||||
|
||||
pair_oxdna2_dh.cpp: Debye-Hueckel electrostatic interaction between backbone
|
||||
sites
|
||||
sites
|
||||
|
||||
** Fixes provided by this package:
|
||||
|
||||
|
|
|
|||
|
|
@ -683,7 +683,7 @@ void *PairLJCutTholeLong::extract(const char *str, int &dim)
|
|||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 6;
|
||||
dim = 2;
|
||||
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str,"scale") == 0) return (void *) scale;
|
||||
|
|
|
|||
|
|
@ -414,7 +414,7 @@ double PairThole::single(int i, int j, int itype, int jtype,
|
|||
|
||||
void *PairThole::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 4;
|
||||
dim = 2;
|
||||
if (strcmp(str,"scale") == 0) return (void *) scale;
|
||||
if (strcmp(str,"polar") == 0) return (void *) polar;
|
||||
if (strcmp(str,"thole") == 0) return (void *) thole;
|
||||
|
|
|
|||
|
|
@ -65,6 +65,7 @@ ComputeEntropyAtom(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (cutoff <= 0.0) error->all(FLERR,"Illegal compute entropy/atom"
|
||||
" command; cutoff must be positive");
|
||||
|
||||
cutoff2 = 0.;
|
||||
avg_flag = 0;
|
||||
local_flag = 0;
|
||||
|
||||
|
|
@ -137,15 +138,20 @@ void ComputeEntropyAtom::init()
|
|||
if (count > 1 && comm->me == 0)
|
||||
error->warning(FLERR,"More than one compute entropy/atom");
|
||||
|
||||
// need a full neighbor list with neighbors of the ghost atoms
|
||||
|
||||
// Request neighbor list
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
neighbor->requests[irequest]->occasional = 0;
|
||||
neighbor->requests[irequest]->ghost = 1;
|
||||
if (avg_flag) {
|
||||
// need a full neighbor list with neighbors of the ghost atoms
|
||||
neighbor->requests[irequest]->ghost = 1;
|
||||
} else {
|
||||
// need a full neighbor list
|
||||
neighbor->requests[irequest]->ghost = 0;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -58,7 +58,8 @@ static const char cite_fix_bond_react[] =
|
|||
#define BIG 1.0e20
|
||||
#define DELTA 16
|
||||
#define MAXLINE 256
|
||||
#define MAXGUESS 20
|
||||
#define MAXGUESS 20 // max # of guesses allowed by superimpose algorithm
|
||||
#define MAXCONARGS 5 // max # of arguments for any type of constraint
|
||||
|
||||
// various statuses of superimpose algorithm:
|
||||
// ACCEPT: site successfully matched to pre-reacted template
|
||||
|
|
@ -161,10 +162,15 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
|
|||
memory->create(unreacted_mol,nreacts,"bond/react:unreacted_mol");
|
||||
memory->create(reacted_mol,nreacts,"bond/react:reacted_mol");
|
||||
memory->create(fraction,nreacts,"bond/react:fraction");
|
||||
memory->create(max_rxn,nreacts,"bond/react:max_rxn");
|
||||
memory->create(nlocalskips,nreacts,"bond/react:nlocalskips");
|
||||
memory->create(nghostlyskips,nreacts,"bond/react:nghostlyskips");
|
||||
memory->create(seed,nreacts,"bond/react:seed");
|
||||
memory->create(limit_duration,nreacts,"bond/react:limit_duration");
|
||||
memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
|
||||
memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag");
|
||||
memory->create(nconstraints,nreacts,"bond/react:nconstraints");
|
||||
memory->create(constraints,nreacts,MAXCONARGS,"bond/react:constraints");
|
||||
memory->create(iatomtype,nreacts,"bond/react:iatomtype");
|
||||
memory->create(jatomtype,nreacts,"bond/react:jatomtype");
|
||||
memory->create(ibonding,nreacts,"bond/react:ibonding");
|
||||
|
|
@ -179,8 +185,10 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
|
|||
for (int i = 0; i < nreacts; i++) {
|
||||
fraction[i] = 1;
|
||||
seed[i] = 12345;
|
||||
max_rxn[i] = INT_MAX;
|
||||
stabilize_steps_flag[i] = 0;
|
||||
update_edges_flag[i] = 0;
|
||||
nconstraints[i] = 0;
|
||||
// set default limit duration to 60 timesteps
|
||||
limit_duration[i] = 60;
|
||||
reaction_count[i] = 0;
|
||||
|
|
@ -244,6 +252,13 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (seed[rxn] <= 0) error->all(FLERR,"Illegal fix bond/react command: "
|
||||
"probability seed must be positive");
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"max_rxn") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
|
||||
"'max_rxn' has too few arguments");
|
||||
max_rxn[rxn] = force->inumeric(FLERR,arg[iarg+1]);
|
||||
if (max_rxn[rxn] < 0) error->all(FLERR,"Illegal fix bond/react command: "
|
||||
"'max_rxn' cannot be negative");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"stabilize_steps") == 0) {
|
||||
if (stabilization_flag == 0) error->all(FLERR,"Stabilize_steps keyword "
|
||||
"used without stabilization keyword");
|
||||
|
|
@ -403,6 +418,9 @@ FixBondReact::~FixBondReact()
|
|||
memory->destroy(reacted_mol);
|
||||
memory->destroy(fraction);
|
||||
memory->destroy(seed);
|
||||
memory->destroy(max_rxn);
|
||||
memory->destroy(nlocalskips);
|
||||
memory->destroy(nghostlyskips);
|
||||
memory->destroy(limit_duration);
|
||||
memory->destroy(stabilize_steps_flag);
|
||||
memory->destroy(update_edges_flag);
|
||||
|
|
@ -674,6 +692,8 @@ void FixBondReact::post_integrate()
|
|||
reaction_count[i] = 0;
|
||||
local_rxn_count[i] = 0;
|
||||
ghostly_rxn_count[i] = 0;
|
||||
nlocalskips[i] = 0;
|
||||
nghostlyskips[i] = 0;
|
||||
}
|
||||
|
||||
if (nevery_check) {
|
||||
|
|
@ -745,6 +765,7 @@ void FixBondReact::post_integrate()
|
|||
|
||||
int j;
|
||||
for (rxnID = 0; rxnID < nreacts; rxnID++) {
|
||||
if (max_rxn[rxnID] <= reaction_count_total[rxnID]) continue;
|
||||
for (int ii = 0; ii < nall; ii++) {
|
||||
partner[ii] = 0;
|
||||
finalpartner[ii] = 0;
|
||||
|
|
@ -941,6 +962,7 @@ void FixBondReact::far_partner()
|
|||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
domain->minimum_image(delx,dely,delz); // ghost location fix
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
|
||||
if (rsq >= cutsq[rxnID][1] || rsq <= cutsq[rxnID][0]) {
|
||||
|
|
@ -1122,7 +1144,7 @@ void FixBondReact::superimpose_algorithm()
|
|||
}
|
||||
}
|
||||
|
||||
if (status == ACCEPT) { // reaction site found successfully!
|
||||
if (status == ACCEPT && check_constraints()) { // reaction site found successfully!
|
||||
glove_ghostcheck();
|
||||
}
|
||||
hang_catch++;
|
||||
|
|
@ -1142,14 +1164,43 @@ void FixBondReact::superimpose_algorithm()
|
|||
|
||||
MPI_Allreduce(&local_rxn_count[0],&reaction_count[0],nreacts,MPI_INT,MPI_SUM,world);
|
||||
|
||||
for (int i = 0; i < nreacts; i++)
|
||||
reaction_count_total[i] += reaction_count[i];
|
||||
|
||||
// this assumes compute_vector is called from process 0
|
||||
// ...so doesn't bother to bcast ghostly_rxn_count
|
||||
if (me == 0)
|
||||
for (int i = 0; i < nreacts; i++)
|
||||
reaction_count_total[i] += ghostly_rxn_count[i];
|
||||
reaction_count_total[i] += reaction_count[i] + ghostly_rxn_count[i];
|
||||
|
||||
MPI_Bcast(&reaction_count_total[0], nreacts, MPI_INT, 0, world);
|
||||
|
||||
// check if we overstepped our reaction limit
|
||||
for (int i = 0; i < nreacts; i++) {
|
||||
if (reaction_count_total[i] > max_rxn[i]) {
|
||||
// let's randomly choose rxns to skip, unbiasedly from local and ghostly
|
||||
int local_rxncounts[nprocs];
|
||||
int all_localskips[nprocs];
|
||||
MPI_Gather(&local_rxn_count[i],1,MPI_INT,local_rxncounts,1,MPI_INT,0,world);
|
||||
if (me == 0) {
|
||||
int overstep = reaction_count_total[i] - max_rxn[i];
|
||||
int delta_rxn = reaction_count[i] + ghostly_rxn_count[i];
|
||||
int rxn_by_proc[delta_rxn];
|
||||
for (int j = 0; j < delta_rxn; j++)
|
||||
rxn_by_proc[j] = -1; // corresponds to ghostly
|
||||
int itemp = 0;
|
||||
for (int j = 0; j < nprocs; j++)
|
||||
for (int k = 0; k < local_rxn_count[j]; k++)
|
||||
rxn_by_proc[itemp++] = j;
|
||||
std::random_shuffle(&rxn_by_proc[0],&rxn_by_proc[delta_rxn]);
|
||||
for (int j = 0; j < nprocs; j++)
|
||||
all_localskips[j] = 0;
|
||||
nghostlyskips[i] = 0;
|
||||
for (int j = 0; j < overstep; j++) {
|
||||
if (rxn_by_proc[j] == -1) nghostlyskips[i]++;
|
||||
else all_localskips[rxn_by_proc[j]]++;
|
||||
}
|
||||
}
|
||||
reaction_count_total[i] = max_rxn[i];
|
||||
MPI_Scatter(&all_localskips[0],1,MPI_INT,&nlocalskips[i],1,MPI_INT,0,world);
|
||||
MPI_Bcast(&nghostlyskips[i],1,MPI_INT,0,world);
|
||||
}
|
||||
}
|
||||
|
||||
// this updates topology next step
|
||||
next_reneighbor = update->ntimestep;
|
||||
|
|
@ -1526,6 +1577,32 @@ void FixBondReact::ring_check()
|
|||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
evaluate constraints: return 0 if any aren't satisfied
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int FixBondReact::check_constraints()
|
||||
{
|
||||
tagint atom1,atom2;
|
||||
double delx,dely,delz,rsq;
|
||||
|
||||
double **x = atom->x;
|
||||
|
||||
for (int i = 0; i < nconstraints[rxnID]; i++) {
|
||||
if (constraints[rxnID][0] == 0) { // 'distance' type
|
||||
atom1 = atom->map(glove[(int) constraints[rxnID][1]-1][1]);
|
||||
atom2 = atom->map(glove[(int) constraints[rxnID][2]-1][1]);
|
||||
delx = x[atom1][0] - x[atom2][0];
|
||||
dely = x[atom1][1] - x[atom2][1];
|
||||
delz = x[atom1][2] - x[atom2][2];
|
||||
domain->minimum_image(delx,dely,delz); // ghost location fix
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
if (rsq < constraints[rxnID][3] || rsq > constraints[rxnID][4]) return 0;
|
||||
}
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Get xspecials for current molecule templates
|
||||
------------------------------------------------------------------------- */
|
||||
|
|
@ -1945,19 +2022,19 @@ void FixBondReact::glove_ghostcheck()
|
|||
// 'ghosts of another' indication taken from comm->sendlist
|
||||
|
||||
int ghostly = 0;
|
||||
if (comm->style == 0) {
|
||||
for (int i = 0; i < onemol->natoms; i++) {
|
||||
int ilocal = atom->map(glove[i][1]);
|
||||
if (ilocal >= atom->nlocal || localsendlist[ilocal] == 1) {
|
||||
ghostly = 1;
|
||||
break;
|
||||
#if !defined(MPI_STUBS)
|
||||
if (comm->style == 0) {
|
||||
for (int i = 0; i < onemol->natoms; i++) {
|
||||
int ilocal = atom->map(glove[i][1]);
|
||||
if (ilocal >= atom->nlocal || localsendlist[ilocal] == 1) {
|
||||
ghostly = 1;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
#if !defined(MPI_STUBS)
|
||||
} else {
|
||||
ghostly = 1;
|
||||
#endif
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
if (ghostly == 1) {
|
||||
ghostly_mega_glove[0][ghostly_num_mega] = rxnID;
|
||||
|
|
@ -2092,18 +2169,26 @@ void FixBondReact::update_everything()
|
|||
memory->create(update_mega_glove,max_natoms+1,MAX(local_num_mega,global_megasize),"bond/react:update_mega_glove");
|
||||
|
||||
for (int pass = 0; pass < 2; pass++) {
|
||||
|
||||
update_num_mega = 0;
|
||||
int iskip[nreacts];
|
||||
for (int i = 0; i < nreacts; i++) iskip[i] = 0;
|
||||
if (pass == 0) {
|
||||
update_num_mega = local_num_mega;
|
||||
for (int i = 0; i < update_num_mega; i++) {
|
||||
for (int i = 0; i < local_num_mega; i++) {
|
||||
rxnID = local_mega_glove[0][i];
|
||||
// reactions already shuffled from dedup procedure, so can skip first N
|
||||
if (iskip[rxnID]++ < nlocalskips[rxnID]) continue;
|
||||
for (int j = 0; j < max_natoms+1; j++)
|
||||
update_mega_glove[j][i] = local_mega_glove[j][i];
|
||||
update_mega_glove[j][update_num_mega] = local_mega_glove[j][i];
|
||||
update_num_mega++;
|
||||
}
|
||||
} else if (pass == 1) {
|
||||
update_num_mega = global_megasize;
|
||||
for (int i = 0; i < global_megasize; i++) {
|
||||
rxnID = global_mega_glove[0][i];
|
||||
// reactions already shuffled from dedup procedure, so can skip first N
|
||||
if (iskip[rxnID]++ < nghostlyskips[rxnID]) continue;
|
||||
for (int j = 0; j < max_natoms+1; j++)
|
||||
update_mega_glove[j][i] = global_mega_glove[j][i];
|
||||
update_mega_glove[j][update_num_mega] = global_mega_glove[j][i];
|
||||
update_num_mega++;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2623,6 +2708,7 @@ void FixBondReact::read(int myrxn)
|
|||
else if (strstr(line,"equivalences")) sscanf(line,"%d",&nequivalent);
|
||||
else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
|
||||
else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
|
||||
else if (strstr(line,"constraints")) sscanf(line,"%d",&nconstraints[myrxn]);
|
||||
else break;
|
||||
}
|
||||
|
||||
|
|
@ -2654,6 +2740,8 @@ void FixBondReact::read(int myrxn)
|
|||
CustomEdges(line, myrxn);
|
||||
} else if (strcmp(keyword,"DeleteIDs") == 0) {
|
||||
DeleteAtoms(line, myrxn);
|
||||
} else if (strcmp(keyword,"Constraints") == 0) {
|
||||
Constraints(line, myrxn);
|
||||
} else error->one(FLERR,"Unknown section in superimpose file");
|
||||
|
||||
parse_keyword(1,line,keyword);
|
||||
|
|
@ -2729,6 +2817,27 @@ void FixBondReact::DeleteAtoms(char *line, int myrxn)
|
|||
}
|
||||
}
|
||||
|
||||
void FixBondReact::Constraints(char *line, int myrxn)
|
||||
{
|
||||
double tmp[MAXCONARGS];
|
||||
int n = strlen("distance") + 1;
|
||||
char *constraint_type = new char[n];
|
||||
for (int i = 0; i < nconstraints[myrxn]; i++) {
|
||||
readline(line);
|
||||
sscanf(line,"%s",constraint_type);
|
||||
if (strcmp(constraint_type,"distance") == 0) {
|
||||
constraints[myrxn][0] = 0; // 0 = 'distance' ...maybe use another enum eventually
|
||||
sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
|
||||
constraints[myrxn][1] = tmp[0];
|
||||
constraints[myrxn][2] = tmp[1];
|
||||
constraints[myrxn][3] = tmp[2]*tmp[2]; // using square of distance
|
||||
constraints[myrxn][4] = tmp[3]*tmp[3];
|
||||
} else
|
||||
error->one(FLERR,"Illegal constraint type in 'Constraints' section of map file");
|
||||
}
|
||||
delete [] constraint_type;
|
||||
}
|
||||
|
||||
void FixBondReact::open(char *file)
|
||||
{
|
||||
fp = fopen(file,"r");
|
||||
|
|
|
|||
|
|
@ -55,11 +55,14 @@ class FixBondReact : public Fix {
|
|||
int *iatomtype,*jatomtype;
|
||||
int *seed;
|
||||
double **cutsq,*fraction;
|
||||
int *max_rxn,*nlocalskips,*nghostlyskips;
|
||||
tagint lastcheck;
|
||||
int stabilization_flag;
|
||||
int custom_exclude_flag;
|
||||
int *stabilize_steps_flag;
|
||||
int *update_edges_flag;
|
||||
int *nconstraints;
|
||||
double **constraints;
|
||||
int status;
|
||||
int *groupbits;
|
||||
|
||||
|
|
@ -140,6 +143,7 @@ class FixBondReact : public Fix {
|
|||
void Equivalences(char *,int);
|
||||
void CustomEdges(char *,int);
|
||||
void DeleteAtoms(char *,int);
|
||||
void Constraints(char *,int);
|
||||
|
||||
void make_a_guess ();
|
||||
void neighbor_loop();
|
||||
|
|
@ -147,6 +151,7 @@ class FixBondReact : public Fix {
|
|||
void crosscheck_the_neighbor();
|
||||
void inner_crosscheck_loop();
|
||||
void ring_check();
|
||||
int check_constraints();
|
||||
|
||||
void open(char *);
|
||||
void readline(char *);
|
||||
|
|
|
|||
|
|
@ -206,17 +206,17 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
|
|||
dahy = ary-c*hry;
|
||||
dahz = arz-c*hrz;
|
||||
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
|
||||
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
|
||||
|
||||
f4[0] = dahx*rhr;
|
||||
f4[1] = dahy*rhr;
|
||||
f4[2] = dahz*rhr;
|
||||
f4[0] = dahx*rhr*a;
|
||||
f4[1] = dahy*rhr*a;
|
||||
f4[2] = dahz*rhr*a;
|
||||
|
||||
f1[0] = -(f2[0] + f3[0] + f4[0]);
|
||||
f1[1] = -(f2[1] + f3[1] + f4[1]);
|
||||
|
|
@ -225,32 +225,32 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
|
|||
// apply force to each of 4 atoms
|
||||
|
||||
if (newton_bond || i1 < nlocal) {
|
||||
f[i1][0] += f1[0]*a;
|
||||
f[i1][1] += f1[1]*a;
|
||||
f[i1][2] += f1[2]*a;
|
||||
f[i1][0] += f1[0];
|
||||
f[i1][1] += f1[1];
|
||||
f[i1][2] += f1[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i2 < nlocal) {
|
||||
f[i2][0] += f3[0]*a;
|
||||
f[i2][1] += f3[1]*a;
|
||||
f[i2][2] += f3[2]*a;
|
||||
f[i2][0] += f3[0];
|
||||
f[i2][1] += f3[1];
|
||||
f[i2][2] += f3[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i3 < nlocal) {
|
||||
f[i3][0] += f2[0]*a;
|
||||
f[i3][1] += f2[1]*a;
|
||||
f[i3][2] += f2[2]*a;
|
||||
f[i3][0] += f2[0];
|
||||
f[i3][1] += f2[1];
|
||||
f[i3][2] += f2[2];
|
||||
}
|
||||
|
||||
if (newton_bond || i4 < nlocal) {
|
||||
f[i4][0] += f4[0]*a;
|
||||
f[i4][1] += f4[1]*a;
|
||||
f[i4][2] += f4[2]*a;
|
||||
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,f1,f3,f4,
|
||||
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
|
||||
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,
|
||||
-vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
|
|||
|
|
@ -281,6 +281,9 @@ void DumpNetCDF::openfile()
|
|||
|
||||
// get total number of atoms
|
||||
ntotalgr = group->count(igroup);
|
||||
for (int i = 0; i < DUMP_NC_MAX_DIMS; i++) {
|
||||
vector_dim[i] = -1;
|
||||
}
|
||||
|
||||
if (filewriter) {
|
||||
if (append_flag && !multifile && access(filecurrent, F_OK) != -1) {
|
||||
|
|
@ -294,9 +297,6 @@ void DumpNetCDF::openfile()
|
|||
|
||||
// dimensions
|
||||
NCERRX( nc_inq_dimid(ncid, NC_FRAME_STR, &frame_dim), NC_FRAME_STR );
|
||||
NCERRX( nc_inq_dimid(ncid, NC_SPATIAL_STR, &spatial_dim),
|
||||
NC_SPATIAL_STR );
|
||||
NCERRX( nc_inq_dimid(ncid, NC_VOIGT_STR, &Voigt_dim), NC_VOIGT_STR );
|
||||
NCERRX( nc_inq_dimid(ncid, NC_ATOM_STR, &atom_dim), NC_ATOM_STR );
|
||||
NCERRX( nc_inq_dimid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_dim),
|
||||
NC_CELL_SPATIAL_STR );
|
||||
|
|
@ -304,6 +304,26 @@ void DumpNetCDF::openfile()
|
|||
NC_CELL_ANGULAR_STR );
|
||||
NCERRX( nc_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR );
|
||||
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
int dims = perat[i].dims;
|
||||
if (vector_dim[dims] < 0) {
|
||||
char dimstr[1024];
|
||||
if (dims == 3) {
|
||||
strcpy(dimstr, NC_SPATIAL_STR);
|
||||
}
|
||||
else if (dims == 6) {
|
||||
strcpy(dimstr, NC_VOIGT_STR);
|
||||
}
|
||||
else {
|
||||
sprintf(dimstr, "vec%i", dims);
|
||||
}
|
||||
if (dims != 1) {
|
||||
NCERRX( nc_inq_dimid(ncid, dimstr, &vector_dim[dims]),
|
||||
dimstr );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// default variables
|
||||
NCERRX( nc_inq_varid(ncid, NC_SPATIAL_STR, &spatial_var),
|
||||
NC_SPATIAL_STR );
|
||||
|
|
@ -320,7 +340,6 @@ void DumpNetCDF::openfile()
|
|||
NCERRX( nc_inq_varid(ncid, NC_CELL_ANGLES_STR, &cell_angles_var),
|
||||
NC_CELL_ANGLES_STR);
|
||||
|
||||
// variables specified in the input file
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
NCERRX( nc_inq_varid(ncid, perat[i].name, &perat[i].var),
|
||||
perat[i].name );
|
||||
|
|
@ -359,10 +378,6 @@ void DumpNetCDF::openfile()
|
|||
// dimensions
|
||||
NCERRX( nc_def_dim(ncid, NC_FRAME_STR, NC_UNLIMITED, &frame_dim),
|
||||
NC_FRAME_STR );
|
||||
NCERRX( nc_def_dim(ncid, NC_SPATIAL_STR, 3, &spatial_dim),
|
||||
NC_SPATIAL_STR );
|
||||
NCERRX( nc_def_dim(ncid, NC_VOIGT_STR, 6, &Voigt_dim),
|
||||
NC_VOIGT_STR );
|
||||
NCERRX( nc_def_dim(ncid, NC_ATOM_STR, ntotalgr, &atom_dim),
|
||||
NC_ATOM_STR );
|
||||
NCERRX( nc_def_dim(ncid, NC_CELL_SPATIAL_STR, 3, &cell_spatial_dim),
|
||||
|
|
@ -372,13 +387,33 @@ void DumpNetCDF::openfile()
|
|||
NCERRX( nc_def_dim(ncid, NC_LABEL_STR, 10, &label_dim),
|
||||
NC_LABEL_STR );
|
||||
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
int dims = perat[i].dims;
|
||||
if (vector_dim[dims] < 0) {
|
||||
char dimstr[1024];
|
||||
if (dims == 3) {
|
||||
strcpy(dimstr, NC_SPATIAL_STR);
|
||||
}
|
||||
else if (dims == 6) {
|
||||
strcpy(dimstr, NC_VOIGT_STR);
|
||||
}
|
||||
else {
|
||||
sprintf(dimstr, "vec%i", dims);
|
||||
}
|
||||
if (dims != 1) {
|
||||
NCERRX( nc_def_dim(ncid, dimstr, dims, &vector_dim[dims]),
|
||||
dimstr );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// default variables
|
||||
dims[0] = spatial_dim;
|
||||
dims[0] = vector_dim[3];
|
||||
NCERRX( nc_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var),
|
||||
NC_SPATIAL_STR );
|
||||
NCERRX( nc_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims,
|
||||
&cell_spatial_var), NC_CELL_SPATIAL_STR );
|
||||
dims[0] = spatial_dim;
|
||||
dims[0] = vector_dim[3];
|
||||
dims[1] = label_dim;
|
||||
NCERRX( nc_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims,
|
||||
&cell_angular_var), NC_CELL_ANGULAR_STR );
|
||||
|
|
@ -400,7 +435,7 @@ void DumpNetCDF::openfile()
|
|||
// variables specified in the input file
|
||||
dims[0] = frame_dim;
|
||||
dims[1] = atom_dim;
|
||||
dims[2] = spatial_dim;
|
||||
dims[2] = vector_dim[3];
|
||||
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
nc_type xtype;
|
||||
|
|
@ -419,53 +454,27 @@ void DumpNetCDF::openfile()
|
|||
|
||||
if (perat[i].constant) {
|
||||
// this quantity will only be written once
|
||||
if (perat[i].dims == 6) {
|
||||
// this is a tensor in Voigt notation
|
||||
dims[2] = Voigt_dim;
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else if (perat[i].dims == 3) {
|
||||
// this is a vector, we need to store x-, y- and z-coordinates
|
||||
dims[2] = spatial_dim;
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else if (perat[i].dims == 1) {
|
||||
if (perat[i].dims == 1) {
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 1, dims+1,
|
||||
&perat[i].var), perat[i].name );
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else {
|
||||
char errstr[1024];
|
||||
sprintf(errstr, "%i dimensions for '%s'. Not sure how to write "
|
||||
"this to the NetCDF trajectory file.", perat[i].dims,
|
||||
perat[i].name);
|
||||
error->all(FLERR,errstr);
|
||||
// this is a vector
|
||||
dims[1] = vector_dim[perat[i].dims];
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims+1,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
}
|
||||
else {
|
||||
if (perat[i].dims == 6) {
|
||||
// this is a tensor in Voigt notation
|
||||
dims[2] = Voigt_dim;
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else if (perat[i].dims == 3) {
|
||||
// this is a vector, we need to store x-, y- and z-coordinates
|
||||
dims[2] = spatial_dim;
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else if (perat[i].dims == 1) {
|
||||
if (perat[i].dims == 1) {
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 2, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else {
|
||||
char errstr[1024];
|
||||
sprintf(errstr, "%i dimensions for '%s'. Not sure how to write "
|
||||
"this to the NetCDF trajectory file.", perat[i].dims,
|
||||
perat[i].name);
|
||||
error->all(FLERR,errstr);
|
||||
// this is a vector
|
||||
dims[2] = vector_dim[perat[i].dims];
|
||||
NCERRX( nc_def_var(ncid, perat[i].name, xtype, 3, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -72,8 +72,7 @@ class DumpNetCDF : public DumpCustom {
|
|||
int ncid;
|
||||
|
||||
int frame_dim;
|
||||
int spatial_dim;
|
||||
int Voigt_dim;
|
||||
int vector_dim[DUMP_NC_MAX_DIMS];
|
||||
int atom_dim;
|
||||
int cell_spatial_dim;
|
||||
int cell_angular_dim;
|
||||
|
|
|
|||
|
|
@ -278,6 +278,9 @@ void DumpNetCDFMPIIO::openfile()
|
|||
|
||||
// get total number of atoms
|
||||
ntotalgr = group->count(igroup);
|
||||
for (int i = 0; i < DUMP_NC_MPIIO_MAX_DIMS; i++) {
|
||||
vector_dim[i] = -1;
|
||||
}
|
||||
|
||||
if (append_flag && !multifile && access(filecurrent, F_OK) != -1) {
|
||||
// Fixme! Perform checks if dimensions and variables conform with
|
||||
|
|
@ -294,9 +297,6 @@ void DumpNetCDFMPIIO::openfile()
|
|||
|
||||
// dimensions
|
||||
NCERRX( ncmpi_inq_dimid(ncid, NC_FRAME_STR, &frame_dim), NC_FRAME_STR );
|
||||
NCERRX( ncmpi_inq_dimid(ncid, NC_SPATIAL_STR, &spatial_dim),
|
||||
NC_SPATIAL_STR );
|
||||
NCERRX( ncmpi_inq_dimid(ncid, NC_VOIGT_STR, &Voigt_dim), NC_VOIGT_STR );
|
||||
NCERRX( ncmpi_inq_dimid(ncid, NC_ATOM_STR, &atom_dim), NC_ATOM_STR );
|
||||
NCERRX( ncmpi_inq_dimid(ncid, NC_CELL_SPATIAL_STR, &cell_spatial_dim),
|
||||
NC_CELL_SPATIAL_STR );
|
||||
|
|
@ -304,6 +304,26 @@ void DumpNetCDFMPIIO::openfile()
|
|||
NC_CELL_ANGULAR_STR );
|
||||
NCERRX( ncmpi_inq_dimid(ncid, NC_LABEL_STR, &label_dim), NC_LABEL_STR );
|
||||
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
int dims = perat[i].dims;
|
||||
if (vector_dim[dims] < 0) {
|
||||
char dimstr[1024];
|
||||
if (dims == 3) {
|
||||
strcpy(dimstr, NC_SPATIAL_STR);
|
||||
}
|
||||
else if (dims == 6) {
|
||||
strcpy(dimstr, NC_VOIGT_STR);
|
||||
}
|
||||
else {
|
||||
sprintf(dimstr, "vec%i", dims);
|
||||
}
|
||||
if (dims != 1) {
|
||||
NCERRX( ncmpi_inq_dimid(ncid, dimstr, &vector_dim[dims]),
|
||||
dimstr );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// default variables
|
||||
NCERRX( ncmpi_inq_varid(ncid, NC_SPATIAL_STR, &spatial_var),
|
||||
NC_SPATIAL_STR );
|
||||
|
|
@ -358,10 +378,6 @@ void DumpNetCDFMPIIO::openfile()
|
|||
// dimensions
|
||||
NCERRX( ncmpi_def_dim(ncid, NC_FRAME_STR, NC_UNLIMITED, &frame_dim),
|
||||
NC_FRAME_STR );
|
||||
NCERRX( ncmpi_def_dim(ncid, NC_SPATIAL_STR, 3, &spatial_dim),
|
||||
NC_SPATIAL_STR );
|
||||
NCERRX( ncmpi_def_dim(ncid, NC_VOIGT_STR, 6, &Voigt_dim),
|
||||
NC_VOIGT_STR );
|
||||
NCERRX( ncmpi_def_dim(ncid, NC_ATOM_STR, ntotalgr, &atom_dim),
|
||||
NC_ATOM_STR );
|
||||
NCERRX( ncmpi_def_dim(ncid, NC_CELL_SPATIAL_STR, 3, &cell_spatial_dim),
|
||||
|
|
@ -371,13 +387,33 @@ void DumpNetCDFMPIIO::openfile()
|
|||
NCERRX( ncmpi_def_dim(ncid, NC_LABEL_STR, 10, &label_dim),
|
||||
NC_LABEL_STR );
|
||||
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
int dims = perat[i].dims;
|
||||
if (vector_dim[dims] < 0) {
|
||||
char dimstr[1024];
|
||||
if (dims == 3) {
|
||||
strcpy(dimstr, NC_SPATIAL_STR);
|
||||
}
|
||||
else if (dims == 6) {
|
||||
strcpy(dimstr, NC_VOIGT_STR);
|
||||
}
|
||||
else {
|
||||
sprintf(dimstr, "vec%i", dims);
|
||||
}
|
||||
if (dims != 1) {
|
||||
NCERRX( ncmpi_def_dim(ncid, dimstr, dims, &vector_dim[dims]),
|
||||
dimstr );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// default variables
|
||||
dims[0] = spatial_dim;
|
||||
dims[0] = vector_dim[3];
|
||||
NCERRX( ncmpi_def_var(ncid, NC_SPATIAL_STR, NC_CHAR, 1, dims, &spatial_var),
|
||||
NC_SPATIAL_STR );
|
||||
NCERRX( ncmpi_def_var(ncid, NC_CELL_SPATIAL_STR, NC_CHAR, 1, dims,
|
||||
&cell_spatial_var), NC_CELL_SPATIAL_STR );
|
||||
dims[0] = spatial_dim;
|
||||
dims[0] = vector_dim[3];
|
||||
dims[1] = label_dim;
|
||||
NCERRX( ncmpi_def_var(ncid, NC_CELL_ANGULAR_STR, NC_CHAR, 2, dims,
|
||||
&cell_angular_var), NC_CELL_ANGULAR_STR );
|
||||
|
|
@ -399,7 +435,7 @@ void DumpNetCDFMPIIO::openfile()
|
|||
// variables specified in the input file
|
||||
dims[0] = frame_dim;
|
||||
dims[1] = atom_dim;
|
||||
dims[2] = spatial_dim;
|
||||
dims[2] = vector_dim[3];
|
||||
|
||||
for (int i = 0; i < n_perat; i++) {
|
||||
nc_type xtype;
|
||||
|
|
@ -416,28 +452,15 @@ void DumpNetCDFMPIIO::openfile()
|
|||
xtype = NC_FLOAT;
|
||||
}
|
||||
|
||||
if (perat[i].dims == 6) {
|
||||
// this is a tensor in Voigt notation
|
||||
dims[2] = Voigt_dim;
|
||||
NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else if (perat[i].dims == 3) {
|
||||
// this is a vector, we need to store x-, y- and z-coordinates
|
||||
dims[2] = spatial_dim;
|
||||
NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else if (perat[i].dims == 1) {
|
||||
if (perat[i].dims == 1) {
|
||||
NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 2, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
else {
|
||||
char errstr[1024];
|
||||
sprintf(errstr, "%i dimensions for '%s'. Not sure how to write "
|
||||
"this to the NetCDF trajectory file.", perat[i].dims,
|
||||
perat[i].name);
|
||||
error->all(FLERR,errstr);
|
||||
// this is a vector
|
||||
dims[2] = vector_dim[perat[i].dims];
|
||||
NCERRX( ncmpi_def_var(ncid, perat[i].name, xtype, 3, dims,
|
||||
&perat[i].var), perat[i].name );
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -71,8 +71,7 @@ class DumpNetCDFMPIIO : public DumpCustom {
|
|||
int ncid;
|
||||
|
||||
int frame_dim;
|
||||
int spatial_dim;
|
||||
int Voigt_dim;
|
||||
int vector_dim[DUMP_NC_MPIIO_MAX_DIMS];
|
||||
int atom_dim;
|
||||
int cell_spatial_dim;
|
||||
int cell_angular_dim;
|
||||
|
|
|
|||
|
|
@ -239,17 +239,17 @@ void ImproperFourierOMP::add1_thr(const int i1,const int i2,
|
|||
dahy = ary-c*hry;
|
||||
dahz = arz-c*hrz;
|
||||
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
|
||||
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
|
||||
|
||||
f4[0] = dahx*rhr;
|
||||
f4[1] = dahy*rhr;
|
||||
f4[2] = dahz*rhr;
|
||||
f4[0] = dahx*rhr*a;
|
||||
f4[1] = dahy*rhr*a;
|
||||
f4[2] = dahz*rhr*a;
|
||||
|
||||
f1[0] = -(f2[0] + f3[0] + f4[0]);
|
||||
f1[1] = -(f2[1] + f3[1] + f4[1]);
|
||||
|
|
@ -258,30 +258,31 @@ void ImproperFourierOMP::add1_thr(const int i1,const int i2,
|
|||
// apply force to each of 4 atoms
|
||||
|
||||
if (NEWTON_BOND || i1 < nlocal) {
|
||||
f[i1][0] += f1[0]*a;
|
||||
f[i1][1] += f1[1]*a;
|
||||
f[i1][2] += f1[2]*a;
|
||||
f[i1][0] += f1[0];
|
||||
f[i1][1] += f1[1];
|
||||
f[i1][2] += f1[2];
|
||||
}
|
||||
|
||||
if (NEWTON_BOND || i2 < nlocal) {
|
||||
f[i2][0] += f3[0]*a;
|
||||
f[i2][1] += f3[1]*a;
|
||||
f[i2][2] += f3[2]*a;
|
||||
f[i2][0] += f3[0];
|
||||
f[i2][1] += f3[1];
|
||||
f[i2][2] += f3[2];
|
||||
}
|
||||
|
||||
if (NEWTON_BOND || i3 < nlocal) {
|
||||
f[i3][0] += f2[0]*a;
|
||||
f[i3][1] += f2[1]*a;
|
||||
f[i3][2] += f2[2]*a;
|
||||
f[i3][0] += f2[0];
|
||||
f[i3][1] += f2[1];
|
||||
f[i3][2] += f2[2];
|
||||
}
|
||||
|
||||
if (NEWTON_BOND || i4 < nlocal) {
|
||||
f[i4][0] += f4[0]*a;
|
||||
f[i4][1] += f4[1]*a;
|
||||
f[i4][2] += f4[2]*a;
|
||||
f[i4][0] += f4[0];
|
||||
f[i4][1] += f4[1];
|
||||
f[i4][2] += f4[2];
|
||||
}
|
||||
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4,
|
||||
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr);
|
||||
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f2,f4,
|
||||
-vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z,thr);
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -212,17 +212,17 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
|
|||
dahy = ary-c*hry;
|
||||
dahz = arz-c*hrz;
|
||||
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
|
||||
f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
|
||||
f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
|
||||
f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
|
||||
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
|
||||
f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
|
||||
f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
|
||||
f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
|
||||
|
||||
f4[0] = dahx*rhr;
|
||||
f4[1] = dahy*rhr;
|
||||
f4[2] = dahz*rhr;
|
||||
f4[0] = dahx*rhr*a;
|
||||
f4[1] = dahy*rhr*a;
|
||||
f4[2] = dahz*rhr*a;
|
||||
|
||||
f1[0] = -(f2[0] + f3[0] + f4[0]);
|
||||
f1[1] = -(f2[1] + f3[1] + f4[1]);
|
||||
|
|
@ -231,27 +231,27 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
|
|||
// apply force to each of 4 atoms
|
||||
|
||||
if (NEWTON_BOND || i1 < nlocal) {
|
||||
f[i1].x += f1[0]*a;
|
||||
f[i1].y += f1[1]*a;
|
||||
f[i1].z += f1[2]*a;
|
||||
f[i1].x += f1[0];
|
||||
f[i1].y += f1[1];
|
||||
f[i1].z += f1[2];
|
||||
}
|
||||
|
||||
if (NEWTON_BOND || i2 < nlocal) {
|
||||
f[i2].x += f3[0]*a;
|
||||
f[i2].y += f3[1]*a;
|
||||
f[i2].z += f3[2]*a;
|
||||
f[i2].x += f3[0];
|
||||
f[i2].y += f3[1];
|
||||
f[i2].z += f3[2];
|
||||
}
|
||||
|
||||
if (NEWTON_BOND || i3 < nlocal) {
|
||||
f[i3].x += f2[0]*a;
|
||||
f[i3].y += f2[1]*a;
|
||||
f[i3].z += f2[2]*a;
|
||||
f[i3].x += f2[0];
|
||||
f[i3].y += f2[1];
|
||||
f[i3].z += f2[2];
|
||||
}
|
||||
|
||||
if (NEWTON_BOND || i4 < nlocal) {
|
||||
f[i4].x += f4[0]*a;
|
||||
f[i4].y += f4[1]*a;
|
||||
f[i4].z += f4[2]*a;
|
||||
f[i4].x += f4[0];
|
||||
f[i4].y += f4[1];
|
||||
f[i4].z += f4[2];
|
||||
}
|
||||
|
||||
if (EVFLAG) {
|
||||
|
|
@ -270,7 +270,7 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
|
|||
vb3y = x[i4].y - x[i3].y;
|
||||
vb3z = x[i4].z - x[i3].z;
|
||||
|
||||
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4,
|
||||
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f2,f4,
|
||||
vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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/switch3/coulgauss/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] = 0.5*(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 / (2.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(2.0*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 mm3/switch3/coulgauss/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.
|
||||
|
||||
*/
|
||||
|
|
@ -725,6 +725,24 @@ void CreateAtoms::add_lattice()
|
|||
domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxhi[2],
|
||||
xmin,ymin,zmin,xmax,ymax,zmax);
|
||||
|
||||
// narrow down min/max further by extent of the region, if possible
|
||||
|
||||
if (domain->regions[nregion]->bboxflag) {
|
||||
const double rxmin = domain->regions[nregion]->extent_xlo;
|
||||
const double rxmax = domain->regions[nregion]->extent_xhi;
|
||||
const double rymin = domain->regions[nregion]->extent_ylo;
|
||||
const double rymax = domain->regions[nregion]->extent_yhi;
|
||||
const double rzmin = domain->regions[nregion]->extent_zlo;
|
||||
const double rzmax = domain->regions[nregion]->extent_zhi;
|
||||
|
||||
if (rxmin > xmin) xmin = (rxmin > xmax) ? xmax : rxmin;
|
||||
if (rxmax < xmax) xmax = (rxmax < xmin) ? xmin : rxmax;
|
||||
if (rymin > ymin) ymin = (rymin > ymax) ? ymax : rymin;
|
||||
if (rymax < ymax) ymax = (rymax < ymin) ? ymin : rymax;
|
||||
if (rzmin > zmin) zmin = (rzmin > zmax) ? zmax : rzmin;
|
||||
if (rzmax < zmax) zmax = (rzmax < zmin) ? zmin : rzmax;
|
||||
}
|
||||
|
||||
// ilo:ihi,jlo:jhi,klo:khi = loop bounds for lattice overlap of my subbox
|
||||
// overlap = any part of a unit cell (face,edge,pt) in common with my subbox
|
||||
// in lattice space, subbox is a tilted box
|
||||
|
|
@ -751,15 +769,32 @@ void CreateAtoms::add_lattice()
|
|||
// convert lattice coords to box coords
|
||||
// add atom or molecule (on each basis point) if it meets all criteria
|
||||
|
||||
double **basis = domain->lattice->basis;
|
||||
double x[3],lamda[3];
|
||||
double *coord;
|
||||
const double * const * const basis = domain->lattice->basis;
|
||||
|
||||
// rough estimate of total time used for create atoms.
|
||||
// one inner loop takes about 25ns on a typical desktop CPU core in 2019
|
||||
double testimate = 2.5e-8/3600.0; // convert seconds to hours
|
||||
testimate *= static_cast<double>(khi-klo+1);
|
||||
testimate *= static_cast<double>(jhi-jlo+1);
|
||||
testimate *= static_cast<double>(ihi-ilo+1);
|
||||
testimate *= static_cast<double>(nbasis);
|
||||
double maxestimate = 0.0;
|
||||
MPI_Reduce(&testimate,&maxestimate,1,MPI_DOUBLE,MPI_MAX,0,world);
|
||||
|
||||
if ((comm->me == 0) && (maxestimate > 0.01)) {
|
||||
if (screen) fprintf(screen,"WARNING: create_atoms will take "
|
||||
"approx. %.2f hours to complete\n",maxestimate);
|
||||
if (logfile) fprintf(logfile,"WARNING: create_atoms will take "
|
||||
"approx. %.2f hours to complete\n",maxestimate);
|
||||
}
|
||||
|
||||
int i,j,k,m;
|
||||
for (k = klo; k <= khi; k++)
|
||||
for (j = jlo; j <= jhi; j++)
|
||||
for (i = ilo; i <= ihi; i++)
|
||||
for (k = klo; k <= khi; k++) {
|
||||
for (j = jlo; j <= jhi; j++) {
|
||||
for (i = ilo; i <= ihi; i++) {
|
||||
for (m = 0; m < nbasis; m++) {
|
||||
double *coord;
|
||||
double x[3],lamda[3];
|
||||
|
||||
x[0] = i + basis[m][0];
|
||||
x[1] = j + basis[m][1];
|
||||
|
|
@ -794,8 +829,12 @@ void CreateAtoms::add_lattice()
|
|||
if (mode == ATOM) atom->avec->create_atom(basistype[m],x);
|
||||
else add_molecule(x);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
add a randomly rotated molecule with its center at center
|
||||
if quat_user set, perform requested rotation
|
||||
|
|
|
|||
Loading…
Reference in New Issue