Merge pull request lammps#4241 from gsalkuin/develop

Add new fix to compute force and torque due to electric potential
ohenrich · Jan 15, 2025 · ffa4765 · ffa4765
2 parents bbb7d86 + 89370ef
commit ffa4765
Show file tree

Hide file tree

Showing 12 changed files with 1,021 additions and 3 deletions.
diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst
@@ -58,6 +58,7 @@ OPT.
    * :doc:`dt/reset (k) <fix_dt_reset>`
    * :doc:`edpd/source <fix_dpd_source>`
    * :doc:`efield (k) <fix_efield>`
+   * :doc:`efield/lepton <fix_efield_lepton>`
    * :doc:`efield/tip4p <fix_efield>`
    * :doc:`ehex <fix_ehex>`
    * :doc:`electrode/conp (i) <fix_electrode>`

diff --git a/doc/src/fix.rst b/doc/src/fix.rst
@@ -237,6 +237,7 @@ accelerated styles exist.
 * :doc:`dt/reset <fix_dt_reset>` - reset the timestep based on velocity, forces
 * :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
 * :doc:`efield <fix_efield>` - impose electric field on system
+* :doc:`efield/lepton <fix_efield_lepton>` - impose electric field on system using a Lepton expression for the potential
 * :doc:`efield/tip4p <fix_efield>` - impose electric field on system with TIP4P molecules
 * :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm
 * :doc:`electrode/conp <fix_electrode>` - impose electric potential

diff --git a/doc/src/fix_efield.rst b/doc/src/fix_efield.rst
@@ -45,8 +45,9 @@ Description
 
 Add a force :math:`\vec{F} = q\vec{E}` to each charged atom in the group due to an
 external electric field being applied to the system.  If the system
-contains point-dipoles, also add a torque on the dipoles due to the
-external electric field.
+contains point-dipoles, also add a torque :math:`\vec{T} = \vec{p} \times \vec{E}` on the dipoles due to the
+external electric field. This fix does not compute the dipole force :math:`\vec{F} = (\vec{p} \cdot \nabla) \vec{E}`,
+and the :doc:`fix efield/lepton <fix_efield_lepton>` command should be used instead.
 
 .. versionadded:: 28Mar2023
 
@@ -68,6 +69,7 @@ For point-dipoles, equal-style variables can be used, but atom-style
 variables are not currently supported, since they imply a spatial
 gradient in the electric field which means additional terms with
 gradients of the field are required for the force and torque on dipoles.
+The :doc:`fix efield/lepton <fix_efield_lepton>` command should be used instead.
 
 Equal-style variables can specify formulas with various mathematical
 functions, and include :doc:`thermo_style <thermo_style>` command
@@ -229,7 +231,7 @@ Fix style *efield/tip4p* can only be used with tip4p pair styles.
 Related commands
 """"""""""""""""
 
-:doc:`fix addforce <fix_addforce>`
+:doc:`fix addforce <fix_addforce>`, :doc:`fix efield/lepton <fix_efield_lepton>`
 
 Default
 """""""

diff --git a/doc/src/fix_efield_lepton.rst b/doc/src/fix_efield_lepton.rst
@@ -0,0 +1,143 @@
+.. index:: fix efield/lepton
+
+fix efield/lepton command
+=========================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   fix ID group-ID efield/lepton V ...
+
+* ID, group-ID are documented in the :doc:`fix <fix>` command
+* style = *efield/lepton*
+* V = electric potential (electric field * distance units)
+* V must be a Lepton expression (see below)
+* zero or more keyword/value pairs may be appended to args
+* keyword = *region* or *step*
+
+  .. parsed-literal::
+
+       *region* value = region-ID
+         region-ID = ID of region atoms must be in to have effect
+       *step* value = h
+         h = step size for numerical differentiation (distance units)
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   fix ex all efield/lepton "-E*x; E=1"
+   fix dexx all efield/lepton "-0.5*x^2" step 1
+   fix yukawa all efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A=1; B=1" step 1e-6
+   fix infp all efield/lepton "-abs(x)" step 1
+
+   variable th equal 2*PI*ramp(0,1)
+   fix erot all efield/lepton "-(x*cos(v_th)+y*sin(v_th))"
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Add an electric potential :math:`V` that applies to a group of charged atoms a force :math:`\vec{F} = q \vec{E}`,
+and to dipoles a force :math:`\vec{F} = (\vec{p} \cdot \nabla) \vec{E}` and torque :math:`\vec{T} = \vec{p} \times \vec{E}`,
+where :math:`\vec{E} = - \nabla V`. The fix also evaluates the electrostatic energy (:math:`U_{q} = q V` and :math:`U_{p} = - \vec{p} \cdot \vec{E}`)
+due to this potential when the :doc:`fix_modify energy yes <fix_modify>` command is specified (see below).
+
+.. note::
+
+   This command should be used instead of :doc:`fix efield <fix_efield>` if you want to impose a non-uniform electric field on a system with dipoles
+   since the latter does not include the dipole force term. If you only have charges or if the electric field gradient is negligible,
+   :doc:`fix efield <fix_efield>` should be used since it is faster.
+
+The `Lepton library <https://simtk.org/projects/lepton>`_, that the *efield/lepton* fix style interfaces with, evaluates
+the expression string at run time to compute the energy, forces, and torques. It creates an analytical representation
+of :math:`V` and :math:`\vec{E}`, while the gradient force is computed using a central difference scheme
+
+.. math::
+
+   \vec{F} = \frac{|\vec{p}|}{2h} \left[ \vec{E}(\vec{x} + h \hat{p}) - \vec{E}(\vec{x} - h \hat{p}) \right] .
+
+The Lepton expression must be either enclosed in quotes or must not contain any whitespace so that LAMMPS
+recognizes it as a single keyword. More on valid Lepton expressions below. The final Lepton expression must
+be a function of only :math:`x, y, z`, which refer to the current *unwrapped* coordinates of the atoms to ensure continuity.
+Special care must be taken when using this fix with periodic boundary conditions or box-changing commands.
+
+----------
+
+.. include:: lepton_expression.rst
+
+----------
+
+If the *region* keyword is used, the atom must also be in the specified
+geometric :doc:`region <region>` in order to be affected by the potential.
+
+The *step* keyword is required when :doc:`atom_style dipole <atom_style>` is used and the electric field is non-uniform.
+
+----------
+
+Restart, fix_modify, output, run start/stop, minimize info
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+No information about this fix is written to :doc:`binary restart files
+<restart>`.
+
+The :doc:`fix_modify <fix_modify>` *energy* option is supported by this
+fix to add the potential energy defined above to the global potential energy
+of the system as part of :doc:`thermodynamic output <thermo_style>`.
+The default setting for this fix is :doc:`fix_modify energy no <fix_modify>`.
+
+The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
+fix to add the contribution due to the added ***forces*** on charges and dipoles
+to both the global pressure and per-atom stress of the system via the
+:doc:`compute pressure <compute_pressure>` and :doc:`compute stress/atom
+<compute_stress_atom>` commands. The former can be accessed by
+:doc:`thermodynamic output <thermo_style>`. The default setting for
+this fix is :doc:`fix_modify virial no <fix_modify>`.
+
+The :doc:`fix_modify <fix_modify>` *respa* option is supported by this
+fix. This allows to set at which level of the :doc:`r-RESPA <run_style>`
+integrator the fix adding its forces. Default is the outermost level.
+
+This fix computes a global scalar and a global 3-vector of forces,
+which can be accessed by various :doc:`output commands <Howto_output>`.
+The scalar is the potential energy discussed above.
+The vector is the total force added to the group of atoms.
+The scalar and vector values calculated by this fix are "extensive".
+
+This fix cannot be used with the *start/stop* keywords of
+the :doc:`run <run>` command.
+
+The forces due to this fix are imposed during an energy minimization,
+invoked by the :doc:`minimize <minimize>` command. You should not
+specify force components with a variable that has time-dependence for
+use with a minimizer, since the minimizer increments the timestep as
+the iteration count during the minimization.
+
+.. note::
+
+   If you want the electric potential energy to be included in the
+   total potential energy of the system (the quantity being minimized),
+   you MUST enable the :doc:`fix_modify <fix_modify>` *energy* option for this fix.
+
+----------
+
+Restrictions
+""""""""""""
+
+Fix style *efield/lepton* is part of the LEPTON package. It is only enabled if LAMMPS was built with that package.
+See the :doc:`Build package <Build_package>` page for more info.
+
+
+Related commands
+""""""""""""""""
+
+:doc:`fix efield <fix_efield>`
+
+Default
+"""""""
+
+none
diff --git a/examples/LEPTON/in.efield-lepton b/examples/LEPTON/in.efield-lepton
@@ -0,0 +1,50 @@
+# Point dipoles in a 3d box with an external potential (ignoring dipolar interactions)
+
+units           lj
+atom_style      hybrid sphere dipole
+dimension       3
+boundary        s s s
+region          box block -2 2 -2 2 -2 2
+
+create_box      1 box
+create_atoms    1 random 100 12345 NULL
+
+# need both mass settings due to hybrid atom style
+mass            1 1.0
+set             group all mass 1.0
+set             group all diameter 0.1
+
+set             group all dipole/random 98934 0.01
+pair_style      none
+comm_modify     cutoff 3.0
+
+velocity        all create 0.0 87287 mom yes rot yes
+
+fix             1 all nve/sphere update dipole
+
+###############################################################################################################
+## Yukawa potential
+#fix            2 all efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A = 0.1; B = 5" step 1e-8
+
+## Gradually increasing uniform field 
+#variable       E equal ramp(0,1)
+#fix            2 all efield/lepton "-v_E*(x+y+z)"
+
+## Linear gradient field
+fix             2 all efield/lepton "-1/6*x^3" step 1e-6
+
+fix_modify      2 energy yes
+
+###############################################################################################################
+
+timestep        1e-3
+
+compute         erot all erotate/sphere
+variable        etotal equal "ke + c_erot + pe" # thermo etotal doesn't include erot
+thermo_style    custom step temp ke c_erot pe v_etotal
+thermo          500
+thermo_modify   norm no
+
+#dump           1 all custom 500 dump.dipole id x y z diameter mux muy muz fx fy fz tqx tqy tqz
+
+run             10000
diff --git a/examples/LEPTON/log.13Jan25.efield-lepton.g++.1 b/examples/LEPTON/log.13Jan25.efield-lepton.g++.1
@@ -0,0 +1,115 @@
+LAMMPS (19 Nov 2024 - Development - patch_19Nov2024-283-g742c869534-modified)
+  using 1 OpenMP thread(s) per MPI task
+# Point dipoles in a 3d box with an external potential (ignoring dipolar interactions)
+
+units           lj
+atom_style      hybrid sphere dipole
+WARNING: Atom style hybrid defines both, per-type and per-atom masses; both must be set, but only per-atom masses will be used (src/atom_vec_hybrid.cpp:132)
+dimension       3
+boundary        s s s
+region          box block -2 2 -2 2 -2 2
+
+create_box      1 box
+Created orthogonal box = (-2 -2 -2) to (2 2 2)
+  1 by 1 by 1 MPI processor grid
+create_atoms    1 random 100 12345 NULL
+Created 100 atoms
+  using lattice units in orthogonal box = (-2.0004 -2.0004 -2.0004) to (2.0004 2.0004 2.0004)
+  create_atoms CPU = 0.000 seconds
+
+# need both mass settings due to hybrid atom style
+mass            1 1.0
+set             group all mass 1.0
+Setting atom values ...
+  100 settings made for mass
+set             group all diameter 0.1
+Setting atom values ...
+  100 settings made for diameter
+
+set             group all dipole/random 98934 0.01
+Setting atom values ...
+  100 settings made for dipole/random
+pair_style      none
+comm_modify     cutoff 3.0
+
+velocity        all create 0.0 87287 mom yes rot yes
+
+fix             1 all nve/sphere update dipole
+
+###############################################################################################################
+## Yukawa potential
+#fix            2 all efield/lepton "A*exp(-B*r)/r; r=abs(sqrt(x^2+y^2+z^2)); A = 0.1; B = 5" step 1e-8
+
+## Gradually increasing uniform field
+#variable       E equal ramp(0,1)
+#fix            2 all efield/lepton "-v_E*(x+y+z)"
+
+## Linear gradient field
+fix             2 all efield/lepton "-1/6*x^3" step 1e-6
+
+fix_modify      2 energy yes
+
+###############################################################################################################
+
+timestep        1e-3
+
+compute         erot all erotate/sphere
+variable        etotal equal "ke + c_erot + pe" # thermo etotal doesn't include erot
+thermo_style    custom step temp ke c_erot pe v_etotal
+thermo          500
+thermo_modify   norm no
+
+#dump           1 all custom 500 dump.dipole id x y z diameter mux muy muz fx fy fz tqx tqy tqz
+
+run             10000
+WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (src/atom.cpp:2442)
+Per MPI rank memory allocation (min/avg/max) = 4.273 | 4.273 | 4.273 Mbytes
+   Step          Temp          KinEng         c_erot         PotEng        v_etotal   
+         0   0              0              0              0.036419797    0.036419797  
+       500   3.7159175e-06  0.00055181374  0.44262618    -0.40675701     0.036420985  
+      1000   1.2808438e-05  0.0019020531   0.24499116    -0.21047295     0.036420259  
+      1500   2.8343769e-05  0.0042090498   0.26504485    -0.2328336      0.036420307  
+      2000   4.8796894e-05  0.0072463388   0.30953526    -0.28036098     0.036420618  
+      2500   7.8933715e-05  0.011721657    0.2015076     -0.17680909     0.036420173  
+      3000   0.00011381678  0.016901791    0.31002163    -0.29050294     0.036420476  
+      3500   0.00015650339  0.023240753    0.27837968    -0.26520001     0.036420418  
+      4000   0.00020429109  0.030337227    0.26201101    -0.25592795     0.036420289  
+      4500   0.00026362339  0.039148074    0.29769952    -0.3004271      0.036420499  
+      5000   0.00033328941  0.049493478    0.21642442    -0.22949776     0.036420131  
+      5500   0.00040914224  0.060757622    0.28422322    -0.30856047     0.036420377  
+      6000   0.00049425119  0.073396302    0.31767       -0.35464572     0.03642058   
+      6500   0.00058508892  0.086885704    0.29079532    -0.34126075     0.036420276  
+      7000   0.00069845073  0.10371993     0.25776048    -0.32506015     0.036420262  
+      7500   0.0008215656   0.12200249     0.27033777    -0.35591972     0.036420539  
+      8000   0.00095528125  0.14185927     0.33943527    -0.44487406     0.036420479  
+      8500   0.0011052502   0.16412965     0.26727165    -0.39498109     0.036420218  
+      9000   0.0012738298   0.18916373     0.31082058    -0.46356382     0.036420485  
+      9500   0.001464197    0.21743325     0.25669856    -0.43771158     0.036420224  
+     10000   0.0016627654   0.24692067     0.36273185    -0.57323194     0.036420578  
+Loop time of 0.84714 on 1 procs for 10000 steps with 100 atoms
+
+Performance: 1019901.911 tau/day, 11804.420 timesteps/s, 1.180 Matom-step/s
+62.3% CPU use with 1 MPI tasks x 1 OpenMP threads
+
+MPI task timing breakdown:
+Section |  min time  |  avg time  |  max time  |%varavg| %total
+---------------------------------------------------------------
+Pair    | 0          | 0          | 0          |   0.0 |  0.00
+Neigh   | 9.21e-07   | 9.21e-07   | 9.21e-07   |   0.0 |  0.00
+Comm    | 0.00094138 | 0.00094138 | 0.00094138 |   0.0 |  0.11
+Output  | 0.0001983  | 0.0001983  | 0.0001983  |   0.0 |  0.02
+Modify  | 0.84105    | 0.84105    | 0.84105    |   0.0 | 99.28
+Other   |            | 0.004946   |            |       |  0.58
+
+Nlocal:            100 ave         100 max         100 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Nghost:              0 ave           0 max           0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+Neighs:              0 ave           0 max           0 min
+Histogram: 1 0 0 0 0 0 0 0 0 0
+
+Total # of neighbors = 0
+Ave neighs/atom = 0
+Neighbor list builds = 5
+Dangerous builds = 0
+Total wall time: 0:00:00