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

Merge pull request #2 from akohlmey/bosonic-pimd-langevin

Browse Source 
Updates to LAMMPS pull request 4479
This commit is contained in:
yotamfe
2025-02-23 20:15:19 +02:00
committed by GitHub
parent df360f0d55 c7e288eb8c
commit aee34a92d3
12 changed files with 650 additions and 560 deletions
Download Patch File Download Diff File Expand all files Collapse all files

50
doc/src/fix_pimd.rst
View File

@ -1,5 +1,5 @@
.. index:: fix pimd/langevin
.. index:: fix pimd/nvt
.. index:: fix pimd/nvt
fix pimd/langevin command
=========================
@ -20,35 +20,37 @@ Syntax
* keywords for style *pimd/nvt*
.. parsed-literal::
*keywords* = *method* or *fmass* or *sp* or *temp* or *nhc*
*method* value = *pimd* or *nmpimd* or *cmd*
*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* = *method* or *fmass* or *sp* or *temp* or *nhc*
*method* value = *pimd* or *nmpimd* or *cmd*
*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*
.. parsed-literal::
*keywords* = *method* or *integrator* or *ensemble* or *fmmode* or *fmass* or *scale* or *temp* or *thermostat* or *tau* or *iso* or *aniso* or *barostat* or *taup* or *fixcom* or *lj*
*method* value = *nmpimd* (default) or *pimd*
*integrator* value = *obabo* or *baoab*
*ensemble* value = *nvt* or *nve* or *nph* or *npt*
*fmmode* value = *physical* or *normal*
*fmass* value = scaling factor on mass
*temp* value = temperature (temperature unit)
*keywords* = *method* or *integrator* or *ensemble* or *fmmode* or *fmass* or *scale* or *temp* or *thermostat* or *tau* or *iso* or *aniso* or *barostat* or *taup* or *fixcom* or *lj*
*method* value = *nmpimd* (default) or *pimd*
*integrator* value = *obabo* or *baoab*
*ensemble* value = *nvt* or *nve* or *nph* or *npt*
*fmmode* value = *physical* or *normal*
*fmass* value = scaling factor on mass
*temp* value = temperature (temperature unit)
temperature = target temperature of the thermostat
*thermostat* values = style seed
*thermostat* values = style seed
style value = *PILE_L*
seed = random number generator seed
*tau* value = thermostat damping parameter (time unit)
*scale* value = scaling factor of the damping times of non-centroid modes of PILE_L thermostat
*iso* or *aniso* values = pressure (pressure unit)
*tau* value = thermostat damping parameter (time unit)
*scale* value = scaling factor of the damping times of non-centroid modes of PILE_L thermostat
*iso* or *aniso* values = pressure (pressure unit)
pressure = scalar external pressure of the barostat
*barostat* value = *BZP* or *MTTK*
*taup* value = barostat damping parameter (time unit)
*fixcom* value = *yes* or *no*
*lj* values = epsilon sigma mass planck mvv2e
*barostat* value = *BZP* or *MTTK*
*taup* value = barostat 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)
@ -250,7 +252,7 @@ system, each running on one of the 4 partitions of processors. Each
replica (shown in green) owns one quasi-bead in each ring.
.. image:: JPG/pimd.jpg
:align: center
:align: center
To run a PIMD simulation with M quasi-beads in each ring polymer using
N MPI tasks for each partition's domain-decomposition, you would use P
@ -419,7 +421,7 @@ Default
The keyword defaults for fix *pimd/nvt* are method = pimd, fmass = 1.0, sp
= 1.0, temp = 300.0, and nhc = 2.
The keyord defaults for fix *pimd/langevin* are integrator = obabo, method = nmpimd, ensemble = nvt, fmmode = physical, fmass = 1.0,
The keyword defaults for fix *pimd/langevin* are integrator = obabo, method = nmpimd, ensemble = nvt, fmmode = physical, fmass = 1.0,
scale = 1, temp = 298.15, thermostat = PILE_L, tau = 1.0, iso = 1.0, taup = 1.0, barostat = BZP, fixcom = yes, and lj = 1 for all its arguments.
----------

165
doc/src/fix_pimdb.rst
View File

@ -1,5 +1,5 @@
.. index:: fix pimdb/langevin
.. index:: fix pimdb/nvt
.. index:: fix pimdb/nvt
fix pimdb/langevin command
==========================
@ -20,34 +20,36 @@ Syntax
* keywords for style *pimdb/nvt*
.. 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* = *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 *pimdb/langevin*
.. 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)
*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
*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
*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*
*esynch* value = *yes* or *no*
Examples
""""""""
@ -60,28 +62,48 @@ Examples
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 *pimdb/nvt* and fix *pimdb/langevin* perform simulations of bosons, including exchange effects.
The *pimdb* commands share syntax with the equivalent *pimd* commands. The user is referred to the documentation of the *pimd* fix for a
detailed syntax description and additional, general capabilities of the commands.
The major differences from fix *pimd* in terms of capabilities are:
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 *pimdb/nvt* and fix
*pimdb/langevin* perform simulations of bosons, including exchange
effects. The *pimdb* commands share syntax with the equivalent *pimd*
commands. The user is referred to the documentation of the *pimd* fix
for a detailed syntax description and additional, general capabilities
of the commands. The major differences from fix *pimd* in terms of
capabilities are:
* Fix *pimdb/nvt* the only supports the "pimd" and "nmpimd" methods. Fix *pimdb/langevin* 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 *pimdb/langevin* currently does not support *ensemble* other than *nve*, *nvt*. The barostat related keywords *iso*, *aniso*, *barostat*, *taup* are not supported.
* Fix *pimdb/langevin* 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.
* Fix *pimdb/nvt* only supports the "pimd" and "nmpimd" methods. Fix
*pimdb/langevin* 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 *pimdb/langevin* currently does not support *ensemble* other than
*nve*, *nvt*. The barostat related keywords *iso*, *aniso*,
*barostat*, *taup* are not supported.
* Fix *pimdb/langevin* 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>`:
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{\bf 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:
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::
@ -90,25 +112,34 @@ distinguishable particles. The sum is over all permutations :math:`\sigma`. Reca
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:`-\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>`.
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>`.
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 *pimdb/nvt* computes a global 4-vector, which can be accessed by
various :doc:`output commands <Howto_output>`. The quantities in
@ -121,11 +152,12 @@ the global vector are:
#. 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 *pimdb/nvt* are "extensive", except for the
temperature, which is "intensive".
The vector values calculated by fix *pimdb/nvt* are "extensive", except
for the temperature, which is "intensive".
Fix *pimdb/langevin* computes a global 6-vector, which
can be accessed by various :doc:`output commands <Howto_output>`. The quantities in the global vector are:
Fix *pimdb/langevin* 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,
@ -134,15 +166,22 @@ can be accessed by various :doc:`output commands <Howto_output>`. The quantities
#. 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 *pimdb/langevin* are "extensive".
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
*pimdb/langevin* are "extensive".
For both *pimdb/nvt* and *pimdb/langevin*, 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 dimentionality) is equally divided over log files.
For both *pimdb/nvt* and *pimdb/langevin*, 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
""""""""""""
@ -156,11 +195,13 @@ The restrictions of :doc:`fix pimd <fix_pimd>` apply.
Default
"""""""
The keyword defaults for fix *pimdb/nvt* are method = pimd, fmass = 1.0, sp
= 1.0, temp = 300.0, and nhc = 2.
The keyword defaults for fix *pimdb/nvt* are method = pimd, fmass = 1.0,
sp = 1.0, temp = 300.0, and nhc = 2.
The keyord defaults for fix *pimdb/langevin* 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.
The keyword defaults for fix *pimdb/langevin* 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.
----------
@ -170,19 +211,19 @@ temp = 298.15, thermostat = PILE_L, tau = 1.0, fixcom = yes, esynch = yes, and l
.. _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)
**(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)
**(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)
**(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)
**(Higer)** J. Higer, Y. M. Y. Feldman, and B. Hirshberg, "Periodic Boundary Conditions for Bosonic Path Integral Molecular Dynamics," arXiv:2501.17618 (2025)
.. _HermanBB:

16
doc/utils/sphinx-config/false_positives.txt
View File

@ -8,6 +8,7 @@ abi
abo
Abramyan
absTol
Acad
Acc
Accelerys
Accelrys
@ -338,6 +339,7 @@ bodyflag
bodyforce
bodystyle
Bogaerts
Bogoliubov
Bogusz
Bohrs
boltz
@ -357,6 +359,8 @@ boostostat
boostostatting
Boresch
borophene
bosonic
bosons
Botero
Botu
Bouguet
@ -1079,6 +1083,7 @@ estretch
esu
esub
esw
esynch
et
etag
etap
@ -1142,6 +1147,7 @@ fdotr
fdt
fe
Fehlberg
Feldman
Fellinger
femtosecond
femtoseconds
@ -1151,6 +1157,8 @@ Fennell
fep
FEP
fermi
fermion
fermions
Fermionic
Ferrand
fexternal
@ -1191,6 +1199,7 @@ Finnis
Fiorin
fitpod
fivebody
fixcom
fixID
fj
Fji
@ -1211,6 +1220,7 @@ fmag
fmass
fmatch
fmm
fmmode
fmt
fmtlib
fmx
@ -1484,6 +1494,7 @@ hgrid
hhmrr
Hibbs
Higdon
Higer
hiID
Hijazi
Hilger
@ -2608,6 +2619,7 @@ nn
nnodes
npits
npj
nmpimd
nO
Nocedal
nocite
@ -2736,6 +2748,7 @@ nylo
nz
Nz
nzlo
obabo
ochre
ocl
octahedral
@ -2935,6 +2948,7 @@ Pieniazek
Pieter
pIm
pimd
pimdb
Piola
pIp
pipelining
@ -3284,6 +3298,7 @@ RiRj
Risi
Rix
Riy
Rizzi
rj
Rj
Rjinner
@ -3717,6 +3732,7 @@ Tanmoy
Tartakovsky
taskset
taubi
taup
Tavenner
taylor
tb

516
src/REPLICA/bosonic_exchange.cpp
View File

@ -24,339 +24,357 @@
------------------------------------------------------------------------- */
#include "bosonic_exchange.h"
#include "domain.h"
#include "universe.h"
#include "memory.h"
#include "error.h"
#include <cassert>
#include "memory.h"
#include "universe.h"
using namespace LAMMPS_NS;
BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic, bool beta_convention) :
Pointers(lmp),
nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic), physical_beta_convention(beta_convention){
memory->create(temp_nbosons_array, nbosons + 1, "BosonicExchange: temp_nbosons_array");
memory->create(E_kn, (nbosons * (nbosons + 1) / 2), "BosonicExchange: E_kn");
memory->create(V, nbosons + 1, "BosonicExchange: V");
memory->create(V_backwards, nbosons + 1, "BosonicExchange: V_backwards");
memory->create(connection_probabilities, nbosons * nbosons, "BosonicExchange: connection probabilities");
BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic,
bool beta_convention) :
Pointers(lmp), nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic),
physical_beta_convention(beta_convention), x(nullptr), x_prev(nullptr), x_next(nullptr),
E_kn(nullptr), V(nullptr), V_backwards(nullptr), connection_probabilities(nullptr),
temp_nbosons_array(nullptr)
{
memory->create(temp_nbosons_array, nbosons + 1, "BosonicExchange: temp_nbosons_array");
memory->create(E_kn, (nbosons * (nbosons + 1) / 2), "BosonicExchange: E_kn");
memory->create(V, nbosons + 1, "BosonicExchange: V");
memory->create(V_backwards, nbosons + 1, "BosonicExchange: V_backwards");
memory->create(connection_probabilities, nbosons * nbosons,
"BosonicExchange: connection probabilities");
}
void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next,
double beta, double spring_constant) {
this->x = x;
this->x_prev = x_prev;
this->x_next = x_next;
this->beta = beta;
this->spring_constant = spring_constant;
void BosonicExchange::prepare_with_coordinates(const double *x, const double *x_prev,
const double *x_next, double beta,
double spring_constant)
{
this->x = x;
this->x_prev = x_prev;
this->x_next = x_next;
this->beta = beta;
this->spring_constant = spring_constant;
if (bead_num == 0 || bead_num == np - 1) {
// exterior beads
evaluate_cycle_energies();
Evaluate_VBn();
Evaluate_V_backwards();
evaluate_connection_probabilities();
}
if (bead_num == 0 || bead_num == np - 1) {
// exterior beads
evaluate_cycle_energies();
Evaluate_VBn();
Evaluate_V_backwards();
evaluate_connection_probabilities();
}
}
/* ---------------------------------------------------------------------- */
BosonicExchange::~BosonicExchange() {
memory->destroy(connection_probabilities);
memory->destroy(V_backwards);
memory->destroy(V);
memory->destroy(E_kn);
memory->destroy(temp_nbosons_array);
BosonicExchange::~BosonicExchange()
{
memory->destroy(connection_probabilities);
memory->destroy(V_backwards);
memory->destroy(V);
memory->destroy(E_kn);
memory->destroy(temp_nbosons_array);
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::diff_two_beads(const double* x1, int l1, const double* x2, int l2,
double diff[3]) const {
l1 = l1 % nbosons;
l2 = l2 % nbosons;
double delx2 = x2[3 * l2 + 0] - x1[3 * l1 + 0];
double dely2 = x2[3 * l2 + 1] - x1[3 * l1 + 1];
double delz2 = x2[3 * l2 + 2] - x1[3 * l1 + 2];
if (apply_minimum_image) {
domain->minimum_image(delx2, dely2, delz2);
}
void BosonicExchange::diff_two_beads(const double *x1, int l1, const double *x2, int l2,
double diff[3]) const
{
l1 = l1 % nbosons;
l2 = l2 % nbosons;
double delx2 = x2[3 * l2 + 0] - x1[3 * l1 + 0];
double dely2 = x2[3 * l2 + 1] - x1[3 * l1 + 1];
double delz2 = x2[3 * l2 + 2] - x1[3 * l1 + 2];
if (apply_minimum_image) { domain->minimum_image(delx2, dely2, delz2); }
diff[0] = delx2;
diff[1] = dely2;
diff[2] = delz2;
diff[0] = delx2;
diff[1] = dely2;
diff[2] = delz2;
}
/* ---------------------------------------------------------------------- */
double BosonicExchange::distance_squared_two_beads(const double* x1, int l1, const double* x2, int l2) const {
double diff[3];
diff_two_beads(x1, l1, x2, l2, diff);
return diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
double BosonicExchange::distance_squared_two_beads(const double *x1, int l1, const double *x2,
int l2) const
{
double diff[3];
diff_two_beads(x1, l1, x2, l2, diff);
return diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::evaluate_cycle_energies()
{
const double* x_first_bead;
const double* x_last_bead;
const double *x_first_bead;
const double *x_last_bead;
if (bead_num == 0) {
x_first_bead = x;
x_last_bead = x_prev;
} else {
x_first_bead = x_next;
x_last_bead = x;
}
for (int i = 0; i < nbosons; i++) {
temp_nbosons_array[i] = distance_squared_two_beads(x_first_bead, i, x_last_bead, i);
}
for (int v = 0; v < nbosons; v++) {
set_Enk(v + 1, 1,
0.5 * spring_constant * (temp_nbosons_array[v]));
for (int u = v - 1; u >= 0; u--) {
double val = get_Enk(v + 1, v - u) +
0.5 * spring_constant * (
// connect u to u+1
+ distance_squared_two_beads(x_last_bead, u, x_first_bead, u + 1)
// break cycle [u+1,v]
- distance_squared_two_beads(x_first_bead, u + 1, x_last_bead, v)
// close cycle from v to u
+ distance_squared_two_beads(x_first_bead, u, x_last_bead, v));
set_Enk(v + 1, v - u + 1, val);
}
if (bead_num == 0) {
x_first_bead = x;
x_last_bead = x_prev;
} else {
x_first_bead = x_next;
x_last_bead = x;
}
for (int i = 0; i < nbosons; i++) {
temp_nbosons_array[i] = distance_squared_two_beads(x_first_bead, i, x_last_bead, i);
}
for (int v = 0; v < nbosons; v++) {
set_Enk(v + 1, 1, 0.5 * spring_constant * (temp_nbosons_array[v]));
for (int u = v - 1; u >= 0; u--) {
double val = get_Enk(v + 1, v - u) +
0.5 * spring_constant *
(
// connect u to u+1
+distance_squared_two_beads(x_last_bead, u, x_first_bead, u + 1)
// break cycle [u+1,v]
- distance_squared_two_beads(x_first_bead, u + 1, x_last_bead, v)
// close cycle from v to u
+ distance_squared_two_beads(x_first_bead, u, x_last_bead, v));
set_Enk(v + 1, v - u + 1, val);
}
}
}
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_Enk(int m, int k) const {
int end_of_m = m * (m + 1) / 2;
return E_kn[end_of_m - k];
double BosonicExchange::get_Enk(int m, int k) const
{
int end_of_m = m * (m + 1) / 2;
return E_kn[end_of_m - k];
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::set_Enk(int m, int k, double val) {
int end_of_m = m * (m + 1) / 2;
E_kn[end_of_m - k] = val;
void BosonicExchange::set_Enk(int m, int k, double val)
{
int end_of_m = m * (m + 1) / 2;
E_kn[end_of_m - k] = val;
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::Evaluate_VBn()
{
V[0] = 0.0;
V[0] = 0.0;
for (int m = 1; m < nbosons + 1; m++) {
double Elongest = std::numeric_limits<double>::max();
for (int m = 1; m < nbosons + 1; m++) {
double Elongest = std::numeric_limits<double>::max();
for (int k = m; k > 0; k--) {
double val = get_Enk(m,k) + V[m-k];
Elongest = std::min(Elongest, val); // Numerical stability (Xiong & Xiong, doi:10.1103/PhysRevE.106.025309)
temp_nbosons_array[k - 1] = val;
}
double sig_denom = 0.0;
for (int k = m; k > 0; k--) {
sig_denom += exp(-beta * (temp_nbosons_array[k - 1] - Elongest));
}
V[m] = Elongest - (1.0 / beta) * log(sig_denom / (double)m);
if (!std::isfinite(V[m])) {
error->universe_one(
FLERR,
fmt::format("Invalid sig_denom {} with Elongest {} in bosonic exchange potential",
sig_denom, Elongest));
}
for (int k = m; k > 0; k--) {
double val = get_Enk(m, k) + V[m - k];
Elongest = std::min(
Elongest, val); // Numerical stability (Xiong & Xiong, doi:10.1103/PhysRevE.106.025309)
temp_nbosons_array[k - 1] = val;
}
double sig_denom = 0.0;
for (int k = m; k > 0; k--) {
sig_denom += exp(-beta * (temp_nbosons_array[k - 1] - Elongest));
}
V[m] = Elongest - (1.0 / beta) * log(sig_denom / (double) m);
if (!std::isfinite(V[m])) {
error->universe_one(
FLERR,
fmt::format("Invalid sig_denom {} with Elongest {} in bosonic exchange potential",
sig_denom, Elongest));
}
}
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::Evaluate_V_backwards() {
V_backwards[nbosons] = 0.0;
void BosonicExchange::Evaluate_V_backwards()
{
V_backwards[nbosons] = 0.0;
for (int l = nbosons - 1; l > 0; l--) {
double Elongest = std::numeric_limits<double>::max();
for (int p = l; p < nbosons; p++) {
double val = get_Enk(p + 1, p - l + 1) + V_backwards[p + 1];
Elongest = std::min(Elongest, val);
temp_nbosons_array[p] = val;
}
double sig_denom = 0.0;
for (int p = l; p < nbosons; p++) {
sig_denom += 1.0 / (p + 1) * exp(-beta *
(temp_nbosons_array[p]
- Elongest)
);
}
V_backwards[l] = Elongest - log(sig_denom) / beta;
if (!std::isfinite(V_backwards[l])) {
error->universe_one(
FLERR,
fmt::format("Invalid sig_denom {} with Elongest {} in bosonic exchange potential backwards",
sig_denom, Elongest));
}
for (int l = nbosons - 1; l > 0; l--) {
double Elongest = std::numeric_limits<double>::max();
for (int p = l; p < nbosons; p++) {
double val = get_Enk(p + 1, p - l + 1) + V_backwards[p + 1];
Elongest = std::min(Elongest, val);
temp_nbosons_array[p] = val;
}
V_backwards[0] = V[nbosons];
double sig_denom = 0.0;
for (int p = l; p < nbosons; p++) {
sig_denom += 1.0 / (p + 1) * exp(-beta * (temp_nbosons_array[p] - Elongest));
}
V_backwards[l] = Elongest - log(sig_denom) / beta;
if (!std::isfinite(V_backwards[l])) {
error->universe_one(
FLERR,
fmt::format(
"Invalid sig_denom {} with Elongest {} in bosonic exchange potential backwards",
sig_denom, Elongest));
}
}
V_backwards[0] = V[nbosons];
}
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_potential() const {
return V[nbosons];
double BosonicExchange::get_potential() const
{
return V[nbosons];
}
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_interior_bead_spring_energy() const {
double spring_energy_for_bead = 0.;
for (int i = 0; i < nbosons; i++) {
spring_energy_for_bead += 0.5 * spring_constant * distance_squared_two_beads(x, i, x_prev, i);
}
return spring_energy_for_bead;
double BosonicExchange::get_interior_bead_spring_energy() const
{
double spring_energy_for_bead = 0.;
for (int i = 0; i < nbosons; i++) {
spring_energy_for_bead += 0.5 * spring_constant * distance_squared_two_beads(x, i, x_prev, i);
}
return spring_energy_for_bead;
}
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_bead_spring_energy() const {
double spring_energy_for_bead = (bead_num == 0 ? get_potential() : get_interior_bead_spring_energy());
return spring_energy_for_bead;
double BosonicExchange::get_bead_spring_energy() const
{
double spring_energy_for_bead =
(bead_num == 0 ? get_potential() : get_interior_bead_spring_energy());
return spring_energy_for_bead;
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force(double** f) const {
if (bead_num == np - 1) {
spring_force_last_bead(f);
} else if (bead_num == 0) {
spring_force_first_bead(f);
} else {
spring_force_interior_bead(f);
}
void BosonicExchange::spring_force(double **f) const
{
if (bead_num == np - 1) {
spring_force_last_bead(f);
} else if (bead_num == 0) {
spring_force_first_bead(f);
} else {
spring_force_interior_bead(f);
}
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::evaluate_connection_probabilities() {
for (int l = 0; l < nbosons - 1; l++) {
double direct_link_probability = 1.0 - (exp(-beta *
(V[l + 1] + V_backwards[l + 1] -
V[nbosons])));
connection_probabilities[nbosons * l + (l + 1)] = direct_link_probability;
}
for (int u = 0; u < nbosons; u++) {
for (int l = u; l < nbosons; l++) {
double close_cycle_probability = 1.0 / (l + 1) *
exp(-beta * (V[u] + get_Enk(l + 1, l - u + 1) + V_backwards[l + 1]
- V[nbosons]));
connection_probabilities[nbosons * l + u] = close_cycle_probability;
}
void BosonicExchange::evaluate_connection_probabilities()
{
for (int l = 0; l < nbosons - 1; l++) {
double direct_link_probability =
1.0 - (exp(-beta * (V[l + 1] + V_backwards[l + 1] - V[nbosons])));
connection_probabilities[nbosons * l + (l + 1)] = direct_link_probability;
}
for (int u = 0; u < nbosons; u++) {
for (int l = u; l < nbosons; l++) {
double close_cycle_probability = 1.0 / (l + 1) *
exp(-beta * (V[u] + get_Enk(l + 1, l - u + 1) + V_backwards[l + 1] - V[nbosons]));
connection_probabilities[nbosons * l + u] = close_cycle_probability;
}
}
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force_last_bead(double** f) const {
const double* x_first_bead = x_next;
const double* x_last_bead = x;
void BosonicExchange::spring_force_last_bead(double **f) const
{
const double *x_first_bead = x_next;
const double *x_last_bead = x;
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
double sum_z = 0.0;
for (int next_l = 0; next_l <= l + 1 && next_l < nbosons; next_l++) {
double diff_next[3];
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
double sum_z = 0.0;
for (int next_l = 0; next_l <= l + 1 && next_l < nbosons; next_l++) {
double diff_next[3];
diff_two_beads(x_last_bead, l, x_first_bead, next_l, diff_next);
diff_two_beads(x_last_bead, l, x_first_bead, next_l, diff_next);
double prob = connection_probabilities[nbosons * l + next_l];
double prob = connection_probabilities[nbosons * l + next_l];
sum_x += prob * diff_next[0];
sum_y += prob * diff_next[1];
sum_z += prob * diff_next[2];
}
double diff_prev[3];
diff_two_beads(x_last_bead, l, x_prev, l, diff_prev);
sum_x += diff_prev[0];
sum_y += diff_prev[1];
sum_z += diff_prev[2];
f[l][0] += sum_x * spring_constant;
f[l][1] += sum_y * spring_constant;
f[l][2] += sum_z * spring_constant;
sum_x += prob * diff_next[0];
sum_y += prob * diff_next[1];
sum_z += prob * diff_next[2];
}
double diff_prev[3];
diff_two_beads(x_last_bead, l, x_prev, l, diff_prev);
sum_x += diff_prev[0];
sum_y += diff_prev[1];
sum_z += diff_prev[2];
f[l][0] += sum_x * spring_constant;
f[l][1] += sum_y * spring_constant;
f[l][2] += sum_z * spring_constant;
}
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force_first_bead(double** f) const {
const double* x_first_bead = x;
const double* x_last_bead = x_prev;
void BosonicExchange::spring_force_first_bead(double **f) const
{
const double *x_first_bead = x;
const double *x_last_bead = x_prev;
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
double sum_z = 0.0;
for (int prev_l = std::max(0, l - 1); prev_l < nbosons; prev_l++) {
double diff_prev[3];
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
double sum_z = 0.0;
for (int prev_l = std::max(0, l - 1); prev_l < nbosons; prev_l++) {
double diff_prev[3];
diff_two_beads(x_first_bead, l, x_last_bead, prev_l, diff_prev);
diff_two_beads(x_first_bead, l, x_last_bead, prev_l, diff_prev);
double prob = connection_probabilities[nbosons * prev_l + l];
double prob = connection_probabilities[nbosons * prev_l + l];
sum_x += prob * diff_prev[0];
sum_y += prob * diff_prev[1];
sum_z += prob * diff_prev[2];
}
double diff_next[3];
diff_two_beads(x_first_bead, l, x_next, l, diff_next);
sum_x += diff_next[0];
sum_y += diff_next[1];
sum_z += diff_next[2];
f[l][0] += sum_x * spring_constant;
f[l][1] += sum_y * spring_constant;
f[l][2] += sum_z * spring_constant;
sum_x += prob * diff_prev[0];
sum_y += prob * diff_prev[1];
sum_z += prob * diff_prev[2];
}
double diff_next[3];
diff_two_beads(x_first_bead, l, x_next, l, diff_next);
sum_x += diff_next[0];
sum_y += diff_next[1];
sum_z += diff_next[2];
f[l][0] += sum_x * spring_constant;
f[l][1] += sum_y * spring_constant;
f[l][2] += sum_z * spring_constant;
}
}
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force_interior_bead(double **f) const {
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
double sum_z = 0.0;
void BosonicExchange::spring_force_interior_bead(double **f) const
{
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
double sum_z = 0.0;
double diff_prev[3];
diff_two_beads(x, l, x_prev, l, diff_prev);
sum_x += diff_prev[0];
sum_y += diff_prev[1];
sum_z += diff_prev[2];
double diff_prev[3];
diff_two_beads(x, l, x_prev, l, diff_prev);
sum_x += diff_prev[0];
sum_y += diff_prev[1];
sum_z += diff_prev[2];
double diff_next[3];
diff_two_beads(x, l, x_next, l, diff_next);
sum_x += diff_next[0];
sum_y += diff_next[1];
sum_z += diff_next[2];
double diff_next[3];
diff_two_beads(x, l, x_next, l, diff_next);
sum_x += diff_next[0];
sum_y += diff_next[1];
sum_z += diff_next[2];
f[l][0] += sum_x * spring_constant;
f[l][1] += sum_y * spring_constant;
f[l][2] += sum_z * spring_constant;
}
f[l][0] += sum_x * spring_constant;
f[l][1] += sum_y * spring_constant;
f[l][2] += sum_z * spring_constant;
}
}
/* ---------------------------------------------------------------------- */
@ -367,7 +385,8 @@ double BosonicExchange::prim_estimator()
// Adds the contribution of an interior spring, which is the same as for distinguishable particles.
if (bead_num != 0) {
return convention_correction * (0.5 * domain->dimension * nbosons / beta - get_interior_bead_spring_energy());
return convention_correction *
(0.5 * domain->dimension * nbosons / beta - get_interior_bead_spring_energy());
}
// For the first bead, add the contribution of the exterior bead (See Equations 4-5 in the
@ -375,28 +394,27 @@ double BosonicExchange::prim_estimator()
temp_nbosons_array[0] = 0.0;
for (int m = 1; m < nbosons + 1; ++m) {
double sig = 0.0;
temp_nbosons_array[m] = 0.0;
double sig = 0.0;
temp_nbosons_array[m] = 0.0;
double Elongest = std::numeric_limits<double>::max();
double Elongest = std::numeric_limits<double>::max();
// Numerical stability (Xiong & Xiong, doi:10.1103/PhysRevE.106.025309)
for (int k = m; k > 0; k--) {
Elongest = std::min(Elongest, get_Enk(m, k) + V[m - k]);
}
// Numerical stability (Xiong & Xiong, doi:10.1103/PhysRevE.106.025309)
for (int k = m; k > 0; k--) { Elongest = std::min(Elongest, get_Enk(m, k) + V[m - k]); }
for (int k = m; k > 0; --k) {
double E_kn_val = get_Enk(m, k);
for (int k = m; k > 0; --k) {
double E_kn_val = get_Enk(m, k);
sig += (temp_nbosons_array[m - k] - E_kn_val) * exp(-beta * (E_kn_val + V[m - k] - Elongest));
}
sig += (temp_nbosons_array[m - k] - E_kn_val) * exp(-beta * (E_kn_val + V[m - k] - Elongest));
}
double sig_denom_m = m * exp(-beta * (V[m] - Elongest));
double sig_denom_m = m * exp(-beta * (V[m] - Elongest));
temp_nbosons_array[m] = sig / sig_denom_m;
temp_nbosons_array[m] = sig / sig_denom_m;
}
return convention_correction * (0.5 * domain->dimension * nbosons / beta + temp_nbosons_array[nbosons]);
return convention_correction *
(0.5 * domain->dimension * nbosons / beta + temp_nbosons_array[nbosons]);
}
/* ---------------------------------------------------------------------- */

106
src/REPLICA/bosonic_exchange.h
View File

@ -18,66 +18,66 @@
namespace LAMMPS_NS {
class BosonicExchange : protected Pointers {
public:
BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic, bool beta_convention);
~BosonicExchange();
class BosonicExchange : protected Pointers {
public:
BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic,
bool beta_convention);
~BosonicExchange();
void prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next,
double beta, double spring_constant);
void prepare_with_coordinates(const double *x, const double *x_prev, const double *x_next,
double beta, double spring_constant);
double get_potential() const;
double get_bead_spring_energy() const;
double get_potential() const;
double get_bead_spring_energy() const;
void spring_force(double** f) const;
void spring_force(double **f) const;
double prim_estimator();
double prim_estimator();
private:
void evaluate_cycle_energies();
void diff_two_beads(const double* x1, int l1, const double* x2, int l2, double diff[3]) const;
double get_interior_bead_spring_energy() const;
double distance_squared_two_beads(const double* x1, int l1, const double* x2, int l2) const;
double get_Enk(int m, int k) const;
void set_Enk(int m, int k, double val);
void evaluate_connection_probabilities();
void spring_force_last_bead(double** f) const;
void spring_force_first_bead(double** f) const;
void spring_force_interior_bead(double** f) const;
void Evaluate_VBn();
void Evaluate_V_backwards();
private:
void evaluate_cycle_energies();
void diff_two_beads(const double *x1, int l1, const double *x2, int l2, double diff[3]) const;
double get_interior_bead_spring_energy() const;
double distance_squared_two_beads(const double *x1, int l1, const double *x2, int l2) const;
double get_Enk(int m, int k) const;
void set_Enk(int m, int k, double val);
void evaluate_connection_probabilities();
void spring_force_last_bead(double **f) const;
void spring_force_first_bead(double **f) const;
void spring_force_interior_bead(double **f) const;
void Evaluate_VBn();
void Evaluate_V_backwards();
const int nbosons;
const int np;
const int bead_num;
const bool apply_minimum_image;
// In the "reduced-beta convention" [e.g. in J. Chem. Phys. 133, 124104 (2010); also J. Chem. Phys. 74, 4078-4095 (1981)],
// the Boltzmann exponents have the form exp[-(beta/P)H], where H is the classical Hamiltonian of the
// ring polymers. This results in a canonical distribution at P times the physical temperature.
// In contrast, the "physical-beta convention" [e.g. in J. Chem. Phys. 99, 2796-2808 (1993)] uses weights of the form exp(-beta*H),
// such that the temperature of the canonical ensemble coincides with the physical temperature.
// Notably, the classical Hamiltonians of the two conventions differ, with the spring constant
// in the reduced-beta convention being P times larger than that in the physical-beta convention. Additionally, the reduced-beta convention
// lacks a 1/P prefactor in front of the external potential. The Hamiltonians of the two conventions are related through
// H_physical = H_reduced / P. Note however that the expressions for the various estimators are unaffected by this choice,
// so as the algorithm for bosonic exchange. The code below was designed to be compatible with both conventions,
// and the choice of convention only affects a single calculation within it.
// Setting the following boolian variable to false amounts to adopting the reduced-beta convention.
const bool physical_beta_convention;
const int nbosons;
const int np;
const int bead_num;
const bool apply_minimum_image;
// In the "reduced-beta convention" [e.g. in J. Chem. Phys. 133, 124104 (2010); also J. Chem. Phys. 74, 4078-4095 (1981)],
// the Boltzmann exponents have the form exp[-(beta/P)H], where H is the classical Hamiltonian of the
// ring polymers. This results in a canonical distribution at P times the physical temperature.
// In contrast, the "physical-beta convention" [e.g. in J. Chem. Phys. 99, 2796-2808 (1993)] uses weights of the form exp(-beta*H),
// such that the temperature of the canonical ensemble coincides with the physical temperature.
// Notably, the classical Hamiltonians of the two conventions differ, with the spring constant
// in the reduced-beta convention being P times larger than that in the physical-beta convention. Additionally, the reduced-beta convention
// lacks a 1/P prefactor in front of the external potential. The Hamiltonians of the two conventions are related through
// H_physical = H_reduced / P. Note however that the expressions for the various estimators are unaffected by this choice,
// so as the algorithm for bosonic exchange. The code below was designed to be compatible with both conventions,
// and the choice of convention only affects a single calculation within it.
// Setting the following boolian variable to false amounts to adopting the reduced-beta convention.
const bool physical_beta_convention;
double spring_constant;
double beta;
const double* x;
const double* x_prev;
const double* x_next;
double spring_constant;
double beta;
const double *x;
const double *x_prev;
const double *x_next;
double* E_kn;
double* V;
double* V_backwards;
double* connection_probabilities;
double* temp_nbosons_array;
};
}
double *E_kn;
double *V;
double *V_backwards;
double *connection_probabilities;
double *temp_nbosons_array;
};
} // namespace LAMMPS_NS
#endif

10
src/REPLICA/fix_pimd_langevin.cpp
View File

@ -53,8 +53,12 @@ using MathConst::MY_SQRT2;
using MathConst::THIRD;
using MathSpecial::powint;
static std::map<int, std::string> Barostats{{MTTK, "MTTK"}, {BZP, "BZP"}};
static std::map<int, std::string> Ensembles{{NVE, "NVE"}, {NVT, "NVT"}, {NPH, "NPH"}, {NPT, "NPT"}};
static std::map<int, std::string> Barostats{{FixPIMDLangevin::MTTK, "MTTK"},
{FixPIMDLangevin::BZP, "BZP"}};
static std::map<int, std::string> Ensembles{{FixPIMDLangevin::NVE, "NVE"},
{FixPIMDLangevin::NVT, "NVT"},
{FixPIMDLangevin::NPH, "NPH"},
{FixPIMDLangevin::NPT, "NPT"}};
/* ---------------------------------------------------------------------- */
@ -240,7 +244,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[i + 1], "no") == 0)
removecomflag = 0;
} else if (strcmp(arg[i], "") != 0) {
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix {}", arg[i], style));
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix {}", arg[i], style));
}
}

18
src/REPLICA/fix_pimd_langevin.h
View File

@ -22,15 +22,6 @@ FixStyle(pimd/langevin,FixPIMDLangevin);
#include "fix.h"
enum { PIMD, NMPIMD };
enum { PHYSICAL, NORMAL };
enum { BAOAB, OBABO };
enum { ISO, ANISO, TRICLINIC };
enum { PILE_L };
enum { MTTK, BZP };
enum { NVE, NVT, NPH, NPT };
enum { SINGLE_PROC, MULTI_PROC };
namespace LAMMPS_NS {
class FixPIMDLangevin : public Fix {
@ -38,6 +29,15 @@ class FixPIMDLangevin : public Fix {
FixPIMDLangevin(class LAMMPS *, int, char **);
~FixPIMDLangevin() override;
enum { PIMD, NMPIMD };
enum { PHYSICAL, NORMAL };
enum { BAOAB, OBABO };
enum { ISO, ANISO, TRICLINIC };
enum { PILE_L };
enum { MTTK, BZP };
enum { NVE, NVT, NPH, NPT };
enum { SINGLE_PROC, MULTI_PROC };
int setmask() override;
void init() override;

5
src/REPLICA/fix_pimd_nvt.h
View File

@ -25,13 +25,13 @@ FixStyle(pimd/nvt,FixPIMDNVT);
namespace LAMMPS_NS {
enum { PIMD, NMPIMD, CMD };
class FixPIMDNVT : public Fix {
public:
FixPIMDNVT(class LAMMPS *, int, char **);
~FixPIMDNVT() override;
enum { PIMD, NMPIMD, CMD };
int setmask() override;
void init() override;
@ -55,7 +55,6 @@ class FixPIMDNVT : public Fix {
void unpack_forward_comm(int, int, double *) override;
protected:
int method;
int np;
double inverse_np;

217
src/REPLICA/fix_pimdb_langevin.cpp
View File

@ -22,165 +22,172 @@
Version 1.0
------------------------------------------------------------------------- */
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <chrono>
#include <fstream>
#include "fix_pimdb_langevin.h"
#include "universe.h"
#include "comm.h"
#include "force.h"
#include "bosonic_exchange.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#include <algorithm>
#include "force.h"
#include "memory.h"
#include "universe.h"
#include "update.h"
#include <cmath>
#include <cstdlib>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) :
FixPIMDLangevin(lmp, narg, filtered_args = filter_args(narg, arg)), nbosons(atom->nlocal),
bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, false)
FixPIMDLangevin(lmp, narg, filtered_args = filter_args(narg, arg)), filtered_narg(narg),
nbosons(atom->nlocal)
{
synch_energies = true;
// Loop over the arguments with i++ instead of i+=2 like the parent,
// because the parent's loop has i++ inside some if blocks.
for (int i = 3; i < narg - 1; i++) {
if ((strcmp(arg[i], "method") == 0) && (strcmp(arg[i+1], "pimd") != 0)) {
error->universe_all(FLERR, "Method not supported in fix pimdb/langevin; only method PIMD");
}
else if (strcmp(arg[i], "scale") == 0) {
error->universe_all(FLERR, "The scale parameter of the PILE_L thermostat is not supported for pimdb, and should be removed.");
}
else if ((strcmp(arg[i], "iso") == 0) || (strcmp(arg[i], "aniso") == 0) || (strcmp(arg[i], "barostat") == 0) || (strcmp(arg[i], "taup") == 0)) {
error->universe_all(FLERR, "Barostat parameters are not available for pimdb.");
}
else if (strcmp(arg[i], "esynch") == 0) {
if (strcmp(arg[i + 1], "yes") == 0) {
synch_energies = true;
}
else if (strcmp(arg[i + 1], "no") == 0) {
synch_energies = false;
}
else {
error->universe_all(FLERR, "The esynch parameter can only receive yes or no!");
}
}
bosonic_exchange = new BosonicExchange(lmp, atom->nlocal, np, universe->me, true, false);
synch_energies = true;
// Loop over the arguments with i++ instead of i+=2 like the parent,
// because the parent's loop has i++ inside some if blocks.
for (int i = 3; i < narg - 1; i++) {
if ((strcmp(arg[i], "method") == 0) && (strcmp(arg[i + 1], "pimd") != 0)) {
error->universe_all(FLERR, "Method not supported in fix pimdb/langevin; only method PIMD");
} else if (strcmp(arg[i], "scale") == 0) {
error->universe_all(FLERR,
"The scale parameter of the PILE_L thermostat is not supported for "
"pimdb, and should be removed.");
} else if ((strcmp(arg[i], "iso") == 0) || (strcmp(arg[i], "aniso") == 0) ||
(strcmp(arg[i], "barostat") == 0) || (strcmp(arg[i], "taup") == 0)) {
error->universe_all(FLERR, "Barostat parameters are not available for pimdb.");
} else if (strcmp(arg[i], "esynch") == 0) {
if (strcmp(arg[i + 1], "yes") == 0) {
synch_energies = true;
} else if (strcmp(arg[i + 1], "no") == 0) {
synch_energies = false;
} else {
error->universe_all(FLERR, "The esynch parameter can only receive yes or no!");
}
}
}
if (fmmode != PHYSICAL) {
error->universe_all(FLERR, "The only available fmmode for pimdb is physical, please remove the fmmode keyword.");
}
if (ensemble != NVE && ensemble != NVT) {
error->universe_all(FLERR, "The only available ensembles for pimdb are nve and nvt, please choose one of these ensembles.");
}
if (fmmode != PHYSICAL) {
error->universe_all(
FLERR,
"The only available fmmode for pimdb is physical, please remove the fmmode keyword.");
}
if (ensemble != NVE && ensemble != NVT) {
error->universe_all(FLERR,
"The only available ensembles for pimdb are nve and nvt, please choose one "
"of these ensembles.");
}
method = PIMD;
method = PIMD;
size_vector = 6;
size_vector = 6;
memory->create(f_tag_order, nbosons, 3, "FixPIMDBLangevin:f_tag_order");
memory->create(f_tag_order, nbosons, 3, "FixPIMDBLangevin:f_tag_order");
}
/* ---------------------------------------------------------------------- */
FixPIMDBLangevin::~FixPIMDBLangevin() {
memory->destroy(f_tag_order);
memory->destroy(filtered_args);
FixPIMDBLangevin::~FixPIMDBLangevin()
{
memory->destroy(f_tag_order);
for (int i = 0; i < filtered_narg; ++i) delete[] filtered_args[i];
memory->destroy(filtered_args);
delete bosonic_exchange;
}
/* ---------------------------------------------------------------------- */
char** FixPIMDBLangevin::filter_args(int narg, char **arg)
char **FixPIMDBLangevin::filter_args(int narg, char **arg)
{
char** filtered_args = new char*[narg];
for (int i = 0; i < narg; i++) {
if (strcmp(arg[i], "esynch") == 0) {
filtered_args[i] = "";
}
else {
filtered_args[i] = arg[i];
}
filtered_narg = narg;
char **filtered_args = new char *[narg];
for (int i = 0; i < narg; i++) {
if (strcmp(arg[i], "esynch") == 0) {
filtered_args[i] = utils::strdup("");
} else {
filtered_args[i] = utils::strdup(arg[i]);
}
return filtered_args;
}
return filtered_args;
}
/* ---------------------------------------------------------------------- */
void FixPIMDBLangevin::prepare_coordinates()
{
inter_replica_comm(atom->x);
double ff = fbond * atom->mass[atom->type[0]];
int nlocal = atom->nlocal;
double* me_bead_positions = *(atom->x);
double* last_bead_positions = &bufsortedall[x_last * nlocal][0];
double* next_bead_positions = &bufsortedall[x_next * nlocal][0];
inter_replica_comm(atom->x);
double ff = fbond * atom->mass[atom->type[0]];
int nlocal = atom->nlocal;
double *me_bead_positions = *(atom->x);
double *last_bead_positions = &bufsortedall[x_last * nlocal][0];
double *next_bead_positions = &bufsortedall[x_next * nlocal][0];
bosonic_exchange.prepare_with_coordinates(me_bead_positions,
last_bead_positions, next_bead_positions,
beta_np, ff);
bosonic_exchange->prepare_with_coordinates(me_bead_positions, last_bead_positions,
next_bead_positions, beta_np, ff);
}
/* ---------------------------------------------------------------------- */
void FixPIMDBLangevin::spring_force() {
void FixPIMDBLangevin::spring_force()
{
for (int i = 0; i < nbosons; i++) {
f_tag_order[i][0] = 0.0;
f_tag_order[i][1] = 0.0;
f_tag_order[i][2] = 0.0;
}
bosonic_exchange.spring_force(f_tag_order);
for (int i = 0; i < nbosons; i++) {
f_tag_order[i][0] = 0.0;
f_tag_order[i][1] = 0.0;
f_tag_order[i][2] = 0.0;
}
bosonic_exchange->spring_force(f_tag_order);
double** f = atom->f;
tagint* tag = atom->tag;
for (int i = 0; i < nbosons; i++) {
f[i][0] += f_tag_order[tag[i] - 1][0];
f[i][1] += f_tag_order[tag[i] - 1][1];
f[i][2] += f_tag_order[tag[i] - 1][2];
}
double **f = atom->f;
tagint *tag = atom->tag;
for (int i = 0; i < nbosons; i++) {
f[i][0] += f_tag_order[tag[i] - 1][0];
f[i][1] += f_tag_order[tag[i] - 1][1];
f[i][2] += f_tag_order[tag[i] - 1][2];
}
}
/* ---------------------------------------------------------------------- */
void FixPIMDBLangevin::compute_spring_energy() {
se_bead = bosonic_exchange.get_bead_spring_energy();
void FixPIMDBLangevin::compute_spring_energy()
{
se_bead = bosonic_exchange->get_bead_spring_energy();
if (synch_energies) {
MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
total_spring_energy /= universe->procs_per_world[universe->iworld];
}
else {
total_spring_energy = 0;
}
if (synch_energies) {
MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
total_spring_energy /= universe->procs_per_world[universe->iworld];
} else {
total_spring_energy = 0;
}
}
/* ---------------------------------------------------------------------- */
void FixPIMDBLangevin::compute_t_prim()
{
if (synch_energies) {
double prim = bosonic_exchange.prim_estimator();
MPI_Allreduce(&prim, &t_prim, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
}
else {
t_prim = bosonic_exchange.prim_estimator();
}
if (synch_energies) {
double prim = bosonic_exchange->prim_estimator();
MPI_Allreduce(&prim, &t_prim, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
} else {
t_prim = bosonic_exchange->prim_estimator();
}
}
/* ---------------------------------------------------------------------- */
double FixPIMDBLangevin::compute_vector(int n)
{
if (0 <= n && n < 6) {
return FixPIMDLangevin::compute_vector(n);
}
else {
error->universe_all(FLERR, "Fix only has 6 outputs!");
}
return 0.0;
if (0 <= n && n < 6) {
return FixPIMDLangevin::compute_vector(n);
} else {
error->universe_all(FLERR, "Fix only has 6 outputs!");
}
return 0.0;
}

43
src/REPLICA/fix_pimdb_langevin.h
View File

@ -12,43 +12,44 @@
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
FixStyle(pimdb/langevin,FixPIMDBLangevin);
// clang-format off
FixStyle(pimdb/langevin, FixPIMDBLangevin);
// clang-format on
#else
#ifndef FIX_PIMDB_LANGEVIN_H
#define FIX_PIMDB_LANGEVIN_H
#include "fix_pimd_langevin.h"
#include "bosonic_exchange.h"
namespace LAMMPS_NS {
class FixPIMDBLangevin : public FixPIMDLangevin {
public:
FixPIMDBLangevin(class LAMMPS *, int, char **);
~FixPIMDBLangevin();
public:
FixPIMDBLangevin(class LAMMPS *, int, char **);
~FixPIMDBLangevin();
double compute_vector(int) override;
void compute_spring_energy() override;
void compute_t_prim() override;
double compute_vector(int) override;
void compute_spring_energy() override;
void compute_t_prim() override;
char** filtered_args;
char **filtered_args;
int filtered_narg;
protected:
void prepare_coordinates() override;
void spring_force() override;
protected:
void prepare_coordinates() override;
void spring_force() override;
private:
const int nbosons;
bool synch_energies;
BosonicExchange bosonic_exchange;
double** f_tag_order;
char** filter_args(int, char **); // for hold memory of filtered arguments when calling the parent constructor
private:
const int nbosons;
bool synch_energies;
class BosonicExchange *bosonic_exchange;
double **f_tag_order;
char **filter_args(
int, char **); // for hold memory of filtered arguments when calling the parent constructor
};
}
} // namespace LAMMPS_NS
#endif
#endif

47
src/REPLICA/fix_pimdb_nvt.cpp
View File

@ -22,6 +22,9 @@
------------------------------------------------------------------------- */
#include "fix_pimdb_nvt.h"
#include "bosonic_exchange.h"
#include "atom.h"
#include "error.h"
#include "force.h"
@ -31,22 +34,25 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) :
FixPIMDNVT(lmp, narg, arg),
bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, true)
FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) : FixPIMDNVT(lmp, narg, arg)
{
virial = 0.;
prim = 0.;
spring_energy = 0.;
bosonic_exchange = new BosonicExchange(lmp, atom->nlocal, np, universe->me, true, true);
virial = 0.0;
prim = 0.0;
spring_energy = 0.0;
size_vector = 4;
if (method != PIMD && method != NMPIMD) {
error->universe_all(FLERR, "Method not supported in fix pimdb/nvt; only methods PIMD and NMPIMD");
error->universe_all(FLERR,
"Method not supported in fix pimdb/nvt; only methods PIMD and NMPIMD");
}
}
/* ---------------------------------------------------------------------- */
FixPIMDBNVT::~FixPIMDBNVT() {
FixPIMDBNVT::~FixPIMDBNVT()
{
delete bosonic_exchange;
}
/* ---------------------------------------------------------------------- */
@ -54,8 +60,8 @@ FixPIMDBNVT::~FixPIMDBNVT() {
void FixPIMDBNVT::pre_spring_force_estimators()
{
FixPIMDNVT::pre_spring_force_estimators();
spring_energy = bosonic_exchange.get_bead_spring_energy();
prim = bosonic_exchange.prim_estimator();
spring_energy = bosonic_exchange->get_bead_spring_energy();
prim = bosonic_exchange->prim_estimator();
}
/* ---------------------------------------------------------------------- */
@ -67,30 +73,27 @@ void FixPIMDBNVT::prepare_coordinates()
double *xlast = buf_beads[x_last];
double *xnext = buf_beads[x_next];
double ff = fbond * atom->mass[atom->type[0]];
bosonic_exchange.prepare_with_coordinates(*x, xlast, xnext, beta, -ff);
bosonic_exchange->prepare_with_coordinates(*x, xlast, xnext, beta, -ff);
}
void FixPIMDBNVT::spring_force()
{
double **f = atom->f;
bosonic_exchange.spring_force(f);
bosonic_exchange->spring_force(f);
}
/* ---------------------------------------------------------------------- */
double FixPIMDBNVT::compute_vector(int n)
{
if (0 <= n && n < 3) {
return FixPIMDNVT::compute_vector(n);
}
if (0 <= n && n < 3) { return FixPIMDNVT::compute_vector(n); }
if (n == 3) {
return prim;
}
else {
error->universe_all(FLERR, "Fix only has 4 outputs!");
}
if (n == 3) {
return prim;
} else {
error->universe_all(FLERR, "Fix only has 4 outputs!");
}
return 0.0;
return 0.0;
}

17
src/REPLICA/fix_pimdb_nvt.h
View File

@ -21,24 +21,23 @@ FixStyle(pimdb/nvt,FixPIMDBNVT);
#define FIX_PIMDB_NVT_H
#include "fix_pimd_nvt.h"
#include "bosonic_exchange.h"
namespace LAMMPS_NS {
class FixPIMDBNVT : public FixPIMDNVT {
public:
FixPIMDBNVT(class LAMMPS *, int, char **);
~FixPIMDBNVT();
double compute_vector(int) override;
FixPIMDBNVT(class LAMMPS *, int, char **);
~FixPIMDBNVT();
double compute_vector(int) override;
protected:
void prepare_coordinates() override;
void spring_force() override;
void pre_spring_force_estimators() override;
void prepare_coordinates() override;
void spring_force() override;
void pre_spring_force_estimators() override;
private:
BosonicExchange bosonic_exchange;
double prim;
class BosonicExchange *bosonic_exchange;
double prim;
};
} // namespace LAMMPS_NS