[pull] develop from lammps:develop

doc/src/fix_nonaffine_displacement.rst

@@ -112,14 +112,16 @@ The reference state is saved to :doc:`binary restart files <restart>`.
 None of the :doc:`fix_modify <fix_modify>` options are relevant to this
 fix.
 
-This fix computes a peratom array with 3 columns, which can be accessed
-by indices 1-3 using any command that uses per-atom values from a fix
-as input.
+This fix computes a peratom array with either 3 or 9 columns, which can
+be accessed by indices 1-9 using any command that uses per-atom values
+from a fix as input.
 
 For the *integrated* style, the three columns are the nonaffine
 displacements in the x, y, and z directions. For the *d2min* style,
-the three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`, the
-volumetric strain, and the deviatoric strain.
+the first three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`,
+the volumetric strain, and the deviatoric strain. The following 6
+columns are the xx, yy, zz, xy, xz, and yz components of the calculated
+strain tensor.
 
 Restrictions
 """"""""""""

doc/src/fix_qeq.rst

@@ -42,7 +42,9 @@ Syntax
 
   .. parsed-literal::
 
-       *alpha* value = Slater type orbital exponent (qeq/slater only)
+       *alpha* value = Slater type orbital exponent (qeq/slater only). Can be followed by optional arguments:
+         *wolf* value = width of taper to terminate Coulomb integrals for the Wolf summation (default value is zero)
+         *dsf* value = width of taper to terminate Coulomb integrals for the Fennell-Gezelter summation (default value is zero)
        *cdamp* value = damping parameter for Coulomb interactions (qeq/ctip only)
        *maxrepeat* value = number of equilibration cycles allowed to ensure no atoms cross charge bounds (qeq/ctip only)
        *qdamp* value = damping factor for damped dynamics charge solver (qeq/dynamic and qeq/fire only)
@@ -57,6 +59,10 @@ Examples
    fix 1 all qeq/point 1 10 1.0e-6 200 param.qeq1
    fix 1 qeq qeq/shielded 1 8 1.0e-6 100 param.qeq2
    fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2
+   fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 wolf
+   fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 wolf 2.0
+   fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 dsf
+   fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 dsf 2.0
    fix 1 all qeq/ctip 1 12 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10
    fix 1 qeq qeq/dynamic 1 12 1.0e-3 100 my_qeq
    fix 1 all qeq/fire 1 10 1.0e-3 100 my_qeq qdamp 0.2 qstep 0.1
@@ -249,6 +255,24 @@ larger sizes, and *qeq/fire* is faster than *qeq/dynamic*\ .
    arbitrary choices of these parameters.  We do not develop these QEq
    parameters.  See the examples/qeq directory for some examples.
 
+.. versionadded:: TBD
+
+In previous versions of LAMMPS, the real-space summations of Coulomb
+interactions were done by replacing *1/r* using a damped potential
+*erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of
+decay. However, any finite value of *alpha* leads to a jump at the
+cutoff, which interferes with equilibration if atoms move across the
+cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) <Wolf4>`
+(*wolf*) and its extension by :ref:`(Fennell and Gezelter) <Fennell3>`
+(*dsf*) solve this problem. An extension was implemented to specify the
+width of taper (see :ref:`(Mei et al.) <Mei2>`) to smoothly terminate the
+Coulomb integrals at the cutoff. This is done by specifying the optional
+arguments *wolf* and *dsf* with the value representing the width of
+taper that smoothly terminates the Coulomb integrals. For example, if
+the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are
+smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For
+backward compatibility, the default taper width is zero.
+
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
@@ -314,3 +338,15 @@ Physical Chemistry, 105, 9396-9049 (2001)
 .. _Shan:
 
 **(QEq/Fire)** T.-R. Shan, A. P. Thompson, S. J. Plimpton, in preparation
+
+.. _Wolf4:
+
+**(Wolf)** D. Wolf, P. Keblinski, S. R. Phillpot, J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999).
+
+.. _Fennell3:
+
+**(Fennell)** J. Fennell, J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006).
+
+.. _Mei2:
+
+**(Mei)** J. Mei, J. W. Davenport, G. W. Fernando, Phys. Rev. B 43, 4653 (1991).

doc/src/fix_rheo.rst

@@ -34,10 +34,10 @@ Syntax
          optional args = *exclude/type* or *scale/cross/type*
            *exclude/type* values = *types*
              *types* = list of types
-           *scale/cross/type* values = *shiftscale* *cmin* *wmin*
+           *scale/cross/type* values = *shiftscale* *cmin* *rmin*
              *shiftscale* = fraction of shifting in normal direction to preserve (unitless)
              *cmin* = minimum color function value required for scaling (unitless)
-             *wmin* = minimum local same-type support required for any shifting (unitless)
+             *rmin* = minimum local same-type weighted distance required for any shifting (unitless)
        *rho/sum* density evolution performed by a kernel summation
          values = none
          optional args = *self/mass*
@@ -123,14 +123,14 @@ zero implies there is no shifting in the normal direction and a value of
 *scaleshift* of one implies no change in behavior. This scaling is only applied
 to atoms with a color function value greater than *cmin*. To handle scenarios
 of a small inclusion of one fluid type (e.g. a single atom) inside another,
-the degree of same-type support is calculated
+a weighted distance of same-type support is calculated
 
 .. math::
-   W_{i,\mathrm{same}} = \sum_{j} W_{ij} \delta_{ij}
+   R_{i,\mathrm{same}} = \sum_{j} r_{ij} W_{ij} \delta_{ij}
 
 where :math:`\delta_{ij}` is zero if atoms :math:`i` and :math:`j` have different
-types but unity otherwise. If :math:`W_{i,\mathrm{same}}` is ever less than the
-specified value of *wmin*, shifting is turned off for particle :math:`i`
+types but unity otherwise. If :math:`R_{i,\mathrm{same}}` is ever less than the
+specified value of *rmin*, shifting is turned off for particle :math:`i`
 
 In systems with free surfaces (atom-vacuum), the *surface/detection* keyword
 can classify the location of particles as being within the bulk fluid, on a

doc/src/pair_coul.rst

@@ -146,6 +146,12 @@ Examples
    pair_style coul/streitz 12.0 wolf 0.30
    pair_coeff * * AlO.streitz Al O
 
+   pair_style coul/streitz 12.0 wolf 0.3         # default taper width is zero
+   pair_style coul/streitz 12.0 wolf 0.3 2.0
+   pair_style coul/streitz 12.0 dsf 0.3          # default taper width is zero
+   pair_style coul/streitz 12.0 dsf 0.3 2.0
+   pair_coeff * * GaN.streitz Ga N
+
    pair_style tip4p/cut 1 2 7 8 0.15 12.0
    pair_coeff * *
 
@@ -302,17 +308,49 @@ command doc page.  Alternatively *qfile* can be replaced by
 "coul/streitz", in which case the fix will extract QEq parameters from
 the coul/streitz pair style itself.
 
-See the examples/strietz directory for an example input script that
+See the ``examples/streitz`` directory for an example input script that
 uses the Streitz-Mintmire potential.  The potentials directory has the
-AlO.eam.alloy and AlO.streitz potential files used by the example.
+``AlO.eam.alloy`` and ``AlO.streitz`` potential files used by the example.
 
-Note that the Streiz-Mintmire potential is generally used for oxides,
+Note that the Streitz-Mintmire potential is generally used for oxides,
 but there is no conceptual problem with extending it to nitrides and
-carbides (such as SiC, TiN).  Pair coul/strietz used by itself or with
+carbides (such as SiC, TiN).  Pair coul/streitz used by itself or with
 any other pair style such as EAM, MEAM, Tersoff, or LJ in
 hybrid/overlay mode.  To do this, you would need to provide a
 Streitz-Mintmire parameterization for the material being modeled.
 
+.. versionchanged:: TBD
+
+In previous versions of LAMMPS, the real-space summations of Coulomb
+interactions were done by replacing *1/r* using a damped potential
+*erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of
+decay. However, any finite value of *alpha* leads to a jump at the
+cutoff, which interferes with equilibration if atoms move across the
+cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) <Wolf1>`
+(*wolf*) and its extension by :ref:`(Fennell and Gezelter) <Fennell1>`
+(*dsf*) solve this problem. An extension was implemented to specify the
+width of taper (see :ref:`(Mei et al.) <Mei1>`) to smoothly terminate
+the Coulomb integrals at the cutoff. This is done by specifying the
+optional arguments *wolf* and *dsf* with the value representing the
+width of taper that smoothly terminates the Coulomb integrals. For
+example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb
+integrals are smoothly rescaled from their actual value at r=6 A to zero
+at r=8 A. For backward compatibility, the default taper width is zero.
+
+An implementation of the Streitz-Mintmire potential for GaN due to
+:ref:`(Groger and Fikar) <Groger1>` can be found in the examples/streitz
+directory. The electrostatic parameters of Ga and N are stored in file
+GaN.streitz and the short-range tersoff/mod potential in the file
+GaN.streitz+tersoff.mod. The total potential must be specified as:
+
+.. code-block:: LAMMPS
+
+   pair_style hybrid/overlay tersoff/mod coul/streitz 12.0 dsf 0.3 2.0
+   pair_coeff * * tersoff/mod GaN.streitz+tersoff.mod Ga N
+   pair_coeff * * coul/streitz GaN.streitz Ga N
+
+where the last three parameters specify the method of calculation of Coulomb interactions (*ewald*, *wolf* or *dsf*), the value of *alpha*, and the width of taper to smoothly terminate the Coulomb integrals.
+
 ----------
 
 Pair style *coul/exclude* computes Coulombic interactions like *coul/cut*
@@ -477,3 +515,11 @@ J Chemical Physics, 162, 054709 (2025).
 
 **(Jorgensen)** Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
 Phys, 79, 926 (1983).
+
+.. _Mei1:
+
+**(Mei)** J. Mei, J. W. Davenport, G. W. Fernando, Phys. Rev. B 43, 4653 (1991).
+
+.. _Groger1:
+
+**(Groger)** R. Groger, J. Fikar, Acta Mater. (2026) in review.

doc/src/set.rst

@@ -27,12 +27,12 @@ Syntax
   or *density* or *density/disc* or *diameter* or *dihedral* or *dipole*
   or *dipole/random* or *dpd/theta* or *edpd/cv* or *edpd/temp* or
   *epsilon* or *image* or *improper* or *length* or *mass* or *mol* or
-  *omega* or *quat* or *quat/random* or *radius/electron* or *shape* or
-  *smd/contact/radius* or *smd/mass/density* or *sph/cv* or *sph/e* or
-  *sph/rho* or *spin/atom* or *spin/atom/random* or *spin/electron* or
-  *temperature* or *theta* or *theta/random* or *tri* or *type* or
-  *type/fraction* or *type/ratio* or *type/subset* or *volume* or *vx*
-  or *vy* or *vz* or *x* or *y* or *z* or *i_name* or *d_name* or
+  *omega* or *quat* or *quat/random* or *radius/electron* or *rheo/rho* or
+  *rheo/status* or *shape* or *smd/contact/radius* or *smd/mass/density* or
+  *sph/cv* or *sph/e* or *sph/rho* or *spin/atom* or *spin/atom/random* or
+  *spin/electron* or *temperature* or *theta* or *theta/random* or *tri* or
+  *type* or *type/fraction* or *type/ratio* or *type/subset* or *volume* or
+  *vx* or *vy* or *vz* or *x* or *y* or *z* or *i_name* or *d_name* or
   *i2_name* or *d2_name*
 
   .. parsed-literal::
@@ -97,6 +97,10 @@ Syntax
        *radius/electron* value = eradius
          eradius = electron radius (or fixed-core radius) (distance units)
          value can be an atom-style variable (see below)
+       *rheo/rho* value = density of RHEO particles (mass/distance\^3)
+         value can be an atom-style variable (see below)
+       *rheo/status* value = status or phase of RHEO particles (unitless)
+         value can be an atom-style variable (see below)
        *shape* values = Sx Sy Sz
          Sx,Sy,Sz = 3 diameters of ellipsoid (distance units)
          any of Sx,Sy,Sz can be an atom-style variable (see below)
@@ -512,6 +516,10 @@ use of an atom-style variable.
 Keyword *radius/electron* uses the specified value to set the radius
 of electrons or fixed cores.
 
+Keywords *rheo/rho* and *rheo/status* set the density and the status of
+rheo particles. In particular, one can only set the phase in the status
+as described by the :doc:`RHEO howto page <Howto_rheo>`.
+
 Keyword *shape* sets the size and shape of the selected atoms.  The
 particles must be ellipsoids as defined by the :doc:`atom_style
 ellipsoid <atom_style>` command.  The *Sx*, *Sy*, *Sz* settings are

doc/utils/sphinx-config/false_positives.txt

@@ -1236,6 +1236,7 @@ Fickian
 fieldname
 figshare
 Fij
+Fikar
 filelink
 filename
 Filename
@@ -1464,6 +1465,7 @@ Grigera
 Grimme
 grmask
 Grmask
+Groger
 gromacs
 Gromacs
 gromos
@@ -3788,10 +3790,8 @@ Straub
 strcmp
 streitz
 Streitz
-Streiz
 strerror
 strided
-strietz
 stringstreams
 strmatch
 strncmp

examples/streitz/GaN.streitz

@@ -0,0 +1 @@
+../../potentials/GaN.streitz

examples/streitz/GaN.streitz+tersoff.mod

@@ -0,0 +1 @@
+../../potentials/GaN.streitz+tersoff.mod

examples/streitz/in.gan

@@ -0,0 +1,73 @@
+variable a0 equal 3.0  # trial values
+variable c0 equal 5.0
+variable u0 equal 0.38
+variable Efu_inf equal -0.9183986628325442  # energy of a pair of isolated Ga and N ions
+
+units metal
+dimension 3
+boundary p p p
+atom_style charge
+atom_modify map yes
+
+lattice custom ${a0} &
+   a1  1.0 0.0               0.0 &
+   a2 -0.5 0.866025403784439 0.0 &
+   a3  0.0 0.0               $(v_c0/v_a0) &
+   spacing 1 0.866025403784439 $(v_c0/v_a0) &
+   basis 0.33333333 0.66666667 0 &
+   basis 0.66666667 0.33333333 0.5 &
+   basis 0.33333333 0.66666667 ${u0} &
+   basis 0.66666667 0.33333333 $(v_u0+0.5)
+region box prism 0 1 0 1 0 1 -0.5 0 0 units lattice
+create_box 2 box
+create_atoms 2 box basis 1 1 basis 2 1 basis 3 2 basis 4 2
+mass * 1.0
+
+set type 1 charge 1.0
+set type 2 charge -1.0
+group g1 type 1
+group g2 type 2
+
+compute eat all pe/atom
+compute q1all g1 property/atom q
+compute q1 g1 reduce ave c_q1all
+
+pair_style hybrid/overlay tersoff/mod coul/streitz 12.0 dsf 0.3 2.0
+pair_coeff * * tersoff/mod GaN.streitz+tersoff.mod Ga N
+pair_coeff * * coul/streitz GaN.streitz Ga N
+
+fix fq all qeq/slater 1 12.0 1e-6 1000 coul/streitz alpha 0.3 dsf 2.0
+
+fix 1 all box/relax aniso 0.0
+
+#dump 1 all cfg 100 dump_*.cfg mass type xs ys zs c_eat q
+#dump_modify 1 element Ga N
+
+thermo 100
+thermo_style custom step pe fmax press vol lx ly lz
+include scmin.mod  # self-consistent minimization & charge equilibration
+
+variable Efu_nocorr equal 2*pe/count(all)
+variable Efu equal ${Efu_nocorr}-${Efu_inf}
+variable a equal (lx+ly*2/sqrt(3))/2
+variable c equal lz
+variable p equal press/1e4
+variable q1 equal c_q1
+
+print '--------------------------------------------------------------------------------------------'
+print 'Energy (no correction): E = ${Efu_nocorr} eV/fu'
+print 'Energy: E-Einf = ${Efu} eV/fu'
+print 'Pressure: p = $p GPa'
+print 'Energy of a pair of isolated ions: Einf = ${Efu_inf} eV/fu'
+print 'Lattice parameters: a = $a A'
+print '                    c = $c A'
+print 'Cationic charge: q = ${q1}*e'
+print '--------------------------------------------------------------------------------------------'
+
+#include polariz.mod
+#include polariz_diff.mod
+#include polariz_wz.mod  # no relaxation, fixed cell shape
+#include polariz_wz_relax.mod  # with relaxation - transforms into different structures?
+
+
+

examples/streitz/scmin.mod

@@ -0,0 +1,37 @@
+# QEq cannot be run together with minimize, because each new calculation of charges alters the energy
+# hypersurface on which we are looking for a minimum. Minimization frequently gets stuck with the
+# message "linesearch alpha is zero" before reaching zero pressure.
+#
+# Below is a self-consistency loop that minimizes the potential energy while distributing atomic
+# charges, which replaces the minimize command. It works as follows:
+#   1. For given atomic positions, run QEq to compute charges. The initial charges are q0 and the
+#      final ones q.
+#   2. Fix charges (q) and run energy minimization with changing box shape to get new atomic
+#      coordinates. This gives zero pressure.
+#   3. If max|q-q0| >= 0.01, set q0 := q and go to 1.
+#      If max|q-q0| < 0.01, we have equilibrium atomic positions and charges => exit.
+#
+# Local variables end with _ to avoid interference with the main script.
+
+variable dq_ atom abs(q-f_q0_)
+compute dqmax_ all reduce max v_dq_
+variable dqmax_ equal c_dqmax_
+
+variable iter loop 100   # number of iterations to equilibrate charges & get zero pressure
+label self_consistency
+  fix q0_ all store/state 0 q
+  fix fq all qeq/slater 1 12.0 1e-6 1000 coul/streitz alpha 0.3 dsf 2.0
+  run 0	      	      	      	# calculate charges -> nonzero pressure
+
+  unfix fq
+  minimize 0 1e-6 10000 10000 	# fix charges -> zero pressure
+
+  print "***** self-consistent loop: iter=${iter}, dqmax=${dqmax_} *****"
+  if "${dqmax_} < 0.01" then "jump SELF break"
+  next iter
+jump SELF self_consistency
+
+label break
+unfix q0_
+uncompute dqmax_
+variable iter delete

potentials/GaN.streitz

@@ -0,0 +1,3 @@
+# chi (eV), J (eV), gamma (1/ang), zeta (1/ang), Z (e)
+Ga -1.25 14.51 0 1.2333912356832242 0
+N   3.38  6.91 0 0.8509319105529678 0