From 0448651a9087d2d28c5171876c5b8ce2948769eb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:05:53 -0500 Subject: [PATCH 1/7] reformatting and correcting doc page issues and appy spelling fixes and false positives --- doc/src/fix_pimd.rst | 50 +++--- doc/src/fix_pimdb.rst | 165 ++++++++++++-------- doc/utils/sphinx-config/false_positives.txt | 16 ++ 3 files changed, 145 insertions(+), 86 deletions(-) diff --git a/doc/src/fix_pimd.rst b/doc/src/fix_pimd.rst index a054489f8c..fb9e4ad938 100644 --- a/doc/src/fix_pimd.rst +++ b/doc/src/fix_pimd.rst @@ -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. ---------- diff --git a/doc/src/fix_pimdb.rst b/doc/src/fix_pimdb.rst index ed77f934ea..3a2dc2492c 100644 --- a/doc/src/fix_pimdb.rst +++ b/doc/src/fix_pimdb.rst @@ -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 ` 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 ` 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) `: +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) +`: .. 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) `: +: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) `: .. 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) `. -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) `. -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) `. +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) `. 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) `. 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) `. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -The use of :doc:`binary restart files ` and :doc:`fix_modify ` is the same as in :doc:`fix pimd `. +The use of :doc:`binary restart files ` and :doc:`fix_modify +` is the same as in :doc:`fix pimd `. Fix *pimdb/nvt* computes a global 4-vector, which can be accessed by various :doc:`output commands `. 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) `. -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 `. The quantities in the global vector are: +Fix *pimdb/langevin* computes a global 6-vector, which can be accessed +by various :doc:`output commands `. 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 `. The quantities #. primitive kinetic energy estimator :ref:`(Hirshberg1) ` #. virial energy estimator :ref:`(Herman) ` (see the justification in the supporting information of :ref:`(Hirshberg2) `). -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 ` 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: diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index f94d528cfe..fb580eafac 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -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 From f29a433fa9c4978e1cb238b220f3deb743856436 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:12:36 -0500 Subject: [PATCH 2/7] apply clang-format and follow LAMMPS programming conventions more closely --- src/REPLICA/bosonic_exchange.cpp | 516 ++++++++++++++++--------------- src/REPLICA/bosonic_exchange.h | 106 +++---- 2 files changed, 320 insertions(+), 302 deletions(-) diff --git a/src/REPLICA/bosonic_exchange.cpp b/src/REPLICA/bosonic_exchange.cpp index 396c551509..53bfddd262 100644 --- a/src/REPLICA/bosonic_exchange.cpp +++ b/src/REPLICA/bosonic_exchange.cpp @@ -24,339 +24,357 @@ ------------------------------------------------------------------------- */ #include "bosonic_exchange.h" + #include "domain.h" -#include "universe.h" -#include "memory.h" #include "error.h" -#include +#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_pref(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::max(); + for (int m = 1; m < nbosons + 1; m++) { + double Elongest = std::numeric_limits::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::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::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::max(); + double Elongest = std::numeric_limits::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]); } /* ---------------------------------------------------------------------- */ diff --git a/src/REPLICA/bosonic_exchange.h b/src/REPLICA/bosonic_exchange.h index 4a4b55d880..58e2148d24 100644 --- a/src/REPLICA/bosonic_exchange.h +++ b/src/REPLICA/bosonic_exchange.h @@ -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 From 35fca290fce1b9502fb0099dfe1511b87bc8548e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:17:48 -0500 Subject: [PATCH 3/7] enumerators must be defined within the scope of the class, not globally --- src/REPLICA/bosonic_exchange.cpp | 2 +- src/REPLICA/fix_pimd_langevin.cpp | 10 +++++++--- src/REPLICA/fix_pimd_langevin.h | 18 +++++++++--------- src/REPLICA/fix_pimd_nvt.h | 5 ++--- 4 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/REPLICA/bosonic_exchange.cpp b/src/REPLICA/bosonic_exchange.cpp index 53bfddd262..4491369013 100644 --- a/src/REPLICA/bosonic_exchange.cpp +++ b/src/REPLICA/bosonic_exchange.cpp @@ -35,7 +35,7 @@ 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), x(nullptr), x_pref(nullptr), x_next(nullptr), + 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) { diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index a8c1ee9fe7..711ce4989c 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -53,8 +53,12 @@ using MathConst::MY_SQRT2; using MathConst::THIRD; using MathSpecial::powint; -static std::map Barostats{{MTTK, "MTTK"}, {BZP, "BZP"}}; -static std::map Ensembles{{NVE, "NVE"}, {NVT, "NVT"}, {NPH, "NPH"}, {NPT, "NPT"}}; +static std::map Barostats{{FixPIMDLangevin::MTTK, "MTTK"}, + {FixPIMDLangevin::BZP, "BZP"}}; +static std::map 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)); } } diff --git a/src/REPLICA/fix_pimd_langevin.h b/src/REPLICA/fix_pimd_langevin.h index f2447597f6..4f1bac2532 100644 --- a/src/REPLICA/fix_pimd_langevin.h +++ b/src/REPLICA/fix_pimd_langevin.h @@ -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; diff --git a/src/REPLICA/fix_pimd_nvt.h b/src/REPLICA/fix_pimd_nvt.h index 5b3ac1cde2..a0bff062ab 100644 --- a/src/REPLICA/fix_pimd_nvt.h +++ b/src/REPLICA/fix_pimd_nvt.h @@ -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; From 16b2988106c289ddc5741b5dd291a378d02313aa Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:32:44 -0500 Subject: [PATCH 4/7] make filter_args() fully C++ compatible --- src/REPLICA/fix_pimdb_langevin.cpp | 32 ++++++++++++++++-------------- src/REPLICA/fix_pimdb_langevin.h | 1 + 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/src/REPLICA/fix_pimdb_langevin.cpp b/src/REPLICA/fix_pimdb_langevin.cpp index d0e6ebec08..c750598813 100644 --- a/src/REPLICA/fix_pimdb_langevin.cpp +++ b/src/REPLICA/fix_pimdb_langevin.cpp @@ -43,8 +43,8 @@ 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), bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, false) { synch_energies = true; // Loop over the arguments with i++ instead of i+=2 like the parent, @@ -88,25 +88,27 @@ FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -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); } /* ---------------------------------------------------------------------- */ -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; } /* ---------------------------------------------------------------------- */ diff --git a/src/REPLICA/fix_pimdb_langevin.h b/src/REPLICA/fix_pimdb_langevin.h index fd7cace7f8..5034fd8673 100644 --- a/src/REPLICA/fix_pimdb_langevin.h +++ b/src/REPLICA/fix_pimdb_langevin.h @@ -35,6 +35,7 @@ public: void compute_t_prim() override; char** filtered_args; + int filtered_narg; protected: void prepare_coordinates() override; From f1ef94aadef132a06e025e262f4e9f75ea7c5e53 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:33:26 -0500 Subject: [PATCH 5/7] apply clang-format --- src/REPLICA/fix_pimdb_langevin.cpp | 180 ++++++++++++++--------------- 1 file changed, 90 insertions(+), 90 deletions(-) diff --git a/src/REPLICA/fix_pimdb_langevin.cpp b/src/REPLICA/fix_pimdb_langevin.cpp index c750598813..748fb90cdf 100644 --- a/src/REPLICA/fix_pimdb_langevin.cpp +++ b/src/REPLICA/fix_pimdb_langevin.cpp @@ -22,21 +22,21 @@ Version 1.0 ------------------------------------------------------------------------- */ -#include -#include -#include -#include -#include #include "fix_pimdb_langevin.h" -#include "universe.h" -#include "comm.h" -#include "force.h" + #include "atom.h" +#include "comm.h" #include "domain.h" -#include "update.h" -#include "memory.h" #include "error.h" -#include +#include "force.h" +#include "memory.h" +#include "universe.h" +#include "update.h" + +#include +#include +#include + using namespace LAMMPS_NS; using namespace FixConst; @@ -46,44 +46,46 @@ FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) : FixPIMDLangevin(lmp, narg, filtered_args = filter_args(narg, arg)), filtered_narg(narg), nbosons(atom->nlocal), bosonic_exchange(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!"); - } - } + 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"); } /* ---------------------------------------------------------------------- */ @@ -115,74 +117,72 @@ char **FixPIMDBLangevin::filter_args(int narg, char **arg) 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; } From f4275ae44cc90a374107ab6a687eec353557fb8c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:41:33 -0500 Subject: [PATCH 6/7] PIMPLify access to BosonicExchange class --- src/REPLICA/fix_pimdb_langevin.cpp | 19 ++++++++++++------- src/REPLICA/fix_pimdb_langevin.h | 3 +-- src/REPLICA/fix_pimdb_nvt.cpp | 18 +++++++++++------- src/REPLICA/fix_pimdb_nvt.h | 3 +-- 4 files changed, 25 insertions(+), 18 deletions(-) diff --git a/src/REPLICA/fix_pimdb_langevin.cpp b/src/REPLICA/fix_pimdb_langevin.cpp index 748fb90cdf..c2d7e38f2d 100644 --- a/src/REPLICA/fix_pimdb_langevin.cpp +++ b/src/REPLICA/fix_pimdb_langevin.cpp @@ -24,6 +24,8 @@ #include "fix_pimdb_langevin.h" +#include "bosonic_exchange.h" + #include "atom.h" #include "comm.h" #include "domain.h" @@ -44,9 +46,11 @@ using namespace FixConst; FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) : FixPIMDLangevin(lmp, narg, filtered_args = filter_args(narg, arg)), filtered_narg(narg), - nbosons(atom->nlocal), bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, false) + nbosons(atom->nlocal) { + 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++) { @@ -95,6 +99,7 @@ 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; } /* ---------------------------------------------------------------------- */ @@ -124,8 +129,8 @@ void FixPIMDBLangevin::prepare_coordinates() 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); } /* ---------------------------------------------------------------------- */ @@ -138,7 +143,7 @@ void FixPIMDBLangevin::spring_force() f_tag_order[i][1] = 0.0; f_tag_order[i][2] = 0.0; } - bosonic_exchange.spring_force(f_tag_order); + bosonic_exchange->spring_force(f_tag_order); double **f = atom->f; tagint *tag = atom->tag; @@ -153,7 +158,7 @@ void FixPIMDBLangevin::spring_force() void FixPIMDBLangevin::compute_spring_energy() { - se_bead = bosonic_exchange.get_bead_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); @@ -168,10 +173,10 @@ void FixPIMDBLangevin::compute_spring_energy() void FixPIMDBLangevin::compute_t_prim() { if (synch_energies) { - double prim = bosonic_exchange.prim_estimator(); + 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(); + t_prim = bosonic_exchange->prim_estimator(); } } diff --git a/src/REPLICA/fix_pimdb_langevin.h b/src/REPLICA/fix_pimdb_langevin.h index 5034fd8673..2917a045d1 100644 --- a/src/REPLICA/fix_pimdb_langevin.h +++ b/src/REPLICA/fix_pimdb_langevin.h @@ -21,7 +21,6 @@ FixStyle(pimdb/langevin,FixPIMDBLangevin); #define FIX_PIMDB_LANGEVIN_H #include "fix_pimd_langevin.h" -#include "bosonic_exchange.h" namespace LAMMPS_NS { @@ -44,9 +43,9 @@ protected: 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 + class BosonicExchange *bosonic_exchange; }; } diff --git a/src/REPLICA/fix_pimdb_nvt.cpp b/src/REPLICA/fix_pimdb_nvt.cpp index 8c5b1eb107..43510ddcbc 100644 --- a/src/REPLICA/fix_pimdb_nvt.cpp +++ b/src/REPLICA/fix_pimdb_nvt.cpp @@ -22,6 +22,9 @@ ------------------------------------------------------------------------- */ #include "fix_pimdb_nvt.h" + +#include "bosonic_exchange.h" + #include "atom.h" #include "error.h" #include "force.h" @@ -31,10 +34,10 @@ 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) + { + bosonic_exchange = new BosonicExchange(lmp, atom->nlocal, np, universe->me, true, true); virial = 0.; prim = 0.; spring_energy = 0.; @@ -47,6 +50,7 @@ FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ FixPIMDBNVT::~FixPIMDBNVT() { + delete bosonic_exchange; } /* ---------------------------------------------------------------------- */ @@ -54,8 +58,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,14 +71,14 @@ 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); } /* ---------------------------------------------------------------------- */ diff --git a/src/REPLICA/fix_pimdb_nvt.h b/src/REPLICA/fix_pimdb_nvt.h index 677dd4c7fb..b0757089ee 100644 --- a/src/REPLICA/fix_pimdb_nvt.h +++ b/src/REPLICA/fix_pimdb_nvt.h @@ -21,7 +21,6 @@ FixStyle(pimdb/nvt,FixPIMDBNVT); #define FIX_PIMDB_NVT_H #include "fix_pimd_nvt.h" -#include "bosonic_exchange.h" namespace LAMMPS_NS { @@ -37,7 +36,7 @@ class FixPIMDBNVT : public FixPIMDNVT { void pre_spring_force_estimators() override; private: - BosonicExchange bosonic_exchange; + class BosonicExchange *bosonic_exchange; double prim; }; From c7e288eb8c169cc56713e908925b5526276b7bb9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 23 Feb 2025 12:46:34 -0500 Subject: [PATCH 7/7] apply clang-format --- src/REPLICA/fix_pimdb_langevin.h | 41 ++++++++++++++++---------------- src/REPLICA/fix_pimdb_nvt.cpp | 29 +++++++++++----------- src/REPLICA/fix_pimdb_nvt.h | 16 ++++++------- 3 files changed, 43 insertions(+), 43 deletions(-) diff --git a/src/REPLICA/fix_pimdb_langevin.h b/src/REPLICA/fix_pimdb_langevin.h index 2917a045d1..a1efe94e5f 100644 --- a/src/REPLICA/fix_pimdb_langevin.h +++ b/src/REPLICA/fix_pimdb_langevin.h @@ -12,9 +12,9 @@ ------------------------------------------------------------------------- */ #ifdef FIX_CLASS - -FixStyle(pimdb/langevin,FixPIMDBLangevin); - +// clang-format off +FixStyle(pimdb/langevin, FixPIMDBLangevin); +// clang-format on #else #ifndef FIX_PIMDB_LANGEVIN_H @@ -25,30 +25,31 @@ FixStyle(pimdb/langevin,FixPIMDBLangevin); 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; - int filtered_narg; + 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; - 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 diff --git a/src/REPLICA/fix_pimdb_nvt.cpp b/src/REPLICA/fix_pimdb_nvt.cpp index 43510ddcbc..431b0d3a70 100644 --- a/src/REPLICA/fix_pimdb_nvt.cpp +++ b/src/REPLICA/fix_pimdb_nvt.cpp @@ -38,18 +38,20 @@ FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) : FixPIMDNVT(lmp, na { bosonic_exchange = new BosonicExchange(lmp, atom->nlocal, np, universe->me, true, true); - virial = 0.; - prim = 0.; - spring_energy = 0.; + 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; } @@ -85,16 +87,13 @@ void FixPIMDBNVT::spring_force() 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; } diff --git a/src/REPLICA/fix_pimdb_nvt.h b/src/REPLICA/fix_pimdb_nvt.h index b0757089ee..72db4cbe7f 100644 --- a/src/REPLICA/fix_pimdb_nvt.h +++ b/src/REPLICA/fix_pimdb_nvt.h @@ -26,18 +26,18 @@ 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: - class BosonicExchange *bosonic_exchange; - double prim; + class BosonicExchange *bosonic_exchange; + double prim; }; } // namespace LAMMPS_NS