Merge branch 'develop' into collected-small-changes

This commit is contained in:
Axel Kohlmeyer
2022-04-29 05:51:05 -04:00
26 changed files with 3194 additions and 44 deletions

View File

@ -2,7 +2,7 @@
DOXYFILE_ENCODING = UTF-8
PROJECT_NAME = "LAMMPS Programmer's Guide"
PROJECT_NUMBER = "24 August 2020"
PROJECT_NUMBER = "4 May 2022"
PROJECT_BRIEF = "Documentation of the LAMMPS library interface and Python wrapper"
PROJECT_LOGO = lammps-logo.png
CREATE_SUBDIRS = NO
@ -437,6 +437,8 @@ INPUT = @LAMMPS_SOURCE_DIR@/utils.cpp \
@LAMMPS_SOURCE_DIR@/math_eigen.h \
@LAMMPS_SOURCE_DIR@/platform.h \
@LAMMPS_SOURCE_DIR@/platform.cpp \
@LAMMPS_SOURCE_DIR@/math_special.h \
@LAMMPS_SOURCE_DIR@/math_special.cpp \
# The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
# directories that are symbolic links (a Unix file system feature) are excluded

View File

@ -245,6 +245,8 @@ OPT.
* :doc:`resquared (go) <pair_resquared>`
* :doc:`saip/metal <pair_saip_metal>`
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>`
* :doc:`smatb <pair_smatb>`
* :doc:`smatb/single <pair_smatb>`
* :doc:`smd/hertz <pair_smd_hertz>`
* :doc:`smd/tlsph <pair_smd_tlsph>`
* :doc:`smd/tri_surface <pair_smd_triangulated_surface>`

View File

@ -246,6 +246,44 @@ Customized standard functions
---------------------------
Special Math functions
----------------------
The ``MathSpecial`` namespace implements a selection of custom and optimized
mathematical functions for a variety of applications.
.. doxygenfunction:: factorial
:project: progguide
.. doxygenfunction:: exp2_x86
:project: progguide
.. doxygenfunction:: fm_exp
:project: progguide
.. doxygenfunction:: my_erfcx
:project: progguide
.. doxygenfunction:: expmsq
:project: progguide
.. doxygenfunction:: square
:project: progguide
.. doxygenfunction:: cube
:project: progguide
.. doxygenfunction:: powsign
:project: progguide
.. doxygenfunction:: powint
:project: progguide
.. doxygenfunction:: powsinxx
:project: progguide
---------------------------
Tokenizer classes
-----------------

View File

@ -250,9 +250,9 @@ on` comments around that block.
Error or warning messages and explanations (preferred)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. versionchanged:: 27Apr2022
.. versionchanged:: 4May2022
Starting with LAMMPS version 27 April 2022 the LAMMPS developers have
Starting with LAMMPS version 4 May 2022 the LAMMPS developers have
agreed on a new policy for error and warning messages.
Previously, all error and warning strings were supposed to be listed in

View File

@ -2578,18 +2578,20 @@ SMTBQ package
**Contents:**
A pair style which implements a Second Moment Tight Binding model with
QEq charge equilibration (SMTBQ) potential for the description of
ionocovalent bonds in oxides.
Pair styles which implement Second Moment Tight Binding models.
One with QEq charge equilibration (SMTBQ) for the description of
ionocovalent bonds in oxides, and two more as plain SMATB models.
**Authors:** Nicolas Salles, Emile Maras, Olivier Politano, and Robert
Tetot (LAAS-CNRS, France).
**Authors:** SMTBQ: Nicolas Salles, Emile Maras, Olivier Politano, and Robert
Tetot (LAAS-CNRS, France);
SMATB: Daniele Rapetti (Politecnico di Torino)
**Supporting info:**
* src/SMTBQ: filenames -> commands
* src/SMTBQ/README
* :doc:`pair_style smtbq <pair_smtbq>`
* :doc:`pair_style smatb <pair_smatb>`, :doc:`pair_style smatb/single <pair_smatb>`
* examples/PACKAGES/smtbq
----------

View File

@ -429,8 +429,8 @@ whether an extra library is needed to build and use the package:
- n/a
- no
* - :ref:`SMTBQ <PKG-SMTBQ>`
- second moment tight binding potential
- :doc:`pair_style smtbq <pair_smtbq>`
- second moment tight binding potentials
- :doc:`pair_style smtbq <pair_smtbq>` :doc:`pair_style smatb <pair_smatb>`
- PACKAGES/smtbq
- no
* - :ref:`SPH <PKG-SPH>`

131
doc/src/pair_smatb.rst Normal file
View File

@ -0,0 +1,131 @@
.. index:: pair_style smatb
.. index:: pair_style smatb/single
pair_style smatb command
=========================
pair_style smatb/single command
===============================
Syntax
""""""
.. code-block:: LAMMPS
pair_style style args
* style = *smatb*
* args = none
.. parsed-literal::
*smatb*
Examples
""""""""
.. code-block:: LAMMPS
pair_style smatb
pair_coeff 1 1 2.88 10.35 4.178 0.210 1.818 4.07293506 4.9883063257983666
Description
"""""""""""
The *smatb* styles compute the Second Moment Approximation to the Tight Binding
:ref:`(Cyrot) <Cyrot>`, :ref:`(Gupta) <Gupta>`, :ref:`(Rosato) <Rosato>`,
given by
.. math::
E_{i} = \sum_{j,R_{ij}\leq R_{c}} \alpha(R_{ij}) - \sqrt{\sum_{j,R_{ij}\leq R_{c}}\Xi^2(R_{ij})}
:math:`R_{ij}` is the distance between the atom :math:`i` and :math:`j`.
And the two functions :math:`\alpha\left(r\right)` and :math:`\Xi\left(r\right)` are:
.. math::
\alpha\left(r\right)=\left\lbrace\begin{array}{ll}
A e^{-p \left(\frac{r}{R_{0}}-1\right)} & r < R_{sc}\\
a_3\left(r-R_{c}\right)^3+a_4\left(r-R_{c}\right)^4
+a_5\left(r-R_{c}\right)^5& R_{sc} < r < R_{c}
\end{array}
\right.
.. math::
\Xi\left(r\right)=\left\lbrace\begin{array}{ll}
\xi e^{-q \left(\frac{r}{R_{0}}-1\right)} & r < R_{sc}\\
x_3\left(r-R_{c}\right)^3+x_4\left(r-R_{c}\right)^4
+x_5\left(r-R_{c}\right)^5& R_{sc} < r < R_{c}
\end{array}
\right.
The polynomial coefficients :math:`a_3`, :math:`a_4`, :math:`a_5`,
:math:`x_3`, :math:`x_4`, :math:`x_5` are computed by LAMMPS: the two
exponential terms and their first and second derivatives are smoothly
reduced to zero, from the inner cutoff :math:`R_{sc}` to the outer
cutoff :math:`R_{c}`.
Coefficients
""""""""""""
The following coefficients must be defined for each pair of atoms types via the
:doc:`pair_coeff <pair_coeff>` command as in the examples above, or in the data
file or restart files read by the :doc:`read_data <read_data>` or
:doc:`read_restart <read_restart>` commands, or by mixing as described below:
* :math:`R_{0}` (distance units)
* :math:`p` (dimensionless)
* :math:`q` (dimensionless)
* :math:`A` (energy units)
* :math:`\xi` (energy units)
* :math:`R_{cs}` (distance units)
* :math:`R_{c}` (distance units)
Note that: :math:`R_{0}` is the nearest neighbor distance, usually coincides
with the diameter of the atoms
See the :doc:`run_style <run_style>` command for details.
----------
Mixing info
"""""""""""
For atom type pairs I,J and I != J the coefficients are not automatically mixed.
----------
Restrictions
""""""""""""
This pair style is part of the SMTBQ package and is only enabled
if LAMMPS is built with that package. See the :doc:`Build package <Build_package>` page for more info.
These pair potentials require the :doc:`newton <newton>` setting to be "on" for pair interactions.
Related commands
""""""""""""""""
* :doc:`pair_coeff <pair_coeff>`
Default
"""""""
none
----------
.. _Cyrot:
**(Cyrot)** Cyrot-Lackmann and Ducastelle, Phys Rev. B, 4, 2406-2412 (1971).
.. _Gupta:
**(Gupta)** Gupta ,Phys Rev. B, 23, 6265-6270 (1981).
.. _Rosato:
**(Rosato)** Rosato and Guillope and Legrand, Philosophical Magazine A, 59.2, 321-336 (1989).

View File

@ -323,6 +323,8 @@ accelerated styles exist.
* :doc:`resquared <pair_resquared>` - Everaers RE-Squared ellipsoidal potential
* :doc:`saip/metal <pair_saip_metal>` - interlayer potential for hetero-junctions formed with hexagonal 2D materials and metal surfaces
* :doc:`sdpd/taitwater/isothermal <pair_sdpd_taitwater_isothermal>` - smoothed dissipative particle dynamics for water at isothermal conditions
* :doc:`smatb <pair_smatb>` - Second Moment Approximation to the Tight Binding
* :doc:`smatb/single <pair_smatb>` - Second Moment Approximation to the Tight Binding for single-element systems
* :doc:`smd/hertz <pair_smd_hertz>` -
* :doc:`smd/tlsph <pair_smd_tlsph>` -
* :doc:`smd/tri_surface <pair_smd_triangulated_surface>` -

View File

@ -693,6 +693,7 @@ DFT
dftb
dh
dhex
di
dia
diag
diagonalization
@ -784,6 +785,7 @@ dtheta
dtshrink
du
dU
Ducastelle
Dudarev
Duin
Dullweber
@ -890,6 +892,7 @@ Embt
emi
emol
eN
endian
energetics
energyCorr
eng
@ -1067,6 +1070,7 @@ fld
floralwhite
Florez
flv
fm
fmackay
fmag
fmass
@ -1246,6 +1250,7 @@ Gubbins
Guenole
Guericke
gui
Guillope
Gumbsch
Gunsteren
Gunzenmuller
@ -1722,6 +1727,7 @@ lebedeva
Lebedeva
Lebold
Lechman
Legrand
Lehoucq
Leimkuhler
Leite
@ -2636,6 +2642,7 @@ Polarizable
polarizables
polarizer
Politano
Politecnico
polyA
polybond
polydisperse
@ -2794,6 +2801,7 @@ Randisi
randomizations
rann
RANN
Rapetti
Raphson
Rappe
Ravelo
@ -2939,6 +2947,7 @@ Rohart
Ronchetti
Ronevich
Rosati
Rosato
Rosenberger
Rossky
rosybrown
@ -3115,6 +3124,7 @@ smallbig
smallint
Smallint
smallsmall
smatb
smd
SMD
smi
@ -3383,6 +3393,7 @@ toolchain
topologies
Toporov
Torder
Torino
torsions
Tosi
Toukmaji

View File

@ -0,0 +1,48 @@
# AgCu Pancake : energy should be around -90.16 eV
34 atoms
2 atom types
0 30 xlo xhi
0 30 ylo yhi
0 30 zlo zhi
Masses
1 108 # Ag
2 64 # Cu
Atoms # atomic
1 1 11.8677744 17.4748811 16.8202155
2 1 14.4591543 12.6388264 17.3769114
3 1 19.2905996 14.8698601 15.9074284
4 1 13.6418392 11.2583912 15.0376329
5 1 16.5295136 18.8875825 16.3408808
6 1 11.4394217 12.6680604 13.934792
7 1 17.1772792 13.4579369 17.0971284
8 1 11.7477198 12.5836832 16.6835448
9 1 12.3254647 15.127665 12.9151285
10 1 11.5595413 17.5592601 14.0713777
11 1 18.9820568 14.9536515 13.1587506
12 1 14.5653354 15.1189885 11.12255
13 1 14.5797485 17.250847 17.5049696
14 1 13.9305528 12.7833817 12.6633235
15 1 12.8538576 14.9830633 17.6286053
16 1 16.7229337 16.4525301 12.4627193
17 1 14.0512024 17.3954338 12.7914031
18 1 18.2824041 17.2316517 14.6999911
19 1 16.2210303 18.9719703 13.5921161
20 1 17.2517138 16.3080477 17.1761679
21 1 16.0269821 11.057668 13.3712031
22 1 16.3354169 10.9735982 16.1199203
23 1 13.9149278 18.8223918 15.2392839
24 1 18.1136886 12.5568867 14.5759053
25 1 16.6484505 13.6019749 12.3834541
26 1 11.0468778 15.1307677 15.4471491
27 1 15.4347551 14.8811145 18.8773723
28 2 15.1502129 14.9589199 16.3394186
29 2 13.2837251 13.9076768 15.1588849
30 2 13.4425035 16.3005712 15.2145264
31 2 15.7534979 16.8960911 14.9738557
32 2 15.4967806 13.024304 14.8837423
33 2 14.8498337 15.0410649 13.6605697
34 2 17.0232044 14.8712576 14.7690776

View File

@ -0,0 +1,38 @@
# -*- lammps -*-
units metal
atom_style atomic
boundary p p p
read_data AgCuPancake.data
pair_style smatb
# NN p q a qsi cutOff_Start cutOff_End
pair_coeff 1 1 2.89 10.85 3.18 0.1031 1.1895 4.08707719 5.0056268338740553
pair_coeff 1 2 2.725 10.70 2.805 0.0977 1.2275 4.08707719 4.4340500673763259
pair_coeff 2 2 2.56 10.55 2.43 0.0894 1.2799 3.62038672 4.4340500673763259
neighbor 8.0 bin
neigh_modify every 1 delay 0 check yes
thermo 1
minimize 1.0e-8 1.0e-10 10000 100000
velocity all create 600.0 761341 rot yes mom yes
fix 1 all nve
thermo 10
timestep 0.005
#dump 1 all atom 50 dump.smatb
#dump 2 all image 10 image.*.jpg element element &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 2 pad 3 element Ag Cu
#dump 3 all movie 10 movie.mpg element element &
# axes yes 0.8 0.02 view 60 -30
#dump_modify 3 pad 3 element Ag Cu
run 10000

View File

@ -0,0 +1,27 @@
# -*- lammps -*-
units metal
atom_style atomic
boundary p p p
lattice fcc 4.0782
region myreg block 0 8 0 8 0 8
create_box 1 myreg
create_atoms 1 box
mass 1 196.96655 # Au
pair_style smatb/single
pair_coeff 1 1 2.88 10.35 4.178 0.210 1.818 4.07293506 4.9883063257983666
neighbor 8.0 bin
neigh_modify every 1 delay 0 check yes
thermo 1
fix boxmin all box/relax iso 1.0
minimize 1.0e-8 1.0e-10 10000 100000
unfix boxmin
minimize 1.0e-8 1.0e-10 10000 100000

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,160 @@
LAMMPS (27 Oct 2021)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
# -*- lammps -*-
units metal
atom_style atomic
boundary p p p
lattice fcc 4.0782
Lattice spacing in x,y,z = 4.0782000 4.0782000 4.0782000
region myreg block 0 8 0 8 0 8
create_box 1 myreg
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (32.625600 32.625600 32.625600)
1 by 1 by 1 MPI processor grid
create_atoms 1 box
Created 2048 atoms
using lattice units in orthogonal box = (0.0000000 0.0000000 0.0000000) to (32.625600 32.625600 32.625600)
create_atoms CPU = 0.001 seconds
mass 1 196.96655 # Au
pair_style smatb/single
pair_coeff 1 1 2.88 10.35 4.178 0.210 1.818 4.07293506 4.9883063257983666
neighbor 8.0 bin
neigh_modify every 1 delay 0 check yes
thermo 1
fix boxmin all box/relax iso 1.0
minimize 1.0e-8 1.0e-10 10000 100000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12.988306
ghost atom cutoff = 12.988306
binsize = 6.4941532, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair smatb/single, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.601 | 6.601 | 6.601 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 0 -7800.9629 0 -7800.9629 -17598.853 34727.66
1 0 -7801.0757 0 -7801.0757 -17102.698 34717.243
2 0 -7801.1852 0 -7801.1852 -16605.672 34706.828
3 0 -7801.2915 0 -7801.2915 -16107.773 34696.415
4 0 -7801.3946 0 -7801.3946 -15609 34686.004
5 0 -7801.4944 0 -7801.4944 -15109.353 34675.595
6 0 -7801.5909 0 -7801.5909 -14608.829 34665.188
7 0 -7801.6841 0 -7801.6841 -14107.429 34654.783
8 0 -7801.7741 0 -7801.7741 -13605.15 34644.38
9 0 -7801.8608 0 -7801.8608 -13101.992 34633.98
10 0 -7801.9442 0 -7801.9442 -12597.953 34623.581
11 0 -7802.0243 0 -7802.0243 -12093.033 34613.185
12 0 -7802.1011 0 -7802.1011 -11587.23 34602.79
13 0 -7802.1746 0 -7802.1746 -11080.543 34592.398
14 0 -7802.2448 0 -7802.2448 -10572.902 34582.008
15 0 -7802.3117 0 -7802.3117 -10064.258 34571.62
16 0 -7802.3753 0 -7802.3753 -9554.6096 34561.234
17 0 -7802.4356 0 -7802.4356 -9043.9555 34550.85
18 0 -7802.4925 0 -7802.4925 -8532.2942 34540.468
19 0 -7802.5462 0 -7802.5462 -8019.6245 34530.088
20 0 -7802.5964 0 -7802.5964 -7505.945 34519.711
21 0 -7802.6434 0 -7802.6434 -6991.2543 34509.335
22 0 -7802.687 0 -7802.687 -6475.5513 34498.961
23 0 -7802.7272 0 -7802.7272 -5958.8344 34488.59
24 0 -7802.7641 0 -7802.7641 -5441.1024 34478.221
25 0 -7802.7977 0 -7802.7977 -4922.354 34467.853
26 0 -7802.8278 0 -7802.8278 -4402.5878 34457.488
27 0 -7802.8546 0 -7802.8546 -3881.8024 34447.125
28 0 -7802.878 0 -7802.878 -3359.9966 34436.764
29 0 -7802.8981 0 -7802.8981 -2837.1689 34426.405
30 0 -7802.9147 0 -7802.9147 -2313.3181 34416.048
31 0 -7802.928 0 -7802.928 -1788.4427 34405.693
32 0 -7802.9378 0 -7802.9378 -1262.5414 34395.34
33 0 -7802.9443 0 -7802.9443 -735.61295 34384.99
34 0 -7802.9473 0 -7802.9473 -207.6559 34374.641
35 0 -7802.9476 0 -7802.9476 0.90227419 34370.559
36 0 -7802.9476 0 -7802.9476 0.99992446 34370.557
Loop time of 0.142744 on 1 procs for 36 steps with 2048 atoms
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-7800.9628521055 -7802.94781441221 -7802.94781442797
Force two-norm initial, final = 1144.4464 4.8784902e-06
Force max component initial, final = 1144.4464 4.8784902e-06
Final line search alpha, max atom move = 0.015845171 7.7300512e-08
Iterations, force evaluations = 36 38
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.13461 | 0.13461 | 0.13461 | 0.0 | 94.30
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.0012706 | 0.0012706 | 0.0012706 | 0.0 | 0.89
Output | 0.00066993 | 0.00066993 | 0.00066993 | 0.0 | 0.47
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.006191 | | | 4.34
Nlocal: 2048.00 ave 2048 max 2048 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 10147.0 ave 10147 max 10147 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 567296.0 ave 567296 max 567296 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 567296
Ave neighs/atom = 277.00000
Neighbor list builds = 0
Dangerous builds = 0
unfix boxmin
minimize 1.0e-8 1.0e-10 10000 100000
Per MPI rank memory allocation (min/avg/max) = 6.601 | 6.601 | 6.601 Mbytes
Step Temp E_pair E_mol TotEng Press
36 0 -7802.9476 0 -7802.9476 0.99992446
37 0 -7802.9476 0 -7802.9476 0.99992446
Loop time of 0.0105782 on 1 procs for 1 steps with 2048 atoms
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-7802.94759154184 -7802.94759154184 -7802.94759154185
Force two-norm initial, final = 4.7040841e-12 1.3779243e-12
Force max component initial, final = 1.1096422e-13 4.1164848e-14
Final line search alpha, max atom move = 1.0000000 4.1164848e-14
Iterations, force evaluations = 1 2
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.010394 | 0.010394 | 0.010394 | 0.0 | 98.25
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 8.6608e-05 | 8.6608e-05 | 8.6608e-05 | 0.0 | 0.82
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 9.804e-05 | | | 0.93
Nlocal: 2048.00 ave 2048 max 2048 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 10147.0 ave 10147 max 10147 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 567296.0 ave 567296 max 567296 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 567296
Ave neighs/atom = 277.00000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00

4
src/.gitignore vendored
View File

@ -1537,6 +1537,10 @@
/pair_mgpt.h
/pair_morse_smooth_linear.cpp
/pair_morse_smooth_linear.h
/pair_smatb.cpp
/pair_smatb.h
/pair_smatb_single.cpp
/pair_smatb_single.h
/pair_smtbq.cpp
/pair_smtbq.h
/pair_vashishta*.cpp

View File

@ -319,8 +319,7 @@ void PairEAM::compute(int eflag, int vflag)
}
if (eflag) evdwl = scale[itype][jtype]*phi;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
}
}

View File

@ -1,16 +0,0 @@
This package implements the Second Moment Tight Binding - QEq (SMTB-Q)
potential for the description of ionocovalent bonds in oxides.
Authors: Nicolas Salles, Emile Maras, Olivier Politano, Robert Tetot
at ICB, Universite de Bourgogne and ICMMO, Universite Paris-Sud.
Contact emails: lammps@u-bourgogne.fr, nsalles33@gmail.com
This package is occasionally maintained.
See the doc page for the pair_style smtbq command to get started.
There are potential files for this potential in the potentials dir.
There are example scripts for using this package in
examples/PACKAGES/smtbq.

558
src/SMTBQ/pair_smatb.cpp Normal file
View File

@ -0,0 +1,558 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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: Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#include "pair_smatb.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairSMATB::PairSMATB(LAMMPS *_lmp) :
Pair(_lmp), nmax(0), on_eb(nullptr), r0(nullptr), p(nullptr), A(nullptr), q(nullptr),
QSI(nullptr), cutOffStart(nullptr), cutOffEnd(nullptr), cutOffEnd2(nullptr), a3(nullptr),
a4(nullptr), a5(nullptr), x3(nullptr), x4(nullptr), x5(nullptr)
{
single_enable = 0; // 1 if single() routine exists
restartinfo = 1; // 1 if pair style writes restart info
respa_enable = 0; // 1 if inner/middle/outer rRESPA routines
one_coeff = 0; // 1 if allows only one coeff * * call
manybody_flag = 1; // 1 if a manybody potential
no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
writedata = 1; // 1 if writes coeffs to data file
ghostneigh = 0; // 1 if pair style needs neighbors of ghosts
// set comm size needed by this Pair
comm_forward = 1;
comm_reverse = 1;
}
/* ---------------------------------------------------------------------- */
PairSMATB::~PairSMATB()
{
if (copymode) { return; }
memory->destroy(on_eb);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(r0);
memory->destroy(p);
memory->destroy(A);
memory->destroy(q);
memory->destroy(QSI);
memory->destroy(cutOffStart);
memory->destroy(cutOffEnd);
memory->destroy(cutOffEnd2);
memory->destroy(a3);
memory->destroy(a4);
memory->destroy(a5);
memory->destroy(x5);
memory->destroy(x4);
memory->destroy(x3);
}
}
/* ---------------------------------------------------------------------- */
void PairSMATB::compute(int eflag, int vflag)
{
int i, j, ii, jj, jnum, itype, jtype;
double xtmp, ytmp, ztmp, del[3], fpair;
double dijsq, dij;
double espo, aexpp, qsiexpq, eb_i, Fb, Fr;
double polyval, polyval2, polyval3, polyval4, polyval5;
if (eflag || vflag) {
ev_setup(eflag, vflag);
eng_vdwl = 0;
} else {
evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
}
// grow on_eb array if necessary
if (atom->nmax > nmax) {
nmax = atom->nmax;
memory->grow(on_eb, nmax, "pair_smatb:on_eb");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair;
// zero out on_eb
memset(on_eb, 0, nall * sizeof(double));
int inum = list->inum;
int *ilist = list->ilist;
int *jlist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
// FIRST LOOP: CALCULATES the squared bonding energy and accumulate it in on_eb for each atom
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; ++jj) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
del[0] = xtmp - x[j][0];
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if (dijsq < cutOffEnd2[itype][jtype]) {
dij = sqrt(dijsq);
if (dij < cutOffStart[itype][jtype]) {
qsiexpq = (QSI[itype][jtype] * QSI[itype][jtype]) *
exp(2.0 * q[itype][jtype] * (1.0 - dij / r0[itype][jtype]));
} else {
polyval = dij - cutOffEnd[itype][jtype];
polyval3 = polyval * polyval * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
qsiexpq = x5[itype][jtype] * polyval5 + x4[itype][jtype] * polyval4 +
x3[itype][jtype] * polyval3;
qsiexpq = qsiexpq * qsiexpq;
}
on_eb[i] += qsiexpq;
if (newton_pair) on_eb[j] += qsiexpq;
}
}
}
// communicate the squared bonding energy between the various bins
if (newton_pair) comm->reverse_comm(this);
// Support Loop: take the square root of the bonding energy and
// accumulate it in the energy accumulator if needed the store the
// reciprocal in on_eb in order to not do it in the SECOND LOOP
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (i < nlocal) {
eb_i = sqrt(on_eb[i]);
if (eb_i != 0.0) {
on_eb[i] = 1.0 / eb_i;
} else {
on_eb[i] = 0.0;
}
// if needed the bonding energy is accumulated:
if (eflag_either) {
if (eflag_atom) { eatom[i] -= eb_i; }
if (eflag_global) { eng_vdwl -= eb_i; }
}
}
}
// this communication stores the denominators in the ghosts atoms,
// this is needed because of how forces are calculated
comm->forward_comm(this);
// SECOND LOOP: given on_eb[i] calculates forces and energies
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
del[0] = xtmp - x[j][0];
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if (dijsq < cutOffEnd2[itype][jtype]) {
dij = sqrt(dijsq);
if (dij < cutOffStart[itype][jtype]) {
espo = 1.0 - dij / r0[itype][jtype];
aexpp = exp(p[itype][jtype] * espo) * A[itype][jtype];
Fr = (2.0 * aexpp) * (p[itype][jtype] / r0[itype][jtype]);
qsiexpq = (QSI[itype][jtype] * QSI[itype][jtype]) * exp(2.0 * q[itype][jtype] * espo);
Fb = -qsiexpq * q[itype][jtype] / r0[itype][jtype];
} else {
polyval = dij - cutOffEnd[itype][jtype];
polyval2 = polyval * polyval;
polyval3 = polyval2 * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
aexpp = a5[itype][jtype] * polyval5 + a4[itype][jtype] * polyval4 +
a3[itype][jtype] * polyval3;
Fr = -2.0 *
(5.0 * a5[itype][jtype] * polyval4 + 4.0 * a4[itype][jtype] * polyval3 +
3.0 * a3[itype][jtype] * polyval2);
qsiexpq = x5[itype][jtype] * polyval5 + x4[itype][jtype] * polyval4 +
x3[itype][jtype] * polyval3;
Fb = ((5.0 * x5[itype][jtype] * polyval4 + 4.0 * x4[itype][jtype] * polyval3 +
3.0 * x3[itype][jtype] * polyval2)) *
qsiexpq;
}
// calculates the module of the pair energy between i and j
fpair = (Fb * (on_eb[i] + on_eb[j]) + Fr) / dij;
f[i][0] += del[0] * fpair;
f[i][1] += del[1] * fpair;
f[i][2] += del[2] * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= del[0] * fpair;
f[j][1] -= del[1] * fpair;
f[j][2] -= del[2] * fpair;
}
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, 2.0 * aexpp, 0.0, fpair, del[0], del[1], del[2]);
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSMATB::settings(int narg, char **)
{
if (narg > 0) error->all(FLERR, "Illegal pair_style command: smatb accepts no options");
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSMATB::allocate()
{
const int np1 = atom->ntypes + 1;
memory->create(setflag, np1, np1, "pair_smatb:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq, np1, np1, "pair_smatb:cutsq");
memory->create(r0, np1, np1, "pair_smatb:r0");
memory->create(p, np1, np1, "pair_smatb:p");
memory->create(A, np1, np1, "pair_smatb:A");
memory->create(q, np1, np1, "pair_smatb:q");
memory->create(QSI, np1, np1, "pair_smatb:QSI");
memory->create(cutOffStart, np1, np1, "pair_smatb:cutOffStart");
memory->create(cutOffEnd, np1, np1, "pair_smatb:cutOffEnd");
memory->create(cutOffEnd2, np1, np1, "pair_smatb:cutOffEnd2");
memory->create(a3, np1, np1, "pair_smatb:a1");
memory->create(a4, np1, np1, "pair_smatb:a2");
memory->create(a5, np1, np1, "pair_smatb:a5");
memory->create(x3, np1, np1, "pair_smatb:x1");
memory->create(x4, np1, np1, "pair_smatb:x2");
memory->create(x5, np1, np1, "pair_smatb:x3");
allocated = 1;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSMATB::coeff(int narg, char **arg)
{
if (!allocated) { allocate(); }
if (narg != 9) utils::missing_cmd_args(FLERR, "pair_style smatb", error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double myr0 = utils::numeric(FLERR, arg[2], false, lmp);
double myp = utils::numeric(FLERR, arg[3], false, lmp);
double myq = utils::numeric(FLERR, arg[4], false, lmp);
double myA = utils::numeric(FLERR, arg[5], false, lmp);
double myQSI = utils::numeric(FLERR, arg[6], false, lmp);
double mycutOffStart = utils::numeric(FLERR, arg[7], false, lmp);
double mycutOffEnd = utils::numeric(FLERR, arg[8], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
r0[i][j] = myr0;
p[i][j] = myp;
A[i][j] = myA;
q[i][j] = myq;
QSI[i][j] = myQSI;
cutOffStart[i][j] = mycutOffStart;
cutOffEnd[i][j] = mycutOffEnd;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ------------------------------------------------------------------------ */
void PairSMATB::init_style()
{
if (force->newton_pair == 0) error->all(FLERR, "Pair style smatb requires newton pair on");
neighbor->add_request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSMATB::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
///@todo implement smatb mixing rules
cutOffStart[i][j] = MIN(cutOffStart[i][i], cutOffStart[j][j]);
cutOffEnd[i][j] = MAX(cutOffEnd[i][i], cutOffEnd[j][j]);
error->all(FLERR, "All pair coeffs are not set");
}
double es = cutOffEnd[i][j] - cutOffStart[i][j];
double es2 = es * es;
double es3 = es2 * es;
// variables for poly for p and A
double expp = A[i][j] * exp(p[i][j] * (1. - cutOffStart[i][j] / r0[i][j]));
double ap = -1. / es3;
double bp = p[i][j] / (r0[i][j] * es2);
double cp = -(p[i][j] * p[i][j]) / (es * r0[i][j] * r0[i][j]);
a5[i][j] = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
a4[i][j] = expp * (15. * ap + 7. * bp + cp) / es;
a3[i][j] = expp * (20. * ap + 8. * bp + cp) / 2.;
// variables for poly for q and qsi
double expq = QSI[i][j] * exp(q[i][j] * (1. - cutOffStart[i][j] / r0[i][j]));
double aq = -1 / es3;
double bq = q[i][j] / (es2 * r0[i][j]);
double cq = -(q[i][j] * q[i][j]) / (es * r0[i][j] * r0[i][j]);
x5[i][j] = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
x4[i][j] = expq * (15. * aq + 7. * bq + cq) / es;
x3[i][j] = expq * (20. * aq + 8. * bq + cq) / 2.;
cutOffEnd2[i][j] = cutOffEnd[i][j] * cutOffEnd[i][j];
if (i != j) {
setflag[j][i] = 1;
cutOffEnd2[j][i] = cutOffEnd2[i][j];
r0[j][i] = r0[i][j];
p[j][i] = p[i][j];
q[j][i] = q[i][j];
A[j][i] = A[i][j];
QSI[j][i] = QSI[i][j];
cutOffStart[j][i] = cutOffStart[i][j];
cutOffEnd[j][i] = cutOffEnd[i][j];
a3[j][i] = a3[i][j];
a4[j][i] = a4[i][j];
a5[j][i] = a5[i][j];
x3[j][i] = x3[i][j];
x4[j][i] = x4[i][j];
x5[j][i] = x5[i][j];
}
return cutOffEnd[i][j];
}
/* ---------------------------------------------------------------------- */
int PairSMATB::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i, j, m;
m = 0;
for (i = 0; i < n; ++i) {
j = list[i];
buf[m++] = on_eb[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairSMATB::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) { on_eb[i] = buf[m++]; }
}
/* ---------------------------------------------------------------------- */
int PairSMATB::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) { buf[m++] = on_eb[i]; }
return m;
}
/* ---------------------------------------------------------------------- */
void PairSMATB::unpack_reverse_comm(int n, int *list, double *buf)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
on_eb[j] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
void PairSMATB::write_restart_settings(FILE *fp)
{
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag, sizeof(int), 1, fp);
}
/* ---------------------------------------------------------------------- */
void PairSMATB::read_restart_settings(FILE *fp)
{
int me = comm->me;
size_t result;
if (me == 0) {
result = fread(&offset_flag, sizeof(int), 1, fp);
result = fread(&mix_flag, sizeof(int), 1, fp);
result = fread(&tail_flag, sizeof(int), 1, fp);
}
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);
}
/* ---------------------------------------------------------------------- */
void PairSMATB::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(&r0[i][j], sizeof(double), 1, fp);
fwrite(&p[i][j], sizeof(double), 1, fp);
fwrite(&q[i][j], sizeof(double), 1, fp);
fwrite(&A[i][j], sizeof(double), 1, fp);
fwrite(&QSI[i][j], sizeof(double), 1, fp);
fwrite(&cutOffStart[i][j], sizeof(double), 1, fp);
fwrite(&cutOffEnd[i][j], sizeof(double), 1, fp);
}
}
}
}
/* ---------------------------------------------------------------------- */
void PairSMATB::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
size_t result;
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) { result = 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) {
utils::sfread(FLERR, &r0[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &p[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &q[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &A[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &QSI[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cutOffStart[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cutOffEnd[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&r0[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&p[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&q[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&A[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&QSI[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffStart[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffEnd[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
/* ---------------------------------------------------------------------- */
void PairSMATB::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++) {
fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, r0[i][i], p[i][i], q[i][i], A[i][i], QSI[i][i],
cutOffStart[i][i], cutOffEnd[i][i]);
}
}
/* ---------------------------------------------------------------------- */
void PairSMATB::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 %g\n", i, j, r0[i][j], p[i][j], q[i][j], A[i][j],
QSI[i][j], cutOffStart[i][j], cutOffEnd[i][j]);
}
}
}

80
src/SMTBQ/pair_smatb.h Normal file
View File

@ -0,0 +1,80 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This pair style is written by Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(smatb,PairSMATB);
// clang-format on
#else
#ifndef LMP_PAIR_SMATB_H
#define LMP_PAIR_SMATB_H
#include "pair.h"
namespace LAMMPS_NS {
class PairSMATB : public Pair {
public:
PairSMATB(class LAMMPS *);
~PairSMATB() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
protected:
virtual void allocate();
// allocated size of per-atom arrays
int nmax;
//allocated to store up calculation values
double *on_eb;
// interaction radius, user-given
double **r0;
// parameters user-given
double **p;
double **A;
double **q;
double **QSI;
//extremes of the cut off, user given
double **cutOffStart;
double **cutOffEnd;
//squared cut off end, calculated
double **cutOffEnd2;
//polynomial for cutoff linking to zero: Ae^p substitution
double **a3;
double **a4;
double **a5;
//polynomial for cutoff linking to zero: QSIe^q substitution
double **x3;
double **x4;
double **x5;
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -0,0 +1,503 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, 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: Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#include "pair_smatb_single.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairSMATBSingle::PairSMATBSingle(LAMMPS *_lmp) :
Pair(_lmp), nmax(0), on_eb(nullptr), r0(0), p(0), A(0), q(0), QSI(0), cutOffStart(0),
cutOffEnd(0), cutOffEnd2(0), a3(0), a4(0), a5(0), x3(0), x4(0), x5(0)
{
single_enable = 0; // 1 if single() routine exists
restartinfo = 1; // 1 if pair style writes restart info
respa_enable = 0; // 1 if inner/middle/outer rRESPA routines
one_coeff = 0; // 1 if allows only one coeff * * call
manybody_flag = 1; // 1 if a manybody potential
no_virial_fdotr_compute = 0; // 1 if does not invoke virial_fdotr_compute()
writedata = 1; // 1 if writes coeffs to data file
ghostneigh = 0; // 1 if pair style needs neighbors of ghosts
// set comm size needed by this Pair
comm_forward = 1;
comm_reverse = 1;
}
/* ---------------------------------------------------------------------- */
PairSMATBSingle::~PairSMATBSingle()
{
if (copymode) { return; }
memory->destroy(on_eb);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
}
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::compute(int eflag, int vflag)
{
int i, j, ii, jj, jnum;
double xtmp, ytmp, ztmp, del[3], fpair;
double dijsq, dij;
double espo, aexpp, qsiexpq, eb_i, Fb, Fr;
double polyval, polyval2, polyval3, polyval4, polyval5;
if (eflag || vflag) {
ev_setup(eflag, vflag);
eng_vdwl = 0;
} else {
evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
}
// grow on_eb array if necessary
if (atom->nmax > nmax) {
nmax = atom->nmax;
memory->grow(on_eb, nmax, "pair_smatb:on_eb");
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair;
// zero out on_eb
memset(on_eb, 0, nall * sizeof(double));
int inum = list->inum;
int *ilist = list->ilist;
int *jlist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
// FIRST LOOP: CALCULATES the squared bonding energy and accumulate it in on_eb for each atom
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; ++jj) {
j = jlist[jj];
j &= NEIGHMASK;
del[0] = xtmp - x[j][0];
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if (dijsq < cutOffEnd2) {
dij = sqrt(dijsq);
if (dij < cutOffStart) {
qsiexpq = (QSI * QSI) * exp(2.0 * q * (1.0 - dij / r0));
} else {
polyval = dij - cutOffEnd;
polyval3 = polyval * polyval * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
qsiexpq = x5 * polyval5 + x4 * polyval4 + x3 * polyval3;
qsiexpq = qsiexpq * qsiexpq;
}
on_eb[i] += qsiexpq;
if (newton_pair) on_eb[j] += qsiexpq;
}
}
}
// communicate the squared bonding energy between the various bins
if (newton_pair) comm->reverse_comm(this);
// Support Loop: take the square root of the bonding energy and
// accumulate it in the energy accumulator if needed the store the
// reciprocal in on_eb in order to not do it in the SECOND LOOP
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
if (i < nlocal) {
eb_i = sqrt(on_eb[i]);
if (eb_i != 0.0) {
on_eb[i] = 1.0 / eb_i;
} else {
on_eb[i] = 0.0;
}
// if needed the bonding energy is accumulated:
if (eflag_either) {
if (eflag_atom) { eatom[i] -= eb_i; }
if (eflag_global) { eng_vdwl -= eb_i; }
}
}
}
// this communication stores the denominators in the ghosts atoms,
// this is needed because of how forces are calculated
comm->forward_comm(this);
// SECOND LOOP: given on_eb[i] calculates forces and energies
for (ii = 0; ii < inum; ++ii) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
del[0] = xtmp - x[j][0];
del[1] = ytmp - x[j][1];
del[2] = ztmp - x[j][2];
dijsq = del[0] * del[0] + del[1] * del[1] + del[2] * del[2];
if (dijsq < cutOffEnd2) {
dij = sqrt(dijsq);
if (dij < cutOffStart) {
espo = 1.0 - dij / r0;
aexpp = exp(p * espo) * A;
Fr = (2.0 * aexpp) * (p / r0);
qsiexpq = (QSI * QSI) * exp(2.0 * q * espo);
Fb = -qsiexpq * q / r0;
} else {
polyval = dij - cutOffEnd;
polyval2 = polyval * polyval;
polyval3 = polyval2 * polyval;
polyval4 = polyval3 * polyval;
polyval5 = polyval4 * polyval;
aexpp = a5 * polyval5 + a4 * polyval4 + a3 * polyval3;
Fr = -2.0 * (5.0 * a5 * polyval4 + 4.0 * a4 * polyval3 + 3.0 * a3 * polyval2);
qsiexpq = x5 * polyval5 + x4 * polyval4 + x3 * polyval3;
Fb = ((5.0 * x5 * polyval4 + 4.0 * x4 * polyval3 + 3.0 * x3 * polyval2)) * qsiexpq;
}
// calculates the module of the pair energy between i and j
fpair = (Fb * (on_eb[i] + on_eb[j]) + Fr) / dij;
f[i][0] += del[0] * fpair;
f[i][1] += del[1] * fpair;
f[i][2] += del[2] * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= del[0] * fpair;
f[j][1] -= del[1] * fpair;
f[j][2] -= del[2] * fpair;
}
if (evflag) {
ev_tally(i, j, nlocal, newton_pair, 2.0 * aexpp, 0.0, fpair, del[0], del[1], del[2]);
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSMATBSingle::settings(int narg, char **)
{
if (narg > 0) error->all(FLERR, "Illegal pair_style command: smatb/single accepts no options");
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSMATBSingle::allocate()
{
int n = atom->ntypes;
int natoms = atom->natoms;
memory->create(setflag, n + 1, n + 1, "pair_smatb: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_smatb:cutsq");
allocated = 1;
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSMATBSingle::coeff(int narg, char **arg)
{
if (!allocated) { allocate(); }
if (narg != 9) utils::missing_cmd_args(FLERR, "pair_style smatb/single", error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
r0 = utils::numeric(FLERR, arg[2], false, lmp);
p = utils::numeric(FLERR, arg[3], false, lmp);
q = utils::numeric(FLERR, arg[4], false, lmp);
A = utils::numeric(FLERR, arg[5], false, lmp);
QSI = utils::numeric(FLERR, arg[6], false, lmp);
cutOffStart = utils::numeric(FLERR, arg[7], false, lmp);
cutOffEnd = utils::numeric(FLERR, arg[8], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ------------------------------------------------------------------------ */
void PairSMATBSingle::init_style()
{
if (force->newton_pair == 0) error->all(FLERR, "Pair style smatb/single requires newton pair on");
neighbor->add_request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSMATBSingle::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
//calculating the polynomial linking to zero
double es = cutOffEnd - cutOffStart;
double es2 = es * es;
double es3 = es2 * es;
//variables for poly for p and A
double expp = A * exp(p * (1. - cutOffStart / r0));
double ap = -1. / es3;
double bp = p / (r0 * es2);
double cp = -(p * p) / (es * r0 * r0);
a5 = expp * (12. * ap + 6. * bp + cp) / (2. * es2);
a4 = expp * (15. * ap + 7. * bp + cp) / es;
a3 = expp * (20. * ap + 8. * bp + cp) / 2.;
//variables for poly for q and qsi
double expq = QSI * exp(q * (1. - cutOffStart / r0));
double aq = -1 / es3;
double bq = q / (es2 * r0);
double cq = -(q * q) / (es * r0 * r0);
x5 = expq * (12. * aq + 6. * bq + cq) / (2. * es2);
x4 = expq * (15. * aq + 7. * bq + cq) / es;
x3 = expq * (20. * aq + 8. * bq + cq) / 2.;
cutOffEnd2 = cutOffEnd * cutOffEnd;
if (i != j) {
setflag[j][i] = 1;
cutOffEnd2 = cutOffEnd2;
r0 = r0;
p = p;
q = q;
A = A;
QSI = QSI;
cutOffStart = cutOffStart;
cutOffEnd = cutOffEnd;
a3 = a3;
a4 = a4;
a5 = a5;
x3 = x3;
x4 = x4;
x5 = x5;
}
return cutOffEnd;
}
/* ---------------------------------------------------------------------- */
int PairSMATBSingle::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i, j, m;
m = 0;
for (i = 0; i < n; ++i) {
j = list[i];
buf[m++] = on_eb[j];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::unpack_forward_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) { on_eb[i] = buf[m++]; }
}
/* ---------------------------------------------------------------------- */
int PairSMATBSingle::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; ++i) { buf[m++] = on_eb[i]; }
return m;
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::unpack_reverse_comm(int n, int *list, double *buf)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
on_eb[j] += buf[m++];
}
}
/* ---------------------------------------------------------------------- */
//write binary data of this simulation:
void PairSMATBSingle::write_restart_settings(FILE *fp)
{
fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag, sizeof(int), 1, fp);
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::read_restart_settings(FILE *fp)
{
int me = comm->me;
size_t result;
if (me == 0) {
result = fread(&offset_flag, sizeof(int), 1, fp);
result = fread(&mix_flag, sizeof(int), 1, fp);
result = fread(&tail_flag, sizeof(int), 1, fp);
}
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);
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::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(&r0, sizeof(double), 1, fp);
fwrite(&p, sizeof(double), 1, fp);
fwrite(&q, sizeof(double), 1, fp);
fwrite(&A, sizeof(double), 1, fp);
fwrite(&QSI, sizeof(double), 1, fp);
fwrite(&cutOffStart, sizeof(double), 1, fp);
fwrite(&cutOffEnd, sizeof(double), 1, fp);
}
}
}
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
size_t result;
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) { result = 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) {
utils::sfread(FLERR, &r0, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &p, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &q, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &A, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &QSI, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cutOffStart, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cutOffEnd, sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&r0, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&p, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&q, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&A, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&QSI, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffStart, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cutOffEnd, 1, MPI_DOUBLE, 0, world);
}
}
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++) {
fprintf(fp, "%d %g %g %g %g %g %g %g\n", i, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
}
}
/* ---------------------------------------------------------------------- */
void PairSMATBSingle::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 %g\n", i, j, r0, p, q, A, QSI, cutOffStart, cutOffEnd);
}
}
}

View File

@ -0,0 +1,80 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
This pair style is written by Daniele Rapetti (iximiel@gmail.com)
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(smatb/single,PairSMATBSingle);
// clang-format on
#else
#ifndef LMP_PAIR_SMATB_SINGLE_H
#define LMP_PAIR_SMATB_SINGLE_H
#include "pair.h"
namespace LAMMPS_NS {
class PairSMATBSingle : public Pair {
public:
PairSMATBSingle(class LAMMPS *);
~PairSMATBSingle() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
protected:
virtual void allocate();
// allocated size of per-atom arrays
int nmax;
//allocated to store up calculation values
double *on_eb{nullptr};
// interaction radius, user-given
double r0;
// parameters user-given
double p;
double A;
double q;
double QSI;
//cut offs, user given
double cutOffStart;
double cutOffEnd;
//squared cut off end, calculated
double cutOffEnd2;
//polynomial for cutoff linking to zero: Ae^p substitution
double a3;
double a4;
double a5;
//polynomial for cutoff linking to zero: QSIe^q substitution
double x3;
double x4;
double x5;
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -706,6 +706,7 @@ static const double fm_exp2_p[] = {
double MathSpecial::exp2_x86(double x)
{
#if defined(__BYTE_ORDER__) && (__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__)
double ipart, fpart, px, qx;
udi_t epart;
@ -726,6 +727,9 @@ double MathSpecial::exp2_x86(double x)
x = 1.0 + 2.0*(px/(qx-px));
return epart.f*x;
#else
return pow(2.0, x);
#endif
}
double MathSpecial::fm_exp(double x)

View File

@ -20,21 +20,61 @@ namespace LAMMPS_NS {
namespace MathSpecial {
// tabulated factorial function
/*! Fast tabulated factorial function
*
* This function looks up pre-computed factorial values for arguments of n = 0
* to a maximum of 167, which is the maximal value representable by a double
* precision floating point number. For other values of n a NaN value is returned.
*
* \param n argument (valid: 0 <= n <= 167)
* \return value of n! as double precision number or NaN */
extern double factorial(const int);
extern double factorial(const int n);
/*! Fast implementation of 2^x without argument checks for little endian CPUs
*
* This function implements an optimized version of pow(2.0, x) that does not
* check for valid arguments and thus may only be used where arguments are well
* behaved. The implementation makes assumptions about the layout of double
* precision floating point numbers in memory and thus will only work on little
* endian CPUs. If little endian cannot be safely detected, the result of
* calling pow(2.0, x) will be returned. This function also is the basis for
* the fast exponential fm_exp(x).
*
* \param x argument
* \return value of 2^x as double precision number */
extern double exp2_x86(double x);
/*! Fast implementation of exp(x) for little endian CPUs
*
* This function implements an optimized version of exp(x) for little endian CPUs.
* It calls the exp2_x86(x) function with a suitable prefactor to x to return exp(x).
* The implementation makes assumptions about the layout of double
* precision floating point numbers in memory and thus will only work on little
* endian CPUs. If little endian cannot be safely detected, the result of
* calling the exp(x) implementation in the standard math library will be returned.
*
* \param x argument
* \return value of e^x as double precision number */
extern double fm_exp(double x);
// support function for scaled error function complement
extern double erfcx_y100(const double y100);
// fast 2**x function without argument checks for little endian CPUs
extern double exp2_x86(double x);
// fast e**x function for little endian CPUs, falls back to libc on other platforms
extern double fm_exp(double x);
// scaled error function complement exp(x*x)*erfc(x) for coul/long styles
/*! Fast scaled error function complement exp(x*x)*erfc(x) for coul/long styles
*
* This is a portable fast implementation of exp(x*x)*erfc(x) that can be used
* in coul/long pair styles as a replacement for the polynomial expansion that
* is/was widely used. Unlike the polynomial expansion, that is only accurate
* at the level of single precision floating point it provides full double precision
* accuracy, but at comparable speed (unlike the erfc() implementation shipped
* with GNU standard math library).
*
* \param x argument
* \return value of e^(x*x)*erfc(x) */
static inline double my_erfcx(const double x)
{
@ -44,7 +84,15 @@ namespace MathSpecial {
return 2.0 * exp(x * x) - erfcx_y100(400.0 / (4.0 - x));
}
// exp(-x*x) for coul/long styles
/*! Fast implementation of exp(-x*x) for little endian CPUs for coul/long styles
*
* This function implements an optimized version of exp(-x*x) based on exp2_x86()
* for use with little endian CPUs. If little endian cannot be safely detected,
* the result of calling the exp(-x*x) implementation in the standard math
* library will be returned.
*
* \param x argument
* \return value of e^(-x*x) as double precision number */
static inline double expmsq(double x)
{
@ -57,18 +105,34 @@ namespace MathSpecial {
#endif
}
// x**2, use instead of pow(x,2.0)
/*! Fast inline version of pow(x, 2.0)
*
* \param x argument
* \return x*x */
static inline double square(const double &x) { return x * x; }
// x**3, use instead of pow(x,3.0)
/*! Fast inline version of pow(x, 3.0)
*
* \param x argument
* \return x*x */
static inline double cube(const double &x) { return x * x * x; }
// return -1.0 for odd n, 1.0 for even n, like pow(-1.0,n)
/* Fast inline version of pow(-1.0, n)
*
* \param n argument (integer)
* \return -1 if n is odd, 1.0 if n is even */
static inline double powsign(const int n) { return (n & 1) ? -1.0 : 1.0; }
// optimized version of pow(x,n) with n being integer
// up to 10x faster than pow(x,y)
/* Fast inline version of pow(x,n) for integer n
*
* This is a version of pow(x,n) optimized for n being integer.
* Speedups of up to 10x faster than pow(x,y) have been measured.
*
* \param n argument (integer)
* \return value of x^n */
static inline double powint(const double &x, const int n)
{
@ -84,7 +148,12 @@ namespace MathSpecial {
return (n > 0) ? yy : 1.0 / yy;
}
// optimized version of (sin(x)/x)**n with n being a _positive_ integer
/* Fast inline version of (sin(x)/x)^n as used by PPPM kspace styles
*
* This is an optimized function to compute (sin(x)/x)^n as frequently used by PPPM.
*
* \param n argument (integer). Expected to be positive.
* \return value of (sin(x)/x)^n */
static inline double powsinxx(const double &x, int n)
{

View File

@ -0,0 +1,95 @@
---
lammps_version: 24 Mar 2022
date_generated: Wed Apr 27 17:35:13 2022
epsilon: 6e-12
skip_tests: single
prerequisites: ! |
pair smatb
pre_commands: ! |
variable units index metal
variable newton_pair delete
variable newton_pair index on
post_commands: ! ""
input_file: in.metal
pair_style: smatb
pair_coeff: ! |
1 1 2.88 10.35 4.178 0.210 1.818 4.07293506 4.9883063257983666
1 2 2.725 10.70 2.805 0.0977 1.2275 4.08707719 4.4340500673763259
2 2 2.56 10.55 2.43 0.0894 1.2799 3.62038672 4.4340500673763259
extract: ! ""
natoms: 32
init_vdwl: -60.168389119061764
init_coul: 0
init_stress: ! |2-
3.0379713999046641e+02 2.8481689618050700e+02 2.7487901025801546e+02 1.2210883607368261e+01 -4.6310157292178644e+00 4.6132589248268313e-01
init_forces: ! |2
1 1.7537039354439259e+00 7.0360956531760914e+00 -8.4184005958845576e-01
2 -9.6588965012743466e-01 -4.2968327220274549e+00 9.9999262128352517e-01
3 -2.2460353121901719e+00 1.0242304966158118e+01 -1.8443560888714017e+00
4 9.9448733303045811e-01 7.5230397435187024e-01 -6.0766236587933242e-01
5 3.2183954634823486e+00 1.5733656263362944e+01 4.7039945398761146e-01
6 2.8052970436237228e+00 4.8969292026494609e+00 2.9489788602236294e-01
7 -5.7345085273481313e+00 8.6530112130150911e+00 7.7447706601282529e-01
8 -2.7028241379012670e+00 -1.2472144397235256e+00 3.5063865351057539e+00
9 -2.2386799966388922e+00 -5.1178927074297116e+00 -3.6574951500303547e+00
10 -4.0842175514437216e+00 -9.2805347046989155e+00 -9.5216421574104011e+00
11 -3.0158013348581383e+00 -9.8825164063184836e+00 4.7066807600425307e+00
12 2.9355737505075075e-01 -9.1277559819075993e+00 1.7503554943451509e+00
13 6.5442807463873454e-01 3.5771666022763327e-01 -2.9268279206132601e+00
14 2.9546102292027019e+00 -2.1076657421595297e+00 2.4027353745102911e+00
15 -1.4362459620281001e+00 -1.4469338088460937e-01 1.2638040704141933e+00
16 1.1853600123947110e+01 -8.2217359454928722e+00 3.5602820971226925e+00
17 2.9760114434025717e+00 9.1510881878371890e+00 2.6527712778949613e+00
18 -3.8448816734858200e+00 -3.8596834411195530e-01 -1.5363106241864222e+00
19 -5.4287609604186988e+00 8.3639132891808646e-02 1.2597633537120057e+00
20 1.8993344143373079e+00 1.7539397493690578e+00 -2.3674883662900768e+00
21 9.3202877569111724e+00 1.0161906623748195e+01 2.5117245084501487e+00
22 -5.4286006790341602e-02 -4.1235834232008373e-01 -1.6165289278166828e-01
23 5.1618239402933046e+00 7.8084801332958875e+00 -7.8617666396584083e+00
24 2.1573958095505019e+00 3.3064758348776766e+00 1.0356666170852131e+00
25 2.7317865068673508e-01 -6.7983594715200391e-01 -8.4962145425034485e-01
26 -6.1973411118185515e+00 -9.9363067654978252e+00 1.1125063212404200e+01
27 8.8836361676639608e-01 -2.5101152574520498e+00 -1.5693011315761736e+00
28 -1.5415779635281726e+00 3.0301797542504927e-01 1.5248141023284689e+00
29 2.7349816106722007e+00 8.8545204263186583e-01 -8.5526109552538948e-02
30 -1.0720483420459732e+00 -4.4498068532010109e+00 -2.2225391974890889e+00
31 1.1773042862086689e+00 -3.7757672138358065e+00 2.8670718827673392e-01
32 -1.0553662576625207e+01 -9.5490168588044995e+00 -4.0724914608207419e+00
run_vdwl: -60.20669005164827
run_coul: 0
run_stress: ! |2-
3.0366990231859160e+02 2.8472030171897950e+02 2.7473728347285567e+02 1.2172424459197687e+01 -4.6332954461814584e+00 4.7608342060898434e-01
run_forces: ! |2
1 1.7287621761096983e+00 7.0355573735520842e+00 -8.4682174175985714e-01
2 -9.6697895447504500e-01 -4.2958415607798441e+00 1.0015306640771833e+00
3 -2.2449153785930140e+00 1.0220034215661133e+01 -1.8504820220633480e+00
4 9.9601791340276336e-01 7.5369431710262291e-01 -6.0681572641646842e-01
5 3.2060009999484147e+00 1.5701015711784372e+01 4.8372081899056063e-01
6 2.8080670120413069e+00 4.8912991510058559e+00 2.8963381062351951e-01
7 -5.7324770842913582e+00 8.6591358821788482e+00 7.6431887086981398e-01
8 -2.7059540524354420e+00 -1.2467200304482609e+00 3.5070033937113965e+00
9 -2.2338690489351816e+00 -5.1254178402905834e+00 -3.6627325643998292e+00
10 -4.0440304555481070e+00 -9.2663471307481604e+00 -9.4950789175699022e+00
11 -3.0095603324750355e+00 -9.8809272965546722e+00 4.6968379350949254e+00
12 2.9955821111103709e-01 -9.1186689644052663e+00 1.7487225514885398e+00
13 6.4891903429679598e-01 3.5328774712317623e-01 -2.9200027226328809e+00
14 2.9573134095573934e+00 -2.0910918248790211e+00 2.3937660891560109e+00
15 -1.4362727166791407e+00 -1.4487584341495263e-01 1.2659041092021559e+00
16 1.1811093106525520e+01 -8.2096829056411114e+00 3.5467098338575100e+00
17 2.9812217984029723e+00 9.1459474011872697e+00 2.6544874139680337e+00
18 -3.8455925628350474e+00 -3.8293266140015725e-01 -1.5378864958810305e+00
19 -5.4246332382008848e+00 8.2734848602692174e-02 1.2542661008733695e+00
20 1.8971229950794164e+00 1.7555326782428182e+00 -2.3672155825569949e+00
21 9.3028027186473548e+00 1.0133810518512346e+01 2.5079668062515421e+00
22 -5.2108622808623904e-02 -4.0963871892230036e-01 -1.5978127490210359e-01
23 5.1748427747828964e+00 7.8093248331208418e+00 -7.8361064692451841e+00
24 2.1561221337678749e+00 3.3060392424596703e+00 1.0375590296124864e+00
25 2.7439938289817345e-01 -6.7913212538661649e-01 -8.4934088300559951e-01
26 -6.1703634942784520e+00 -9.9078856414327738e+00 1.1113100533969693e+01
27 8.8899668170198953e-01 -2.5084240247927481e+00 -1.5752281408111217e+00
28 -1.5419740545035272e+00 3.0574247513776465e-01 1.5244579333305823e+00
29 2.7318838331391233e+00 8.8556581077732033e-01 -8.1366934099631305e-02
30 -1.0707952569409716e+00 -4.4481845742260617e+00 -2.2189658111567914e+00
31 1.1734335189797629e+00 -3.7753076875705340e+00 2.8648052295520054e-01
32 -1.0557032447392665e+01 -9.5476433755557508e+00 -4.0686411315317814e+00
...

View File

@ -0,0 +1,93 @@
---
lammps_version: 24 Mar 2022
date_generated: Wed Apr 27 17:51:25 2022
epsilon: 6e-12
skip_tests: single
prerequisites: ! |
pair smatb
pre_commands: ! |
variable units index metal
variable newton_pair delete
variable newton_pair index on
post_commands: ! ""
input_file: in.metal
pair_style: smatb
pair_coeff: ! |
* * 2.88 10.35 4.178 0.210 1.818 4.07293506 4.9883063257983666
extract: ! ""
natoms: 32
init_vdwl: 15.546009002462199
init_coul: 0
init_stress: ! |2-
7.0981045919895291e+02 6.9464048009855617e+02 6.8185854046738905e+02 1.8940099448628423e+01 -2.3190641785495831e+00 -7.6126579722972734e+00
init_forces: ! |2
1 1.8383788817685609e+00 7.3043414500016803e+00 1.2928059118046571e+00
2 -2.6848648395616506e+00 -7.8218002510265290e+00 5.1835989139783809e+00
3 -4.1262133070389613e+00 2.0356333849891929e+01 -5.0300134853162675e+00
4 2.1640213680130582e+00 -8.2080827254125310e-01 -5.7189855300172958e+00
5 5.1336080958834502e+00 1.4999943174537506e+01 4.1898025714112999e+00
6 4.5017623545159724e+00 3.6176208409156039e+00 -4.4379897612576352e-01
7 -6.2917785937694219e+00 9.6471529013643753e+00 -1.0169471822201503e+00
8 -6.9079266006711242e+00 -6.2411913330296631e+00 4.0631933657348354e+00
9 1.4448974765969274e+00 -6.3335992210787921e+00 -2.7877666693063965e+00
10 -3.2369455571053383e+00 -1.4726480521388170e+01 -1.2343244613879385e+01
11 -4.4279873426372749e+00 -5.8875813372324703e+00 4.2645243184789097e+00
12 2.7168432086644443e+00 -9.2808529346155240e+00 -4.2510938219210220e-01
13 -2.0372275262801081e+00 2.2026130441803748e+00 -4.7893800758259584e+00
14 6.1358565129828193e+00 -2.4215787035807863e+00 7.4237482521334002e+00
15 -3.2151583983209826e+00 2.8722155568519296e+00 -3.1022228268650602e+00
16 1.1834167029249446e+01 -9.9960409045212177e+00 4.0046914221724927e+00
17 4.4438667082387253e+00 1.4015434663357324e+01 1.1980157677253445e-01
18 -1.5049446792950388e+01 -2.9866974100159291e+00 -7.3487669643132998e+00
19 -1.5962776779061105e+01 -7.7159633670829013e+00 4.0578908190472625e+00
20 2.6368639403789449e+00 2.1948034299735144e+00 -3.2117143706849434e+00
21 1.6938751336137905e+01 5.3405439576559601e+00 3.9965661174944200e-01
22 4.6915262682186748e-01 -8.7314686216914801e+00 -1.2388516494108810e+00
23 4.7426470160406025e+00 1.2927761843820045e+01 -8.4044784783918658e+00
24 1.2732138070191686e+01 9.8600435340252197e+00 8.6676713483852375e+00
25 7.4659209593593703e+00 1.1048037573717906e-01 -5.0287486768790597e+00
26 -5.0656394516203669e+00 -1.5673717710206070e+01 1.3291597641889949e+01
27 4.4816496229054108e+00 -4.8760375214009000e-01 3.1769655817027243e-01
28 1.7501136457128652e+00 5.5544282266261300e+00 7.7052630849006061e+00
29 4.4375276260835363e+00 1.0588320818463746e+01 -5.0996447319488158e-01
30 -4.6722349677858155e+00 -3.1632380030448481e+00 -8.9559874574015321e+00
31 -2.1137799706905263e+00 -1.2576277946084598e+01 7.8183097013203859e+00
32 -2.0076186352052531e+01 -6.7271373781221957e+00 -2.4442712859248417e+00
run_vdwl: 15.467597736353799
run_coul: 0
run_stress: ! |2-
7.0958747804689210e+02 6.9441827563947470e+02 6.8163127749154512e+02 1.8878502566994296e+01 -2.3370443631780007e+00 -7.6034095912045760e+00
run_forces: ! |2
1 1.8060069687912315e+00 7.2937921411443272e+00 1.2874043790398257e+00
2 -2.6860686380065211e+00 -7.8171831729036949e+00 5.1804211817416146e+00
3 -4.1171784798928526e+00 2.0287394801872782e+01 -5.0380109293520263e+00
4 2.1656515555309701e+00 -8.1400307560916696e-01 -5.7145340820348807e+00
5 5.1215865715331210e+00 1.4962842274862354e+01 4.1990209437793720e+00
6 4.5053950543110819e+00 3.6131448405817745e+00 -4.5093219067517526e-01
7 -6.2879800842830242e+00 9.6521438182158885e+00 -1.0299290710899491e+00
8 -6.9178598166543104e+00 -6.2386160865726810e+00 4.0623272209295660e+00
9 1.4558477052948295e+00 -6.3449931159549564e+00 -2.7980116422532566e+00
10 -3.1947087252955644e+00 -1.4694132852113814e+01 -1.2296068939837111e+01
11 -4.4213006161890327e+00 -5.8917488565993494e+00 4.2536421369191544e+00
12 2.7297061846000523e+00 -9.2628857408236485e+00 -4.2352406239391233e-01
13 -2.0463406900814625e+00 2.1917055360874724e+00 -4.7773710542900174e+00
14 6.1383962160687062e+00 -2.3802746888087309e+00 7.3930657535424977e+00
15 -3.2106271078248252e+00 2.8691911150262972e+00 -3.0935315047336651e+00
16 1.1775118023237834e+01 -9.9744044265688689e+00 3.9903927678211364e+00
17 4.4470872747424268e+00 1.3995945810173936e+01 1.2365235440776301e-01
18 -1.5026673304996480e+01 -2.9750422082331016e+00 -7.3424362144564483e+00
19 -1.5930045931183475e+01 -7.7015474120118768e+00 4.0406518181177002e+00
20 2.6298049917805058e+00 2.2000986429230047e+00 -3.2110928815587529e+00
21 1.6909713615092866e+01 5.3101658559139251e+00 4.0868602173987068e-01
22 4.7983382632912780e-01 -8.7146552778884860e+00 -1.2319261458206390e+00
23 4.7543321502399483e+00 1.2910146662851659e+01 -8.3714393526052007e+00
24 1.2703010580857182e+01 9.8446813695739301e+00 8.6611807480415717e+00
25 7.4606294163504421e+00 1.1676144511890918e-01 -5.0227835859998491e+00
26 -5.0389376565333173e+00 -1.5622579837196479e+01 1.3265181089097052e+01
27 4.4802107201737105e+00 -4.8768004812534205e-01 3.0515995799100715e-01
28 1.7426613249450380e+00 5.5559456274721297e+00 7.6944148346751184e+00
29 4.4190985096740434e+00 1.0572340858974512e+01 -4.9630774662960464e-01
30 -4.6643835948427661e+00 -3.1602285684767137e+00 -8.9347962499921074e+00
31 -2.1238163435419626e+00 -1.2568441864236553e+01 7.8082109733945497e+00
32 -2.0058169700227520e+01 -6.7278835686694558e+00 -2.4407165275151792e+00
...

View File

@ -0,0 +1,93 @@
---
lammps_version: 24 Mar 2022
date_generated: Wed Apr 27 17:43:01 2022
epsilon: 6e-12
skip_tests: single
prerequisites: ! |
pair smatb/single
pre_commands: ! |
variable units index metal
variable newton_pair delete
variable newton_pair index on
post_commands: ! ""
input_file: in.metal
pair_style: smatb/single
pair_coeff: ! |
* * 2.88 10.35 4.178 0.210 1.818 4.07293506 4.9883063257983666
extract: ! ""
natoms: 32
init_vdwl: 15.546009002462199
init_coul: 0
init_stress: ! |2-
7.0981045919895291e+02 6.9464048009855617e+02 6.8185854046738905e+02 1.8940099448628423e+01 -2.3190641785495831e+00 -7.6126579722972734e+00
init_forces: ! |2
1 1.8383788817685609e+00 7.3043414500016803e+00 1.2928059118046571e+00
2 -2.6848648395616506e+00 -7.8218002510265290e+00 5.1835989139783809e+00
3 -4.1262133070389613e+00 2.0356333849891929e+01 -5.0300134853162675e+00
4 2.1640213680130582e+00 -8.2080827254125310e-01 -5.7189855300172958e+00
5 5.1336080958834502e+00 1.4999943174537506e+01 4.1898025714112999e+00
6 4.5017623545159724e+00 3.6176208409156039e+00 -4.4379897612576352e-01
7 -6.2917785937694219e+00 9.6471529013643753e+00 -1.0169471822201503e+00
8 -6.9079266006711242e+00 -6.2411913330296631e+00 4.0631933657348354e+00
9 1.4448974765969274e+00 -6.3335992210787921e+00 -2.7877666693063965e+00
10 -3.2369455571053383e+00 -1.4726480521388170e+01 -1.2343244613879385e+01
11 -4.4279873426372749e+00 -5.8875813372324703e+00 4.2645243184789097e+00
12 2.7168432086644443e+00 -9.2808529346155240e+00 -4.2510938219210220e-01
13 -2.0372275262801081e+00 2.2026130441803744e+00 -4.7893800758259584e+00
14 6.1358565129828193e+00 -2.4215787035807863e+00 7.4237482521334002e+00
15 -3.2151583983209826e+00 2.8722155568519296e+00 -3.1022228268650602e+00
16 1.1834167029249446e+01 -9.9960409045212177e+00 4.0046914221724927e+00
17 4.4438667082387253e+00 1.4015434663357324e+01 1.1980157677253445e-01
18 -1.5049446792950388e+01 -2.9866974100159291e+00 -7.3487669643132998e+00
19 -1.5962776779061105e+01 -7.7159633670829013e+00 4.0578908190472625e+00
20 2.6368639403789449e+00 2.1948034299735144e+00 -3.2117143706849434e+00
21 1.6938751336137905e+01 5.3405439576559601e+00 3.9965661174944195e-01
22 4.6915262682186748e-01 -8.7314686216914801e+00 -1.2388516494108810e+00
23 4.7426470160406025e+00 1.2927761843820045e+01 -8.4044784783918658e+00
24 1.2732138070191686e+01 9.8600435340252197e+00 8.6676713483852375e+00
25 7.4659209593593703e+00 1.1048037573717906e-01 -5.0287486768790597e+00
26 -5.0656394516203669e+00 -1.5673717710206070e+01 1.3291597641889949e+01
27 4.4816496229054108e+00 -4.8760375214009000e-01 3.1769655817027243e-01
28 1.7501136457128652e+00 5.5544282266261300e+00 7.7052630849006061e+00
29 4.4375276260835363e+00 1.0588320818463746e+01 -5.0996447319488158e-01
30 -4.6722349677858155e+00 -3.1632380030448481e+00 -8.9559874574015321e+00
31 -2.1137799706905263e+00 -1.2576277946084598e+01 7.8183097013203859e+00
32 -2.0076186352052531e+01 -6.7271373781221957e+00 -2.4442712859248417e+00
run_vdwl: 15.467597736353799
run_coul: 0
run_stress: ! |2-
7.0958747804689210e+02 6.9441827563947470e+02 6.8163127749154512e+02 1.8878502566994296e+01 -2.3370443631780007e+00 -7.6034095912045760e+00
run_forces: ! |2
1 1.8060069687912315e+00 7.2937921411443272e+00 1.2874043790398257e+00
2 -2.6860686380065211e+00 -7.8171831729036949e+00 5.1804211817416146e+00
3 -4.1171784798928526e+00 2.0287394801872782e+01 -5.0380109293520263e+00
4 2.1656515555309701e+00 -8.1400307560916696e-01 -5.7145340820348807e+00
5 5.1215865715331210e+00 1.4962842274862354e+01 4.1990209437793720e+00
6 4.5053950543110819e+00 3.6131448405817745e+00 -4.5093219067517526e-01
7 -6.2879800842830242e+00 9.6521438182158885e+00 -1.0299290710899491e+00
8 -6.9178598166543104e+00 -6.2386160865726810e+00 4.0623272209295660e+00
9 1.4558477052948295e+00 -6.3449931159549564e+00 -2.7980116422532566e+00
10 -3.1947087252955644e+00 -1.4694132852113814e+01 -1.2296068939837111e+01
11 -4.4213006161890327e+00 -5.8917488565993494e+00 4.2536421369191544e+00
12 2.7297061846000523e+00 -9.2628857408236485e+00 -4.2352406239391233e-01
13 -2.0463406900814625e+00 2.1917055360874724e+00 -4.7773710542900174e+00
14 6.1383962160687062e+00 -2.3802746888087309e+00 7.3930657535424977e+00
15 -3.2106271078248252e+00 2.8691911150262972e+00 -3.0935315047336651e+00
16 1.1775118023237834e+01 -9.9744044265688672e+00 3.9903927678211364e+00
17 4.4470872747424268e+00 1.3995945810173936e+01 1.2365235440776301e-01
18 -1.5026673304996480e+01 -2.9750422082331016e+00 -7.3424362144564483e+00
19 -1.5930045931183475e+01 -7.7015474120118768e+00 4.0406518181177002e+00
20 2.6298049917805058e+00 2.2000986429230047e+00 -3.2110928815587529e+00
21 1.6909713615092866e+01 5.3101658559139251e+00 4.0868602173987068e-01
22 4.7983382632912785e-01 -8.7146552778884860e+00 -1.2319261458206390e+00
23 4.7543321502399483e+00 1.2910146662851659e+01 -8.3714393526052007e+00
24 1.2703010580857182e+01 9.8446813695739301e+00 8.6611807480415717e+00
25 7.4606294163504421e+00 1.1676144511890918e-01 -5.0227835859998500e+00
26 -5.0389376565333173e+00 -1.5622579837196479e+01 1.3265181089097052e+01
27 4.4802107201737105e+00 -4.8768004812534205e-01 3.0515995799100759e-01
28 1.7426613249450380e+00 5.5559456274721297e+00 7.6944148346751184e+00
29 4.4190985096740434e+00 1.0572340858974512e+01 -4.9630774662960464e-01
30 -4.6643835948427661e+00 -3.1602285684767137e+00 -8.9347962499921074e+00
31 -2.1238163435419626e+00 -1.2568441864236553e+01 7.8082109733945497e+00
32 -2.0058169700227520e+01 -6.7278835686694558e+00 -2.4407165275151792e+00
...