forked from lijiext/lammps
Merge pull request #782 from hheenen/user-mofff-contribution
User mofff contribution
This commit is contained in:
Steve Plimpton
GitHub
commit
160edc9532
|
|
@ -104,7 +104,7 @@ set(OTHER_PACKAGES KIM PYTHON MSCG MPIIO VORONOI POEMS LATTE
|
|||
USER-ATC USER-AWPMD USER-CGDNA USER-MESO
|
||||
USER-CGSDK USER-COLVARS USER-DIFFRACTION USER-DPD USER-DRUDE USER-EFF
|
||||
USER-FEP USER-H5MD USER-LB USER-MANIFOLD USER-MEAMC USER-MGPT USER-MISC
|
||||
USER-MOLFILE USER-NETCDF USER-PHONON USER-QTB USER-REAXC USER-SMD
|
||||
USER-MOFFF USER-MOLFILE USER-NETCDF USER-PHONON USER-QTB USER-REAXC USER-SMD
|
||||
USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM)
|
||||
set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU)
|
||||
foreach(PKG ${DEFAULT_PACKAGES})
|
||||
|
|
|
|||
Binary file not shown.
|
After Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 5.2 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 4.5 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 6.9 KiB |
|
|
@ -0,0 +1,9 @@
|
|||
\documentclass[12pt]{article}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\begin{eqnarray*}
|
||||
E = A e^{-\kappa r} - \frac{C}{r^6} \cdot \frac{1}{1 + D r^{14}} \qquad r < r_c \\
|
||||
\end{eqnarray*}
|
||||
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 7.1 KiB |
|
|
@ -138,6 +138,7 @@ Package, Description, Doc page, Example, Library
|
|||
"USER-MESO"_#USER-MESO, mesoscale DPD models, "pair_style edpd"_pair_meso.html, USER/meso, -
|
||||
"USER-MGPT"_#USER-MGPT, fast MGPT multi-ion potentials, "pair_style mgpt"_pair_mgpt.html, USER/mgpt, -
|
||||
"USER-MISC"_#USER-MISC, single-file contributions, USER-MISC/README, USER/misc, -
|
||||
"USER-MOFFF"_#USER-MOFFF, styles for "MOF-FF"_MOFplus force field, "pair_style buck6d/coul/gauss"_pair_buck6d_coul_gauss.html, USER/mofff, -
|
||||
"USER-MOLFILE"_#USER-MOLFILE, "VMD"_vmd_home molfile plug-ins,"dump molfile"_dump_molfile.html, -, ext
|
||||
"USER-NETCDF"_#USER-NETCDF, dump output via NetCDF,"dump netcdf"_dump_netcdf.html, -, ext
|
||||
"USER-OMP"_#USER-OMP, OpenMP-enabled styles,"Section 5.3.4"_accelerate_omp.html, "Benchmarks"_http://lammps.sandia.gov/bench.html, -
|
||||
|
|
@ -2259,6 +2260,44 @@ http://lammps.sandia.gov/movies.html#mesodpd :ul
|
|||
|
||||
:line
|
||||
|
||||
USER-MOFFF package :link(USER-MOFFF),h4
|
||||
|
||||
[Contents:]
|
||||
|
||||
Pair, angle and improper styles needed to employ the MOF-FF
|
||||
force field by Schmid and coworkers with LAMMPS.
|
||||
MOF-FF is a first principles derived force field with the primary aim
|
||||
to simulate MOFs and related porous framework materials, using spherical
|
||||
Gaussian charges. It is described in S. Bureekaew et al., Phys. Stat. Sol. B
|
||||
2013, 250, 1128-1141.
|
||||
For the usage of MOF-FF see the example in the example directory as
|
||||
well as the "MOF+"_MOFplus website.
|
||||
|
||||
:link(MOFplus,https://www.mofplus.org/content/show/MOF-FF)
|
||||
|
||||
[Author:] Hendrik Heenen (Technical U of Munich),
|
||||
Rochus Schmid (Ruhr-University Bochum).
|
||||
|
||||
[Install or un-install:]
|
||||
|
||||
make yes-user-mofff
|
||||
make machine :pre
|
||||
|
||||
make no-user-mofff
|
||||
make machine :pre
|
||||
|
||||
[Supporting info:]
|
||||
|
||||
src/USER-MOFFF: filenames -> commands
|
||||
src/USER-MOFFF/README
|
||||
"pair_style buck6d/coul/gauss"_pair_buck6d_coul_gauss.html
|
||||
"angle_style class2"_angle_class2.html
|
||||
"angle_style cosine/buck6d"_angle_cosine_buck6d.html
|
||||
"improper_style inversion/harmonic"_improper_inversion_harmonic.html
|
||||
examples/USER/mofff :ul
|
||||
|
||||
:line
|
||||
|
||||
USER-MOLFILE package :link(USER-MOLFILE),h4
|
||||
|
||||
[Contents:]
|
||||
|
|
|
|||
|
|
@ -9,6 +9,7 @@
|
|||
angle_style class2 command :h3
|
||||
angle_style class2/omp command :h3
|
||||
angle_style class2/kk command :h3
|
||||
angle_style class2/p6 command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
|
|
@ -102,11 +103,29 @@ more instructions on how to use the accelerated styles effectively.
|
|||
|
||||
:line
|
||||
|
||||
The {class2/p6} angle style uses the {class2} potential expanded to sixth order:
|
||||
|
||||
:c,image(Eqs/angle_class2_p6.jpg)
|
||||
|
||||
In this expanded term 6 coefficients for the Ea formula need to be set:
|
||||
|
||||
theta0 (degrees)
|
||||
K2 (energy/radian^2)
|
||||
K3 (energy/radian^3)
|
||||
K4 (energy/radian^4)
|
||||
K5 (energy/radian^5)
|
||||
K6 (energy/radian^6) :ul
|
||||
|
||||
The bond-bond and bond-angle terms remain unchanged.
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
This angle style can only be used if LAMMPS was built with the CLASS2
|
||||
package. See the "Making LAMMPS"_Section_start.html#start_3 section
|
||||
for more info on packages.
|
||||
package. For the {class2/p6} style LAMMPS needs to be built with the
|
||||
USER-MOFFF package. See the "Making LAMMPS"_Section_start.html#start_3
|
||||
section for more info on packages.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,65 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
angle_style cosine/buck6d command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
angle_style cosine/buck6d :pre
|
||||
|
||||
[Examples:]
|
||||
|
||||
angle_style cosine/buck6d
|
||||
angle_coeff 1 cosine/buck6d 1.978350 4 180.000000 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
The {cosine/buck6d} angle style uses the potential
|
||||
|
||||
:c,image(Eqs/angle_cosine_buck6d.jpg)
|
||||
|
||||
where K is the energy constant, n is the periodic multiplicity and
|
||||
Theta0 is the equilibrium angle.
|
||||
|
||||
The 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 in the following order:
|
||||
|
||||
K (energy)
|
||||
n
|
||||
Theta0 (degrees) :ul
|
||||
|
||||
Theta0 is specified in degrees, but LAMMPS converts it to radians
|
||||
internally.
|
||||
|
||||
Additional to the cosine term the {cosine/buck6d} angle style computes
|
||||
the short range (vdW) interaction belonging to the
|
||||
"pair_buck6d"_pair_buck6d_coul_gauss.html between the end atoms of
|
||||
the angle. For this reason this angle style only works in combination
|
||||
with the "pair_buck6d"_pair_buck6d_coul_gauss.html styles and needs
|
||||
the "special_bonds"_special_bonds.html 1-3 interactions to be weighted
|
||||
0.0 to prevent double counting.
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
{cosine/buck6d} can only be used in combination with the
|
||||
"pair_buck6d"_pair_buck6d_coul_gauss.html style and with a
|
||||
"special_bonds"_special_bonds.html 0.0 weighting of 1-3 interactions.
|
||||
|
||||
This angle style can only be used if LAMMPS was built with the
|
||||
USER-MOFFF package. See the "Making
|
||||
LAMMPS"_Section_start.html#start_3 section for more info on packages.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"angle_coeff"_angle_coeff.html
|
||||
|
||||
[Default:] none
|
||||
|
|
@ -8,6 +8,7 @@ Angle Styles :h1
|
|||
angle_charmm
|
||||
angle_class2
|
||||
angle_cosine
|
||||
angle_cosine_buck6d
|
||||
angle_cosine_delta
|
||||
angle_cosine_periodic
|
||||
angle_cosine_shift
|
||||
|
|
|
|||
|
|
@ -0,0 +1,65 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
improper_style inversion/harmonic command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
improper_style inversion/harmonic :pre
|
||||
|
||||
[Examples:]
|
||||
|
||||
improper_style inversion/harmonic
|
||||
improper_coeff 1 18.776340 0.000000 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
The {inversion/harmonic} improper style follows the Wilson-Decius
|
||||
out-of-plane angle definition and uses an harmonic potential:
|
||||
|
||||
:c,image(Eqs/improper_inversion_harmonic.jpg)
|
||||
|
||||
where K is the force constant and omega is the angle evaluated for
|
||||
all three axis-plane combinations centered around the atom I. For
|
||||
the IL axis and the IJK plane omega looks as follows:
|
||||
|
||||
:c,image(Eqs/umbrella.jpg)
|
||||
|
||||
Note that the {inversion/harmonic} angle term evaluation differs to
|
||||
the "improper_umbrella"_improper_umbrella.html due to the cyclic
|
||||
evaluation of all possible angles omega.
|
||||
|
||||
The following coefficients must be defined for each improper type via
|
||||
the "improper_coeff"_improper_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)
|
||||
omega0 (degrees) :ul
|
||||
|
||||
If omega0 = 0 the potential term has a minimum for the planar
|
||||
structure. Otherwise it has two minima at +/- omega0, with a barrier
|
||||
in between.
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
This improper style can only be used if LAMMPS was built with the
|
||||
USER-MOFFF package. See the "Making
|
||||
LAMMPS"_Section_start.html#start_3 section for more info on packages.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"improper_coeff"_improper_coeff.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
|
|
@ -12,6 +12,7 @@ Improper Styles :h1
|
|||
improper_fourier
|
||||
improper_harmonic
|
||||
improper_hybrid
|
||||
improper_inversion_harmonic
|
||||
improper_none
|
||||
improper_ring
|
||||
improper_umbrella
|
||||
|
|
|
|||
|
|
@ -0,0 +1,138 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
pair_style buck6d/coul/gauss/dsf :h3
|
||||
pair_style buck6d/coul/gauss/long :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
pair_style style args :pre
|
||||
|
||||
style = {buck6d/coul/gauss/dsf} or {buck6d/coul/gauss/long}
|
||||
args = list of arguments for a particular style :ul
|
||||
{buck6d/coul/gauss/dsf} args = smooth cutoff (cutoff2)
|
||||
smooth = smoothing onset within Buckingham cutoff (ratio)
|
||||
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
|
||||
cutoff2 = global cutoff for Coulombic (optional) (distance units)
|
||||
{buck6d/coul/gauss/long} args = smooth smooth2 cutoff (cutoff2)
|
||||
smooth = smoothing onset within Buckingham cutoff (ratio)
|
||||
smooth2 = smoothing onset within Coulombic cutoff (ratio)
|
||||
cutoff = global cutoff for Buckingham (and Coulombic if only 1 arg) (distance units)
|
||||
cutoff2 = global cutoff for Coulombic (optional) (distance units) :pre
|
||||
|
||||
[Examples:]
|
||||
|
||||
pair_style buck6d/coul/gauss/dsf 0.9000 12.0000
|
||||
pair_coeff 1 1 1030. 3.061 457.179 4.521 0.608 :pre
|
||||
|
||||
pair_style buck6d/coul/gauss/long 0.9000 1.0000 12.0000
|
||||
pair_coeff 1 1 1030. 3.061 457.179 4.521 0.608 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
The {buck6d/coul/gauss} styles evaluate vdW and Coulomb
|
||||
interactions following the MOF-FF force field after
|
||||
"(Schmid)"_#Schmid. The vdW term of the {buck6d} styles
|
||||
computes a dispersion damped Buckingham potential:
|
||||
|
||||
:c,image(Eqs/pair_buck6d.jpg)
|
||||
|
||||
where A and C are a force constant, kappa is an ionic-pair dependent
|
||||
reciprocal length parameter, D is a dispersion correction parameter,
|
||||
and the cutoff Rc truncates the interaction distance.
|
||||
The first term in the potential corresponds to the Buckingham
|
||||
repulsion term and the second term to the dispersion attraction with
|
||||
a damping correction analog to the Grimme correction used in DFT.
|
||||
The latter corrects for artifacts occurring at short distances which
|
||||
become an issue for soft vdW potentials.
|
||||
|
||||
The {buck6d} styles include a smoothing function which is invoked
|
||||
according to the global smooting parameter within the specified
|
||||
cutoff. Hereby a parameter of i.e. 0.9 invokes the smoothing
|
||||
within 90% of the cutoff. No smoothing is applied at a value
|
||||
of 1.0. For the {gauss/dsf} style this smoothing is only applicable
|
||||
for the dispersion damped Buckingham potential. For the {gauss/long}
|
||||
styles the smoothing function can also be invoked for the real
|
||||
space coulomb interactions which enforce continous energies and
|
||||
forces at the cutoff.
|
||||
|
||||
Both styles {buck6d/coul/gauss/dsf} and {buck6d/coul/gauss/long}
|
||||
evaluate a Coulomb potential using spherical Gaussian type charge
|
||||
distributions which effectively dampen electrostatic ineractions
|
||||
for high charges at close distances. The electrostatic potential
|
||||
is thus evaluated as:
|
||||
|
||||
:c,image(Eqs/pair_coul_gauss.jpg)
|
||||
|
||||
where C is an energy-conversion constant, 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, alpha is
|
||||
ion pair dependent damping parameter and erf() is the error-function.
|
||||
The cutoff Rc truncates the interaction distance.
|
||||
|
||||
The style {buck6d/coul/gauss/dsf} computes the Coulomb interaction
|
||||
via the damped shifted force model described in "(Fennell)"_#Fennell
|
||||
approximating an Ewald sum similar to the "pair coul/dsf"_pair_coul.html
|
||||
styles. In {buck6d/coul/gauss/long} an additional damping factor is
|
||||
applied to the Coulombic term so it can be used in conjunction with the
|
||||
"kspace_style"_kspace_style.html command and its {ewald} or {pppm}
|
||||
options. The Coulombic cutoff in this case separates the real and
|
||||
reciprocal space evaluation of the Ewald sum.
|
||||
|
||||
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:
|
||||
|
||||
A (energy units)
|
||||
rho (distance^-1 units)
|
||||
C (energy-distance^6 units)
|
||||
D (distance^14 units)
|
||||
alpha (distance^-1 units)
|
||||
cutoff (distance units) :ul
|
||||
|
||||
The second coefficient, rho, must be greater than zero. The latter
|
||||
coefficient is optional. If not specified, the global vdW cutoff
|
||||
is used.
|
||||
|
||||
:line
|
||||
|
||||
[Mixing, shift, table, tail correction, restart, rRESPA info]:
|
||||
|
||||
These pair styles do not support mixing. Thus, coefficients for all
|
||||
I,J pairs must be specified explicitly.
|
||||
|
||||
These styles do not support the "pair_modify"_pair_modify.html shift
|
||||
option for the energy. Instead the smoothing function should be applied
|
||||
by setting the global smoothing parameter to a value < 1.0.
|
||||
|
||||
These styles write their information to "binary restart
|
||||
files"_restart.html, so pair_style and pair_coeff commands do not need
|
||||
to be specified in an input script that reads a restart file.
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
These styles are part of the USER-MOFFF package. They are only enabled
|
||||
if LAMMPS was built with that package. See the
|
||||
"Making LAMMPS"_Section_start.html#start_3 section for more info.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"pair_coeff"_pair_coeff.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:link(Schmid)
|
||||
[(Schmid)] S. Bureekaew, S. Amirjalayer, M. Tafipolsky, C. Spickermann, T.K. Roy and R. Schmid, Phys. Status Solidi B, 6, 1128 (2013).
|
||||
:link(Fennell)
|
||||
[(Fennell)] C. J. Fennell, J. D. Gezelter, J Chem Phys, 124, 234104 (2006).
|
||||
|
|
@ -17,6 +17,7 @@ Pair Styles :h1
|
|||
pair_brownian
|
||||
pair_buck
|
||||
pair_buck_long
|
||||
pair_buck6d_coul_gauss
|
||||
pair_charmm
|
||||
pair_class2
|
||||
pair_colloid
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,93 @@
|
|||
clear
|
||||
units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
read_data hkust1.data
|
||||
|
||||
neighbor 2.0 bin
|
||||
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD ------------------------------
|
||||
|
||||
pair_style buck6d/coul/gauss/dsf 0.9000 12.0000
|
||||
|
||||
pair_coeff 1 1 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW
|
||||
pair_coeff 1 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 1 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 1 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 1 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 1 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 2 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 2 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 2 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 2 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 2 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 3 3 3680 4 32.805 0.10690769 0.9771554 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 3 4 14157.243 3.1914894 489.18197 2.5231391 0.45538909 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 3 5 6320.6076 3.6144578 103.50278 0.44181916 0.75109952 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 3 6 6157.8178 3.4883721 124.7792 0.72632262 0.73006542 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 4 4 54464 2.6548673 5679.4311 33.208515 0.34105936 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 4 5 24315.863 2.9411765 1371.5617 7.9168726 0.42457748 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 4 6 23689.598 2.8571429 1590.0769 11.87959 0.42066711 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 5 5 10856 3.2967033 308.7755 1.6022517 0.63272774 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 5 6 10576.399 3.1914894 365.45138 2.5231391 0.61999948 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 6 6 10304 3.0927835 429.89352 3.9170177 0.60800971 # buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
|
||||
bond_style hybrid class2 morse
|
||||
|
||||
bond_coeff 5 morse 50.000000 1.451345 1.914000 # morse->(cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 4 class2 2.536000 75.465060 -192.435903 286.248406 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2)|CuPW
|
||||
bond_coeff 3 class2 1.094000 390.634200 -996.117210 1481.724350 # mm3->(c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
bond_coeff 6 class2 1.278000 585.591600 -1493.258580 2221.222138 # mm3->(c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 1 class2 1.394000 509.335200 -1298.804760 1931.972080 # mm3->(c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
bond_coeff 2 class2 1.485000 360.635220 -919.619811 1367.934469 # mm3->(c3_c1o2@co2,c3_c3@ph)|CuPW
|
||||
|
||||
angle_style hybrid class2/p6 cosine/buck6d
|
||||
|
||||
angle_coeff 2 class2/p6 117.711000 57.408120 -46.049402 10.553745 -7.558563 13.610890 # mm3->(c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
angle_coeff 2 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 2 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 126.814000 13.740540 -11.021849 2.526022 -1.809130 3.257744 # mm3->(c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
angle_coeff 6 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 7 class2/p6 123.490000 111.075360 -89.098091 20.419778 -14.624589 26.334856 # mm3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 7 class2/p6 bb 14.244120 1.278000 1.278000
|
||||
angle_coeff 7 class2/p6 ba 3.309240 3.309240 1.278000 1.278000
|
||||
angle_coeff 1 class2/p6 127.050000 53.307540 -42.760159 9.799907 -7.018666 12.638684 # mm3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
angle_coeff 1 class2/p6 bb 71.796120 1.394000 1.394000
|
||||
angle_coeff 1 class2/p6 ba 6.762360 6.762360 1.394000 1.394000
|
||||
angle_coeff 4 class2/p6 84.336000 29.351520 -23.544055 5.395900 -3.864529 6.958951 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 4 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 4 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 3 class2/p6 120.350000 36.185820 -29.026127 6.652298 -4.764358 8.579296 # mm3->(c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
angle_coeff 3 class2/p6 bb 93.378120 1.394000 1.094000
|
||||
angle_coeff 3 class2/p6 ba -25.179000 53.523360 1.394000 1.094000
|
||||
angle_coeff 8 class2/p6 115.098000 79.493700 -63.765149 14.613896 -10.466432 18.847160 # mm3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 8 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 8 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 5 cosine/buck6d 1.978350 4 180.000000 # fourier->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
|
||||
dihedral_style opls
|
||||
|
||||
dihedral_coeff 3 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 9 0.000000 4.528000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 8 0.000000 0.000000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 5 0.000000 1.741000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
dihedral_coeff 2 0.000000 6.316000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
dihedral_coeff 1 0.000000 4.379000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
dihedral_coeff 6 0.000000 0.000000 0.000000 0.609000 # cos4->(o2_c1cu1@co2,cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
dihedral_coeff 4 0.000000 0.000000 0.000000 0.000000 # cos3->(h1_c1@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 10 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 7 0.000000 0.000000 0.000000 0.000000 # cos3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
|
||||
improper_style inversion/harmonic
|
||||
|
||||
improper_coeff 1 18.776340 0.000000 # harm->(c3_c3@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c1o2@co2)|CuPW
|
||||
improper_coeff 3 41.005800 0.000000 # harm->(c3_c1o2@co2,c3_c3@ph,o2_c1cu1@co2,o2_c1cu1@co2)|CuPW
|
||||
improper_coeff 2 4.100580 0.000000 # harm->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
|
||||
special_bonds lj 0.0 0.0 1.0 coul 1.0 1.0 1.0
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD END --------------------------
|
||||
|
||||
run 0
|
||||
|
|
@ -0,0 +1,94 @@
|
|||
clear
|
||||
units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
read_data hkust1.data
|
||||
|
||||
neighbor 2.0 bin
|
||||
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD ------------------------------
|
||||
kspace_style ewald 1.0e-6
|
||||
|
||||
pair_style buck6d/coul/gauss/long 0.9 0.9 12.0000
|
||||
|
||||
pair_coeff 1 1 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW
|
||||
pair_coeff 1 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 1 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 1 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 1 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 1 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 2 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 2 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 2 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 2 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 2 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 3 3 3680 4 32.805 0.10690769 0.9771554 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 3 4 14157.243 3.1914894 489.18197 2.5231391 0.45538909 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 3 5 6320.6076 3.6144578 103.50278 0.44181916 0.75109952 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 3 6 6157.8178 3.4883721 124.7792 0.72632262 0.73006542 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 4 4 54464 2.6548673 5679.4311 33.208515 0.34105936 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 4 5 24315.863 2.9411765 1371.5617 7.9168726 0.42457748 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 4 6 23689.598 2.8571429 1590.0769 11.87959 0.42066711 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 5 5 10856 3.2967033 308.7755 1.6022517 0.63272774 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 5 6 10576.399 3.1914894 365.45138 2.5231391 0.61999948 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 6 6 10304 3.0927835 429.89352 3.9170177 0.60800971 # buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
|
||||
bond_style hybrid class2 morse
|
||||
|
||||
bond_coeff 5 morse 50.000000 1.451345 1.914000 # morse->(cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 4 class2 2.536000 75.465060 -192.435903 286.248406 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2)|CuPW
|
||||
bond_coeff 3 class2 1.094000 390.634200 -996.117210 1481.724350 # mm3->(c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
bond_coeff 6 class2 1.278000 585.591600 -1493.258580 2221.222138 # mm3->(c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 1 class2 1.394000 509.335200 -1298.804760 1931.972080 # mm3->(c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
bond_coeff 2 class2 1.485000 360.635220 -919.619811 1367.934469 # mm3->(c3_c1o2@co2,c3_c3@ph)|CuPW
|
||||
|
||||
angle_style hybrid class2/p6 cosine/buck6d
|
||||
|
||||
angle_coeff 2 class2/p6 117.711000 57.408120 -46.049402 10.553745 -7.558563 13.610890 # mm3->(c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
angle_coeff 2 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 2 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 126.814000 13.740540 -11.021849 2.526022 -1.809130 3.257744 # mm3->(c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
angle_coeff 6 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 7 class2/p6 123.490000 111.075360 -89.098091 20.419778 -14.624589 26.334856 # mm3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 7 class2/p6 bb 14.244120 1.278000 1.278000
|
||||
angle_coeff 7 class2/p6 ba 3.309240 3.309240 1.278000 1.278000
|
||||
angle_coeff 1 class2/p6 127.050000 53.307540 -42.760159 9.799907 -7.018666 12.638684 # mm3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
angle_coeff 1 class2/p6 bb 71.796120 1.394000 1.394000
|
||||
angle_coeff 1 class2/p6 ba 6.762360 6.762360 1.394000 1.394000
|
||||
angle_coeff 4 class2/p6 84.336000 29.351520 -23.544055 5.395900 -3.864529 6.958951 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 4 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 4 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 3 class2/p6 120.350000 36.185820 -29.026127 6.652298 -4.764358 8.579296 # mm3->(c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
angle_coeff 3 class2/p6 bb 93.378120 1.394000 1.094000
|
||||
angle_coeff 3 class2/p6 ba -25.179000 53.523360 1.394000 1.094000
|
||||
angle_coeff 8 class2/p6 115.098000 79.493700 -63.765149 14.613896 -10.466432 18.847160 # mm3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 8 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 8 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 5 cosine/buck6d 1.978350 4 180.000000 #1.0 # fourier->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
|
||||
dihedral_style opls
|
||||
|
||||
dihedral_coeff 3 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 9 0.000000 4.528000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 8 0.000000 0.000000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 5 0.000000 1.741000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
dihedral_coeff 2 0.000000 6.316000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
dihedral_coeff 1 0.000000 4.379000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
dihedral_coeff 6 0.000000 0.000000 0.000000 0.609000 # cos4->(o2_c1cu1@co2,cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
dihedral_coeff 4 0.000000 0.000000 0.000000 0.000000 # cos3->(h1_c1@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 10 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 7 0.000000 0.000000 0.000000 0.000000 # cos3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
|
||||
improper_style inversion/harmonic
|
||||
|
||||
improper_coeff 1 18.776340 0.000000 # harm->(c3_c3@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c1o2@co2)|CuPW
|
||||
improper_coeff 3 41.005800 0.000000 # harm->(c3_c1o2@co2,c3_c3@ph,o2_c1cu1@co2,o2_c1cu1@co2)|CuPW
|
||||
improper_coeff 2 4.100580 0.000000 # harm->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
|
||||
special_bonds lj 0.0 0.0 1.0 coul 1.0 1.0 1.0
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD END --------------------------
|
||||
|
||||
run 0
|
||||
|
|
@ -0,0 +1,166 @@
|
|||
LAMMPS (17 Jan 2018)
|
||||
clear
|
||||
units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
read_data hkust1.data
|
||||
triclinic box = (0 0 0) to (26.4408 26.4408 26.4408) with tilt (0 0 0)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
624 atoms
|
||||
scanning bonds ...
|
||||
5 = max bonds/atom
|
||||
scanning angles ...
|
||||
10 = max angles/atom
|
||||
scanning dihedrals ...
|
||||
32 = max dihedrals/atom
|
||||
scanning impropers ...
|
||||
2 = max impropers/atom
|
||||
reading bonds ...
|
||||
792 bonds
|
||||
reading angles ...
|
||||
1536 angles
|
||||
reading dihedrals ...
|
||||
2688 dihedrals
|
||||
reading impropers ...
|
||||
288 impropers
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
20 = max # of 1-4 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
neighbor 2.0 bin
|
||||
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD ------------------------------
|
||||
|
||||
pair_style buck6d/coul/gauss/dsf 0.9000 12.0000
|
||||
|
||||
pair_coeff 1 1 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW
|
||||
pair_coeff 1 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 1 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 1 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 1 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 1 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 2 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 2 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 2 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 2 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 2 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 3 3 3680 4 32.805 0.10690769 0.9771554 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 3 4 14157.243 3.1914894 489.18197 2.5231391 0.45538909 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 3 5 6320.6076 3.6144578 103.50278 0.44181916 0.75109952 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 3 6 6157.8178 3.4883721 124.7792 0.72632262 0.73006542 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 4 4 54464 2.6548673 5679.4311 33.208515 0.34105936 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 4 5 24315.863 2.9411765 1371.5617 7.9168726 0.42457748 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 4 6 23689.598 2.8571429 1590.0769 11.87959 0.42066711 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 5 5 10856 3.2967033 308.7755 1.6022517 0.63272774 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 5 6 10576.399 3.1914894 365.45138 2.5231391 0.61999948 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 6 6 10304 3.0927835 429.89352 3.9170177 0.60800971 # buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
|
||||
bond_style hybrid class2 morse
|
||||
|
||||
bond_coeff 5 morse 50.000000 1.451345 1.914000 # morse->(cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 4 class2 2.536000 75.465060 -192.435903 286.248406 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2)|CuPW
|
||||
bond_coeff 3 class2 1.094000 390.634200 -996.117210 1481.724350 # mm3->(c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
bond_coeff 6 class2 1.278000 585.591600 -1493.258580 2221.222138 # mm3->(c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 1 class2 1.394000 509.335200 -1298.804760 1931.972080 # mm3->(c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
bond_coeff 2 class2 1.485000 360.635220 -919.619811 1367.934469 # mm3->(c3_c1o2@co2,c3_c3@ph)|CuPW
|
||||
|
||||
angle_style hybrid class2/p6 cosine/buck6d
|
||||
|
||||
angle_coeff 2 class2/p6 117.711000 57.408120 -46.049402 10.553745 -7.558563 13.610890 # mm3->(c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
angle_coeff 2 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 2 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 126.814000 13.740540 -11.021849 2.526022 -1.809130 3.257744 # mm3->(c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
angle_coeff 6 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 7 class2/p6 123.490000 111.075360 -89.098091 20.419778 -14.624589 26.334856 # mm3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 7 class2/p6 bb 14.244120 1.278000 1.278000
|
||||
angle_coeff 7 class2/p6 ba 3.309240 3.309240 1.278000 1.278000
|
||||
angle_coeff 1 class2/p6 127.050000 53.307540 -42.760159 9.799907 -7.018666 12.638684 # mm3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
angle_coeff 1 class2/p6 bb 71.796120 1.394000 1.394000
|
||||
angle_coeff 1 class2/p6 ba 6.762360 6.762360 1.394000 1.394000
|
||||
angle_coeff 4 class2/p6 84.336000 29.351520 -23.544055 5.395900 -3.864529 6.958951 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 4 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 4 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 3 class2/p6 120.350000 36.185820 -29.026127 6.652298 -4.764358 8.579296 # mm3->(c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
angle_coeff 3 class2/p6 bb 93.378120 1.394000 1.094000
|
||||
angle_coeff 3 class2/p6 ba -25.179000 53.523360 1.394000 1.094000
|
||||
angle_coeff 8 class2/p6 115.098000 79.493700 -63.765149 14.613896 -10.466432 18.847160 # mm3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 8 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 8 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 5 cosine/buck6d 1.978350 4 180.000000 # fourier->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
|
||||
dihedral_style opls
|
||||
|
||||
dihedral_coeff 3 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 9 0.000000 4.528000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 8 0.000000 0.000000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 5 0.000000 1.741000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
dihedral_coeff 2 0.000000 6.316000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
dihedral_coeff 1 0.000000 4.379000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
dihedral_coeff 6 0.000000 0.000000 0.000000 0.609000 # cos4->(o2_c1cu1@co2,cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
dihedral_coeff 4 0.000000 0.000000 0.000000 0.000000 # cos3->(h1_c1@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 10 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 7 0.000000 0.000000 0.000000 0.000000 # cos3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
|
||||
improper_style inversion/harmonic
|
||||
|
||||
improper_coeff 1 18.776340 0.000000 # harm->(c3_c3@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c1o2@co2)|CuPW
|
||||
improper_coeff 3 41.005800 0.000000 # harm->(c3_c1o2@co2,c3_c3@ph,o2_c1cu1@co2,o2_c1cu1@co2)|CuPW
|
||||
improper_coeff 2 4.100580 0.000000 # harm->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
|
||||
special_bonds lj 0.0 0.0 1.0 coul 1.0 1.0 1.0
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD END --------------------------
|
||||
|
||||
run 0
|
||||
WARNING: No fixes defined, atoms won't move (../verlet.cpp:55)
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 4 4 4
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair buck6d/coul/gauss/dsf, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/newton/tri
|
||||
stencil: half/bin/3d/newton/tri
|
||||
bin: standard
|
||||
WARNING: Inconsistent image flags (../domain.cpp:786)
|
||||
Per MPI rank memory allocation (min/avg/max) = 21.23 | 21.23 | 21.23 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -11833.81 343.7619 -11490.048 -5.8508834
|
||||
Loop time of 9.53674e-07 on 1 procs for 0 steps with 624 atoms
|
||||
|
||||
0.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Bond | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 9.537e-07 | | |100.00
|
||||
|
||||
Nlocal: 624 ave 624 max 624 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 4464 ave 4464 max 4464 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 115368 ave 115368 max 115368 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 115368
|
||||
Ave neighs/atom = 184.885
|
||||
Ave special neighs/atom = 7.46154
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -0,0 +1,166 @@
|
|||
LAMMPS (17 Jan 2018)
|
||||
clear
|
||||
units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
read_data hkust1.data
|
||||
triclinic box = (0 0 0) to (26.4408 26.4408 26.4408) with tilt (0 0 0)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
624 atoms
|
||||
scanning bonds ...
|
||||
5 = max bonds/atom
|
||||
scanning angles ...
|
||||
10 = max angles/atom
|
||||
scanning dihedrals ...
|
||||
32 = max dihedrals/atom
|
||||
scanning impropers ...
|
||||
2 = max impropers/atom
|
||||
reading bonds ...
|
||||
792 bonds
|
||||
reading angles ...
|
||||
1536 angles
|
||||
reading dihedrals ...
|
||||
2688 dihedrals
|
||||
reading impropers ...
|
||||
288 impropers
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
20 = max # of 1-4 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
neighbor 2.0 bin
|
||||
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD ------------------------------
|
||||
|
||||
pair_style buck6d/coul/gauss/dsf 0.9000 12.0000
|
||||
|
||||
pair_coeff 1 1 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW
|
||||
pair_coeff 1 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 1 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 1 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 1 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 1 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 2 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 2 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 2 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 2 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 2 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 3 3 3680 4 32.805 0.10690769 0.9771554 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 3 4 14157.243 3.1914894 489.18197 2.5231391 0.45538909 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 3 5 6320.6076 3.6144578 103.50278 0.44181916 0.75109952 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 3 6 6157.8178 3.4883721 124.7792 0.72632262 0.73006542 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 4 4 54464 2.6548673 5679.4311 33.208515 0.34105936 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 4 5 24315.863 2.9411765 1371.5617 7.9168726 0.42457748 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 4 6 23689.598 2.8571429 1590.0769 11.87959 0.42066711 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 5 5 10856 3.2967033 308.7755 1.6022517 0.63272774 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 5 6 10576.399 3.1914894 365.45138 2.5231391 0.61999948 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 6 6 10304 3.0927835 429.89352 3.9170177 0.60800971 # buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
|
||||
bond_style hybrid class2 morse
|
||||
|
||||
bond_coeff 5 morse 50.000000 1.451345 1.914000 # morse->(cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 4 class2 2.536000 75.465060 -192.435903 286.248406 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2)|CuPW
|
||||
bond_coeff 3 class2 1.094000 390.634200 -996.117210 1481.724350 # mm3->(c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
bond_coeff 6 class2 1.278000 585.591600 -1493.258580 2221.222138 # mm3->(c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 1 class2 1.394000 509.335200 -1298.804760 1931.972080 # mm3->(c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
bond_coeff 2 class2 1.485000 360.635220 -919.619811 1367.934469 # mm3->(c3_c1o2@co2,c3_c3@ph)|CuPW
|
||||
|
||||
angle_style hybrid class2/p6 cosine/buck6d
|
||||
|
||||
angle_coeff 2 class2/p6 117.711000 57.408120 -46.049402 10.553745 -7.558563 13.610890 # mm3->(c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
angle_coeff 2 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 2 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 126.814000 13.740540 -11.021849 2.526022 -1.809130 3.257744 # mm3->(c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
angle_coeff 6 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 7 class2/p6 123.490000 111.075360 -89.098091 20.419778 -14.624589 26.334856 # mm3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 7 class2/p6 bb 14.244120 1.278000 1.278000
|
||||
angle_coeff 7 class2/p6 ba 3.309240 3.309240 1.278000 1.278000
|
||||
angle_coeff 1 class2/p6 127.050000 53.307540 -42.760159 9.799907 -7.018666 12.638684 # mm3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
angle_coeff 1 class2/p6 bb 71.796120 1.394000 1.394000
|
||||
angle_coeff 1 class2/p6 ba 6.762360 6.762360 1.394000 1.394000
|
||||
angle_coeff 4 class2/p6 84.336000 29.351520 -23.544055 5.395900 -3.864529 6.958951 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 4 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 4 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 3 class2/p6 120.350000 36.185820 -29.026127 6.652298 -4.764358 8.579296 # mm3->(c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
angle_coeff 3 class2/p6 bb 93.378120 1.394000 1.094000
|
||||
angle_coeff 3 class2/p6 ba -25.179000 53.523360 1.394000 1.094000
|
||||
angle_coeff 8 class2/p6 115.098000 79.493700 -63.765149 14.613896 -10.466432 18.847160 # mm3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 8 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 8 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 5 cosine/buck6d 1.978350 4 180.000000 # fourier->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
|
||||
dihedral_style opls
|
||||
|
||||
dihedral_coeff 3 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 9 0.000000 4.528000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 8 0.000000 0.000000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 5 0.000000 1.741000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
dihedral_coeff 2 0.000000 6.316000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
dihedral_coeff 1 0.000000 4.379000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
dihedral_coeff 6 0.000000 0.000000 0.000000 0.609000 # cos4->(o2_c1cu1@co2,cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
dihedral_coeff 4 0.000000 0.000000 0.000000 0.000000 # cos3->(h1_c1@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 10 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 7 0.000000 0.000000 0.000000 0.000000 # cos3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
|
||||
improper_style inversion/harmonic
|
||||
|
||||
improper_coeff 1 18.776340 0.000000 # harm->(c3_c3@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c1o2@co2)|CuPW
|
||||
improper_coeff 3 41.005800 0.000000 # harm->(c3_c1o2@co2,c3_c3@ph,o2_c1cu1@co2,o2_c1cu1@co2)|CuPW
|
||||
improper_coeff 2 4.100580 0.000000 # harm->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
|
||||
special_bonds lj 0.0 0.0 1.0 coul 1.0 1.0 1.0
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD END --------------------------
|
||||
|
||||
run 0
|
||||
WARNING: No fixes defined, atoms won't move (../verlet.cpp:55)
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 4 4 4
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair buck6d/coul/gauss/dsf, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/newton/tri
|
||||
stencil: half/bin/3d/newton/tri
|
||||
bin: standard
|
||||
WARNING: Inconsistent image flags (../domain.cpp:786)
|
||||
Per MPI rank memory allocation (min/avg/max) = 20.68 | 20.68 | 20.68 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -11833.81 343.7619 -11490.048 -5.8508834
|
||||
Loop time of 2.20537e-06 on 4 procs for 0 steps with 624 atoms
|
||||
|
||||
0.0% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Bond | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 2.205e-06 | | |100.00
|
||||
|
||||
Nlocal: 156 ave 156 max 156 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 2718 ave 2718 max 2718 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 28842 ave 28870 max 28814 min
|
||||
Histogram: 1 1 0 0 0 0 0 0 1 1
|
||||
|
||||
Total # of neighbors = 115368
|
||||
Ave neighs/atom = 184.885
|
||||
Ave special neighs/atom = 7.46154
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -0,0 +1,175 @@
|
|||
LAMMPS (17 Jan 2018)
|
||||
clear
|
||||
units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
read_data hkust1.data
|
||||
triclinic box = (0 0 0) to (26.4408 26.4408 26.4408) with tilt (0 0 0)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
624 atoms
|
||||
scanning bonds ...
|
||||
5 = max bonds/atom
|
||||
scanning angles ...
|
||||
10 = max angles/atom
|
||||
scanning dihedrals ...
|
||||
32 = max dihedrals/atom
|
||||
scanning impropers ...
|
||||
2 = max impropers/atom
|
||||
reading bonds ...
|
||||
792 bonds
|
||||
reading angles ...
|
||||
1536 angles
|
||||
reading dihedrals ...
|
||||
2688 dihedrals
|
||||
reading impropers ...
|
||||
288 impropers
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
20 = max # of 1-4 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
neighbor 2.0 bin
|
||||
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD ------------------------------
|
||||
kspace_style ewald 1.0e-6
|
||||
|
||||
pair_style buck6d/coul/gauss/long 0.9 0.9 12.0000
|
||||
|
||||
pair_coeff 1 1 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW
|
||||
pair_coeff 1 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 1 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 1 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 1 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 1 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 2 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 2 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 2 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 2 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 2 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 3 3 3680 4 32.805 0.10690769 0.9771554 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 3 4 14157.243 3.1914894 489.18197 2.5231391 0.45538909 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 3 5 6320.6076 3.6144578 103.50278 0.44181916 0.75109952 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 3 6 6157.8178 3.4883721 124.7792 0.72632262 0.73006542 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 4 4 54464 2.6548673 5679.4311 33.208515 0.34105936 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 4 5 24315.863 2.9411765 1371.5617 7.9168726 0.42457748 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 4 6 23689.598 2.8571429 1590.0769 11.87959 0.42066711 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 5 5 10856 3.2967033 308.7755 1.6022517 0.63272774 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 5 6 10576.399 3.1914894 365.45138 2.5231391 0.61999948 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 6 6 10304 3.0927835 429.89352 3.9170177 0.60800971 # buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
|
||||
bond_style hybrid class2 morse
|
||||
|
||||
bond_coeff 5 morse 50.000000 1.451345 1.914000 # morse->(cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 4 class2 2.536000 75.465060 -192.435903 286.248406 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2)|CuPW
|
||||
bond_coeff 3 class2 1.094000 390.634200 -996.117210 1481.724350 # mm3->(c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
bond_coeff 6 class2 1.278000 585.591600 -1493.258580 2221.222138 # mm3->(c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 1 class2 1.394000 509.335200 -1298.804760 1931.972080 # mm3->(c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
bond_coeff 2 class2 1.485000 360.635220 -919.619811 1367.934469 # mm3->(c3_c1o2@co2,c3_c3@ph)|CuPW
|
||||
|
||||
angle_style hybrid class2/p6 cosine/buck6d
|
||||
|
||||
angle_coeff 2 class2/p6 117.711000 57.408120 -46.049402 10.553745 -7.558563 13.610890 # mm3->(c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
angle_coeff 2 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 2 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 126.814000 13.740540 -11.021849 2.526022 -1.809130 3.257744 # mm3->(c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
angle_coeff 6 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 7 class2/p6 123.490000 111.075360 -89.098091 20.419778 -14.624589 26.334856 # mm3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 7 class2/p6 bb 14.244120 1.278000 1.278000
|
||||
angle_coeff 7 class2/p6 ba 3.309240 3.309240 1.278000 1.278000
|
||||
angle_coeff 1 class2/p6 127.050000 53.307540 -42.760159 9.799907 -7.018666 12.638684 # mm3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
angle_coeff 1 class2/p6 bb 71.796120 1.394000 1.394000
|
||||
angle_coeff 1 class2/p6 ba 6.762360 6.762360 1.394000 1.394000
|
||||
angle_coeff 4 class2/p6 84.336000 29.351520 -23.544055 5.395900 -3.864529 6.958951 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 4 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 4 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 3 class2/p6 120.350000 36.185820 -29.026127 6.652298 -4.764358 8.579296 # mm3->(c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
angle_coeff 3 class2/p6 bb 93.378120 1.394000 1.094000
|
||||
angle_coeff 3 class2/p6 ba -25.179000 53.523360 1.394000 1.094000
|
||||
angle_coeff 8 class2/p6 115.098000 79.493700 -63.765149 14.613896 -10.466432 18.847160 # mm3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 8 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 8 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 5 cosine/buck6d 1.978350 4 180.000000 #1.0 # fourier->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
|
||||
dihedral_style opls
|
||||
|
||||
dihedral_coeff 3 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 9 0.000000 4.528000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 8 0.000000 0.000000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 5 0.000000 1.741000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
dihedral_coeff 2 0.000000 6.316000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
dihedral_coeff 1 0.000000 4.379000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
dihedral_coeff 6 0.000000 0.000000 0.000000 0.609000 # cos4->(o2_c1cu1@co2,cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
dihedral_coeff 4 0.000000 0.000000 0.000000 0.000000 # cos3->(h1_c1@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 10 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 7 0.000000 0.000000 0.000000 0.000000 # cos3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
|
||||
improper_style inversion/harmonic
|
||||
|
||||
improper_coeff 1 18.776340 0.000000 # harm->(c3_c3@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c1o2@co2)|CuPW
|
||||
improper_coeff 3 41.005800 0.000000 # harm->(c3_c1o2@co2,c3_c3@ph,o2_c1cu1@co2,o2_c1cu1@co2)|CuPW
|
||||
improper_coeff 2 4.100580 0.000000 # harm->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
|
||||
special_bonds lj 0.0 0.0 1.0 coul 1.0 1.0 1.0
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD END --------------------------
|
||||
|
||||
run 0
|
||||
WARNING: No fixes defined, atoms won't move (../verlet.cpp:55)
|
||||
Ewald initialization ...
|
||||
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
|
||||
G vector (1/distance) = 0.267593
|
||||
estimated absolute RMS force accuracy = 0.000333665
|
||||
estimated relative force accuracy = 1.00482e-06
|
||||
KSpace vectors: actual max1d max3d = 1054 8 2456
|
||||
kxmax kymax kzmax = 8 8 8
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 4 4 4
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair buck6d/coul/gauss/long, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/newton/tri
|
||||
stencil: half/bin/3d/newton/tri
|
||||
bin: standard
|
||||
WARNING: Inconsistent image flags (../domain.cpp:786)
|
||||
Per MPI rank memory allocation (min/avg/max) = 34.64 | 34.64 | 34.64 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -16541.109 343.7619 -16197.347 -629.64956
|
||||
Loop time of 9.53674e-07 on 1 procs for 0 steps with 624 atoms
|
||||
|
||||
0.0% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Bond | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Kspace | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 9.537e-07 | | |100.00
|
||||
|
||||
Nlocal: 624 ave 624 max 624 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 4464 ave 4464 max 4464 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 115368 ave 115368 max 115368 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 115368
|
||||
Ave neighs/atom = 184.885
|
||||
Ave special neighs/atom = 7.46154
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -0,0 +1,175 @@
|
|||
LAMMPS (17 Jan 2018)
|
||||
clear
|
||||
units real
|
||||
boundary p p p
|
||||
atom_style full
|
||||
read_data hkust1.data
|
||||
triclinic box = (0 0 0) to (26.4408 26.4408 26.4408) with tilt (0 0 0)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
624 atoms
|
||||
scanning bonds ...
|
||||
5 = max bonds/atom
|
||||
scanning angles ...
|
||||
10 = max angles/atom
|
||||
scanning dihedrals ...
|
||||
32 = max dihedrals/atom
|
||||
scanning impropers ...
|
||||
2 = max impropers/atom
|
||||
reading bonds ...
|
||||
792 bonds
|
||||
reading angles ...
|
||||
1536 angles
|
||||
reading dihedrals ...
|
||||
2688 dihedrals
|
||||
reading impropers ...
|
||||
288 impropers
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
20 = max # of 1-4 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
neighbor 2.0 bin
|
||||
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD ------------------------------
|
||||
kspace_style ewald 1.0e-6
|
||||
|
||||
pair_style buck6d/coul/gauss/long 0.9 0.9 12.0000
|
||||
|
||||
pair_coeff 1 1 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW
|
||||
pair_coeff 1 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 1 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 1 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 1 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 1 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c3@ph)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 2 2 10304 3.0612245 457.17971 4.5218516 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene
|
||||
pair_coeff 2 3 6157.8178 3.4682081 129.19572 0.78772886 0.73006542 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 2 4 23689.598 2.8436019 1636.052 12.696549 0.42066711 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 2 5 10576.399 3.1746032 377.27092 2.7176691 0.61999948 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 2 6 10304 3.0769231 443.36171 4.2093581 0.60800971 # buck6d->(c3_c2h1@ph)|benzene/gaussian->(c3_c2h1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 3 3 3680 4 32.805 0.10690769 0.9771554 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene
|
||||
pair_coeff 3 4 14157.243 3.1914894 489.18197 2.5231391 0.45538909 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 3 5 6320.6076 3.6144578 103.50278 0.44181916 0.75109952 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 3 6 6157.8178 3.4883721 124.7792 0.72632262 0.73006542 # buck6d->(h1_c1@ph)|benzene/gaussian->(h1_c1@ph)|benzene <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 4 4 54464 2.6548673 5679.4311 33.208515 0.34105936 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW
|
||||
pair_coeff 4 5 24315.863 2.9411765 1371.5617 7.9168726 0.42457748 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 4 6 23689.598 2.8571429 1590.0769 11.87959 0.42066711 # buck6d->(cu5_cu1o4@cu2)|CuPW/gaussian->(cu5_cu1o4@cu2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 5 5 10856 3.2967033 308.7755 1.6022517 0.63272774 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW
|
||||
pair_coeff 5 6 10576.399 3.1914894 365.45138 2.5231391 0.61999948 # buck6d->(o2_c1cu1@co2)|CuPW/gaussian->(o2_c1cu1@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
pair_coeff 6 6 10304 3.0927835 429.89352 3.9170177 0.60800971 # buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW <--> buck6d->(c3_c1o2@co2)|CuPW/gaussian->(c3_c1o2@co2)|CuPW
|
||||
|
||||
bond_style hybrid class2 morse
|
||||
|
||||
bond_coeff 5 morse 50.000000 1.451345 1.914000 # morse->(cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 4 class2 2.536000 75.465060 -192.435903 286.248406 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2)|CuPW
|
||||
bond_coeff 3 class2 1.094000 390.634200 -996.117210 1481.724350 # mm3->(c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
bond_coeff 6 class2 1.278000 585.591600 -1493.258580 2221.222138 # mm3->(c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
bond_coeff 1 class2 1.394000 509.335200 -1298.804760 1931.972080 # mm3->(c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
bond_coeff 2 class2 1.485000 360.635220 -919.619811 1367.934469 # mm3->(c3_c1o2@co2,c3_c3@ph)|CuPW
|
||||
|
||||
angle_style hybrid class2/p6 cosine/buck6d
|
||||
|
||||
angle_coeff 2 class2/p6 117.711000 57.408120 -46.049402 10.553745 -7.558563 13.610890 # mm3->(c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
angle_coeff 2 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 2 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 126.814000 13.740540 -11.021849 2.526022 -1.809130 3.257744 # mm3->(c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
angle_coeff 6 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 6 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 7 class2/p6 123.490000 111.075360 -89.098091 20.419778 -14.624589 26.334856 # mm3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 7 class2/p6 bb 14.244120 1.278000 1.278000
|
||||
angle_coeff 7 class2/p6 ba 3.309240 3.309240 1.278000 1.278000
|
||||
angle_coeff 1 class2/p6 127.050000 53.307540 -42.760159 9.799907 -7.018666 12.638684 # mm3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
angle_coeff 1 class2/p6 bb 71.796120 1.394000 1.394000
|
||||
angle_coeff 1 class2/p6 ba 6.762360 6.762360 1.394000 1.394000
|
||||
angle_coeff 4 class2/p6 84.336000 29.351520 -23.544055 5.395900 -3.864529 6.958951 # mm3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 4 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 4 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 3 class2/p6 120.350000 36.185820 -29.026127 6.652298 -4.764358 8.579296 # mm3->(c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
angle_coeff 3 class2/p6 bb 93.378120 1.394000 1.094000
|
||||
angle_coeff 3 class2/p6 ba -25.179000 53.523360 1.394000 1.094000
|
||||
angle_coeff 8 class2/p6 115.098000 79.493700 -63.765149 14.613896 -10.466432 18.847160 # mm3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2)|CuPW
|
||||
angle_coeff 8 class2/p6 bb 0.0 1.0 1.0
|
||||
angle_coeff 8 class2/p6 ba 0.0 0.0 1.0 1.0
|
||||
angle_coeff 5 cosine/buck6d 1.978350 4 180.000000 #1.0 # fourier->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
|
||||
dihedral_style opls
|
||||
|
||||
dihedral_coeff 3 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 9 0.000000 4.528000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 8 0.000000 0.000000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 5 0.000000 1.741000 0.000000 0.000000 # cos3->(o2_c1cu1@co2,c3_c1o2@co2,c3_c3@ph,c3_c2h1@ph)|CuPW
|
||||
dihedral_coeff 2 0.000000 6.316000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
dihedral_coeff 1 0.000000 4.379000 0.000000 0.000000 # cos3->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph)|benzene
|
||||
dihedral_coeff 6 0.000000 0.000000 0.000000 0.609000 # cos4->(o2_c1cu1@co2,cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2)|CuPW
|
||||
dihedral_coeff 4 0.000000 0.000000 0.000000 0.000000 # cos3->(h1_c1@ph,c3_c2h1@ph,c3_c3@ph,c3_c1o2@co2)|CuPW
|
||||
dihedral_coeff 10 0.000000 0.000000 0.000000 0.000000 # cos3->(c3_c3@ph,c3_c1o2@co2,o2_c1cu1@co2,cu5_cu1o4@cu2)|CuPW
|
||||
dihedral_coeff 7 0.000000 0.000000 0.000000 0.000000 # cos3->(cu5_cu1o4@cu2,cu5_cu1o4@cu2,o2_c1cu1@co2,c3_c1o2@co2)|CuPW
|
||||
|
||||
improper_style inversion/harmonic
|
||||
|
||||
improper_coeff 1 18.776340 0.000000 # harm->(c3_c3@ph,c3_c2h1@ph,c3_c2h1@ph,c3_c1o2@co2)|CuPW
|
||||
improper_coeff 3 41.005800 0.000000 # harm->(c3_c1o2@co2,c3_c3@ph,o2_c1cu1@co2,o2_c1cu1@co2)|CuPW
|
||||
improper_coeff 2 4.100580 0.000000 # harm->(c3_c2h1@ph,c3_c2h1@ph,c3_c2h1@ph,h1_c1@ph)|benzene
|
||||
|
||||
special_bonds lj 0.0 0.0 1.0 coul 1.0 1.0 1.0
|
||||
5 = max # of 1-2 neighbors
|
||||
8 = max # of 1-3 neighbors
|
||||
17 = max # of special neighbors
|
||||
|
||||
# ------------------------ MOF-FF FORCE FIELD END --------------------------
|
||||
|
||||
run 0
|
||||
WARNING: No fixes defined, atoms won't move (../verlet.cpp:55)
|
||||
Ewald initialization ...
|
||||
WARNING: Using 12-bit tables for long-range coulomb (../kspace.cpp:321)
|
||||
G vector (1/distance) = 0.267593
|
||||
estimated absolute RMS force accuracy = 0.000333665
|
||||
estimated relative force accuracy = 1.00482e-06
|
||||
KSpace vectors: actual max1d max3d = 1054 8 2456
|
||||
kxmax kymax kzmax = 8 8 8
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 4 4 4
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair buck6d/coul/gauss/long, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/newton/tri
|
||||
stencil: half/bin/3d/newton/tri
|
||||
bin: standard
|
||||
WARNING: Inconsistent image flags (../domain.cpp:786)
|
||||
Per MPI rank memory allocation (min/avg/max) = 34.1 | 34.1 | 34.1 Mbytes
|
||||
Step Temp E_pair E_mol TotEng Press
|
||||
0 0 -16541.109 343.7619 -16197.347 -629.64956
|
||||
Loop time of 6.13928e-06 on 4 procs for 0 steps with 624 atoms
|
||||
|
||||
0.0% CPU use with 4 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Bond | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Kspace | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Other | | 6.139e-06 | | |100.00
|
||||
|
||||
Nlocal: 156 ave 156 max 156 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 2718 ave 2718 max 2718 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 28842 ave 28870 max 28814 min
|
||||
Histogram: 1 1 0 0 0 0 0 0 1 1
|
||||
|
||||
Total # of neighbors = 115368
|
||||
Ave neighs/atom = 184.885
|
||||
Ave special neighs/atom = 7.46154
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:00
|
||||
|
|
@ -134,6 +134,10 @@
|
|||
/angle_charmm.h
|
||||
/angle_class2.cpp
|
||||
/angle_class2.h
|
||||
/angle_class2_p6.cpp
|
||||
/angle_class2_p6.h
|
||||
/angle_cosine_buck6d.cpp
|
||||
/angle_cosine_buck6d.h
|
||||
/angle_cosine.cpp
|
||||
/angle_cosine.h
|
||||
/angle_cosine_delta.cpp
|
||||
|
|
@ -604,6 +608,8 @@
|
|||
/improper_harmonic.h
|
||||
/improper_hybrid.cpp
|
||||
/improper_hybrid.h
|
||||
/improper_inversion_harmonic.cpp
|
||||
/improper_inversion_harmonic.h
|
||||
/improper_ring.cpp
|
||||
/improper_ring.h
|
||||
/improper_umbrella.cpp
|
||||
|
|
@ -649,6 +655,10 @@
|
|||
/pair_buck_coul.h
|
||||
/pair_buck_long_coul_long.cpp
|
||||
/pair_buck_long_coul_long.h
|
||||
/pair_buck6d_coul_gauss_dsf.cpp
|
||||
/pair_buck6d_coul_gauss_dsf.h
|
||||
/pair_buck6d_coul_gauss_long.cpp
|
||||
/pair_buck6d_coul_gauss_long.h
|
||||
/pair_cdeam.cpp
|
||||
/pair_cdeam.h
|
||||
/pair_cg_cmm.cpp
|
||||
|
|
|
|||
|
|
@ -0,0 +1,32 @@
|
|||
This Package implements pair, angle and improper styles needed to employ
|
||||
the MOF-FF force field by Schmid and coworkers with LAMMPS.
|
||||
MOF-FF is a first principles derived force field with the primary aim
|
||||
to simulate MOFs and related porous framework materials, using spherical
|
||||
Gaussian charges. It is described in S. Bureekaew et al., Phys. Stat. Sol. B
|
||||
2013, 250, 1128-1141.
|
||||
For the usage of MOF-FF see the example in the example directory as
|
||||
well as the "MOF+" website (https://www.mofplus.org/content/show/MOF-FF).
|
||||
The package provides the following features:
|
||||
|
||||
* a dispersion damped Buckhingham potential with spherical Gaussian type
|
||||
charges (dsf and long-range treatment of charges)
|
||||
|
||||
* a modified angle/class2 including 6th order polynomial
|
||||
|
||||
* a modified angle/cosine style which adds a dispersion damped Buckhingham
|
||||
1-3 interaction analog to the dihedral/charmm style
|
||||
|
||||
* an improper style following the Wilson-Decius definition of the
|
||||
out-of-plane angle
|
||||
|
||||
See the file doc/drude_tutorial.html for getting started.
|
||||
See the doc pages for "pair_style buck6d/coul/gauss", "anlge_style class2",
|
||||
"angle_style cosine/buck6d", and "improper_style inversion/harmonic"
|
||||
commands to get started. Also see the above mentioned website and
|
||||
literature for further documentation about the force field.
|
||||
|
||||
There are example scripts for using this force field in examples/USER/mofff.
|
||||
|
||||
The creators of this package are Hendrik Heenen (hendrik.heenen at mytum.de)
|
||||
and Rochus Schmid (rochus.schmid at rub.de). Contact them directly if you
|
||||
have questions.
|
||||
|
|
@ -0,0 +1,487 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Hendrik Heenen (Technical University of Munich)
|
||||
and Rochus Schmid (Ruhr-Universitaet Bochum)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include "angle_class2_p6.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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleClass2P6::AngleClass2P6(LAMMPS *lmp) : Angle(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleClass2P6::~AngleClass2P6()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(setflag_a);
|
||||
memory->destroy(setflag_bb);
|
||||
memory->destroy(setflag_ba);
|
||||
|
||||
memory->destroy(theta0);
|
||||
memory->destroy(k2);
|
||||
memory->destroy(k3);
|
||||
memory->destroy(k4);
|
||||
memory->destroy(k5);
|
||||
memory->destroy(k6);
|
||||
|
||||
memory->destroy(bb_k);
|
||||
memory->destroy(bb_r1);
|
||||
memory->destroy(bb_r2);
|
||||
|
||||
memory->destroy(ba_k1);
|
||||
memory->destroy(ba_k2);
|
||||
memory->destroy(ba_r1);
|
||||
memory->destroy(ba_r2);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void AngleClass2P6::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,dtheta5,dtheta6,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;
|
||||
dtheta5 = dtheta4*dtheta;
|
||||
dtheta6 = dtheta5*dtheta;
|
||||
|
||||
de_angle = 2.0*k2[type]*dtheta + 3.0*k3[type]*dtheta2 +
|
||||
4.0*k4[type]*dtheta3 + 5.0*k5[type]*dtheta4 +
|
||||
6.0*k6[type]*dtheta5;
|
||||
|
||||
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;
|
||||
|
||||
if (eflag) eangle = k2[type]*dtheta2 + k3[type]*dtheta3 + k4[type]*dtheta4
|
||||
+ k5[type]*dtheta5 + k6[type]*dtheta6;
|
||||
|
||||
// force & energy for bond-bond term
|
||||
|
||||
dr1 = r1 - bb_r1[type];
|
||||
dr2 = r2 - bb_r2[type];
|
||||
tk1 = bb_k[type] * dr1;
|
||||
tk2 = bb_k[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 += bb_k[type]*dr1*dr2;
|
||||
|
||||
// force & energy for bond-angle term
|
||||
|
||||
aa1 = s * dr1 * ba_k1[type];
|
||||
aa2 = s * dr2 * ba_k2[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 = ba_k1[type] * dtheta / r1;
|
||||
b2 = ba_k2[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 += ba_k1[type]*dr1*dtheta + ba_k2[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 AngleClass2P6::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
|
||||
memory->create(theta0,n+1,"angle:theta0");
|
||||
memory->create(k2,n+1,"angle:k2");
|
||||
memory->create(k3,n+1,"angle:k3");
|
||||
memory->create(k4,n+1,"angle:k4");
|
||||
memory->create(k5,n+1,"angle:k5");
|
||||
memory->create(k6,n+1,"angle:k6");
|
||||
|
||||
memory->create(bb_k,n+1,"angle:bb_k");
|
||||
memory->create(bb_r1,n+1,"angle:bb_r1");
|
||||
memory->create(bb_r2,n+1,"angle:bb_r2");
|
||||
|
||||
memory->create(ba_k1,n+1,"angle:ba_k1");
|
||||
memory->create(ba_k2,n+1,"angle:ba_k2");
|
||||
memory->create(ba_r1,n+1,"angle:ba_r1");
|
||||
memory->create(ba_r2,n+1,"angle:ba_r2");
|
||||
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
memory->create(setflag_a,n+1,"angle:setflag_a");
|
||||
memory->create(setflag_bb,n+1,"angle:setflag_bb");
|
||||
memory->create(setflag_ba,n+1,"angle:setflag_ba");
|
||||
for (int i = 1; i <= n; i++)
|
||||
setflag[i] = setflag_a[i] = setflag_bb[i] = setflag_ba[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more types
|
||||
arg1 = "bb" -> BondBond coeffs
|
||||
arg1 = "ba" -> BondAngle coeffs
|
||||
else -> Angle coeffs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleClass2P6::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 2) 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;
|
||||
|
||||
if (strcmp(arg[1],"bb") == 0) {
|
||||
if (narg != 5) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
|
||||
double bb_k_one = force->numeric(FLERR,arg[2]);
|
||||
double bb_r1_one = force->numeric(FLERR,arg[3]);
|
||||
double bb_r2_one = force->numeric(FLERR,arg[4]);
|
||||
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
bb_k[i] = bb_k_one;
|
||||
bb_r1[i] = bb_r1_one;
|
||||
bb_r2[i] = bb_r2_one;
|
||||
setflag_bb[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
} else if (strcmp(arg[1],"ba") == 0) {
|
||||
if (narg != 6) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
|
||||
double ba_k1_one = force->numeric(FLERR,arg[2]);
|
||||
double ba_k2_one = force->numeric(FLERR,arg[3]);
|
||||
double ba_r1_one = force->numeric(FLERR,arg[4]);
|
||||
double ba_r2_one = force->numeric(FLERR,arg[5]);
|
||||
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
ba_k1[i] = ba_k1_one;
|
||||
ba_k2[i] = ba_k2_one;
|
||||
ba_r1[i] = ba_r1_one;
|
||||
ba_r2[i] = ba_r2_one;
|
||||
setflag_ba[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
} else {
|
||||
if (narg != 7) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
|
||||
double theta0_one = force->numeric(FLERR,arg[1]);
|
||||
double k2_one = force->numeric(FLERR,arg[2]);
|
||||
double k3_one = force->numeric(FLERR,arg[3]);
|
||||
double k4_one = force->numeric(FLERR,arg[4]);
|
||||
double k5_one = force->numeric(FLERR,arg[5]);
|
||||
double k6_one = force->numeric(FLERR,arg[6]);
|
||||
|
||||
// convert theta0 from degrees to radians
|
||||
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
theta0[i] = theta0_one/180.0 * MY_PI;
|
||||
k2[i] = k2_one;
|
||||
k3[i] = k3_one;
|
||||
k4[i] = k4_one;
|
||||
k5[i] = k5_one;
|
||||
k6[i] = k6_one;
|
||||
setflag_a[i] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
|
||||
for (int i = ilo; i <= ihi; i++)
|
||||
if (setflag_a[i] == 1 && setflag_bb[i] == 1 && setflag_ba[i] == 1)
|
||||
setflag[i] = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleClass2P6::equilibrium_angle(int i)
|
||||
{
|
||||
return theta0[i];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes out coeffs to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleClass2P6::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k2[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k3[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k4[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k5[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k6[1],sizeof(double),atom->nangletypes,fp);
|
||||
|
||||
fwrite(&bb_k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&bb_r1[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&bb_r2[1],sizeof(double),atom->nangletypes,fp);
|
||||
|
||||
fwrite(&ba_k1[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&ba_k2[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&ba_r1[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&ba_r2[1],sizeof(double),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads coeffs from restart file, bcasts them
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleClass2P6::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&theta0[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&k2[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&k3[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&k4[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&k5[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&k6[1],sizeof(double),atom->nangletypes,fp);
|
||||
|
||||
fread(&bb_k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&bb_r1[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&bb_r2[1],sizeof(double),atom->nangletypes,fp);
|
||||
|
||||
fread(&ba_k1[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&ba_k2[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&ba_r1[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&ba_r2[1],sizeof(double),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k3[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k4[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k5[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k6[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
|
||||
MPI_Bcast(&bb_k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&bb_r1[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&bb_r2[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
|
||||
MPI_Bcast(&ba_k1[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&ba_k2[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&ba_r1[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&ba_r2[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 AngleClass2P6::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->nangletypes; i++)
|
||||
fprintf(fp,"%d %g %g %g %g\n",
|
||||
i,theta0[i]/MY_PI*180.0,k2[i],k3[i],k4[i],k5[i],k6[i]);
|
||||
|
||||
fprintf(fp,"\nBondBond Coeffs\n\n");
|
||||
for (int i = 1; i <= atom->nangletypes; i++)
|
||||
fprintf(fp,"%d %g %g %g\n",i,bb_k[i],bb_r1[i],bb_r2[i]);
|
||||
|
||||
fprintf(fp,"\nBondAngle Coeffs\n\n");
|
||||
for (int i = 1; i <= atom->nangletypes; i++)
|
||||
fprintf(fp,"%d %g %g %g %g\n",i,ba_k1[i],ba_k2[i],ba_r1[i],ba_r2[i]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleClass2P6::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 dtheta5 = dtheta4*dtheta;
|
||||
double dtheta6 = dtheta5*dtheta;
|
||||
|
||||
double energy = k2[type]*dtheta2 + k3[type]*dtheta3 + k4[type]*dtheta4
|
||||
+ k5[type]*dtheta5 + k6[type]*dtheta6;
|
||||
|
||||
double dr1 = r1 - bb_r1[type];
|
||||
double dr2 = r2 - bb_r2[type];
|
||||
energy += bb_k[type]*dr1*dr2;
|
||||
|
||||
energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta;
|
||||
return energy;
|
||||
}
|
||||
|
|
@ -0,0 +1,60 @@
|
|||
/* -*- 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(class2/p6,AngleClass2P6)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_ANGLE_CLASS2_P6_H
|
||||
#define LMP_ANGLE_CLASS2_P6_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include "angle.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class AngleClass2P6 : public Angle {
|
||||
public:
|
||||
AngleClass2P6(class LAMMPS *);
|
||||
virtual ~AngleClass2P6();
|
||||
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,*k3,*k4,*k5,*k6;
|
||||
double *bb_k,*bb_r1,*bb_r2;
|
||||
double *ba_k1,*ba_k2,*ba_r1,*ba_r2;
|
||||
int *setflag_a,*setflag_bb,*setflag_ba;
|
||||
|
||||
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,386 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Hendrik Heenen (Technical University of Munich)
|
||||
and Rochus Schmid (Ruhr-Universitaet Bochum)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include "angle_cosine_buck6d.h"
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "math_const.h"
|
||||
#include "math_special.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using namespace MathSpecial;
|
||||
|
||||
#define SMALL 0.001
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCosineBuck6d::AngleCosineBuck6d(LAMMPS *lmp) : Angle(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCosineBuck6d::~AngleCosineBuck6d()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(k);
|
||||
memory->destroy(multiplicity);
|
||||
memory->destroy(th0);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosineBuck6d::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,i1,i2,i3,n,type,itype,jtype;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
|
||||
double tk;
|
||||
|
||||
// extra lj variables
|
||||
double delx3,dely3,delz3,rsq3,r3;
|
||||
double rexp,r32inv,r36inv,r314inv,forcebuck6d,fpair;
|
||||
double term1,term2,term3,term4,term5,ebuck6d,evdwl;
|
||||
double rcu,rqu,sme,smf;
|
||||
|
||||
eangle = evdwl = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 0;
|
||||
|
||||
// insure pair->ev_tally() will use 1-3 virial contribution
|
||||
|
||||
if (vflag_global == 2)
|
||||
force->pair->vflag_either = force->pair->vflag_global = 1;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int **anglelist = neighbor->anglelist;
|
||||
int nanglelist = neighbor->nanglelist;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_pair = force->newton_pair;
|
||||
int newton_bond = force->newton_bond;
|
||||
int *atomtype = atom->type;
|
||||
|
||||
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);
|
||||
|
||||
// c = cosine of angle
|
||||
|
||||
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
|
||||
|
||||
// explicit lj-contribution
|
||||
|
||||
itype = atomtype[i1];
|
||||
jtype = atomtype[i3];
|
||||
|
||||
delx3 = x[i1][0] - x[i3][0];
|
||||
dely3 = x[i1][1] - x[i3][1];
|
||||
delz3 = x[i1][2] - x[i3][2];
|
||||
rsq3 = delx3*delx3 + dely3*dely3 + delz3*delz3;
|
||||
|
||||
if (rsq3 < cut_ljsq[itype][jtype]) {
|
||||
r3 = sqrt(rsq3);
|
||||
r32inv = 1.0/rsq3;
|
||||
r36inv = r32inv*r32inv*r32inv;
|
||||
r314inv = r36inv*r36inv*r32inv;
|
||||
rexp = exp(-r3*buck6d2[itype][jtype]);
|
||||
term1 = buck6d3[itype][jtype]*r36inv;
|
||||
term2 = buck6d4[itype][jtype]*r314inv;
|
||||
term3 = term2*term2;
|
||||
term4 = 1.0/(1.0 + term2);
|
||||
term5 = 1.0/(1.0 + 2.0*term2 + term3);
|
||||
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r3*rexp;
|
||||
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
|
||||
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
|
||||
|
||||
// smoothing term
|
||||
if (rsq3 > rsmooth_sq[itype][jtype]) {
|
||||
rcu = r3*rsq3;
|
||||
rqu = rsq3*rsq3;
|
||||
sme = c5[itype][jtype]*rqu*r3 + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
|
||||
c2[itype][jtype]*rsq3 + c1[itype][jtype]*r3 + c0[itype][jtype];
|
||||
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
|
||||
3.0*c3[itype][jtype]*rsq3 + 2.0*c2[itype][jtype]*r3 + c1[itype][jtype];
|
||||
forcebuck6d = forcebuck6d*sme + ebuck6d*smf;
|
||||
ebuck6d *= sme;
|
||||
}
|
||||
} else forcebuck6d = 0.0;
|
||||
|
||||
// add forces of additional LJ interaction
|
||||
|
||||
fpair = forcebuck6d * r32inv;
|
||||
if (newton_pair || i1 < nlocal) {
|
||||
f[i1][0] += delx3*fpair;
|
||||
f[i1][1] += dely3*fpair;
|
||||
f[i1][2] += delz3*fpair;
|
||||
}
|
||||
if (newton_pair || i3 < nlocal) {
|
||||
f[i3][0] -= delx3*fpair;
|
||||
f[i3][1] -= dely3*fpair;
|
||||
f[i3][2] -= delz3*fpair;
|
||||
}
|
||||
|
||||
evdwl = 0.0;
|
||||
if (eflag) {
|
||||
if (rsq3 < cut_ljsq[itype][jtype]) {
|
||||
evdwl = ebuck6d - offset[itype][jtype];
|
||||
}
|
||||
}
|
||||
|
||||
//update pair energy and velocities
|
||||
|
||||
if (evflag) force->pair->ev_tally(i1,i3,nlocal,newton_pair,
|
||||
evdwl,0.0,fpair,delx3,dely3,delz3);
|
||||
|
||||
tk = multiplicity[type]*acos(c)-th0[type];
|
||||
|
||||
if (eflag) eangle = k[type]*(1.0+cos(tk));
|
||||
|
||||
a = k[type]*multiplicity[type]*sin(tk)*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;
|
||||
|
||||
// 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 AngleCosineBuck6d::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(multiplicity,n+1,"angle:multiplicity");
|
||||
memory->create(th0,n+1,"angle:th0");
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more types
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosineBuck6d::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 4) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
force->bounds(FLERR,arg[0],atom->nangletypes,ilo,ihi);
|
||||
|
||||
double c_one = force->numeric(FLERR,arg[1]);
|
||||
int n_one = force->inumeric(FLERR,arg[2]);
|
||||
int th0_one = force->numeric(FLERR,arg[3]);
|
||||
if (n_one <= 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
k[i] = c_one;
|
||||
multiplicity[i] = n_one;
|
||||
// transform offset angle to radians
|
||||
th0[i] = th0_one/180.0 * MY_PI;
|
||||
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
check special_bond settings are valid and initialize vdwl parameters
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosineBuck6d::init_style()
|
||||
{
|
||||
// set local ptrs to buck6d 13 arrays setup by Pair
|
||||
int itmp;
|
||||
if (force->pair == NULL)
|
||||
error->all(FLERR,"Angle cosine/buck6d is incompatible with Pair style");
|
||||
cut_ljsq = (double **) force->pair->extract("cut_ljsq",itmp);
|
||||
buck6d1 = (double **) force->pair->extract("buck6d1",itmp);
|
||||
buck6d2 = (double **) force->pair->extract("buck6d2",itmp);
|
||||
buck6d3 = (double **) force->pair->extract("buck6d3",itmp);
|
||||
buck6d4 = (double **) force->pair->extract("buck6d4",itmp);
|
||||
rsmooth_sq = (double **) force->pair->extract("rsmooth_sq",itmp);
|
||||
c0 = (double **) force->pair->extract("c0",itmp);
|
||||
c1 = (double **) force->pair->extract("c1",itmp);
|
||||
c2 = (double **) force->pair->extract("c2",itmp);
|
||||
c3 = (double **) force->pair->extract("c3",itmp);
|
||||
c4 = (double **) force->pair->extract("c4",itmp);
|
||||
c5 = (double **) force->pair->extract("c5",itmp);
|
||||
offset = (double **) force->pair->extract("offset",itmp);
|
||||
if (!buck6d1 || !buck6d2 || !buck6d3 || !buck6d4 || !c0 || !c1 || !c2)
|
||||
error->all(FLERR,"Angle cosine/buck6d is incompatible with Pair style");
|
||||
|
||||
// special bonds must be x 0 x or double counting
|
||||
if (force->special_lj[2] != 0.0)
|
||||
error->all(FLERR,"Angle style requires special_bonds lj = x,0,x;"
|
||||
" otherwise buck6d 1-3 interaction are counted twice");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleCosineBuck6d::equilibrium_angle(int i)
|
||||
{
|
||||
return MY_PI;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes out coeffs to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosineBuck6d::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&multiplicity[1],sizeof(int),atom->nangletypes,fp);
|
||||
fwrite(&th0[1],sizeof(double),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads coeffs from restart file, bcasts them
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosineBuck6d::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fread(&multiplicity[1],sizeof(int),atom->nangletypes,fp);
|
||||
fread(&th0[1],sizeof(double),atom->nangletypes,fp);
|
||||
}
|
||||
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&multiplicity[1],atom->nangletypes,MPI_INT,0,world);
|
||||
MPI_Bcast(&th0[1],atom->nangletypes,MPI_INT,0,world);
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleCosineBuck6d::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->nangletypes; i++) {
|
||||
fprintf(fp,"%d %g %d %d\n",i,k[i],multiplicity[i],th0[i]);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleCosineBuck6d::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 tk = multiplicity[type]*acos(c)-th0[type];
|
||||
|
||||
return k[type]*(1.0+cos(tk));
|
||||
}
|
||||
|
|
@ -0,0 +1,62 @@
|
|||
/* -*- 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(cosine/buck6d, AngleCosineBuck6d)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_ANGLE_COSINE_BUCK6D_H
|
||||
#define LMP_ANGLE_COSINE_BUCK6D_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include "angle.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class AngleCosineBuck6d : public Angle {
|
||||
public:
|
||||
AngleCosineBuck6d(class LAMMPS *);
|
||||
virtual ~AngleCosineBuck6d();
|
||||
virtual void compute(int, int);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double equilibrium_angle(int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_data(FILE *);
|
||||
double single(int, int, int, int);
|
||||
|
||||
protected:
|
||||
double *k,*th0;
|
||||
double *eps,*d0;
|
||||
double **buck6d1,**buck6d2,**buck6d3,**buck6d4,**cut_ljsq;
|
||||
double **c0,**c1,**c2,**c3,**c4,**c5,**rsmooth_sq,**offset;
|
||||
int *multiplicity;
|
||||
|
||||
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,337 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Hendrik Heenen (Technical University of Munich)
|
||||
and Rochus Schmid (Ruhr-Universitaet Bochum)
|
||||
[ based on improper_fourier.cpp Loukas D. Peristeras (Scienomics SARL) ]
|
||||
[ based on improper_umbrella.cpp Tod A Pascal (Caltech) ]
|
||||
[ abbreviated from and verified via DLPOLY2.0 ]
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <mpi.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include "string.h"
|
||||
#include "improper_inversion_harmonic.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define TOLERANCE 0.05
|
||||
#define SMALL 0.001
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ImproperInversionHarmonic::ImproperInversionHarmonic(LAMMPS *lmp) : Improper(lmp)
|
||||
{
|
||||
writedata = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ImproperInversionHarmonic::~ImproperInversionHarmonic()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(kw);
|
||||
memory->destroy(w0);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,i4,n,type;
|
||||
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z;
|
||||
double rrvb1,rrvb2,rrvb3,rr2vb1,rr2vb2,rr2vb3;
|
||||
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
int **improperlist = neighbor->improperlist;
|
||||
int nimproperlist = neighbor->nimproperlist;
|
||||
|
||||
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];
|
||||
|
||||
// 1st bond - IJ
|
||||
|
||||
vb1x = x[i2][0] - x[i1][0];
|
||||
vb1y = x[i2][1] - x[i1][1];
|
||||
vb1z = x[i2][2] - x[i1][2];
|
||||
rrvb1 = 1.0/sqrt(vb1x*vb1x+vb1y*vb1y+vb1z*vb1z);
|
||||
rr2vb1 = rrvb1*rrvb1;
|
||||
|
||||
// 2nd bond - IK
|
||||
|
||||
vb2x = x[i3][0] - x[i1][0];
|
||||
vb2y = x[i3][1] - x[i1][1];
|
||||
vb2z = x[i3][2] - x[i1][2];
|
||||
rrvb2 = 1.0/sqrt(vb2x*vb2x+vb2y*vb2y+vb2z*vb2z);
|
||||
rr2vb2 = rrvb2*rrvb2;
|
||||
|
||||
// 3rd bond - IL
|
||||
|
||||
vb3x = x[i4][0] - x[i1][0];
|
||||
vb3y = x[i4][1] - x[i1][1];
|
||||
vb3z = x[i4][2] - x[i1][2];
|
||||
rrvb3 = 1.0/sqrt(vb3x*vb3x+vb3y*vb3y+vb3z*vb3z);
|
||||
rr2vb3 = rrvb3*rrvb3;
|
||||
|
||||
// compute all three inversion angles
|
||||
invang(i1,i2,i3,i4, type,evflag,eflag,
|
||||
vb3x, vb3y, vb3z, rrvb3, rr2vb3,
|
||||
vb2x, vb2y, vb2z, rrvb2, rr2vb2,
|
||||
vb1x, vb1y, vb1z, rrvb1, rr2vb1);
|
||||
invang(i1,i3,i4,i2, type,evflag,eflag,
|
||||
vb1x, vb1y, vb1z, rrvb1, rr2vb1,
|
||||
vb3x, vb3y, vb3z, rrvb3, rr2vb3,
|
||||
vb2x, vb2y, vb2z, rrvb2, rr2vb2);
|
||||
invang(i1,i4,i2,i3, type,evflag,eflag,
|
||||
vb2x, vb2y, vb2z, rrvb2, rr2vb2,
|
||||
vb1x, vb1y, vb1z, rrvb1, rr2vb1,
|
||||
vb3x, vb3y, vb3z, rrvb3, rr2vb3);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute inversion angles + energy and forces
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::invang(const int &i1,const int &i2,
|
||||
const int &i3,const int &i4,
|
||||
const int &type,const int &evflag,const int &eflag,
|
||||
const double &vb1x, const double &vb1y, const double &vb1z,
|
||||
const double &rrvb1, const double &rr2vb1,
|
||||
const double &vb2x, const double &vb2y, const double &vb2z,
|
||||
const double &rrvb2, const double &rr2vb2,
|
||||
const double &vb3x, const double &vb3y, const double &vb3z,
|
||||
const double &rrvb3, const double &rr2vb3)
|
||||
{
|
||||
double eimproper,f1[3],f2[3],f3[3],f4[3];
|
||||
double omega,cosomega,domega,gomega,rjk,rjl;
|
||||
double upx,upy,upz,upn,rup,umx,umy,umz,umn,rum,wwr;
|
||||
double rucb,rudb,rvcb,rvdb,rupupn,rumumn;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
eimproper = 0.0;
|
||||
|
||||
// scalar products of IJ*IK and IJ*IL
|
||||
rjk=vb3x*vb2x+vb3y*vb2y+vb3z*vb2z;
|
||||
rjl=vb1x*vb3x+vb1y*vb3y+vb1z*vb3z;
|
||||
|
||||
// unit-vector: IK+IL
|
||||
upx=vb2x*rrvb2+vb1x*rrvb1;
|
||||
upy=vb2y*rrvb2+vb1y*rrvb1;
|
||||
upz=vb2z*rrvb2+vb1z*rrvb1;
|
||||
upn=1.0/sqrt(upx*upx+upy*upy+upz*upz);
|
||||
upx=upx*upn;
|
||||
upy=upy*upn;
|
||||
upz=upz*upn;
|
||||
rup=vb3x*upx+vb3y*upy+vb3z*upz;
|
||||
|
||||
// unit-vector: IK-IL
|
||||
umx=vb2x*rrvb2-vb1x*rrvb1;
|
||||
umy=vb2y*rrvb2-vb1y*rrvb1;
|
||||
umz=vb2z*rrvb2-vb1z*rrvb1;
|
||||
umn=1.0/sqrt(umx*umx+umy*umy+umz*umz);
|
||||
umx=umx*umn;
|
||||
umy=umy*umn;
|
||||
umz=umz*umn;
|
||||
rum=vb3x*umx+vb3y*umy+vb3z*umz;
|
||||
|
||||
// angle theta
|
||||
wwr=sqrt(rup*rup+rum*rum);
|
||||
|
||||
cosomega=wwr*rrvb3;
|
||||
|
||||
if (cosomega > 1.0) cosomega = 1.0;
|
||||
|
||||
omega=acos(cosomega);
|
||||
|
||||
domega = acos(cosomega)-w0[type];
|
||||
if (eflag) eimproper = kw[type]*(domega*domega);
|
||||
|
||||
// kw[type] is divided by 3 -> threefold contribution
|
||||
gomega=0.0;
|
||||
if (omega*omega > 1.0e-24) gomega=2.0*kw[type]*(domega)/(sin(omega));
|
||||
|
||||
// projection IK and IL on unit vectors and contribution on IK and IL
|
||||
rucb = rjk-rup*(vb2x*upx+vb2y*upy+vb2z*upz);
|
||||
rudb = rjl-rup*(vb1x*upx+vb1y*upy+vb1z*upz);
|
||||
rvcb = rjk-rum*(vb2x*umx+vb2y*umy+vb2z*umz);
|
||||
rvdb = rjl-rum*(vb1x*umx+vb1y*umy+vb1z*umz);
|
||||
|
||||
rupupn = rup*upn;
|
||||
rumumn = rum*umn;
|
||||
|
||||
// force contributions of angle
|
||||
f2[0]=gomega*(-cosomega*vb3x*rr2vb3+rrvb3*(rup*upx+rum*umx)/wwr);
|
||||
f2[1]=gomega*(-cosomega*vb3y*rr2vb3+rrvb3*(rup*upy+rum*umy)/wwr);
|
||||
f2[2]=gomega*(-cosomega*vb3z*rr2vb3+rrvb3*(rup*upz+rum*umz)/wwr);
|
||||
|
||||
f3[0]=gomega*rrvb3*(rupupn*rrvb2*(vb3x-rup*upx-rucb*vb2x*rr2vb2) +
|
||||
rumumn*rrvb2*(vb3x-rum*umx-rvcb*vb2x*rr2vb2))/wwr;
|
||||
f3[1]=gomega*rrvb3*(rupupn*rrvb2*(vb3y-rup*upy-rucb*vb2y*rr2vb2) +
|
||||
rumumn*rrvb2*(vb3y-rum*umy-rvcb*vb2y*rr2vb2))/wwr;
|
||||
f3[2]=gomega*rrvb3*(rupupn*rrvb2*(vb3z-rup*upz-rucb*vb2z*rr2vb2) +
|
||||
rumumn*rrvb2*(vb3z-rum*umz-rvcb*vb2z*rr2vb2))/wwr;
|
||||
|
||||
f4[0]=gomega*rrvb3*(rupupn*rrvb1*(vb3x-rup*upx-rudb*vb1x*rr2vb1) -
|
||||
rumumn*rrvb1*(vb3x-rum*umx-rvdb*vb1x*rr2vb1))/wwr;
|
||||
f4[1]=gomega*rrvb3*(rupupn*rrvb1*(vb3y-rup*upy-rudb*vb1y*rr2vb1) -
|
||||
rumumn*rrvb1*(vb3y-rum*umy-rvdb*vb1y*rr2vb1))/wwr;
|
||||
f4[2]=gomega*rrvb3*(rupupn*rrvb1*(vb3z-rup*upz-rudb*vb1z*rr2vb1) -
|
||||
rumumn*rrvb1*(vb3z-rum*umz-rvdb*vb1z*rr2vb1))/wwr;
|
||||
|
||||
f1[0] = -(f2[0] + f3[0] + f4[0]);
|
||||
f1[1] = -(f2[1] + f3[1] + f4[1]);
|
||||
f1[2] = -(f2[2] + f3[2] + f4[2]);
|
||||
|
||||
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];
|
||||
}
|
||||
|
||||
double rb1x, rb1y, rb1z, rb2x, rb2y, rb2z, rb3x, rb3y, rb3z;
|
||||
if (evflag)
|
||||
|
||||
rb3x = vb1x - vb2x;
|
||||
rb3y = vb1y - vb2y;
|
||||
rb3z = vb1z - vb2z;
|
||||
|
||||
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f2,f3,f4,
|
||||
vb3x,vb3y,vb3z,
|
||||
vb2x,vb2y,vb2z,
|
||||
rb3x,rb3y,rb3z);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nimpropertypes;
|
||||
|
||||
memory->create(kw,n+1,"improper:kw");
|
||||
memory->create(w0,n+1,"improper:w0");
|
||||
|
||||
memory->create(setflag,n+1,"improper:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one type
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::coeff(int narg, char **arg)
|
||||
{
|
||||
|
||||
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 w_one = force->numeric(FLERR,arg[2]);
|
||||
|
||||
// convert w0 from degrees to radians
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
kw[i] = k_one/3.0; // parameter division due to 3 vector averaging
|
||||
w0[i] = w_one/180.0 * MY_PI;
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for improper coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes out coeffs to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&kw[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
fwrite(&w0[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads coeffs from restart file, bcasts them
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
fread(&kw[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
fread(&w0[1],sizeof(double),atom->nimpropertypes,fp);
|
||||
}
|
||||
MPI_Bcast(&kw[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&w0[1],atom->nimpropertypes,MPI_DOUBLE,0,world);
|
||||
|
||||
for (int i = 1; i <= atom->nimpropertypes; i++) setflag[i] = 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ImproperInversionHarmonic::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->nimpropertypes; i++)
|
||||
fprintf(fp,"%d %g %g\n",i,kw[i],w0[i]/MY_PI*180.0);
|
||||
}
|
||||
|
|
@ -0,0 +1,54 @@
|
|||
/* -*- 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(inversion/harmonic,ImproperInversionHarmonic)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_IMPROPER_INVERSION_HARMONIC_H
|
||||
#define LMP_IMPROPER_INVERSION_HARMONIC_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include "improper.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ImproperInversionHarmonic : public Improper {
|
||||
public:
|
||||
ImproperInversionHarmonic(class LAMMPS *);
|
||||
virtual ~ImproperInversionHarmonic();
|
||||
virtual void compute(int, int);
|
||||
void coeff(int, char **);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_data(FILE *);
|
||||
|
||||
protected:
|
||||
double *kw, *w0;
|
||||
void invang(const int &i1,const int &i2,const int &i3,const int &i4,
|
||||
const int &type,const int &evflag,const int &eflag,
|
||||
const double &vb1x, const double &vb1y, const double &vb1z,
|
||||
const double &rrvb1, const double &rr2vb1,
|
||||
const double &vb2x, const double &vb2y, const double &vb2z,
|
||||
const double &rrvb2, const double &rr2vb2,
|
||||
const double &vb3x, const double &vb3y, const double &vb3z,
|
||||
const double &rrvb3, const double &rr2vb3);
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
|
@ -0,0 +1,590 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Hendrik Heenen (Technical University of Munich)
|
||||
and Rochus Schmid (Ruhr-Universitaet Bochum)
|
||||
References: Bureekaew and Schmid, Phys. Status Solidi B 250, 1128 (2013)
|
||||
Fennell and Gezelter, JCP 124, 234104 (2006)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "pair_buck6d_coul_gauss_dsf.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "math_const.h"
|
||||
#include "error.h"
|
||||
#include "math_special.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBuck6dCoulGaussDSF::PairBuck6dCoulGaussDSF(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
single_enable = 1;
|
||||
writedata = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBuck6dCoulGaussDSF::~PairBuck6dCoulGaussDSF()
|
||||
{
|
||||
if (!copymode) {
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(cut_lj);
|
||||
memory->destroy(cut_ljsq);
|
||||
memory->destroy(alpha_ij);
|
||||
memory->destroy(f_shift_ij);
|
||||
memory->destroy(e_shift_ij);
|
||||
memory->destroy(buck6d1);
|
||||
memory->destroy(buck6d2);
|
||||
memory->destroy(buck6d3);
|
||||
memory->destroy(buck6d4);
|
||||
memory->destroy(c0);
|
||||
memory->destroy(c1);
|
||||
memory->destroy(c2);
|
||||
memory->destroy(c3);
|
||||
memory->destroy(c4);
|
||||
memory->destroy(c5);
|
||||
memory->destroy(rsmooth_sq);
|
||||
memory->destroy(offset);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,r2inv,r6inv,r14inv,rexp,forcecoul,forcebuck6d,factor_coul,factor_lj;
|
||||
double term1,term2,term3,term4,term5;
|
||||
double rcu,rqu,sme,smf,ebuck6d;
|
||||
double prefactor,erfcc,erfcd,t,arg;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
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;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
r14inv = r6inv*r6inv*r2inv;
|
||||
rexp = exp(-r*buck6d2[itype][jtype]);
|
||||
term1 = buck6d3[itype][jtype]*r6inv;
|
||||
term2 = buck6d4[itype][jtype]*r14inv;
|
||||
term3 = term2*term2;
|
||||
term4 = 1.0/(1.0 + term2);
|
||||
term5 = 1.0/(1.0 + 2.0*term2 + term3);
|
||||
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
|
||||
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
|
||||
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
|
||||
|
||||
// smoothing term
|
||||
if (rsq > rsmooth_sq[itype][jtype]) {
|
||||
rcu = r*rsq;
|
||||
rqu = rsq*rsq;
|
||||
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
|
||||
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
|
||||
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
|
||||
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
|
||||
// forcebuck6d is -dE/dr*r
|
||||
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r; //RS was here: changed this from +E*smf to -E*smf*r
|
||||
ebuck6d *= sme;
|
||||
}
|
||||
} else forcebuck6d = 0.0;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = qqrd2e*qtmp*q[j]/r;
|
||||
|
||||
//below: parametrization for erfcc = erf(alpha_ij[itype][jtype]*r);
|
||||
arg = alpha_ij[itype][jtype]*r;
|
||||
erfcd = MathSpecial::expmsq(arg);
|
||||
erfcc = 1 - (MathSpecial::my_erfcx(arg) * erfcd);
|
||||
|
||||
forcecoul = prefactor * ((erfcc/r) - (2.0/MY_PIS*alpha_ij[itype][jtype]*erfcd) +
|
||||
r*f_shift_ij[itype][jtype]) * r;
|
||||
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
fpair = (forcecoul + factor_lj*forcebuck6d) * r2inv;
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = ebuck6d - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
ecoul = prefactor * (erfcc - r*e_shift_ij[itype][jtype] -
|
||||
rsq*f_shift_ij[itype][jtype]);
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
}
|
||||
|
||||
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 PairBuck6dCoulGaussDSF::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(alpha_ij,n+1,n+1,"pair:alpha_ij");
|
||||
memory->create(f_shift_ij,n+1,n+1,"pair:f_shift_ij");
|
||||
memory->create(e_shift_ij,n+1,n+1,"pair:e_shift_ij");
|
||||
memory->create(buck6d1,n+1,n+1,"pair:buck6d1");
|
||||
memory->create(buck6d2,n+1,n+1,"pair:buck6d2");
|
||||
memory->create(buck6d3,n+1,n+1,"pair:buck6d3");
|
||||
memory->create(buck6d4,n+1,n+1,"pair:buck6d4");
|
||||
memory->create(c0,n+1,n+1,"pair:c0");
|
||||
memory->create(c1,n+1,n+1,"pair:c1");
|
||||
memory->create(c2,n+1,n+1,"pair:c2");
|
||||
memory->create(c3,n+1,n+1,"pair:c3");
|
||||
memory->create(c4,n+1,n+1,"pair:c4");
|
||||
memory->create(c5,n+1,n+1,"pair:c5");
|
||||
memory->create(rsmooth_sq,n+1,n+1,"pair:rsmooth_sq");
|
||||
memory->create(offset,n+1,n+1,"pair:offset");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 2 || narg > 3) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
vdwl_smooth = force->numeric(FLERR,arg[0]);
|
||||
|
||||
cut_lj_global = force->numeric(FLERR,arg[1]);
|
||||
if (narg == 2) cut_coul = cut_lj_global;
|
||||
else cut_coul = force->numeric(FLERR,arg[2]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i+1; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j])
|
||||
cut_lj[i][j] = cut_lj_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 7 || narg > 8)
|
||||
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 buck6d1_one = force->numeric(FLERR,arg[2]);
|
||||
double buck6d2_one = force->numeric(FLERR,arg[3]);
|
||||
double buck6d3_one = force->numeric(FLERR,arg[4]);
|
||||
double buck6d4_one = force->numeric(FLERR,arg[5]);
|
||||
double alpha_one = force->numeric(FLERR,arg[6]);
|
||||
|
||||
double cut_lj_one = cut_lj_global;
|
||||
if (narg == 8) cut_lj_one = force->numeric(FLERR,arg[7]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
buck6d1[i][j] = buck6d1_one;
|
||||
buck6d2[i][j] = buck6d2_one;
|
||||
buck6d3[i][j] = buck6d3_one;
|
||||
buck6d4[i][j] = buck6d4_one;
|
||||
alpha_ij[i][j] = alpha_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 PairBuck6dCoulGaussDSF::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style buck6d/coul/dsf requires atom attribute q");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
cut_coulsq = cut_coul * cut_coul;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairBuck6dCoulGaussDSF::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
double cut = MAX(cut_lj[i][j],cut_coul);
|
||||
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
|
||||
|
||||
//calculation of smoothing coefficients c0-c5
|
||||
c0[i][j] = c1[i][j] = c2[i][j] = c3[i][j] = c4[i][j] = c5[i][j] = 0.0;
|
||||
rsmooth_sq[i][j] = cut_ljsq[i][j];
|
||||
if (vdwl_smooth < 1.0) {
|
||||
double rsm = vdwl_smooth * cut_lj[i][j];
|
||||
double rsm_sq = rsm * rsm;
|
||||
double denom = pow((cut_lj[i][j]-rsm),5.0);
|
||||
c0[i][j] = cut_lj[i][j]*cut_ljsq[i][j]*(cut_ljsq[i][j]-
|
||||
5.0*cut_lj[i][j]*rsm+10.0*rsm_sq)/denom;
|
||||
c1[i][j] = -30.0*(cut_ljsq[i][j]*rsm_sq)/denom;
|
||||
c2[i][j] = 30.0*(cut_ljsq[i][j]*rsm + cut_lj[i][j]*rsm_sq)/denom;
|
||||
c3[i][j] = -10.0*(cut_ljsq[i][j] + 4.0*cut_lj[i][j]*rsm + rsm_sq)/denom;
|
||||
c4[i][j] = 15.0*(cut_lj[i][j]+rsm)/denom;
|
||||
c5[i][j] = -6.0/denom;
|
||||
rsmooth_sq[i][j] = rsm_sq;
|
||||
}
|
||||
|
||||
// if offset_flag, shift is only invoked if there is not already smoothing
|
||||
if (offset_flag && vdwl_smooth >= 1.0) {
|
||||
double term1 = buck6d3[i][j]/pow(cut_lj[i][j],6.0);
|
||||
double term4 = 1.0/(1.0 + (buck6d4[i][j]/pow(cut_lj[i][j],14.0)));
|
||||
double rexp = exp(-cut_lj[i][j]*buck6d2[i][j]);
|
||||
offset[i][j] = buck6d1[i][j]*rexp - term1*term4;
|
||||
} else offset[i][j] = 0.0;
|
||||
|
||||
double erfcd_c = exp(-alpha_ij[i][j]*alpha_ij[i][j]*cut_coul*cut_coul);
|
||||
double erfcc_c = erf(alpha_ij[i][j]*cut_coul);
|
||||
f_shift_ij[i][j] = -erfcc_c/cut_coulsq + 2.0/MY_PIS*alpha_ij[i][j]*erfcd_c/cut_coul;
|
||||
e_shift_ij[i][j] = erfcc_c/cut_coul - f_shift_ij[i][j]*cut_coul;
|
||||
|
||||
cut_ljsq[j][i] = cut_ljsq[i][j];
|
||||
alpha_ij[j][i] = alpha_ij[i][j];
|
||||
f_shift_ij[j][i] = f_shift_ij[i][j];
|
||||
e_shift_ij[j][i] = e_shift_ij[i][j];
|
||||
buck6d1[j][i] = buck6d1[i][j];
|
||||
buck6d2[j][i] = buck6d2[i][j];
|
||||
buck6d3[j][i] = buck6d3[i][j];
|
||||
buck6d4[j][i] = buck6d4[i][j];
|
||||
c0[j][i] = c0[i][j];
|
||||
c1[j][i] = c1[i][j];
|
||||
c2[j][i] = c2[i][j];
|
||||
c3[j][i] = c3[i][j];
|
||||
c4[j][i] = c4[i][j];
|
||||
c5[j][i] = c5[i][j];
|
||||
rsmooth_sq[j][i] = rsmooth_sq[i][j];
|
||||
offset[j][i] = offset[i][j];
|
||||
|
||||
return cut;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::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(&buck6d1[i][j],sizeof(double),1,fp);
|
||||
fwrite(&buck6d2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&buck6d3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&buck6d4[i][j],sizeof(double),1,fp);
|
||||
fwrite(&alpha_ij[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::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(&buck6d1[i][j],sizeof(double),1,fp);
|
||||
fread(&buck6d2[i][j],sizeof(double),1,fp);
|
||||
fread(&buck6d3[i][j],sizeof(double),1,fp);
|
||||
fread(&buck6d4[i][j],sizeof(double),1,fp);
|
||||
fread(&alpha_ij[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&buck6d1[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&buck6d2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&buck6d3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&buck6d4[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alpha_ij[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 PairBuck6dCoulGaussDSF::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&vdwl_smooth,sizeof(double),1,fp);
|
||||
fwrite(&cut_lj_global,sizeof(double),1,fp);
|
||||
fwrite(&cut_coul,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
fwrite(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&vdwl_smooth,sizeof(double),1,fp);
|
||||
fread(&cut_lj_global,sizeof(double),1,fp);
|
||||
fread(&cut_coul,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
fread(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&vdwl_smooth,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g %g %g %g %g\n",i,
|
||||
buck6d1[i][i],buck6d2[i][i],buck6d3[i][i],
|
||||
buck6d4[i][i],alpha_ij[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussDSF::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 %g %g\n",i,j,
|
||||
buck6d1[i][j],buck6d2[i][j],buck6d3[i][j],
|
||||
buck6d4[i][j],alpha_ij[i][j],cut_lj[i][j]);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairBuck6dCoulGaussDSF::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r,r2inv,r6inv,r14inv,rexp,term1,term2,term3,term4,term5;
|
||||
double rcu,rqu,sme,smf;
|
||||
double erfcc,erfcd,prefactor,arg;
|
||||
double forcecoul,forcebuck6d,ebuck6d,phicoul,phibuck6d;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
r14inv = r6inv*r6inv*r2inv;
|
||||
rexp = exp(-r*buck6d2[itype][jtype]);
|
||||
term1 = buck6d3[itype][jtype]*r6inv;
|
||||
term2 = buck6d4[itype][jtype]*r14inv;
|
||||
term3 = term2*term2;
|
||||
term4 = 1.0/(1.0 + term2);
|
||||
term5 = 1.0/(1.0 + 2.0*term2 + term3);
|
||||
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
|
||||
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
|
||||
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
|
||||
// smoothing term
|
||||
if (rsq > rsmooth_sq[itype][jtype]) {
|
||||
rcu = r*rsq;
|
||||
rqu = rsq*rsq;
|
||||
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
|
||||
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
|
||||
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
|
||||
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
|
||||
// forcebuck6d is -dE/dr*r
|
||||
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r; //RS was here: changed this from +E*smf to -E*smf*r
|
||||
ebuck6d *= sme;
|
||||
}
|
||||
} else forcebuck6d = 0.0;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
prefactor = factor_coul * force->qqrd2e * atom->q[i] * atom->q[j] / r;
|
||||
arg = alpha_ij[itype][jtype]*r;
|
||||
erfcd = MathSpecial::expmsq(arg);
|
||||
erfcc = 1 - (MathSpecial::my_erfcx(arg) * erfcd);
|
||||
forcecoul = prefactor * ((erfcc/r) - (2.0/MY_PIS*alpha_ij[itype][jtype]*erfcd) +
|
||||
r*f_shift_ij[itype][jtype]) * r;
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
fforce = (forcecoul + factor_lj*forcebuck6d) * r2inv;
|
||||
|
||||
double eng = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
phibuck6d = ebuck6d - offset[itype][jtype];
|
||||
eng += factor_lj*phibuck6d;
|
||||
}
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
phicoul = prefactor * (erfcc - r*e_shift_ij[itype][jtype] -
|
||||
rsq*f_shift_ij[itype][jtype]);
|
||||
eng += phicoul;
|
||||
}
|
||||
|
||||
return eng;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairBuck6dCoulGaussDSF::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str,"cut_ljsq") == 0) return (void *) cut_ljsq;
|
||||
if (strcmp(str,"buck6d1") == 0) return (void *) buck6d1;
|
||||
if (strcmp(str,"buck6d2") == 0) return (void *) buck6d2;
|
||||
if (strcmp(str,"buck6d3") == 0) return (void *) buck6d3;
|
||||
if (strcmp(str,"buck6d4") == 0) return (void *) buck6d4;
|
||||
if (strcmp(str,"rsmooth_sq") == 0) return (void *) rsmooth_sq;
|
||||
if (strcmp(str,"c0") == 0) return (void *) c0;
|
||||
if (strcmp(str,"c1") == 0) return (void *) c1;
|
||||
if (strcmp(str,"c2") == 0) return (void *) c2;
|
||||
if (strcmp(str,"c3") == 0) return (void *) c3;
|
||||
if (strcmp(str,"c4") == 0) return (void *) c4;
|
||||
if (strcmp(str,"c5") == 0) return (void *) c5;
|
||||
if (strcmp(str,"offset") == 0) return (void *) offset;
|
||||
|
||||
if (strcmp(str,"cut_coul") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_coul;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,79 @@
|
|||
/* -*- 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(buck6d/coul/gauss/dsf,PairBuck6dCoulGaussDSF)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_BUCK6D_COUL_GAUSS_DSF_H
|
||||
#define LMP_PAIR_BUCK6D_COUL_GAUSS_DSF_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairBuck6dCoulGaussDSF : public Pair {
|
||||
public:
|
||||
PairBuck6dCoulGaussDSF(class LAMMPS *);
|
||||
virtual ~PairBuck6dCoulGaussDSF();
|
||||
virtual void compute(int, int);
|
||||
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 *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
double single(int, int, int, int, double, double, double, double &);
|
||||
void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
double cut_lj_global;
|
||||
double **cut_lj,**cut_ljsq;
|
||||
double **alpha_ij;
|
||||
double **buck6d1,**buck6d2,**buck6d3,**buck6d4,**offset;
|
||||
double **f_shift_ij,**e_shift_ij;
|
||||
double cut_coul,cut_coulsq;
|
||||
double vdwl_smooth;
|
||||
double **c0,**c1,**c2,**c3,**c4,**c5,**rsmooth_sq;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style lj/cut/coul/dsf requires atom attribute q
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,646 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
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: Hendrik Heenen (Technical University of Munich)
|
||||
and Rochus Schmid (Ruhr-Universitaet Bochum)
|
||||
References: Bureekaew and Schmid, Phys. Status Solidi B 250, 1128 (2013)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "pair_buck6d_coul_gauss_long.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "math_const.h"
|
||||
#include "error.h"
|
||||
#include "math_special.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBuck6dCoulGaussLong::PairBuck6dCoulGaussLong(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
ewaldflag = pppmflag = 1;
|
||||
single_enable = 1;
|
||||
writedata = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairBuck6dCoulGaussLong::~PairBuck6dCoulGaussLong()
|
||||
{
|
||||
if (!copymode) {
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(cut_lj);
|
||||
memory->destroy(cut_ljsq);
|
||||
memory->destroy(alpha_ij);
|
||||
memory->destroy(buck6d1);
|
||||
memory->destroy(buck6d2);
|
||||
memory->destroy(buck6d3);
|
||||
memory->destroy(buck6d4);
|
||||
memory->destroy(c0);
|
||||
memory->destroy(c1);
|
||||
memory->destroy(c2);
|
||||
memory->destroy(c3);
|
||||
memory->destroy(c4);
|
||||
memory->destroy(c5);
|
||||
memory->destroy(rsmooth_sq);
|
||||
memory->destroy(offset);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double r,rsq,r2inv,r6inv,r14inv,rexp,forcecoul,forcebuck6d,factor_coul,factor_lj;
|
||||
double grij,expm2,erf;
|
||||
double term1,term2,term3,term4,term5;
|
||||
double rcu,rqu,sme,smf,ebuck6d,ealpha;
|
||||
double prefactor,erfa,expa,t,arg,falpha;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
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;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
r14inv = r6inv*r6inv*r2inv;
|
||||
rexp = exp(-r*buck6d2[itype][jtype]);
|
||||
term1 = buck6d3[itype][jtype]*r6inv;
|
||||
term2 = buck6d4[itype][jtype]*r14inv;
|
||||
term3 = term2*term2;
|
||||
term4 = 1.0/(1.0 + term2);
|
||||
term5 = 1.0/(1.0 + 2.0*term2 + term3);
|
||||
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
|
||||
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
|
||||
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
|
||||
|
||||
// smoothing term
|
||||
if (rsq > rsmooth_sq[itype][jtype]) {
|
||||
rcu = r*rsq;
|
||||
rqu = rsq*rsq;
|
||||
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
|
||||
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
|
||||
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
|
||||
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
|
||||
// forcebuck6d is -dE/dr*r
|
||||
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r;
|
||||
ebuck6d *= sme;
|
||||
}
|
||||
} else forcebuck6d = 0.0;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
// long range - real space
|
||||
grij = g_ewald * r;
|
||||
expm2 = MathSpecial::expmsq(grij);
|
||||
erf = 1 - (MathSpecial::my_erfcx(grij) * expm2);
|
||||
|
||||
// gaussian for 1/r alpha_ij contribution
|
||||
arg = alpha_ij[itype][jtype]*r;
|
||||
expa = MathSpecial::expmsq(arg);
|
||||
erfa = 1 - (MathSpecial::my_erfcx(arg) * expa);
|
||||
|
||||
prefactor = qqrd2e*qtmp*q[j]/r;
|
||||
falpha = erfa - EWALD_F*arg*expa;
|
||||
forcecoul = prefactor * (falpha - erf + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*falpha;
|
||||
|
||||
// (q*q/r) * (gauss(alpha_ij) - gauss(alpha_long)
|
||||
ealpha = prefactor * (erfa-erf);
|
||||
// smoothing term - NOTE: ingnored in special_bonds correction
|
||||
// since likely rsmooth_sq_c >> d(special)
|
||||
if (rsq > rsmooth_sq_c) {
|
||||
rcu = r*rsq;
|
||||
rqu = rsq*rsq;
|
||||
sme = c5_c*rqu*r + c4_c*rqu + c3_c*rcu + c2_c*rsq + c1_c*r + c0_c;
|
||||
smf = 5.0*c5_c*rqu + 4.0*c4_c*rcu + 3.0*c3_c*rsq + 2.0*c2_c*r + c1_c;
|
||||
forcecoul = forcecoul*sme - ealpha*smf*r;
|
||||
ealpha *= sme;
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
|
||||
fpair = (forcecoul + factor_lj*forcebuck6d) * r2inv;
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = ebuck6d - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
ecoul = ealpha;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor*erfa;
|
||||
} else ecoul = 0.0;
|
||||
}
|
||||
|
||||
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 PairBuck6dCoulGaussLong::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(alpha_ij,n+1,n+1,"pair:alpha_ij");
|
||||
memory->create(buck6d1,n+1,n+1,"pair:buck6d1");
|
||||
memory->create(buck6d2,n+1,n+1,"pair:buck6d2");
|
||||
memory->create(buck6d3,n+1,n+1,"pair:buck6d3");
|
||||
memory->create(buck6d4,n+1,n+1,"pair:buck6d4");
|
||||
memory->create(c0,n+1,n+1,"pair:c0");
|
||||
memory->create(c1,n+1,n+1,"pair:c1");
|
||||
memory->create(c2,n+1,n+1,"pair:c2");
|
||||
memory->create(c3,n+1,n+1,"pair:c3");
|
||||
memory->create(c4,n+1,n+1,"pair:c4");
|
||||
memory->create(c5,n+1,n+1,"pair:c5");
|
||||
memory->create(rsmooth_sq,n+1,n+1,"pair:rsmooth_sq");
|
||||
memory->create(offset,n+1,n+1,"pair:offset");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 3 || narg > 4) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
vdwl_smooth = force->numeric(FLERR,arg[0]);
|
||||
coul_smooth = force->numeric(FLERR,arg[1]);
|
||||
|
||||
cut_lj_global = force->numeric(FLERR,arg[2]);
|
||||
if (narg == 3) cut_coul = cut_lj_global;
|
||||
else cut_coul = force->numeric(FLERR,arg[3]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i+1; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j])
|
||||
cut_lj[i][j] = cut_lj_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 7 || narg > 8)
|
||||
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 buck6d1_one = force->numeric(FLERR,arg[2]);
|
||||
double buck6d2_one = force->numeric(FLERR,arg[3]);
|
||||
double buck6d3_one = force->numeric(FLERR,arg[4]);
|
||||
double buck6d4_one = force->numeric(FLERR,arg[5]);
|
||||
double alpha_one = force->numeric(FLERR,arg[6]);
|
||||
|
||||
double cut_lj_one = cut_lj_global;
|
||||
if (narg == 8) cut_lj_one = force->numeric(FLERR,arg[7]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
buck6d1[i][j] = buck6d1_one;
|
||||
buck6d2[i][j] = buck6d2_one;
|
||||
buck6d3[i][j] = buck6d3_one;
|
||||
buck6d4[i][j] = buck6d4_one;
|
||||
alpha_ij[i][j] = alpha_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 PairBuck6dCoulGaussLong::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style buck6d/coul/dsf requires atom attribute q");
|
||||
|
||||
// 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;
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
cut_coulsq = cut_coul * cut_coul;
|
||||
|
||||
//calculation of smoothing coefficients c0_c-c5_c for coulomb smoothing
|
||||
c0_c = c1_c = c2_c = c3_c = c4_c = c5_c = 0.0;
|
||||
rsmooth_sq_c = cut_coulsq;
|
||||
if (coul_smooth < 1.0) {
|
||||
double rsm = coul_smooth * cut_coul;
|
||||
double rsm_sq = rsm * rsm;
|
||||
double denom = pow((cut_coul-rsm),5.0);
|
||||
c0_c = cut_coul*cut_coulsq*(cut_coulsq-
|
||||
5.0*cut_coul*rsm+10.0*rsm_sq)/denom;
|
||||
c1_c = -30.0*(cut_coulsq*rsm_sq)/denom;
|
||||
c2_c = 30.0*(cut_coulsq*rsm + cut_coul*rsm_sq)/denom;
|
||||
c3_c = -10.0*(cut_coulsq + 4.0*cut_coul*rsm + rsm_sq)/denom;
|
||||
c4_c = 15.0*(cut_coul+rsm)/denom;
|
||||
c5_c = -6.0/denom;
|
||||
rsmooth_sq_c = rsm_sq;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairBuck6dCoulGaussLong::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
double cut = MAX(cut_lj[i][j],cut_coul);
|
||||
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
|
||||
|
||||
//calculation of smoothing coefficients c0-c5
|
||||
c0[i][j] = c1[i][j] = c2[i][j] = c3[i][j] = c4[i][j] = c5[i][j] = 0.0;
|
||||
rsmooth_sq[i][j] = cut_ljsq[i][j];
|
||||
if (vdwl_smooth < 1.0) {
|
||||
double rsm = vdwl_smooth * cut_lj[i][j];
|
||||
double rsm_sq = rsm * rsm;
|
||||
double denom = pow((cut_lj[i][j]-rsm),5.0);
|
||||
c0[i][j] = cut_lj[i][j]*cut_ljsq[i][j]*(cut_ljsq[i][j]-
|
||||
5.0*cut_lj[i][j]*rsm+10.0*rsm_sq)/denom;
|
||||
c1[i][j] = -30.0*(cut_ljsq[i][j]*rsm_sq)/denom;
|
||||
c2[i][j] = 30.0*(cut_ljsq[i][j]*rsm + cut_lj[i][j]*rsm_sq)/denom;
|
||||
c3[i][j] = -10.0*(cut_ljsq[i][j] + 4.0*cut_lj[i][j]*rsm + rsm_sq)/denom;
|
||||
c4[i][j] = 15.0*(cut_lj[i][j]+rsm)/denom;
|
||||
c5[i][j] = -6.0/denom;
|
||||
rsmooth_sq[i][j] = rsm_sq;
|
||||
}
|
||||
|
||||
// if offset_flag, shift is only invoked if there is not already smoothing
|
||||
if (offset_flag && vdwl_smooth >= 1.0) {
|
||||
double term1 = buck6d3[i][j]/pow(cut_lj[i][j],6.0);
|
||||
double term4 = 1.0/(1.0 + (buck6d4[i][j]/pow(cut_lj[i][j],14.0)));
|
||||
double rexp = exp(-cut_lj[i][j]*buck6d2[i][j]);
|
||||
offset[i][j] = buck6d1[i][j]*rexp - term1*term4;
|
||||
} else offset[i][j] = 0.0;
|
||||
|
||||
cut_ljsq[j][i] = cut_ljsq[i][j];
|
||||
alpha_ij[j][i] = alpha_ij[i][j];
|
||||
buck6d1[j][i] = buck6d1[i][j];
|
||||
buck6d2[j][i] = buck6d2[i][j];
|
||||
buck6d3[j][i] = buck6d3[i][j];
|
||||
buck6d4[j][i] = buck6d4[i][j];
|
||||
c0[j][i] = c0[i][j];
|
||||
c1[j][i] = c1[i][j];
|
||||
c2[j][i] = c2[i][j];
|
||||
c3[j][i] = c3[i][j];
|
||||
c4[j][i] = c4[i][j];
|
||||
c5[j][i] = c5[i][j];
|
||||
rsmooth_sq[j][i] = rsmooth_sq[i][j];
|
||||
offset[j][i] = offset[i][j];
|
||||
|
||||
return cut;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::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(&buck6d1[i][j],sizeof(double),1,fp);
|
||||
fwrite(&buck6d2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&buck6d3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&buck6d4[i][j],sizeof(double),1,fp);
|
||||
fwrite(&alpha_ij[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::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(&buck6d1[i][j],sizeof(double),1,fp);
|
||||
fread(&buck6d2[i][j],sizeof(double),1,fp);
|
||||
fread(&buck6d3[i][j],sizeof(double),1,fp);
|
||||
fread(&buck6d4[i][j],sizeof(double),1,fp);
|
||||
fread(&alpha_ij[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_lj[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&buck6d1[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&buck6d2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&buck6d3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&buck6d4[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&alpha_ij[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 PairBuck6dCoulGaussLong::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&vdwl_smooth,sizeof(double),1,fp);
|
||||
fwrite(&coul_smooth,sizeof(double),1,fp);
|
||||
fwrite(&cut_lj_global,sizeof(double),1,fp);
|
||||
fwrite(&cut_coul,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
fwrite(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&vdwl_smooth,sizeof(double),1,fp);
|
||||
fread(&coul_smooth,sizeof(double),1,fp);
|
||||
fread(&cut_lj_global,sizeof(double),1,fp);
|
||||
fread(&cut_coul,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
fread(&tail_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&vdwl_smooth,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&coul_smooth,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp,"%d %g %g %g %g %g\n",i,
|
||||
buck6d1[i][i],buck6d2[i][i],buck6d3[i][i],
|
||||
buck6d4[i][i],alpha_ij[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairBuck6dCoulGaussLong::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 %g %g\n",i,j,
|
||||
buck6d1[i][j],buck6d2[i][j],buck6d3[i][j],
|
||||
buck6d4[i][j],alpha_ij[i][j],cut_lj[i][j]);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairBuck6dCoulGaussLong::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
{
|
||||
double r,r2inv,r6inv,r14inv,rexp,term1,term2,term3,term4,term5;
|
||||
double rcu,rqu,sme,smf;
|
||||
double erfa,expa,prefactor,arg,falpha,ealpha;
|
||||
double grij,expm2,erf;
|
||||
double forcecoul,forcebuck6d,ebuck6d,phicoul,phibuck6d;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
r14inv = r6inv*r6inv*r2inv;
|
||||
rexp = exp(-r*buck6d2[itype][jtype]);
|
||||
term1 = buck6d3[itype][jtype]*r6inv;
|
||||
term2 = buck6d4[itype][jtype]*r14inv;
|
||||
term3 = term2*term2;
|
||||
term4 = 1.0/(1.0 + term2);
|
||||
term5 = 1.0/(1.0 + 2.0*term2 + term3);
|
||||
forcebuck6d = buck6d1[itype][jtype]*buck6d2[itype][jtype]*r*rexp;
|
||||
forcebuck6d -= term1*(6.0*term4 - term5*14.0*term2);
|
||||
ebuck6d = buck6d1[itype][jtype]*rexp - term1*term4;
|
||||
// smoothing term
|
||||
if (rsq > rsmooth_sq[itype][jtype]) {
|
||||
rcu = r*rsq;
|
||||
rqu = rsq*rsq;
|
||||
sme = c5[itype][jtype]*rqu*r + c4[itype][jtype]*rqu + c3[itype][jtype]*rcu +
|
||||
c2[itype][jtype]*rsq + c1[itype][jtype]*r + c0[itype][jtype];
|
||||
smf = 5.0*c5[itype][jtype]*rqu + 4.0*c4[itype][jtype]*rcu +
|
||||
3.0*c3[itype][jtype]*rsq + 2.0*c2[itype][jtype]*r + c1[itype][jtype];
|
||||
// forcebuck6d is -dE/dr*r
|
||||
forcebuck6d = forcebuck6d*sme - ebuck6d*smf*r; //RS was here: changed this from +E*smf to -E*smf*r
|
||||
ebuck6d *= sme;
|
||||
}
|
||||
} else forcebuck6d = 0.0;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
// long range - real space
|
||||
grij = g_ewald * r;
|
||||
expm2 = MathSpecial::expmsq(grij);
|
||||
erf = 1 - (MathSpecial::my_erfcx(grij) * expm2);
|
||||
|
||||
// gaussian for 1/r alpha_ij contribution
|
||||
arg = alpha_ij[itype][jtype]*r;
|
||||
expa = MathSpecial::expmsq(arg);
|
||||
erfa = 1 - (MathSpecial::my_erfcx(arg) * expa);
|
||||
|
||||
prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
|
||||
falpha = erfa - EWALD_F*arg*expa;
|
||||
forcecoul = prefactor * (falpha - erf + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*falpha;
|
||||
|
||||
ealpha = prefactor * (erfa-erf);
|
||||
// smoothing term
|
||||
if (rsq > rsmooth_sq_c) {
|
||||
rcu = r*rsq;
|
||||
rqu = rsq*rsq;
|
||||
sme = c5_c*rqu*r + c4_c*rqu + c3_c*rcu + c2_c*rsq + c1_c*r + c0_c;
|
||||
smf = 5.0*c5_c*rqu + 4.0*c4_c*rcu + 3.0*c3_c*rsq + 2.0*c2_c*r + c1_c;
|
||||
forcecoul = forcecoul*sme - ealpha*smf*r;
|
||||
ealpha *= sme;
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
|
||||
fforce = (forcecoul + factor_lj*forcebuck6d) * r2inv;
|
||||
|
||||
double eng = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
phibuck6d = ebuck6d - offset[itype][jtype];
|
||||
eng += factor_lj*phibuck6d;
|
||||
}
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
phicoul = ealpha;
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor*erfa;
|
||||
eng += phicoul;
|
||||
}
|
||||
|
||||
return eng;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairBuck6dCoulGaussLong::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str,"cut_ljsq") == 0) return (void *) cut_ljsq;
|
||||
if (strcmp(str,"buck6d1") == 0) return (void *) buck6d1;
|
||||
if (strcmp(str,"buck6d2") == 0) return (void *) buck6d2;
|
||||
if (strcmp(str,"buck6d3") == 0) return (void *) buck6d3;
|
||||
if (strcmp(str,"buck6d4") == 0) return (void *) buck6d4;
|
||||
if (strcmp(str,"rsmooth_sq") == 0) return (void *) rsmooth_sq;
|
||||
if (strcmp(str,"c0") == 0) return (void *) c0;
|
||||
if (strcmp(str,"c1") == 0) return (void *) c1;
|
||||
if (strcmp(str,"c2") == 0) return (void *) c2;
|
||||
if (strcmp(str,"c3") == 0) return (void *) c3;
|
||||
if (strcmp(str,"c4") == 0) return (void *) c4;
|
||||
if (strcmp(str,"c5") == 0) return (void *) c5;
|
||||
if (strcmp(str,"offset") == 0) return (void *) offset;
|
||||
|
||||
if (strcmp(str,"cut_coul") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_coul;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,80 @@
|
|||
/* -*- 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(buck6d/coul/gauss/long,PairBuck6dCoulGaussLong)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_BUCK6D_COUL_GAUSS_LONG_H
|
||||
#define LMP_PAIR_BUCK6D_COUL_GAUSS_LONG_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairBuck6dCoulGaussLong : public Pair {
|
||||
public:
|
||||
PairBuck6dCoulGaussLong(class LAMMPS *);
|
||||
virtual ~PairBuck6dCoulGaussLong();
|
||||
virtual void compute(int, int);
|
||||
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 *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
double single(int, int, int, int, double, double, double, double &);
|
||||
void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
double cut_lj_global;
|
||||
double **cut_lj,**cut_ljsq;
|
||||
double **alpha_ij;
|
||||
double **buck6d1,**buck6d2,**buck6d3,**buck6d4,**offset;
|
||||
double cut_coul, cut_coulsq;
|
||||
double vdwl_smooth, coul_smooth;
|
||||
double c0_c,c1_c,c2_c,c3_c,c4_c,c5_c,rsmooth_sq_c;
|
||||
double **c0,**c1,**c2,**c3,**c4,**c5,**rsmooth_sq;
|
||||
double g_ewald;
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style lj/cut/coul/dsf requires atom attribute q
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue