corrections for documentation of bosonic PIMD fix styles

This commit is contained in:
Axel Kohlmeyer
2025-03-14 09:54:26 -04:00
parent e97807b92e
commit fcb2eee686
4 changed files with 31 additions and 260 deletions

View File

@ -162,8 +162,8 @@ OPT.
* :doc:`phonon <fix_phonon>`
* :doc:`pimd/langevin <fix_pimd>`
* :doc:`pimd/nvt <fix_pimd>`
* :doc:`pimd/langevin/bosonic <fix_pimd_bosonic>`
* :doc:`pimd/nvt/bosonic <fix_pimd_bosonic>`
* :doc:`pimd/langevin/bosonic <fix_pimd>`
* :doc:`pimd/nvt/bosonic <fix_pimd>`
* :doc:`planeforce <fix_planeforce>`
* :doc:`plumed <fix_plumed>`
* :doc:`poems <fix_poems>`

View File

@ -341,8 +341,8 @@ accelerated styles exist.
* :doc:`phonon <fix_phonon>` - calculate dynamical matrix from MD simulations
* :doc:`pimd/langevin <fix_pimd>` - Feynman path-integral molecular dynamics with stochastic thermostat
* :doc:`pimd/nvt <fix_pimd>` - Feynman path-integral molecular dynamics with Nose-Hoover thermostat
* :doc:`pimd/langevin/bosonic <fix_pimd_bosonic>` - Bosonic Feynman path-integral molecular dynamics for with stochastic thermostat
* :doc:`pimd/nvt/bosonic <fix_pimd_bosonic>` - Bosonic Feynman path-integral molecular dynamics with Nose-Hoover thermostat
* :doc:`pimd/langevin/bosonic <fix_pimd>` - Bosonic Feynman path-integral molecular dynamics for with stochastic thermostat
* :doc:`pimd/nvt/bosonic <fix_pimd>` - Bosonic Feynman path-integral molecular dynamics with Nose-Hoover thermostat
* :doc:`planeforce <fix_planeforce>` - constrain atoms to move in a plane
* :doc:`plumed <fix_plumed>` - wrapper on PLUMED free energy library
* :doc:`poems <fix_poems>` - constrain clusters of atoms to move as coupled rigid bodies

View File

@ -9,11 +9,11 @@ fix pimd/langevin command
fix pimd/nvt command
====================
:doc:`fix pimd/langevin/bosonic <fix_pimd_bosonic>` command
===========================================================
fix pimd/langevin/bosonic command
=================================
:doc:`fix pimd/nvt/bosonic <fix_pimd_bosonic>` command
======================================================
fix pimd/nvt/bosonic command
============================
Syntax
""""""
@ -23,7 +23,7 @@ Syntax
fix ID group-ID style keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* style = *pimd/langevin* or *pimd/nvt* = style name of this fix command
* style = *pimd/langevin* or *pimd/nvt* or *pimd/langevin/bosonic* or *pimd/nvt/bosonic* = style name of this fix command
* zero or more keyword/value pairs may be appended
* keywords for style *pimd/nvt*
@ -89,12 +89,13 @@ partition function for the original system to a classical partition
function for a ring-polymer system is exploited, to efficiently sample
configurations from the canonical ensemble :ref:`(Feynman) <Feynman>`.
.. versionadded:: 11Mar2025
.. versionadded:: TBD
Fix *pimd/langevin/bosonic* and *pimd/nvt/bosonic* were added.
Fix *pimd/nvt* and fix *pimd/langevin* simulate *distinguishable* quantum particles.
Simulations of bosons, including exchange effects, are supported with the :doc:`fix pimd/langevin/bosonic <fix_pimd_bosonic>` and :doc:`fix pimd/nvt/bosonic <fix_pimd_bosonic>` commands.
Simulations of bosons, including exchange effects, are supported with the
fix *pimd/langevin/bosonic* and the *pimd/nvt/bosonic* commands.
For distinguishable particles, the isomorphic classical partition function and its components are given
by the following equations:
@ -173,15 +174,17 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics
Mode *pimd* added to fix pimd/langevin.
Fix pimd/langevin supports the *method* values *nmpimd* and *pimd*. The default value is *nmpimd*.
If *method* is *nmpimd*, the normal mode representation is used to integrate the equations of motion.
The exact solution of harmonic oscillator is used to propagate the free ring polymer part of the Hamiltonian.
If *method* is *pimd*, the Cartesian representation is used to integrate the equations of motion.
The harmonic force is added to the total force of the system, and the numerical integrator is used to propagate the Hamiltonian.
Fix pimd/langevin supports the *method* values *nmpimd* and *pimd*. The default
value is *nmpimd*. If *method* is *nmpimd*, the normal mode representation is
used to integrate the equations of motion. The exact solution of harmonic
oscillator is used to propagate the free ring polymer part of the Hamiltonian.
If *method* is *pimd*, the Cartesian representation is used to integrate the
equations of motion. The harmonic force is added to the total force of the
system, and the numerical integrator is used to propagate the Hamiltonian.
The keyword *integrator* specifies the Trotter splitting method used by *fix pimd/langevin*.
See :ref:`(Liu) <Liu>` for a discussion on the OBABO and BAOAB splitting schemes. Typically
either of the two should work fine.
The keyword *integrator* specifies the Trotter splitting method used by *fix
pimd/langevin*. See :ref:`(Liu) <Liu>` for a discussion on the OBABO and BAOAB
splitting schemes. Typically either of the two should work fine.
The keyword *fmass* sets a further scaling factor for the fictitious
masses of beads, which can be used for the Partial Adiabatic CMD
@ -231,8 +234,8 @@ a positive floating-point number.
For pimd simulations, a temperature values should be specified even for nve ensemble. Temperature will make a difference
for nve pimd, since the spring elastic frequency between the beads will be affected by the temperature.
The keyword *thermostat* reads *style* and *seed* of thermostat for fix style *pimd/langevin*. *style* can only
be *PILE_L* (path integral Langevin equation local thermostat, as described in :ref:`Ceriotti <Ceriotti2>`), and *seed* should a positive integer number, which serves as the seed of the pseudo random number generator.
The keyword *thermostat* reads *style* and *seed* of thermostat for fix style *pimd/langevin*.
*style* can only be *PILE_L* (path integral Langevin equation local thermostat, as described in :ref:`Ceriotti <Ceriotti2>`), and *seed* should a positive integer number, which serves as the seed of the pseudo random number generator.
.. note::
@ -418,7 +421,12 @@ LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Fix *pimd/nvt* cannot be used with :doc:`lj units <units>`.
Fix *pimd/langevin* can be used with :doc:`lj units <units>`. See the above part for how to use it.
Fix *pimd/langevin* can be used with :doc:`lj units <units>`.
See the documentation above for how to use it.
Only some combinations of fix styles and their options support
partitions with multiple processors. LAMMPS will stop with an
error if multi-processor partitions are not supported.
A PIMD simulation can be initialized with a single data file read via
the :doc:`read_data <read_data>` command. However, this means all
@ -435,7 +443,7 @@ variable, e.g.
Related commands
""""""""""""""""
:doc:`pimd/nvt/bosonic <fix_pimd_bosonic>`, :doc:`pimd/langevin/bosonic <fix_pimd_bosonic>`
:doc:`fix ipi <fix_ipi>`
Default
"""""""

View File

@ -1,237 +0,0 @@
.. index:: fix pimd/langevin/bosonic
.. index:: fix pimd/nvt/bosonic
fix pimd/langevin/bosonic command
=================================
fix pimd/nvt/bosonic command
============================
Syntax
""""""
.. code-block:: LAMMPS
fix ID group-ID style keyword value ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* style = *pimd/langevin/bosonic* or *pimd/nvt/bosonic* = style name of this fix command
* zero or more keyword/value pairs may be appended
* keywords for style *pimd/nvt/bosonic*
.. parsed-literal::
*keywords* = *method* or *fmass* or *sp* or *temp* or *nhc*
*method* value = *pimd* or *nmpimd*
*fmass* value = scaling factor on mass
*sp* value = scaling factor on Planck constant
*temp* value = temperature (temperature units)
*nhc* value = Nc = number of chains in Nose-Hoover thermostat
* keywords for style *pimd/langevin/bosonic*
.. parsed-literal::
*keywords* = *integrator* or *ensemble* or *fmass* or *temp* or *thermostat* or *tau* or *fixcom* or *lj* or *esych*
*integrator* value = *obabo* or *baoab*
*ensemble* value = *nvt* or *nve*
*fmass* value = scaling factor on mass
*temp* value = temperature (temperature unit)
temperature = target temperature of the thermostat
*thermostat* values = style seed
style value = *PILE_L*
seed = random number generator seed
*tau* value = thermostat damping parameter (time unit)
*fixcom* value = *yes* or *no*
*lj* values = epsilon sigma mass planck mvv2e
epsilon = energy scale for reduced units (energy units)
sigma = length scale for reduced units (length units)
mass = mass scale for reduced units (mass units)
planck = Planck's constant for other unit style
mvv2e = mass * velocity^2 to energy conversion factor for other unit style
*esynch* value = *yes* or *no*
Examples
""""""""
.. code-block:: LAMMPS
fix 1 all pimd/nvt/bosonic method pimd fmass 1.0 sp 1.0 temp 2.0 nhc 4
fix 1 all pimd/langevin/bosonic integrator obabo temp 113.15 thermostat PILE_L 1234 tau 1.0
Example input files are provided in the examples/PACKAGES/pimd_bosonic directory.
Description
"""""""""""
These fix commands are based on the fixes :doc:`pimd/nvt and
pimd/langevin <fix_pimd>` for performing quantum molecular dynamics
simulations based on the Feynman path-integral formalism. The key
difference is that fix *pimd/nvt* and fix *pimd/langevin* simulate
*distinguishable* particles, while fix *pimd/nvt/bosonic* and fix
*pimd/langevin/bosonic* perform simulations of bosons, including exchange
effects. The *bosonic* commands share syntax with the equivalent commands for distinguishable particles.
The user is referred to the documentation of :doc:`these commands <fix_pimd>`
for a detailed syntax description and additional, general capabilities
of the commands. The major differences from fix *pimd/nvt* and fix *pimd/langevin* in terms of
capabilities are:
* Fix *pimd/nvt/bosonic* only supports the "pimd" and "nmpimd" methods. Fix
*pimd/langevin/bosonic* only supports the "pimd" method, which is the default
in this fix. These restrictions are related to the use of normal
modes, which change in bosons. For similar reasons, *fmmode* of
*pimd/langevin* should not be used, and would raise an error if set to
a value other than *physical*.
* Fix *pimd/langevin/bosonic* currently does not support *ensemble* other than
*nve*, *nvt*. The barostat related keywords *iso*, *aniso*,
*barostat*, *taup* are not supported.
* Fix *pimd/langevin/bosonic* also has a keyword not available in fix
*pimd/langevin*: *esynch*, with default *yes*. If set to *no*, some
time consuming synchronization of spring energies and the primitive
kinetic energy estimator between processors is avoided.
The isomorphism between the partition function of :math:`N` bosonic
quantum particles and that of a system of classical ring polymers at
inverse temperature :math:`\beta` is given by :ref:`(Tuckerman)
<book-Tuckerman>`:
.. math::
Z \propto \int d\mathbf{q} \cdot \frac{1}{N!} \sum_\sigma \textrm{exp} [ -\beta \left( E^\sigma + V \right) ].
Here, :math:`V` is the potential between different particles at the same
imaginary time slice, which is the same for bosons and distinguishable
particles. The sum is over all permutations :math:`\sigma`. Recall that
a permutation matches each element :math:`l` in :math:`1, ..., N` to an
element :math:`\sigma(l)` in :math:`1, ..., N` without repetitions. The
energies :math:`E^\sigma` correspond to the linking of ring polymers of
different particles according to the permutations:
.. math::
E^\sigma = \frac{mP}{2\beta^2 \hbar^2} \sum_{\ell=1}^N \sum_{j=1}^P \left(\mathbf{q}_\ell^j - \mathbf{q}_\ell^{j+1}\right)^2,
where :math:`P` is the number of beads and :math:`\mathbf{q}_\ell^{P+1}=\mathbf{q}_{\sigma(\ell)}^1.`
Hirshberg et. al. showed that the ring polymer potential
:math:`-\frac{1}{\beta}\textrm{ln}\left[ \frac{1}{N!} \sum_\sigma e ^ {
-\beta E^\sigma } \right]`, which scales exponentially with :math:`N`,
can be replaced by a potential :math:`V^{[1,N]}` defined through a
recurrence relation :ref:`(Hirshberg1) <Hirshberg>`:
.. math::
e ^ { -\beta V^{[1,N]} } = \frac{1}{N} \sum_{k=1}^N e ^ { -\beta \left( V^{[1,N-k]} + E^{[N-K+1,N]} \right)}.
Here, :math:`E^{[N-K+1,N]}` is the spring energy of the ring polymer
obtained by connecting the beads of particles :math:`N - k + 1, N - k +
2, ..., N` in a cycle. This potential does not include all :math:`N!`
permutations, but samples the same bosonic partition function. The
implemented algorithm in LAMMPS for calculating the potential is the one
developed by Feldman and Hirshberg, which scales like :math:`N^2+PN`
:ref:`(Feldman) <Feldman>`. The forces are calculated as weighted
averages over the representative permutations, through an algorithm that
scales the same as the one for the potential calculation, :math:`N^2+PN`
:ref:`(Feldman) <Feldman>`. The minimum-image convention is employed on
the springs to account for periodic boundary conditions; an elaborate
discussion of the validity of the approximation is available in
:ref:`(Higer) <HigerFeldman>`.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
The use of :doc:`binary restart files <restart>` and :doc:`fix_modify
<fix_modify>` is the same as in :doc:`fix pimd <fix_pimd>`.
Fix *pimd/nvt/bosonic* computes a global 4-vector, which can be accessed by
various :doc:`output commands <Howto_output>`. The quantities in
the global vector are:
#. the total spring energy of the quasi-beads,
#. the current temperature of the classical system of ring polymers,
#. the current value of the scalar virial estimator for the kinetic
energy of the quantum system :ref:`(Herman) <HermanBB>` (see the justification in the supporting information of :ref:`(Hirshberg2) <HirshbergInvernizzi>`),
#. the current value of the scalar primitive estimator for the kinetic
energy of the quantum system :ref:`(Hirshberg1) <Hirshberg>`.
The vector values calculated by fix *pimd/nvt/bosonic* are "extensive", except
for the temperature, which is "intensive".
Fix *pimd/langevin/bosonic* computes a global 6-vector, which can be accessed
by various :doc:`output commands <Howto_output>`. The quantities in the
global vector are:
#. kinetic energy of the beads,
#. spring elastic energy of the beads,
#. potential energy of the bead,
#. total energy of all beads (conserved if *ensemble* is *nve*) if *esynch* is *yes*
#. primitive kinetic energy estimator :ref:`(Hirshberg1) <Hirshberg>`
#. virial energy estimator :ref:`(Herman) <HermanBB>` (see the justification in the supporting information of :ref:`(Hirshberg2) <HirshbergInvernizzi>`).
The first three are different for different log files, and the others
are the same for different log files, except for the primitive kinetic
energy estimator when setting *esynch* to *no*. Then, the primitive
kinetic energy estimator is obtained by summing over all log files.
Also note that when *esynch* is set to *no*, the fourth output gives the
total energy of all beads excluding the spring elastic energy; the total
classical energy can then be obtained by adding the sum of second output
over all log files. All vector values calculated by fix
*pimd/langevin/bosonic* are "extensive".
For both *pimd/nvt/bosonic* and *pimd/langevin/bosonic*, the contribution of the
exterior spring to the primitive estimator is printed to the first log
file. The contribution of the :math:`P-1` interior springs is printed
to the other :math:`P-1` log files. The contribution of the constant
:math:`\frac{PdN}{2 \beta}` (with :math:`d` being the dimensionality) is
equally divided over log files.
Restrictions
""""""""""""
These fixes are part of the REPLICA package. They are only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
The restrictions of :doc:`fix pimd <fix_pimd>` apply.
Related commands
""""""""""""""""
:doc:`pimd/nvt <fix_pimd>`, :doc:`pimd/langevin <fix_pimd>`
Default
"""""""
The keyword defaults for fix *pimd/nvt/bosonic* are method = pimd, fmass = 1.0,
sp = 1.0, temp = 300.0, and nhc = 2.
The keyword defaults for fix *pimd/langevin/bosonic* are integrator = obabo,
method = pimd, ensemble = nvt, fmass = 1.0, temp = 298.15, thermostat =
PILE_L, tau = 1.0, fixcom = yes, esynch = yes, and lj = 1 for all its
arguments.
----------
.. _book-Tuckerman:
**(Tuckerman)** M. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford University Press, 2010)
.. _Hirshberg:
**(Hirshberg1)** B. Hirshberg, V. Rizzi, and M. Parrinello, "Path integral molecular dynamics for bosons," Proc. Natl. Acad. Sci. U. S. A. 116, 21445 (2019)
.. _HirshbergInvernizzi:
**(Hirshberg2)** B. Hirshberg, M. Invernizzi, and M. Parrinello, "Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality," J Chem Phys, 152, 171102 (2020)
.. _Feldman:
**(Feldman)** Y. M. Y. Feldman and B. Hirshberg, "Quadratic scaling bosonic path integral molecular dynamics," J. Chem. Phys. 159, 154107 (2023)
.. _HigerFeldman:
**(Higer)** J. Higer, Y. M. Y. Feldman, and B. Hirshberg, "Periodic Boundary Conditions for Bosonic Path Integral Molecular Dynamics," arXiv:2501.17618 (2025)
.. _HermanBB:
**(Herman)** M. F. Herman, E. J. Bruskin, B. J. Berne, J Chem Phys, 76, 5150 (1982).