ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into dhairyaiitb_develop

Browse Source
This commit is contained in:
Axel Kohlmeyer
2025-03-04 21:29:40 -05:00
parent 7952281f78 16db8e1515
commit 885099cf77
33 changed files with 418 additions and 306 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.github/CODEOWNERS vendored
View File

@ -72,6 +72,8 @@ src/EXTRA-COMMAND/ndx_group.* @akohlmey
src/EXTRA-COMPUTE/compute_stress_mop*.* @RomainVermorel src/EXTRA-COMPUTE/compute_stress_mop*.* @RomainVermorel
src/EXTRA-COMPUTE/compute_born_matrix.* @Bibobu @athomps src/EXTRA-COMPUTE/compute_born_matrix.* @Bibobu @athomps
src/EXTRA-FIX/fix_deform_pressure.* @jtclemm src/EXTRA-FIX/fix_deform_pressure.* @jtclemm
src/EXTRA-PAIR/pair_dispersion_d3.* @soniasolomoni @arthurfl
src/EXTRA-PAIR/d3_parameters.h @soniasolomoni @arthurfl
src/MISC/*_tracker.* @jtclemm src/MISC/*_tracker.* @jtclemm
src/MC/fix_gcmc.* @athomps src/MC/fix_gcmc.* @athomps
src/MC/fix_sgcmc.* @athomps src/MC/fix_sgcmc.* @athomps

2
doc/src/compute_reduce.rst
View File

@ -87,7 +87,7 @@ values in the vector. The *sumsq* option sums the square of the
values in the vector into a global total. The *avesq* setting does values in the vector into a global total. The *avesq* setting does
the same as *sumsq*, then divides the sum of squares by the number of the same as *sumsq*, then divides the sum of squares by the number of
values. The last two options can be useful for calculating the values. The last two options can be useful for calculating the
variance of some quantity (e.g., variance = sumsq :math:`-` ave\ variance of some quantity (e.g., variance = *avesq* :math:`-` *ave*\
:math:`^2`). The *sumabs* option sums the absolute values in the :math:`^2`). The *sumabs* option sums the absolute values in the
vector into a global total. The *aveabs* setting does the same as vector into a global total. The *aveabs* setting does the same as
*sumabs*, then divides the sum of absolute values by the number of *sumabs*, then divides the sum of absolute values by the number of

2
doc/src/fix_press_langevin.rst
View File

@ -54,7 +54,7 @@ the Langevin equation such as:
f_P = & \frac{N k_B T_{target}}{V} + \frac{1}{V d}\sum_{i=1}^{N} \vec r_i \cdot \vec f_i - P_{target} \\ f_P = & \frac{N k_B T_{target}}{V} + \frac{1}{V d}\sum_{i=1}^{N} \vec r_i \cdot \vec f_i - P_{target} \\
Q\ddot{L} + \alpha{}\dot{L} = & f_P + \beta(t)\\ Q\ddot{L} + \alpha{}\dot{L} = & f_P + \beta(t)\\
L^{n+1} = & L^{n} + bdt\dot{L}^{n} \frac{bdt^{2}}{2Q} \\ L^{n+1} = & L^{n} + bdt\dot{L}^{n} + \frac{bdt^{2}}{2Q} f^{n}_{P} + \frac{bdt}{2Q} \beta^{n+1} \\
\dot{L}^{n+1} = & \alpha\dot{L}^{n} + \frac{dt}{2Q}\left(a f^{n}_{P} + f^{n+1}_{P}\right) + \frac{b}{Q}\beta^{n+1} \\ \dot{L}^{n+1} = & \alpha\dot{L}^{n} + \frac{dt}{2Q}\left(a f^{n}_{P} + f^{n+1}_{P}\right) + \frac{b}{Q}\beta^{n+1} \\
a = & \frac{1-\frac{\alpha{}dt}{2Q}}{1+\frac{\alpha{}dt}{2Q}} \\ a = & \frac{1-\frac{\alpha{}dt}{2Q}}{1+\frac{\alpha{}dt}{2Q}} \\
b = & \frac{1}{1+\frac{\alpha{}dt}{2Q}} \\ b = & \frac{1}{1+\frac{\alpha{}dt}{2Q}} \\

21
doc/src/pair_dispersion_d3.rst
View File

@ -10,7 +10,7 @@ Syntax
pair_style dispersion/d3 damping functional cutoff cn_cutoff pair_style dispersion/d3 damping functional cutoff cn_cutoff
* damping = damping function: *zero*, *zerom*, *bj*, or *bjm* * damping = damping function: *original*, *zerom*, *bj*, or *bjm*
* functional = XC functional form: *pbe*, *pbe0*, ... (see list below) * functional = XC functional form: *pbe*, *pbe0*, ... (see list below)
* cutoff = global cutoff (distance units) * cutoff = global cutoff (distance units)
* cn_cutoff = coordination number cutoff (distance units) * cn_cutoff = coordination number cutoff (distance units)
@ -20,7 +20,7 @@ Examples
.. code-block:: LAMMPS .. code-block:: LAMMPS
pair_style dispersion/d3 zero pbe 30.0 20.0 pair_style dispersion/d3 original pbe 30.0 20.0
pair_coeff * * C pair_coeff * * C
Description Description
@ -55,9 +55,16 @@ factor, and :math:`f_n^{damp}` are damping functions.
contributions, according to, for example, the Axilrod-Teller-Muto contributions, according to, for example, the Axilrod-Teller-Muto
model. model.
Available damping functions are the original "zero-damping" .. versionchanged:: TBD
:ref:`(Grimme1) <Grimme1>`, Becke-Johnson damping :ref:`(Grimme2)
<Grimme2>`, and their revised forms :ref:`(Sherrill) <Sherrill>`. renamed *zero* keyword to *original* to avoid conflicts with
:doc:`pair style zero <pair_zero>` when used as :doc:`hybrid
sub-style <pair_hybrid>`.
Available damping functions are the original "zero-damping" (*original*)
:ref:`(Grimme1) <Grimme1>`, Becke-Johnson damping (*bj*) :ref:`(Grimme2)
<Grimme2>`, and their revised forms (*zerom* and *bjm*, respectively)
:ref:`(Sherrill) <Sherrill>`.
Available XC functional scaling factors are listed in the table below, Available XC functional scaling factors are listed in the table below,
and depend on the selected damping function. and depend on the selected damping function.
@ -67,7 +74,7 @@ and depend on the selected damping function.
+==================+================================================================================+ +==================+================================================================================+
| | | | slater-dirac-exchange, b-lyp, b-p, b97-d, revpbe, pbe, pbesol, rpw86-pbe, | | | | | slater-dirac-exchange, b-lyp, b-p, b97-d, revpbe, pbe, pbesol, rpw86-pbe, |
| | | | rpbe, tpss, b3-lyp, pbe0, hse06, revpbe38, pw6b95, tpss0, b2-plyp, pwpb95, | | | | | rpbe, tpss, b3-lyp, pbe0, hse06, revpbe38, pw6b95, tpss0, b2-plyp, pwpb95, |
| | zero | | b2gp-plyp, ptpss, hf, mpwlyp, bpbe, bh-lyp, tpssh, pwb6k, b1b95, bop, o-lyp, | | | original | | b2gp-plyp, ptpss, hf, mpwlyp, bpbe, bh-lyp, tpssh, pwb6k, b1b95, bop, o-lyp, |
| | | | o-pbe, ssb, revssb, otpss, b3pw91, revpbe0, pbe38, mpw1b95, mpwb1k, bmk, | | | | | o-pbe, ssb, revssb, otpss, b3pw91, revpbe0, pbe38, mpw1b95, mpwb1k, bmk, |
| | | | cam-b3lyp, lc-wpbe, m05, m052x, m06l, m06, m062x, m06hf, hcth120 | | | | | cam-b3lyp, lc-wpbe, m05, m052x, m06l, m06, m062x, m06hf, hcth120 |
+------------------+--------------------------------------------------------------------------------+ +------------------+--------------------------------------------------------------------------------+
@ -122,6 +129,8 @@ Style *dispersion/d3* is part of the EXTRA-PAIR package. It is only
enabled if LAMMPS was built with that package. See the :doc:`Build enabled if LAMMPS was built with that package. See the :doc:`Build
package <Build_package>` page for more info. package <Build_package>` page for more info.
The compiled in parameters require the use of :doc:`metal units <units>`.
It is currently *not* possible to calculate three-body dispersion It is currently *not* possible to calculate three-body dispersion
contributions according to, for example, the Axilrod-Teller-Muto model. contributions according to, for example, the Axilrod-Teller-Muto model.

25
doc/src/region.rst
View File

@ -40,7 +40,7 @@ Syntax
*plane* args = px py pz nx ny nz *plane* args = px py pz nx ny nz
px,py,pz = point on the plane (distance units) px,py,pz = point on the plane (distance units)
nx,ny,nz = direction normal to plane (distance units) nx,ny,nz = direction normal to plane (distance units)
px,py,pz can be a variable (see below) px,py,pz,nx,ny,nz can be a variable (see below)
*prism* args = xlo xhi ylo yhi zlo zhi xy xz yz *prism* args = xlo xhi ylo yhi zlo zhi xy xz yz
xlo,xhi,ylo,yhi,zlo,zhi = bounds of untilted prism (distance units) xlo,xhi,ylo,yhi,zlo,zhi = bounds of untilted prism (distance units)
xy = distance to tilt y in x direction (distance units) xy = distance to tilt y in x direction (distance units)
@ -211,9 +211,11 @@ equal-style :doc:`variable <variable>`. Likewise, for style *sphere*
and *ellipsoid* the x-, y-, and z- coordinates of the center of the and *ellipsoid* the x-, y-, and z- coordinates of the center of the
sphere/ellipsoid can be specified as an equal-style variable. And for sphere/ellipsoid can be specified as an equal-style variable. And for
style *cylinder* the two center positions c1 and c2 for the location of style *cylinder* the two center positions c1 and c2 for the location of
the cylinder axes can be specified as a equal-style variable. For style the cylinder axes can be specified as a equal-style variable. For styles
*cone* and *prism* all properties can be defined via equal-style variables. For *block*, *cone*, *prism*, and *plane* all properties can be defined via
style *plane* the point can be defined via equal-style variables. equal-style variables. For style *plane*, the components of the direction
vector normal to plane should be either all constants or all defined by
equal-style variables.
If the value is a variable, it should be specified as v_name, where If the value is a variable, it should be specified as v_name, where
name is the variable name. In this case, the variable will be name is the variable name. In this case, the variable will be
@ -226,6 +228,21 @@ keywords for the simulation box parameters and timestep and elapsed
time. Thus it is easy to specify a time-dependent radius or have time. Thus it is easy to specify a time-dependent radius or have
a time dependent position of the sphere or cylinder region. a time dependent position of the sphere or cylinder region.
.. note::
Whenever a region property, such as a coordinate or an upper/lower
bound, is defined via an equal-style variable, the variable should
not cause any of the region boundaries to move
too far within a single timestep. Otherwise, bad dynamics will occur.
"Too far" means a small fraction of the approximate distance of
closest approach between two particles, which for the case of Lennard-Jones
particles is the distance of the energy minimum while for granular
particles it is their diameter. An example is a rapidly varying direction
vector in region plane since a small change in the normal to plane will
shift the region surface far away from the region point by a large displacement.
Similarly, bad dynamics can also occur for fast changing variables employed
in the move/rotate options.
See the :doc:`Howto tricilinc <Howto_triclinic>` page for a See the :doc:`Howto tricilinc <Howto_triclinic>` page for a
geometric description of triclinic boxes, as defined by LAMMPS, and geometric description of triclinic boxes, as defined by LAMMPS, and
how to transform these parameters to and from other commonly used how to transform these parameters to and from other commonly used

67
examples/ASPHERE/flat_membrane/log.18Oct22.flat_membrane.g++.1 → examples/ASPHERE/flat_membrane/log.3Mar25.flat_membrane.g++.1
View File

@ -1,4 +1,5 @@
LAMMPS (15 Sep 2022) LAMMPS (4 Feb 2025 - Development - patch_4Feb2025-106-g5d02e140d4)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99)
using 1 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# flat membrane demo # flat membrane demo
variable r0 equal 0.97 variable r0 equal 0.97
@ -46,7 +47,7 @@ set type 1 mass 1.0
Setting atom values ... Setting atom values ...
1920 settings made for mass 1920 settings made for mass
set type 1 shape 1 1 1 set type 1 shape 1 0.99 0.99
Setting atom values ... Setting atom values ...
1920 settings made for shape 1920 settings made for shape
@ -99,7 +100,7 @@ CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ... Neighbor list info ...
update: every = 1 steps, delay = 10 steps, check = yes update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000 max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.6 master list distance cutoff = 3.6
ghost atom cutoff = 3.6 ghost atom cutoff = 3.6
@ -113,47 +114,47 @@ Neighbor list info ...
Per MPI rank memory allocation (min/avg/max) = 5.024 | 5.024 | 5.024 Mbytes Per MPI rank memory allocation (min/avg/max) = 5.024 | 5.024 | 5.024 Mbytes
Step Temp Press Pxx Pyy Step Temp Press Pxx Pyy
0 0.23 -0.0073508785 -0.012283389 -0.012234574 0 0.23 -0.0073508785 -0.012283389 -0.012234574
200 0.20903886 -0.0010605951 -0.0011885957 -0.00198842 200 0.20906916 -0.0010610794 -0.0011895359 -0.00198968
400 0.21898026 -0.00069250685 -0.0013217981 -0.00073225707 400 0.21907691 -0.00068750897 -0.0013175111 -0.00072355516
600 0.22689361 -0.00057919328 -0.00076880503 -0.0010242283 600 0.22648739 -0.00057306694 -0.000831971 -0.00096776143
800 0.22983221 -0.00032145682 -0.00051928834 -0.00059337525 800 0.2368957 -0.00028517767 -0.00050638195 -0.00052670236
1000 0.23819392 -0.00027969126 -0.00088082301 -5.2666567e-05 1000 0.22735705 -0.00032638104 -0.00037959812 -0.00071879257
1200 0.22053795 -0.00029571329 -0.0004446455 -0.00035529929 1200 0.22910882 -0.00019457758 -0.00024451315 -0.00027807764
1400 0.22285021 -0.0002690371 -0.00068896571 -3.6258442e-05 1400 0.22754022 -0.00033048317 -0.00010053263 -0.00075173132
1600 0.22687044 2.8599875e-05 -0.00032651798 0.0004056081 1600 0.22503496 -7.593954e-05 5.2989168e-05 -0.00022521685
1800 0.23356905 -2.28742e-05 -0.00027073251 0.00025081131 1800 0.22673577 -0.00017917699 1.0176667e-05 -0.00047893102
2000 0.22499821 8.8230586e-06 -7.5750159e-05 0.0001988705 2000 0.22480858 9.318318e-05 0.00031024343 1.2738253e-05
2200 0.23162995 -9.026855e-05 -0.00025832535 5.4904927e-05 2200 0.22387294 -6.9083174e-05 -4.6170865e-05 4.6682009e-05
2400 0.22920223 0.00016700455 3.5283125e-05 0.00034955857 2400 0.22686231 0.00022776785 0.00028204391 0.00030116273
2600 0.2260299 5.3095557e-05 0.00025691786 0.00013353467 2600 0.23055637 8.2778175e-05 0.00027137529 4.2254076e-05
2800 0.2296401 0.00043234854 0.00058344966 0.00063645193 2800 0.23126717 0.00044863024 0.00071654709 0.00057652233
3000 0.22564577 2.6423111e-05 8.9918406e-05 0.00022146229 3000 0.22234736 0.00031577792 0.00051797227 0.00039087282
Loop time of 6.76659 on 1 procs for 3000 steps with 1920 atoms Loop time of 7.17991 on 1 procs for 3000 steps with 1920 atoms
Performance: 383058.431 tau/day, 443.355 timesteps/s Performance: 361007.348 tau/day, 417.833 timesteps/s, 802.239 katom-step/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads 99.9% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 5.7968 | 5.7968 | 5.7968 | 0.0 | 85.67 Pair | 6.2257 | 6.2257 | 6.2257 | 0.0 | 86.71
Neigh | 0.086077 | 0.086077 | 0.086077 | 0.0 | 1.27 Neigh | 0.11392 | 0.11392 | 0.11392 | 0.0 | 1.59
Comm | 0.034761 | 0.034761 | 0.034761 | 0.0 | 0.51 Comm | 0.036614 | 0.036614 | 0.036614 | 0.0 | 0.51
Output | 0.00038014 | 0.00038014 | 0.00038014 | 0.0 | 0.01 Output | 0.00052672 | 0.00052672 | 0.00052672 | 0.0 | 0.01
Modify | 0.8268 | 0.8268 | 0.8268 | 0.0 | 12.22 Modify | 0.78577 | 0.78577 | 0.78577 | 0.0 | 10.94
Other | | 0.02181 | | | 0.32 Other | | 0.01743 | | | 0.24
Nlocal: 1920 ave 1920 max 1920 min Nlocal: 1920 ave 1920 max 1920 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 772 ave 772 max 772 min Nghost: 771 ave 771 max 771 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 46804 ave 46804 max 46804 min Neighs: 46750 ave 46750 max 46750 min
Histogram: 1 0 0 0 0 0 0 0 0 0 Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 46804 Total # of neighbors = 46750
Ave neighs/atom = 24.377083 Ave neighs/atom = 24.348958
Neighbor list builds = 99 Neighbor list builds = 98
Dangerous builds = 0 Dangerous builds = 0
Total wall time: 0:00:06 Total wall time: 0:00:07

71
examples/ASPHERE/flat_membrane/log.18Oct22.flat_membrane.g++.4 → examples/ASPHERE/flat_membrane/log.3Mar25.flat_membrane.g++.4
View File

@ -1,4 +1,5 @@
LAMMPS (15 Sep 2022) LAMMPS (4 Feb 2025 - Development - patch_4Feb2025-106-g5d02e140d4)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:99)
using 1 OpenMP thread(s) per MPI task using 1 OpenMP thread(s) per MPI task
# flat membrane demo # flat membrane demo
variable r0 equal 0.97 variable r0 equal 0.97
@ -46,7 +47,7 @@ set type 1 mass 1.0
Setting atom values ... Setting atom values ...
1920 settings made for mass 1920 settings made for mass
set type 1 shape 1 1 1 set type 1 shape 1 0.99 0.99
Setting atom values ... Setting atom values ...
1920 settings made for shape 1920 settings made for shape
@ -99,7 +100,7 @@ CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ... Neighbor list info ...
update: every = 1 steps, delay = 10 steps, check = yes update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000 max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.6 master list distance cutoff = 3.6
ghost atom cutoff = 3.6 ghost atom cutoff = 3.6
@ -113,47 +114,47 @@ Neighbor list info ...
Per MPI rank memory allocation (min/avg/max) = 4.182 | 4.794 | 5.472 Mbytes Per MPI rank memory allocation (min/avg/max) = 4.182 | 4.794 | 5.472 Mbytes
Step Temp Press Pxx Pyy Step Temp Press Pxx Pyy
0 0.23 -0.0073508785 -0.012283389 -0.012234574 0 0.23 -0.0073508785 -0.012283389 -0.012234574
200 0.20647718 -0.0012368074 -0.0021167303 -0.0015343502 200 0.21866566 -0.0011213114 -0.0017540296 -0.0017008786
400 0.21648371 -0.00085695085 -0.0015627331 -0.0011177093 400 0.22409469 -0.00055137604 -0.00050366017 -0.0012846672
600 0.22929515 -0.00050218657 -0.0008332 -0.00062622609 600 0.21623646 -0.0004808564 -0.00063496522 -0.00071617302
800 0.22062664 -0.00049172378 -0.000611884 -0.00075089294 800 0.22248339 -0.00066333806 -0.00093812145 -0.00099410261
1000 0.22422425 -0.00039405068 -0.00037600355 -0.00070786572 1000 0.22367907 -0.0003439459 -0.00076259578 -0.00013354375
1200 0.2298767 -0.00025939082 -0.00021616578 -0.00053125505 1200 0.23276206 8.5764334e-06 -0.00036165178 0.00038543803
1400 0.2335927 5.8028332e-05 0.00017530192 -3.1675138e-05 1400 0.23129049 -5.1997966e-06 -0.00040514787 0.00017093646
1600 0.22884878 -0.0001733902 -0.0008056431 0.00014276754 1600 0.22074564 -0.0001604642 -0.00053959031 9.0463391e-05
1800 0.22813498 0.00019873459 0.00051040124 5.8860949e-05 1800 0.2276062 -0.00015655254 -0.00042925633 5.4938379e-05
2000 0.2273166 -3.3595127e-05 0.0001705632 -0.00026498213 2000 0.22469391 -3.6198836e-05 -1.8538296e-06 -0.00010464749
2200 0.2251643 -2.4517311e-05 -4.0618888e-05 1.066658e-05 2200 0.22737515 -1.8967356e-05 0.0001341006 -0.00022739894
2400 0.22460629 -4.5661259e-05 -0.00019144039 -1.6649099e-05 2400 0.22607533 0.00014650382 -2.4871789e-05 0.00058296255
2600 0.23085675 0.00014029405 0.00017983536 0.00017895001 2600 0.24044743 0.00045066449 0.00049153595 0.00078042437
2800 0.22364591 4.2999164e-05 -0.00011000466 0.00024363243 2800 0.23346628 0.00017994019 -5.0897724e-05 0.00046011801
3000 0.23421357 0.00023505702 0.00020752013 0.00053567111 3000 0.22427206 0.00040437022 0.00049580531 0.00046259202
Loop time of 4.68577 on 4 procs for 3000 steps with 1920 atoms Loop time of 5.04539 on 4 procs for 3000 steps with 1920 atoms
Performance: 553164.568 tau/day, 640.237 timesteps/s Performance: 513736.789 tau/day, 594.603 timesteps/s, 1.142 Matom-step/s
95.6% CPU use with 4 MPI tasks x 1 OpenMP threads 99.0% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown: MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total Section | min time | avg time | max time |%varavg| %total
--------------------------------------------------------------- ---------------------------------------------------------------
Pair | 0.00015072 | 1.6029 | 3.8573 | 131.7 | 34.21 Pair | 4.1163 | 4.2634 | 4.4004 | 6.5 | 84.50
Neigh | 0.00055747 | 0.025423 | 0.065858 | 17.0 | 0.54 Neigh | 0.00088605 | 0.032793 | 0.083516 | 19.1 | 0.65
Comm | 0.0052259 | 0.48173 | 1.624 | 96.5 | 10.28 Comm | 0.006704 | 0.20003 | 0.40672 | 42.5 | 3.96
Output | 0.0003894 | 0.023428 | 0.047223 | 15.0 | 0.50 Output | 0.00033812 | 0.0012192 | 0.0020905 | 2.3 | 0.02
Modify | 0.00037337 | 0.2141 | 0.44595 | 46.3 | 4.57 Modify | 0.00035967 | 0.20167 | 0.42535 | 45.0 | 4.00
Other | | 2.338 | | | 49.90 Other | | 0.3462 | | | 6.86
Nlocal: 480 ave 1011 max 0 min Nlocal: 480 ave 1052 max 0 min
Histogram: 2 0 0 0 0 0 0 0 1 1 Histogram: 2 0 0 0 0 0 0 0 1 1
Nghost: 860 ave 1771 max 0 min Nghost: 854.5 ave 1801 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2 Histogram: 2 0 0 0 0 0 0 0 1 1
Neighs: 11697.8 ave 30095 max 0 min Neighs: 11684.2 ave 31145 max 0 min
Histogram: 2 0 0 0 0 1 0 0 0 1 Histogram: 2 0 0 0 0 1 0 0 0 1
Total # of neighbors = 46791 Total # of neighbors = 46737
Ave neighs/atom = 24.370313 Ave neighs/atom = 24.342188
Neighbor list builds = 99 Neighbor list builds = 98
Dangerous builds = 0 Dangerous builds = 0
Total wall time: 0:00:04 Total wall time: 0:00:05

4
src/ASPHERE/pair_ylz.cpp
View File

@ -32,7 +32,7 @@
#include <cmath> #include <cmath>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using MathConst::MY_4PI; using MathConst::MY_PI;
using MathConst::MY_PI2; using MathConst::MY_PI2;
using MathConst::MY_TWOBYSIXTH; using MathConst::MY_TWOBYSIXTH;
@ -487,7 +487,7 @@ double PairYLZ::ylz_analytic(const int i, const int j, double a1[3][3], double a
uA = -energy_well * t1 * cos_t; uA = -energy_well * t1 * cos_t;
U = uA * phi; U = uA * phi;
dUdr = MY_4PI / (rcut - rmin) * (t1) *sin(t) * phi * energy_well; dUdr = MY_PI * zt / (rcut - rmin) * (t1) *sin(t) * phi * energy_well;
dUdphi = uA; dUdphi = uA;
} }

8
src/DPD-BASIC/pair_dpd.cpp
View File

@ -74,6 +74,10 @@ void PairDPD::compute(int eflag, int vflag)
evdwl = 0.0; evdwl = 0.0;
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -266,10 +270,6 @@ void PairDPD::init_style()
error->warning(FLERR, "Pair dpd needs newton pair on for momentum conservation"); error->warning(FLERR, "Pair dpd needs newton pair on for momentum conservation");
neighbor->add_request(this); neighbor->add_request(this);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

9
src/DPD-BASIC/pair_dpd_coul_slater_long.cpp
View File

@ -90,6 +90,10 @@ void PairDPDCoulSlaterLong::compute(int eflag, int vflag)
evdwl = ecoul = 0.0; evdwl = ecoul = 0.0;
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -334,11 +338,6 @@ void PairDPDCoulSlaterLong::init_style()
neighbor->add_request(this); neighbor->add_request(this);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
// ensure use of KSpace long-range solver, set g_ewald // ensure use of KSpace long-range solver, set g_ewald
if (force->kspace == nullptr) if (force->kspace == nullptr)

8
src/DPD-BASIC/pair_dpd_ext.cpp
View File

@ -81,6 +81,10 @@ void PairDPDExt::compute(int eflag, int vflag)
evdwl = 0.0; evdwl = 0.0;
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -325,10 +329,6 @@ void PairDPDExt::init_style()
error->warning(FLERR, "Pair dpd needs newton pair on for momentum conservation"); error->warning(FLERR, "Pair dpd needs newton pair on for momentum conservation");
neighbor->add_request(this); neighbor->add_request(this);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

4
src/DPD-BASIC/pair_dpd_ext_tstat.cpp
View File

@ -54,6 +54,10 @@ void PairDPDExtTstat::compute(int eflag, int vflag)
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
// adjust sigma if target T is changing // adjust sigma if target T is changing
if (t_start != t_stop) { if (t_start != t_stop) {

4
src/DPD-BASIC/pair_dpd_tstat.cpp
View File

@ -48,6 +48,10 @@ void PairDPDTstat::compute(int eflag, int vflag)
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
// adjust sigma if target T is changing // adjust sigma if target T is changing
if (t_start != t_stop) { if (t_start != t_stop) {

2
src/EXTRA-MOLECULE/angle_cosine_delta.cpp
View File

@ -89,8 +89,6 @@ void AngleCosineDelta::compute(int eflag, int vflag)
s = sqrt(1.0 - c*c); s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL; if (s < SMALL) s = SMALL;
s = 1.0/s;
cot = c/s; cot = c/s;
// force & energy // force & energy

84
src/EXTRA-PAIR/pair_dispersion_d3.cpp
View File

@ -29,6 +29,7 @@
#include "memory.h" #include "memory.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neighbor.h" #include "neighbor.h"
#include "update.h"
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
@ -58,6 +59,7 @@ static constexpr double autoang = 0.52917725; // atomic units (Bohr) to Angst
static constexpr double autoev = 27.21140795; // atomic units (Hartree) to eV static constexpr double autoev = 27.21140795; // atomic units (Hartree) to eV
#include "d3_parameters.h" #include "d3_parameters.h"
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
Constructor (Required) Constructor (Required)
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -75,6 +77,7 @@ PairDispersionD3::PairDispersionD3(LAMMPS *lmp) :
one_coeff = 1; one_coeff = 1;
single_enable = 0; single_enable = 0;
dampingCode = 0;
s6 = s8 = s18 = rs6 = rs8 = rs18 = a1 = a2 = alpha = alpha6 = alpha8 = 0.0; s6 = s8 = s18 = rs6 = rs8 = rs18 = a1 = a2 = alpha = alpha6 = alpha8 = 0.0;
} }
@ -130,12 +133,20 @@ void PairDispersionD3::allocate()
void PairDispersionD3::settings(int narg, char **arg) void PairDispersionD3::settings(int narg, char **arg)
{ {
if (narg != 4) error->all(FLERR, "Pair_style dispersion/d3 needs 4 arguments"); if (narg != 4) error->all(FLERR, "Pair style dispersion/d3 needs 4 arguments");
if (strcmp("metal", update->unit_style) != 0)
error->all(FLERR, Error::NOLASTLINE, "Pair style dispersion/d3 requires metal units");
damping_type = arg[0]; std::string damping_type = arg[0];
std::string functional_name = arg[1]; std::string functional_name = arg[1];
std::transform(damping_type.begin(), damping_type.end(), damping_type.begin(), ::tolower); std::transform(damping_type.begin(), damping_type.end(), damping_type.begin(), ::tolower);
std::unordered_map<std::string, int> dampingMap = {
{"original", 1}, {"zero", 1}, {"zerom", 2}, {"bj", 3}, {"bjm", 4}};
if (!dampingMap.count(damping_type))
error->all(FLERR, Error::NOPOINTER, "Unknown damping type {} for pair style dispersion/d3",
damping_type);
dampingCode = dampingMap[damping_type];
rthr = utils::numeric(FLERR, arg[2], false, lmp); rthr = utils::numeric(FLERR, arg[2], false, lmp);
cn_thr = utils::numeric(FLERR, arg[3], false, lmp); cn_thr = utils::numeric(FLERR, arg[3], false, lmp);
@ -154,7 +165,7 @@ int PairDispersionD3::find_atomic_number(std::string &key)
{ {
std::transform(key.begin(), key.end(), key.begin(), ::tolower); std::transform(key.begin(), key.end(), key.begin(), ::tolower);
if (key.length() == 1) key += " "; if (key.length() == 1) key += " ";
key.resize(2); if (key.length() > 2) return -1;
std::vector<std::string> element_table = { std::vector<std::string> element_table = {
"h ", "he", "li", "be", "b ", "c ", "n ", "o ", "f ", "ne", "na", "mg", "al", "si", "h ", "he", "li", "be", "b ", "c ", "n ", "o ", "f ", "ne", "na", "mg", "al", "si",
@ -283,6 +294,8 @@ void PairDispersionD3::coeff(int narg, char **arg)
for (int i = 0; i < ntypes; i++) { for (int i = 0; i < ntypes; i++) {
element = arg[i + 2]; element = arg[i + 2];
atomic_numbers[i] = find_atomic_number(element); atomic_numbers[i] = find_atomic_number(element);
if (atomic_numbers[i] < 0)
error->all(FLERR, Error::NOLASTLINE, "Element {} not supported", element);
} }
int count = 0; int count = 0;
@ -450,11 +463,6 @@ double *PairDispersionD3::get_dC6(int iat, int jat, double cni, double cnj)
void PairDispersionD3::compute(int eflag, int vflag) void PairDispersionD3::compute(int eflag, int vflag)
{ {
std::unordered_map<std::string, int> dampingMap = {
{"zero", 1}, {"zerom", 2}, {"bj", 3}, {"bjm", 4}};
int dampingCode = dampingMap[damping_type];
double evdwl = 0.0; double evdwl = 0.0;
ev_init(eflag, vflag); ev_init(eflag, vflag);
@ -518,7 +526,8 @@ void PairDispersionD3::compute(int eflag, int vflag)
t6 = t8 = e6 = e8 = evdwl = fpair = fpair1 = fpair2 = 0.0; t6 = t8 = e6 = e8 = evdwl = fpair = fpair1 = fpair2 = 0.0;
switch (dampingCode) { switch (dampingCode) {
case 1: { // zero
case 1: { // original
double r0 = r / r0ab[type[i]][type[j]]; double r0 = r / r0ab[type[i]][type[j]];
@ -539,6 +548,7 @@ void PairDispersionD3::compute(int eflag, int vflag)
fpair = fpair1 + fpair2; fpair = fpair1 + fpair2;
fpair *= factor_lj; fpair *= factor_lj;
} break; } break;
case 2: { // zerom case 2: { // zerom
double r0 = r0ab[type[i]][type[j]]; double r0 = r0ab[type[i]][type[j]];
@ -564,6 +574,7 @@ void PairDispersionD3::compute(int eflag, int vflag)
fpair = fpair1 + fpair2; fpair = fpair1 + fpair2;
fpair *= factor_lj; fpair *= factor_lj;
} break; } break;
case 3: { // bj case 3: { // bj
double r0 = sqrt(C8 / C6); double r0 = sqrt(C8 / C6);
@ -584,6 +595,7 @@ void PairDispersionD3::compute(int eflag, int vflag)
fpair = -(tmp6 + tmp8); fpair = -(tmp6 + tmp8);
fpair *= factor_lj; fpair *= factor_lj;
} break; } break;
case 4: { // bjm case 4: { // bjm
double r0 = sqrt(C8 / C6); double r0 = sqrt(C8 / C6);
@ -603,10 +615,15 @@ void PairDispersionD3::compute(int eflag, int vflag)
fpair = -(tmp6 + tmp8); fpair = -(tmp6 + tmp8);
fpair *= factor_lj; fpair *= factor_lj;
} } break;
default: {
// this should not happen with the error check in the init_style function
error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode);
} break;
} }
if (eflag) { evdwl = -(s6 * e6 + s8 * e8) * factor_lj; } if (eflag) evdwl = -(s6 * e6 + s8 * e8) * factor_lj;
double rest = (s6 * e6 + s8 * e8) / C6; double rest = (s6 * e6 + s8 * e8) / C6;
@ -693,12 +710,6 @@ void PairDispersionD3::compute(int eflag, int vflag)
void PairDispersionD3::set_funcpar(std::string &functional_name) void PairDispersionD3::set_funcpar(std::string &functional_name)
{ {
std::unordered_map<std::string, int> dampingMap = {
{"zero", 1}, {"zerom", 2}, {"bj", 3}, {"bjm", 4}};
int dampingCode = dampingMap[damping_type];
switch (dampingCode) { switch (dampingCode) {
case 1: { // zero case 1: { // zero
@ -962,14 +973,10 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
s8 = 1.206; s8 = 1.206;
break; break;
default: default:
error->all(FLERR, "Functional name unknown"); error->all(FLERR, Error::NOLASTLINE,
"Functional {} not supported with original damping function", functional_name);
break; break;
} }
//fprintf(stderr,"s6 : %f\n", s6);
//fprintf(stderr,"s8 : %f\n", s8);
//fprintf(stderr,"rs6 : %f\n", rs6);
//fprintf(stderr,"rs8 : %f\n", rs8);
//fprintf(stderr,"alpha : %f\n", alpha);
} break; } break;
case 2: { // zerom case 2: { // zerom
@ -1024,14 +1031,10 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
rs8 = 0.003160; rs8 = 0.003160;
break; break;
default: default:
error->all(FLERR, "Functional name unknown"); error->all(FLERR, Error::NOLASTLINE,
"Functional {} not supported with zerom damping function", functional_name);
break; break;
} }
//fprintf(stderr,"s6 : %f\n", s6);
//fprintf(stderr,"s8 : %f\n", s8);
//fprintf(stderr,"rs6 : %f\n", rs6);
//fprintf(stderr,"rs8 : %f\n", rs8);
//fprintf(stderr,"alpha : %f\n", alpha);
rs8 = rs8 / autoang; rs8 = rs8 / autoang;
} break; } break;
@ -1334,16 +1337,11 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
a2 = 4.5000; a2 = 4.5000;
break; break;
default: default:
error->all(FLERR, "Functional name unknown"); error->all(FLERR, Error::NOLASTLINE,
"Functional {} not supported with bj damping function", functional_name);
break; break;
} }
//fprintf(stderr,"s6 : %f\n", s6);
//fprintf(stderr,"s8 : %f\n", s8);
//fprintf(stderr,"a1 : %f\n", a1);
//fprintf(stderr,"a2 : %f\n", a2);
//fprintf(stderr,"alpha : %f\n", alpha);
a2 = a2 * autoang; a2 = a2 * autoang;
} break; } break;
@ -1400,21 +1398,17 @@ void PairDispersionD3::set_funcpar(std::string &functional_name)
a2 = 3.593680; a2 = 3.593680;
break; break;
default: default:
error->all(FLERR, "Functional name unknown"); error->all(FLERR, Error::NOLASTLINE,
"Functional {} not supported with bjm damping function", functional_name);
break; break;
} }
//fprintf(stderr,"s6 : %f\n", s6);
//fprintf(stderr,"s8 : %f\n", s8);
//fprintf(stderr,"a1 : %f\n", a1);
//fprintf(stderr,"a2 : %f\n", a2);
//fprintf(stderr,"alpha : %f\n", alpha);
a2 = a2 * autoang; a2 = a2 * autoang;
} break; } break;
default: default:
error->all(FLERR, "Damping type unknown"); // this should not happen with the error check in the init_style function
error->all(FLERR, Error::NOLASTLINE, "Damping code {} unknown", dampingCode);
break; break;
} }
} }
@ -1436,8 +1430,6 @@ double PairDispersionD3::init_one(int i, int j)
void PairDispersionD3::init_style() void PairDispersionD3::init_style()
{ {
if (atom->tag_enable == 0) error->all(FLERR, "Pair style D3 requires atom IDs"); if (atom->tag_enable == 0) error->all(FLERR, "Pair style D3 requires atom IDs");
//if (force->newton_pair == 0)
// error->all(FLERR,"Pair style D3 requires newton pair on");
// need an half neighbor list // need an half neighbor list
neighbor->add_request(this); neighbor->add_request(this);

2
src/EXTRA-PAIR/pair_dispersion_d3.h
View File

@ -48,7 +48,7 @@ class PairDispersionD3 : public Pair {
double rthr; // R^2 distance to cutoff for D3_calculation double rthr; // R^2 distance to cutoff for D3_calculation
double cn_thr; // R^2 distance to cutoff for CN_calculation double cn_thr; // R^2 distance to cutoff for CN_calculation
std::string damping_type; // damping function type int dampingCode;
double s6, s8, s18, rs6, rs8, rs18; // XC parameters double s6, s8, s18, rs6, rs8, rs18; // XC parameters
double a1, a2, alpha, alpha6, alpha8; double a1, a2, alpha, alpha6, alpha8;

7
src/OPENMP/angle_cosine_delta_omp.cpp
View File

@ -16,15 +16,16 @@
Contributing author: Axel Kohlmeyer (Temple U) Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "omp_compat.h"
#include "angle_cosine_delta_omp.h" #include "angle_cosine_delta_omp.h"
#include <cmath>
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "force.h" #include "force.h"
#include "neighbor.h" #include "neighbor.h"
#include <cmath>
#include "omp_compat.h"
#include "suffix.h" #include "suffix.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -128,8 +129,6 @@ void AngleCosineDeltaOMP::eval(int nfrom, int nto, ThrData * const thr)
s = sqrt(1.0 - c*c); s = sqrt(1.0 - c*c);
if (s < SMALL) s = SMALL; if (s < SMALL) s = SMALL;
s = 1.0/s;
cot = c/s; cot = c/s;
// force & energy // force & energy

2
src/OPENMP/fix_rigid_small_omp.cpp
View File

@ -229,7 +229,7 @@ void FixRigidSmallOMP::compute_forces_and_torques()
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp parallel for LMP_DEFAULT_NONE schedule(static) #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
#endif #endif
for (int ibody = 0; ibody < nbody; ibody++) { for (int ibody = 0; ibody < nlocal_body; ibody++) {
double * _noalias const fcm = body[ibody].fcm; double * _noalias const fcm = body[ibody].fcm;
const double mass = body[ibody].mass; const double mass = body[ibody].mass;
fcm[0] += gvec[0]*mass; fcm[0] += gvec[0]*mass;

4
src/OPENMP/pair_dpd_ext_omp.cpp
View File

@ -60,6 +60,10 @@ void PairDPDExtOMP::compute(int eflag, int vflag)
{ {
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
const int nall = atom->nlocal + atom->nghost; const int nall = atom->nlocal + atom->nghost;
const int inum = list->inum; const int inum = list->inum;

4
src/OPENMP/pair_dpd_ext_tstat_omp.cpp
View File

@ -60,6 +60,10 @@ void PairDPDExtTstatOMP::compute(int eflag, int vflag)
{ {
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
const int nall = atom->nlocal + atom->nghost; const int nall = atom->nlocal + atom->nghost;
const int inum = list->inum; const int inum = list->inum;

4
src/OPENMP/pair_dpd_omp.cpp
View File

@ -59,6 +59,10 @@ void PairDPDOMP::compute(int eflag, int vflag)
{ {
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
const int nall = atom->nlocal + atom->nghost; const int nall = atom->nlocal + atom->nghost;
const int inum = list->inum; const int inum = list->inum;

4
src/OPENMP/pair_dpd_tstat_omp.cpp
View File

@ -60,6 +60,10 @@ void PairDPDTstatOMP::compute(int eflag, int vflag)
{ {
ev_init(eflag,vflag); ev_init(eflag,vflag);
// precompute random force scaling factors
for (int i = 0; i < 4; ++i) special_sqrt[i] = sqrt(force->special_lj[i]);
const int nall = atom->nlocal + atom->nghost; const int nall = atom->nlocal + atom->nghost;
const int inum = list->inum; const int inum = list->inum;

12
src/REAXFF/reaxff_control.cpp
View File

@ -57,8 +57,8 @@ namespace ReaxFF {
class control_parser_error : public std::exception { class control_parser_error : public std::exception {
std::string message; std::string message;
public: public:
explicit control_parser_error(const std::string &msg) { explicit control_parser_error(const std::string &format, const std::string &keyword) {
message = msg; message = fmt::format(fmt::runtime(format), keyword);
} }
const char *what() const noexcept override { return message.c_str(); } const char *what() const noexcept override { return message.c_str(); }
}; };
@ -92,11 +92,11 @@ namespace ReaxFF {
auto keyword = values.next_string(); auto keyword = values.next_string();
if (!values.has_next()) if (!values.has_next())
throw control_parser_error( throw control_parser_error("No value(s) for control parameter: {}\n", keyword);
fmt::format("No value(s) for control parameter: {}\n", keyword));
if (inactive_keywords.find(keyword) != inactive_keywords.end()) { if (inactive_keywords.find(keyword) != inactive_keywords.end()) {
error->warning(FLERR,fmt::format("Ignoring inactive control parameter: {}", keyword)); error->warning(FLERR, fmt::format(fmt::runtime("Ignoring inactive control parameter: {}"),
keyword));
} else if (keyword == "nbrhood_cutoff") { } else if (keyword == "nbrhood_cutoff") {
control->bond_cut = values.next_double(); control->bond_cut = values.next_double();
} else if (keyword == "bond_graph_cutoff") { } else if (keyword == "bond_graph_cutoff") {
@ -114,7 +114,7 @@ namespace ReaxFF {
error->warning(FLERR,"Support for writing native trajectories has " error->warning(FLERR,"Support for writing native trajectories has "
"been removed after LAMMPS version 8 April 2021"); "been removed after LAMMPS version 8 April 2021");
} else { } else {
throw control_parser_error(fmt::format("Unknown parameter {} in control file", keyword)); throw control_parser_error("Unknown parameter {} in control file", keyword);
} }
} }
} catch (LAMMPS_NS::EOFException &) { } catch (LAMMPS_NS::EOFException &) {

9
src/RHEO/compute_rheo_vshift.cpp
View File

@ -42,8 +42,8 @@ using namespace MathExtra;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), vshift(nullptr), ct(nullptr), wsame(nullptr), cgradt(nullptr), Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), wsame(nullptr),
fix_rheo(nullptr), rho0(nullptr), list(nullptr), compute_interface(nullptr), ct(nullptr), cgradt(nullptr), shift_type(nullptr), list(nullptr), compute_interface(nullptr),
compute_kernel(nullptr), compute_surface(nullptr) compute_kernel(nullptr), compute_surface(nullptr)
{ {
if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command"); if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command");
@ -316,8 +316,7 @@ void ComputeRHEOVShift::correct_type_interface()
{ {
int i, j, a, ii, jj, jnum, itype, jtype; int i, j, a, ii, jj, jnum, itype, jtype;
int fluidi, fluidj; int fluidi, fluidj;
double xtmp, ytmp, ztmp, rsq, r, rinv; double xtmp, ytmp, ztmp, rsq, r, w;
double w, wp, dr, w0, prefactor;
double imass, jmass, voli, volj, rhoi, rhoj; double imass, jmass, voli, volj, rhoi, rhoj;
double dx[3]; double dx[3];
int dim = domain->dimension; int dim = domain->dimension;
@ -327,14 +326,12 @@ void ComputeRHEOVShift::correct_type_interface()
int *type = atom->type; int *type = atom->type;
int *status = atom->rheo_status; int *status = atom->rheo_status;
int *mask = atom->mask;
double **x = atom->x; double **x = atom->x;
double *rho = atom->rho; double *rho = atom->rho;
double *mass = atom->mass; double *mass = atom->mass;
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
inum = list->inum; inum = list->inum;

70
src/balance.cpp
View File

@ -1,3 +1,4 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -51,8 +52,6 @@ enum { XYZ, SHIFT, BISECTION };
enum { NONE, UNIFORM, USER }; enum { NONE, UNIFORM, USER };
enum { X, Y, Z }; enum { X, Y, Z };
// clang-format off
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
Balance::Balance(LAMMPS *lmp) : Command(lmp) Balance::Balance(LAMMPS *lmp) : Command(lmp)
@ -114,13 +113,13 @@ Balance::~Balance()
void Balance::command(int narg, char **arg) void Balance::command(int narg, char **arg)
{ {
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Balance command before simulation box is defined"); error->all(FLERR, -1, "Balance command before simulation box is defined");
if (comm->me == 0) utils::logmesg(lmp,"Balancing ...\n"); if (comm->me == 0) utils::logmesg(lmp,"Balancing ...\n");
// parse required arguments // parse required arguments
if (narg < 2) error->all(FLERR,"Illegal balance command"); if (narg < 2) utils::missing_cmd_args(FLERR,"balance", error);
thresh = utils::numeric(FLERR,arg[0],false,lmp); thresh = utils::numeric(FLERR,arg[0],false,lmp);
@ -132,16 +131,15 @@ void Balance::command(int narg, char **arg)
int iarg = 1; int iarg = 1;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"x") == 0) { if (strcmp(arg[iarg],"x") == 0) {
if (style != -1 && style != XYZ) if ((style != -1) && (style != XYZ))
error->all(FLERR,"Illegal balance command"); error->all(FLERR, iarg, "Must not combine multiple balancing styles");
style = XYZ; style = XYZ;
if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"balance x", error);
if (strcmp(arg[iarg+1],"uniform") == 0) { if (strcmp(arg[iarg+1],"uniform") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal balance command");
xflag = UNIFORM; xflag = UNIFORM;
iarg += 2; iarg += 2;
} else { } else {
if (1 + procgrid[0]-1 > narg) if (iarg + procgrid[0]-1 > narg) utils::missing_cmd_args(FLERR,"balance x", error);
error->all(FLERR,"Illegal balance command");
xflag = USER; xflag = USER;
delete[] user_xsplit; delete[] user_xsplit;
user_xsplit = new double[procgrid[0]+1]; user_xsplit = new double[procgrid[0]+1];
@ -152,16 +150,15 @@ void Balance::command(int narg, char **arg)
user_xsplit[procgrid[0]] = 1.0; user_xsplit[procgrid[0]] = 1.0;
} }
} else if (strcmp(arg[iarg],"y") == 0) { } else if (strcmp(arg[iarg],"y") == 0) {
if (style != -1 && style != XYZ) if ((style != -1) && (style != XYZ))
error->all(FLERR,"Illegal balance command"); error->all(FLERR, iarg, "Must not combine multiple balancing styles");
style = XYZ; style = XYZ;
if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"balance y", error);
if (strcmp(arg[iarg+1],"uniform") == 0) { if (strcmp(arg[iarg+1],"uniform") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal balance command");
yflag = UNIFORM; yflag = UNIFORM;
iarg += 2; iarg += 2;
} else { } else {
if (1 + procgrid[1]-1 > narg) if (iarg + procgrid[1]-1 > narg) utils::missing_cmd_args(FLERR,"balance y", error);
error->all(FLERR,"Illegal balance command");
yflag = USER; yflag = USER;
delete[] user_ysplit; delete[] user_ysplit;
user_ysplit = new double[procgrid[1]+1]; user_ysplit = new double[procgrid[1]+1];
@ -172,16 +169,15 @@ void Balance::command(int narg, char **arg)
user_ysplit[procgrid[1]] = 1.0; user_ysplit[procgrid[1]] = 1.0;
} }
} else if (strcmp(arg[iarg],"z") == 0) { } else if (strcmp(arg[iarg],"z") == 0) {
if (style != -1 && style != XYZ) if ((style != -1) && (style != XYZ))
error->all(FLERR,"Illegal balance command"); error->all(FLERR, iarg, "Must not combine multiple balancing styles");
style = XYZ; style = XYZ;
if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"balance z", error);
if (strcmp(arg[iarg+1],"uniform") == 0) { if (strcmp(arg[iarg+1],"uniform") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal balance command");
zflag = UNIFORM; zflag = UNIFORM;
iarg += 2; iarg += 2;
} else { } else {
if (1 + procgrid[2]-1 > narg) if (iarg + procgrid[2]-1 > narg) utils::missing_cmd_args(FLERR,"balance z", error);
error->all(FLERR,"Illegal balance command");
zflag = USER; zflag = USER;
delete[] user_zsplit; delete[] user_zsplit;
user_zsplit = new double[procgrid[2]+1]; user_zsplit = new double[procgrid[2]+1];
@ -193,19 +189,21 @@ void Balance::command(int narg, char **arg)
} }
} else if (strcmp(arg[iarg],"shift") == 0) { } else if (strcmp(arg[iarg],"shift") == 0) {
if (style != -1) error->all(FLERR,"Illegal balance command"); if (style != -1) error->all(FLERR, iarg, "Must not combine multiple balancing styles");
if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "balance shift", error); if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "balance shift", error);
style = SHIFT; style = SHIFT;
bstr = arg[iarg+1]; bstr = arg[iarg+1];
if (bstr.size() > BSTR_SIZE) error->all(FLERR,"Illegal balance shift command"); if (bstr.size() > BSTR_SIZE) error->all(FLERR, iarg + 1, "Illegal balance shift command");
nitermax = utils::inumeric(FLERR,arg[iarg+2],false,lmp); nitermax = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
if (nitermax <= 0) error->all(FLERR,"Illegal balance command"); if (nitermax <= 0)
error->all(FLERR,iarg+2,"Illegal balance shift value {}", arg[iarg+2]);
stopthresh = utils::numeric(FLERR,arg[iarg+3],false,lmp); stopthresh = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (stopthresh < 1.0) error->all(FLERR,"Illegal balance command"); if (stopthresh < 1.0)
error->all(FLERR,iarg+3,"Illegal balance stop threshold value {}", arg[iarg+3]);
iarg += 4; iarg += 4;
} else if (strcmp(arg[iarg],"rcb") == 0) { } else if (strcmp(arg[iarg],"rcb") == 0) {
if (style != -1) error->all(FLERR,"Illegal balance command"); if (style != -1) error->all(FLERR,"Must not combine multiple balancing styles");
style = BISECTION; style = BISECTION;
iarg++; iarg++;
@ -221,32 +219,32 @@ void Balance::command(int narg, char **arg)
if (xflag == USER) if (xflag == USER)
for (int i = 1; i <= procgrid[0]; i++) for (int i = 1; i <= procgrid[0]; i++)
if (user_xsplit[i-1] >= user_xsplit[i]) if (user_xsplit[i-1] >= user_xsplit[i])
error->all(FLERR,"Illegal balance command"); error->all(FLERR,"Slices for balance x command must be in ascending order");
if (yflag == USER) if (yflag == USER)
for (int i = 1; i <= procgrid[1]; i++) for (int i = 1; i <= procgrid[1]; i++)
if (user_ysplit[i-1] >= user_ysplit[i]) if (user_ysplit[i-1] >= user_ysplit[i])
error->all(FLERR,"Illegal balance command"); error->all(FLERR,"Slices for balance y command must be in ascending order");
if (zflag == USER) if (zflag == USER)
for (int i = 1; i <= procgrid[2]; i++) for (int i = 1; i <= procgrid[2]; i++)
if (user_zsplit[i-1] >= user_zsplit[i]) if (user_zsplit[i-1] >= user_zsplit[i])
error->all(FLERR,"Illegal balance command"); error->all(FLERR,"Slices for balance z command must be in ascending order");
} }
if (style == SHIFT) { if (style == SHIFT) {
const int blen = bstr.size(); const int blen = bstr.size();
for (int i = 0; i < blen; i++) { for (int i = 0; i < blen; i++) {
if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z') if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z')
error->all(FLERR,"Balance shift string is invalid"); error->all(FLERR,"Balance shift string {} is invalid", bstr);
if (bstr[i] == 'z' && dimension == 2) if (bstr[i] == 'z' && dimension == 2)
error->all(FLERR,"Balance shift string is invalid"); error->all(FLERR,"Balance shift string {} is invalid", bstr);
for (int j = i+1; j < blen; j++) for (int j = i+1; j < blen; j++)
if (bstr[i] == bstr[j]) if (bstr[i] == bstr[j])
error->all(FLERR,"Balance shift string is invalid"); error->all(FLERR,"Balance shift string {} is invalid", bstr);
} }
} }
if (style == BISECTION && comm->style == Comm::BRICK) if (style == BISECTION && comm->style == Comm::BRICK)
error->all(FLERR,"Balance rcb cannot be used with comm_style brick"); error->all(FLERR, Error::ARGZERO, "Balance rcb cannot be used with comm_style brick");
// process remaining optional args // process remaining optional args
@ -467,16 +465,16 @@ void Balance::options(int iarg, int narg, char **arg, int sortflag_default)
nopt = imb->options(narg-iarg,arg+iarg+2); nopt = imb->options(narg-iarg,arg+iarg+2);
imbalances[nimbalance++] = imb; imbalances[nimbalance++] = imb;
} else { } else {
error->all(FLERR,"Unknown (fix) balance weight method: {}", arg[iarg+1]); error->all(FLERR, iarg + 1, "Unknown (fix) balance weight method: {}", arg[iarg+1]);
} }
iarg += 2+nopt; iarg += 2+nopt;
} else if (strcmp(arg[iarg],"sort") == 0) { } else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "balance sort", error); if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "(fix) balance sort", error);
sortflag = utils::logical(FLERR,arg[iarg+1],false,lmp); sortflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"out") == 0) { } else if (strcmp(arg[iarg],"out") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal (fix) balance command"); if (iarg+2 > narg) utils::missing_cmd_args(FLERR,"(fix) balance out", error);
outflag = 1; outflag = 1;
outarg = iarg+1; outarg = iarg+1;
iarg += 2; iarg += 2;
@ -485,7 +483,7 @@ void Balance::options(int iarg, int narg, char **arg, int sortflag_default)
oldrcb = 1; oldrcb = 1;
iarg++; iarg++;
} else error->all(FLERR,"Illegal (fix) balance command"); } else error->all(FLERR, iarg, "Unknown (fix) balance keyword {}", arg[iarg]);
} }
// output file // output file
@ -493,7 +491,7 @@ void Balance::options(int iarg, int narg, char **arg, int sortflag_default)
if (outflag && comm->me == 0) { if (outflag && comm->me == 0) {
fp = fopen(arg[outarg],"w"); fp = fopen(arg[outarg],"w");
if (fp == nullptr) if (fp == nullptr)
error->one(FLERR,"Cannot open (fix) balance output file {}: {}", error->one(FLERR, outarg, "Cannot open (fix) balance output file {}: {}",
arg[outarg], utils::getsyserror()); arg[outarg], utils::getsyserror());
} }
} }

23
src/fix_balance.cpp
View File

@ -42,7 +42,7 @@ enum { SHIFT, BISECTION };
FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) : FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), balance(nullptr), irregular(nullptr) Fix(lmp, narg, arg), balance(nullptr), irregular(nullptr)
{ {
if (narg < 6) error->all(FLERR,"Illegal fix balance command"); if (narg < 6) utils::missing_cmd_args(FLERR, "fix balance", error);
box_change = BOX_CHANGE_DOMAIN; box_change = BOX_CHANGE_DOMAIN;
pre_exchange_migrate = 1; pre_exchange_migrate = 1;
@ -58,7 +58,7 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
int dimension = domain->dimension; int dimension = domain->dimension;
nevery = utils::inumeric(FLERR,arg[3],false,lmp); nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery < 0) error->all(FLERR,"Illegal fix balance command"); if (nevery < 0) error->all(FLERR, 3, "Illegal fix balance nevery value {}", arg[3]);
thresh = utils::numeric(FLERR,arg[4],false,lmp); thresh = utils::numeric(FLERR,arg[4],false,lmp);
reportonly = 0; reportonly = 0;
@ -69,7 +69,7 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[5],"report") == 0) { } else if (strcmp(arg[5],"report") == 0) {
lbstyle = SHIFT; lbstyle = SHIFT;
reportonly = 1; reportonly = 1;
} else error->all(FLERR,"Unknown fix balance style {}", arg[5]); } else error->all(FLERR, 5, "Unknown fix balance style {}", arg[5]);
int iarg = 5; int iarg = 5;
if (lbstyle == SHIFT) { if (lbstyle == SHIFT) {
@ -84,11 +84,14 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
} else { } else {
if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix balance shift", error); if (iarg+4 > narg) utils::missing_cmd_args(FLERR, "fix balance shift", error);
bstr = arg[iarg+1]; bstr = arg[iarg+1];
if (bstr.size() > Balance::BSTR_SIZE) error->all(FLERR,"Illegal fix balance shift command"); if (bstr.size() > Balance::BSTR_SIZE)
error->all(FLERR, iarg + 1, "Illegal fix balance shift argument {}", bstr);
nitermax = utils::inumeric(FLERR,arg[iarg+2],false,lmp); nitermax = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
if (nitermax <= 0) error->all(FLERR,"Illegal fix balance command"); if (nitermax <= 0)
error->all(FLERR, iarg + 2, "Illegal fix balance shift niter argument {}", nitermax);
stopthresh = utils::numeric(FLERR,arg[iarg+3],false,lmp); stopthresh = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (stopthresh < 1.0) error->all(FLERR,"Illegal fix balance command"); if (stopthresh < 1.0)
error->all(FLERR, iarg + 3, "Illegal fix balance stop threshold argument {}", stopthresh);
iarg += 4; iarg += 4;
} }
@ -102,12 +105,12 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
const int blen = bstr.size(); const int blen = bstr.size();
for (int i = 0; i < blen; i++) { for (int i = 0; i < blen; i++) {
if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z') if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z')
error->all(FLERR,"Fix balance shift string is invalid"); error->all(FLERR,"Fix balance shift string {} is invalid", bstr);
if (bstr[i] == 'z' && dimension == 2) if (bstr[i] == 'z' && dimension == 2)
error->all(FLERR,"Fix balance shift string is invalid"); error->all(FLERR,"Fix balance shift string {} is invalid", bstr);
for (int j = i+1; j < blen; j++) for (int j = i+1; j < blen; j++)
if (bstr[i] == bstr[j]) if (bstr[i] == bstr[j])
error->all(FLERR,"Fix balance shift string is invalid"); error->all(FLERR,"Fix balance shift string {} is invalid", bstr);
} }
} }
@ -125,7 +128,7 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
sortflag = balance->sortflag; sortflag = balance->sortflag;
if (balance->varflag && nevery == 0) if (balance->varflag && nevery == 0)
error->all(FLERR,"Fix balance nevery = 0 cannot be used with weight var"); error->all(FLERR, 3, "Fix balance nevery = 0 cannot be used with weight var");
// create instance of Irregular class // create instance of Irregular class

65
src/region_plane.cpp
View File

@ -25,9 +25,11 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) : RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr) Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr),
nxstr(nullptr), nystr(nullptr), nzstr(nullptr)
{ {
xvar = yvar = zvar = 0.0; xvar = yvar = zvar = 0.0;
nxvar = nyvar = nzvar = 0.0;
options(narg - 8, &arg[8]); options(narg - 8, &arg[8]);
@ -61,15 +63,34 @@ RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
zstyle = CONSTANT; zstyle = CONSTANT;
} }
int nxstyle = (utils::strmatch(arg[5], "^v_")) ? 1 : 0;
int nystyle = (utils::strmatch(arg[6], "^v_")) ? 1 : 0;
int nzstyle = (utils::strmatch(arg[7], "^v_")) ? 1 : 0;
if (!nxstyle && !nystyle && !nzstyle) nstyle = CONSTANT;
if (nxstyle && nystyle && nzstyle) nstyle = VARIABLE;
if (nstyle == CONSTANT) {
normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
} else if (nstyle == VARIABLE) {
normal[0] = 0.0;
normal[1] = 0.0;
normal[2] = 0.0;
nxstr = utils::strdup(arg[5] + 2);
nystr = utils::strdup(arg[6] + 2);
nzstr = utils::strdup(arg[7] + 2);
varshape = 1;
} else {
error->all(FLERR, "The components of the normal vector should be either all variables or all constants");
}
if (varshape) { if (varshape) {
variable_check(); variable_check();
RegPlane::shape_update(); RegPlane::shape_update();
} }
normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
// enforce unit normal // enforce unit normal
double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]; double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
@ -95,6 +116,9 @@ RegPlane::~RegPlane()
delete[] xstr; delete[] xstr;
delete[] ystr; delete[] ystr;
delete[] zstr; delete[] zstr;
delete[] nxstr;
delete[] nystr;
delete[] nzstr;
delete[] contact; delete[] contact;
} }
@ -174,6 +198,20 @@ void RegPlane::shape_update()
if (ystyle == VARIABLE) yp = yscale * input->variable->compute_equal(yvar); if (ystyle == VARIABLE) yp = yscale * input->variable->compute_equal(yvar);
if (zstyle == VARIABLE) zp = zscale * input->variable->compute_equal(zvar); if (zstyle == VARIABLE) zp = zscale * input->variable->compute_equal(zvar);
if (nstyle == VARIABLE) {
normal[0] = xscale * input->variable->compute_equal(nxvar);
normal[1] = yscale * input->variable->compute_equal(nyvar);
normal[2] = zscale * input->variable->compute_equal(nzvar);
// enforce unit normal
double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
if (rsq == 0.0) error->all(FLERR, "Illegal region plane normal vector: {} {} {}", normal[0], normal[1], normal[2]);
normal[0] /= sqrt(rsq);
normal[1] /= sqrt(rsq);
normal[2] /= sqrt(rsq);
}
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -202,4 +240,21 @@ void RegPlane::variable_check()
if (!input->variable->equalstyle(zvar)) if (!input->variable->equalstyle(zvar))
error->all(FLERR, "Variable {} for region plane is invalid style", zstr); error->all(FLERR, "Variable {} for region plane is invalid style", zstr);
} }
if (nstyle == VARIABLE) {
nxvar = input->variable->find(nxstr);
if (nxvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nxstr);
if (!input->variable->equalstyle(nxvar))
error->all(FLERR, "Variable {} for region plane is invalid style", nxstr);
nyvar = input->variable->find(nystr);
if (nyvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nystr);
if (!input->variable->equalstyle(nyvar))
error->all(FLERR, "Variable {} for region plane is invalid style", nystr);
nzvar = input->variable->find(nzstr);
if (nzvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nzstr);
if (!input->variable->equalstyle(nzvar))
error->all(FLERR, "Variable {} for region plane is invalid style", nzstr);
}
} }

4
src/region_plane.h
View File

@ -45,6 +45,10 @@ class RegPlane : public Region {
int zstyle, zvar; int zstyle, zvar;
char *xstr, *ystr, *zstr; char *xstr, *ystr, *zstr;
int nstyle;
int nxvar, nyvar, nzvar;
char *nxstr, *nystr, *nzstr;
void variable_check(); void variable_check();
}; };

2
src/utils.cpp
View File

@ -156,7 +156,7 @@ std::string utils::strcompress(const std::string &text)
} }
// remove trailing blank // remove trailing blank
if (output.back() == ' ') output.erase(output.size() - 1, 1); if (!output.empty() && output.back() == ' ') output.erase(output.size() - 1, 1);
return output; return output;
} }

8
unittest/commands/test_simple_commands.cpp
View File

@ -415,19 +415,19 @@ TEST_F(SimpleCommandsTest, Plugin)
if (!bindir) GTEST_SKIP() << "LAMMPS_PLUGIN_DIR not set"; if (!bindir) GTEST_SKIP() << "LAMMPS_PLUGIN_DIR not set";
std::string loadfmt = "plugin load {}/{}plugin.so"; std::string loadfmt = "plugin load {}/{}plugin.so";
::testing::internal::CaptureStdout(); ::testing::internal::CaptureStdout();
lmp->input->one(fmt::format(loadfmt, bindir, "hello")); lmp->input->one(fmt::format(fmt::runtime(loadfmt), bindir, "hello"));
auto text = ::testing::internal::GetCapturedStdout(); auto text = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << text; if (verbose) std::cout << text;
ASSERT_THAT(text, ContainsRegex(".*\n.*Loading plugin: Hello world command.*")); ASSERT_THAT(text, ContainsRegex(".*\n.*Loading plugin: Hello world command.*"));
::testing::internal::CaptureStdout(); ::testing::internal::CaptureStdout();
lmp->input->one(fmt::format(loadfmt, bindir, "xxx")); lmp->input->one(fmt::format(fmt::runtime(loadfmt), bindir, "xxx"));
text = ::testing::internal::GetCapturedStdout(); text = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << text; if (verbose) std::cout << text;
ASSERT_THAT(text, ContainsRegex(".*Open of file .*xxx.* failed.*")); ASSERT_THAT(text, ContainsRegex(".*Open of file .*xxx.* failed.*"));
::testing::internal::CaptureStdout(); ::testing::internal::CaptureStdout();
lmp->input->one(fmt::format(loadfmt, bindir, "nve2")); lmp->input->one(fmt::format(fmt::runtime(loadfmt), bindir, "nve2"));
text = ::testing::internal::GetCapturedStdout(); text = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << text; if (verbose) std::cout << text;
ASSERT_THAT(text, ContainsRegex(".*Loading plugin: NVE2 variant fix style.*")); ASSERT_THAT(text, ContainsRegex(".*Loading plugin: NVE2 variant fix style.*"));
@ -438,7 +438,7 @@ TEST_F(SimpleCommandsTest, Plugin)
ASSERT_THAT(text, ContainsRegex(".*1: command style plugin hello\n.*2: fix style plugin nve2.*")); ASSERT_THAT(text, ContainsRegex(".*1: command style plugin hello\n.*2: fix style plugin nve2.*"));
::testing::internal::CaptureStdout(); ::testing::internal::CaptureStdout();
lmp->input->one(fmt::format(loadfmt, bindir, "hello")); lmp->input->one(fmt::format(fmt::runtime(loadfmt), bindir, "hello"));
text = ::testing::internal::GetCapturedStdout(); text = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << text; if (verbose) std::cout << text;
ASSERT_THAT(text, ContainsRegex(".*Ignoring load of command style hello: " ASSERT_THAT(text, ContainsRegex(".*Ignoring load of command style hello: "

127
unittest/force-styles/tests/angle-cosine_delta.yaml
View File

@ -1,6 +1,7 @@
--- ---
lammps_version: 17 Feb 2022 lammps_version: 4 Feb 2025
date_generated: Fri Mar 18 22:17:51 2022 tags: generated
date_generated: Fri Feb 21 16:17:55 2025
epsilon: 5e-13 epsilon: 5e-13
skip_tests: numdiff skip_tests: numdiff
prerequisites: ! | prerequisites: ! |
@ -20,68 +21,68 @@ extract: ! ""
natoms: 29 natoms: 29
init_energy: 20.65098133017812 init_energy: 20.65098133017812
init_stress: ! |2- init_stress: ! |2-
4.0957903125808471e+01 -2.2499632589286236e+01 -1.8458270536522217e+01 3.4849378298981264e+01 -4.9327943096381759e+00 -2.6124108916945294e+00 4.4473180621356008e+01 -4.3463840015691233e+01 -1.0093406056647736e+00 4.5637584879601796e+01 -1.3907393586961433e+01 -4.2824939510892590e-01
init_forces: ! |2 init_forces: ! |2
1 1.6023996455832325e+01 2.4226311046592137e+00 -1.8589204266697958e+01 1 2.3703974495881418e+01 3.8584484587903516e+00 -1.5929138796293980e+01
2 -9.3711866901463381e-02 -7.5876216223294790e-01 -6.1013759002633705e-01 2 -5.5592010971213179e-01 -4.5011497307751007e+00 -3.6194749630118039e+00
3 -6.0380408211033587e-01 2.6993032298913775e+01 2.7951913775716850e+01 3 -2.8457465117799385e+00 2.6325352833001663e+01 2.5623331548006441e+01
4 -4.3564633502451269e+00 -1.0062417183984081e+01 2.1642186242701338e+00 4 -7.9894054979056088e+00 -1.2258360452072225e+01 6.4168233111021600e+00
5 -1.5386439746171822e+01 -2.3440313238754023e+01 -1.0067874149174445e+01 5 -2.2204738294831444e+01 -2.3698109015673438e+01 -1.1480776479805904e+01
6 9.6873286040901903e+00 6.5740743013867231e-02 -2.6313308270088624e+00 6 2.3347569650102201e+01 -4.7137825565345990e+00 -1.6078168111448818e+01
7 -4.1923658452230708e+00 3.7873664481235689e+00 -9.2297211491112741e-01 7 -8.0686673030183851e+00 7.1975893900877121e+00 -1.7579980765823877e+00
8 -4.1328674241018124e+00 -1.3752082356191019e+01 9.4169257319147519e+00 8 -3.1824458002558198e+00 -1.2503128164663320e+01 4.0602708742130751e+01
9 -5.2849046978682424e+00 5.1551127779597277e+00 -6.0440022931290430e-02 9 -7.7323033517797342e+00 8.0399091320421157e+00 -2.0951400051039776e-01
10 8.5785584536205555e+00 1.9092658049053345e+01 -1.3500832217122007e+01 10 4.9524658729988271e+00 2.0621211528107622e+01 -3.2638235652308239e+01
11 -8.2548956740792931e+00 -3.8384442674446007e+00 4.3703985032254611e+00 11 -1.0127717019112277e+01 -3.1783095109070576e+00 4.4288595135259623e+00
12 -2.2390516462227508e+00 -2.4339015686808576e+00 -5.7125783103934751e+00 12 -5.4777124732518541e+00 -5.3066394246038655e+00 -4.8409352484994574e+00
13 7.0748633959583906e-01 1.7495942830019229e+00 8.1749690284178111e-01 13 1.4412505843466221e+00 3.6458164082141984e+00 -5.3031504624386550e-02
14 5.0293855034635815e+00 -4.3347315249911471e+00 -2.2794877131015205e+00 14 7.6148261062259106e+00 -4.7428356370992075e+00 -3.3084707170534466e+00
15 4.2836129602031709e+00 -2.6803643297135835e+00 7.8342386854077670e+00 15 6.1963584357832016e+00 -3.1286786554578465e+00 9.2319462456089330e+00
16 -1.6475295059196460e+00 4.5975354658004379e+00 4.3681254253354895e-01 16 1.9503516089515505e-01 5.3411834252929946e+00 3.0732566957386886e+00
17 1.8816655220378973e+00 -2.5626545385335930e+00 1.3828524454567255e+00 17 7.3317605541385511e-01 -9.9851802775000564e-01 5.3881749402587964e-01
18 1.6064693296025290e-01 1.6787555240920624e+00 -6.8129331975320220e+00 18 7.8915981928254553e-02 8.2466959163160869e-01 -3.3467760833494964e+00
19 -1.9670041679005905e+00 -2.4504879220906162e+00 2.9176882347102016e+00 19 -9.6626846530114285e-01 -1.2037743703637904e+00 1.4332812196273430e+00
20 1.8063572349403376e+00 7.7173239799855375e-01 3.8952449628218204e+00 20 8.8735248337288830e-01 3.7910477873218174e-01 1.9134948637221534e+00
21 2.3714039975200247e+00 2.9686059791249759e+00 -8.8717809550958702e+00 21 1.7083446765076102e+00 2.1385652661420855e+00 -6.3911757682915322e+00
22 -4.2227815913112092e+00 -2.8790902397071427e+00 3.1362978814995639e+00 22 -3.0420655692218990e+00 -2.0740786847506296e+00 2.2593694686375940e+00
23 1.8513775937911845e+00 -8.9515739417833418e-02 5.7354830735963063e+00 23 1.3337208927142887e+00 -6.4486581391455855e-02 4.1318062996539382e+00
24 -1.4568411370574479e+00 5.6628278324597936e+00 -3.3714304805753512e+00 24 -6.5262931439741845e-01 2.5368088200155672e+00 -1.5103186662621493e+00
25 -1.1685952378028448e+00 -3.9152363566022288e+00 6.4031275969729240e-01 25 -5.2350217841582225e-01 -1.7539304417736190e+00 2.8684450674231776e-01
26 2.6254363748602927e+00 -1.7475914758575648e+00 2.7311177208780588e+00 26 1.1761314928132407e+00 -7.8287837824194817e-01 1.2234741595198315e+00
27 -6.3244441950917718e-01 5.1918732409438757e+00 -1.7685243510997146e+00 27 -1.4360109768190377e-01 1.1788525179864244e+00 -4.0155629532187376e-01
28 -1.3828773692511866e+00 -3.1042973063575481e+00 7.2112563744890223e-02 28 -3.1399234787785879e-01 -7.0485324397341742e-01 1.6373681213693980e-02
29 2.0153217887603638e+00 -2.0875759345863276e+00 1.6964117873548243e+00 29 4.5759344555976256e-01 -4.7399927401300701e-01 3.8518261410817978e-01
run_energy: 20.572797398210106 run_energy: 20.563089678525095
run_stress: ! |2- run_stress: ! |2-
4.0414601145881797e+01 -2.2259720000355955e+01 -1.8154881145525852e+01 3.4356097430708779e+01 -5.1760682620031364e+00 -3.0219001073831251e+00 4.3757504913648717e+01 -4.3136475705003690e+01 -6.2102920864502398e-01 4.4978336549538170e+01 -1.4195605982625462e+01 -9.5958875069195426e-01
run_forces: ! |2 run_forces: ! |2
1 1.5854493147994468e+01 2.5625041903632799e+00 -1.8403707584219106e+01 1 2.3451032879726370e+01 4.0617257197173897e+00 -1.5621894955519736e+01
2 -1.0778147157154638e-01 -8.7954374413252623e-01 -7.1033263283976922e-01 2 -5.7421401726330101e-01 -4.6835780748393674e+00 -3.7823299076024632e+00
3 -3.5903740231368531e-01 2.6747382244473858e+01 2.7776094765308645e+01 3 -2.4306027199459486e+00 2.6072519970009257e+01 2.5364241217921947e+01
4 -4.3667455132582980e+00 -9.9043336720965058e+00 2.2554137013679947e+00 4 -8.0249953558852880e+00 -1.2074108242088608e+01 6.5482299702187063e+00
5 -1.5435306686274892e+01 -2.3368620551346410e+01 -1.0071205010683158e+01 5 -2.2289140460440823e+01 -2.3631032649621293e+01 -1.1494387406884398e+01
6 9.7004722937707939e+00 2.8185574443369887e-02 -2.6469111215932326e+00 6 2.3354057831390289e+01 -4.7586966738793901e+00 -1.6119827865958250e+01
7 -4.2051376166149819e+00 3.8059875922956321e+00 -9.2948254913796591e-01 7 -8.0844064733133774e+00 7.2149529473510370e+00 -1.7660273021312418e+00
8 -4.0811558306056437e+00 -1.3725195219352738e+01 9.4305267102762365e+00 8 -3.1478296151921232e+00 -1.2498537099228518e+01 4.0622223121616166e+01
9 -5.3018040121457650e+00 5.1712375236873847e+00 -5.5170438746634645e-02 9 -7.7511017876278920e+00 8.0605427659852786e+00 -2.0241485109620189e-01
10 8.5884631851414728e+00 1.8977169595422669e+01 -1.3503733846157115e+01 10 4.9915744598623162e+00 2.0482598351331305e+01 -3.2624972375682837e+01
11 -8.2489793835708181e+00 -3.8307795657532537e+00 4.3670355271621792e+00 11 -1.0123311005215982e+01 -3.1688924167926276e+00 4.4252676678697700e+00
12 -2.2854196014066837e+00 -2.5730964210485521e+00 -5.5327026624884477e+00 12 -5.5513668937759473e+00 -5.4355045396603447e+00 -4.6167965020560917e+00
13 6.8857552776933151e-01 1.7099649678302227e+00 8.4476951940005041e-01 13 1.4059367926919679e+00 3.5728783709008667e+00 2.5935291038754471e-03
14 5.0388676050555423e+00 -4.1380632752237698e+00 -2.2619786920581912e+00 14 7.6143970879681024e+00 -4.5037794326053131e+00 -3.2831360392059907e+00
15 4.2713089863647395e+00 -2.6164691252678471e+00 7.6327855841331855e+00 15 6.2058701073495151e+00 -3.0406636127154227e+00 8.9669093364259496e+00
16 -1.6264161996020632e+00 4.5902884407035716e+00 4.2906258286847709e-01 16 2.2902239968198557e-01 5.3178918884056445e+00 3.0490139755759853e+00
17 1.8756029712680284e+00 -2.5566185549983906e+00 1.3795361474068475e+00 17 7.2507676999013526e-01 -9.8831727226989585e-01 5.3330838740480679e-01
18 1.5949074618033343e-01 1.6020315457790089e+00 -6.4757984092122420e+00 18 7.1127024594080268e-02 7.1445956198119154e-01 -2.8880230496978854e+00
19 -1.8652368265982364e+00 -2.3286173150011749e+00 2.7801057275919980e+00 19 -8.3188306575710236e-01 -1.0385320257838888e+00 1.2398389619370178e+00
20 1.7057460804179030e+00 7.2658576922216600e-01 3.6956926816202440e+00 20 7.6075604116302209e-01 3.2407246380269727e-01 1.6481840877608676e+00
21 2.3015994659382928e+00 2.8609989652412211e+00 -8.6002184142674114e+00 21 1.6140194966539338e+00 2.0063050011497334e+00 -6.0309928490009383e+00
22 -4.0843967543541595e+00 -2.7748751066344699e+00 3.0465239349836377e+00 22 -2.8643056040092603e+00 -1.9459468475255317e+00 2.1363729514270524e+00
23 1.7827972884158667e+00 -8.6123858606751036e-02 5.5536944792837737e+00 23 1.2502861073553264e+00 -6.0358153624201771e-02 3.8946198975738859e+00
24 -1.4105671680157181e+00 5.4757757163248080e+00 -3.2567653472769562e+00 24 -5.8408404963238691e-01 2.2673919964618818e+00 -1.3485518584960339e+00
25 -1.1213719131694626e+00 -3.7848791001392499e+00 6.1939772997237075e-01 25 -4.6437462735412272e-01 -1.5672547533259689e+00 2.5645620543256875e-01
26 2.5319390811851807e+00 -1.6908966161855581e+00 2.6373676173045855e+00 26 1.0484586769865096e+00 -7.0013724313591297e-01 1.0920956530634651e+00
27 -6.3028932344593613e-01 5.1333019072249222e+00 -1.7384473964198608e+00 27 -1.3425430939812166e-01 1.0934202317724093e+00 -3.7029815332428351e-01
28 -1.3620336377103808e+00 -3.0692701554272155e+00 6.6735136389115879e-02 28 -2.9014265041938359e-01 -6.5377692091134065e-01 1.4204173363709049e-02
29 1.9923229611563169e+00 -2.0640317517977067e+00 1.6717122600307448e+00 29 4.2439695981750525e-01 -4.3964331086106867e-01 3.5609397996057446e-01
... ...

2
unittest/force-styles/tests/bond-class2.yaml
View File

@ -1,7 +1,7 @@
--- ---
lammps_version: 17 Feb 2022 lammps_version: 17 Feb 2022
date_generated: Fri Mar 18 22:17:50 2022 date_generated: Fri Mar 18 22:17:50 2022
epsilon: 4e-13 epsilon: 5e-13
skip_tests: skip_tests:
prerequisites: ! | prerequisites: ! |
atom full atom full

38
unittest/utils/test_fmtlib.cpp
View File

@ -24,21 +24,21 @@ using ::testing::Eq;
TEST(FmtLib, insert_string) TEST(FmtLib, insert_string)
{ {
const char val[] = "word"; constexpr char val[] = "word";
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word word")); ASSERT_THAT(text, Eq("word word"));
} }
TEST(FmtLib, insert_int) TEST(FmtLib, insert_int)
{ {
const int val = 333; constexpr int val = 333;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word 333")); ASSERT_THAT(text, Eq("word 333"));
} }
TEST(FmtLib, insert_neg_int) TEST(FmtLib, insert_neg_int)
{ {
const int val = -333; constexpr int val = -333;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word -333")); ASSERT_THAT(text, Eq("word -333"));
} }
@ -46,7 +46,7 @@ TEST(FmtLib, insert_neg_int)
TEST(FmtLib, insert_bigint) TEST(FmtLib, insert_bigint)
{ {
#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG) #if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG)
const bigint val = 9945234592L; constexpr bigint val = 9945234592L;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word 9945234592")); ASSERT_THAT(text, Eq("word 9945234592"));
#else #else
@ -57,7 +57,7 @@ TEST(FmtLib, insert_bigint)
TEST(FmtLib, insert_neg_bigint) TEST(FmtLib, insert_neg_bigint)
{ {
#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG) #if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG)
const bigint val = -9945234592L; constexpr bigint val = -9945234592L;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word -9945234592")); ASSERT_THAT(text, Eq("word -9945234592"));
#else #else
@ -68,7 +68,7 @@ TEST(FmtLib, insert_neg_bigint)
TEST(FmtLib, insert_tagint) TEST(FmtLib, insert_tagint)
{ {
#if defined(LAMMPS_BIGBIG) #if defined(LAMMPS_BIGBIG)
const tagint val = 9945234592L; constexpr tagint val = 9945234592L;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word 9945234592")); ASSERT_THAT(text, Eq("word 9945234592"));
#else #else
@ -79,7 +79,7 @@ TEST(FmtLib, insert_tagint)
TEST(FmtLib, insert_neg_tagint) TEST(FmtLib, insert_neg_tagint)
{ {
#if defined(LAMMPS_BIGBIG) #if defined(LAMMPS_BIGBIG)
const tagint val = -9945234592L; constexpr tagint val = -9945234592L;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word -9945234592")); ASSERT_THAT(text, Eq("word -9945234592"));
#else #else
@ -90,7 +90,7 @@ TEST(FmtLib, insert_neg_tagint)
TEST(FmtLib, insert_imageint) TEST(FmtLib, insert_imageint)
{ {
#if defined(LAMMPS_BIGBIG) #if defined(LAMMPS_BIGBIG)
const imageint val = 9945234592L; constexpr imageint val = 9945234592L;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word 9945234592")); ASSERT_THAT(text, Eq("word 9945234592"));
#else #else
@ -101,7 +101,7 @@ TEST(FmtLib, insert_imageint)
TEST(FmtLib, insert_neg_imageint) TEST(FmtLib, insert_neg_imageint)
{ {
#if defined(LAMMPS_BIGBIG) #if defined(LAMMPS_BIGBIG)
const imageint val = -9945234592L; constexpr imageint val = -9945234592L;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word -9945234592")); ASSERT_THAT(text, Eq("word -9945234592"));
#else #else
@ -111,14 +111,14 @@ TEST(FmtLib, insert_neg_imageint)
TEST(FmtLib, insert_double) TEST(FmtLib, insert_double)
{ {
const double val = 1.5; constexpr double val = 1.5;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word 1.5")); ASSERT_THAT(text, Eq("word 1.5"));
} }
TEST(FmtLib, insert_neg_double) TEST(FmtLib, insert_neg_double)
{ {
const double val = -1.5; constexpr double val = -1.5;
auto text = fmt::format("word {}", val); auto text = fmt::format("word {}", val);
ASSERT_THAT(text, Eq("word -1.5")); ASSERT_THAT(text, Eq("word -1.5"));
} }
@ -126,11 +126,23 @@ TEST(FmtLib, insert_neg_double)
TEST(FmtLib, int_for_double) TEST(FmtLib, int_for_double)
{ {
const double val = -1.5; const double val = -1.5;
ASSERT_THROW(auto text = fmt::format("word {:d}", val), std::exception); ASSERT_THROW(auto text = fmt::format(fmt::runtime("word {:d}"), val), std::exception);
} }
TEST(FmtLib, double_for_int) TEST(FmtLib, double_for_int)
{ {
const int val = 15; const int val = 15;
ASSERT_THROW(auto text = fmt::format("word {:g}", val), std::exception); ASSERT_THROW(auto text = fmt::format(fmt::runtime("word {:g}"), val), std::exception);
}
TEST(FmtLib, double_for_string)
{
ASSERT_THROW(auto text = fmt::format(fmt::runtime("word {:g}"), "I am not a number"),
std::exception);
}
TEST(FmtLib, int_for_string)
{
ASSERT_THROW(auto text = fmt::format(fmt::runtime("word {:d}"), "I am not a number"),
std::exception);
} }