diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 6192029a0a..c173216887 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -162,6 +162,8 @@ OPT. * :doc:`phonon ` * :doc:`pimd/langevin ` * :doc:`pimd/nvt ` + * :doc:`pimd/langevin/bosonic ` + * :doc:`pimd/nvt/bosonic ` * :doc:`planeforce ` * :doc:`plumed ` * :doc:`poems ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 49fbb72726..1d508b4ed4 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -341,6 +341,8 @@ accelerated styles exist. * :doc:`phonon ` - calculate dynamical matrix from MD simulations * :doc:`pimd/langevin ` - Feynman path-integral molecular dynamics with stochastic thermostat * :doc:`pimd/nvt ` - Feynman path-integral molecular dynamics with Nose-Hoover thermostat +* :doc:`pimd/langevin/bosonic ` - Bosonic Feynman path-integral molecular dynamics for with stochastic thermostat +* :doc:`pimd/nvt/bosonic ` - Bosonic Feynman path-integral molecular dynamics with Nose-Hoover thermostat * :doc:`planeforce ` - constrain atoms to move in a plane * :doc:`plumed ` - wrapper on PLUMED free energy library * :doc:`poems ` - constrain clusters of atoms to move as coupled rigid bodies diff --git a/doc/src/fix_ave_correlate.rst b/doc/src/fix_ave_correlate.rst index 659e15105c..91f5c701aa 100644 --- a/doc/src/fix_ave_correlate.rst +++ b/doc/src/fix_ave_correlate.rst @@ -32,13 +32,14 @@ Syntax .. parsed-literal:: - *type* arg = *auto* or *upper* or *lower* or *auto/upper* or *auto/lower* or *full* + *type* arg = *auto* or *upper* or *lower* or *auto/upper* or *auto/lower* or *full* or *first* auto = correlate each value with itself upper = correlate each value with each succeeding value lower = correlate each value with each preceding value auto/upper = auto + upper auto/lower = auto + lower full = correlate each value with every other value, including itself = auto + upper + lower + first = correlate each value with the first value *ave* args = *one* or *running* one = zero the correlation accumulation every Nfreq steps running = accumulate correlations continuously @@ -257,6 +258,9 @@ time. * If *type* is set to *full* then each input value is correlated with itself and every other value (i.e., :math:`C_{ij} = V_i V_j` for :math:`\{i,j\} = \{1,N\}`, so :math:`N_\text{pair} = N^2`). +* If *type* is set to *first* then each input value is correlated with + the first input value (i.e., :math:`C_{ij} = V_1 V_j` for + :math:`\{j\} = \{1,N\}`, so :math:`N_\text{pair} = N`). The *ave* keyword determines what happens to the accumulation of correlation samples every :math:`N_\text{freq}` timesteps. If the *ave* setting is *one*, @@ -369,6 +373,8 @@ above. * For *type* = *full*, the :math:`N_\text{pair} = N^2` columns are ordered: :math:`C_{11}, C_{12}, \dotsc, C_{1N}, C_{21}, C_{22}, \dotsc, C_{2N}, C_{31}, \dotsc, C_{3N}, \dotsc, C_{N1}, \dotsc, C_{N,N-1}, C_{NN}` +* For *type* = *first*, the :math:`N_\text{pair} = N` columns are ordered: + :math:`C_{11}, C_{12}, \dotsc, C_{1N}` The array values calculated by this fix are treated as extensive. If you need to divide them by the number of atoms, you must do this in a diff --git a/doc/src/fix_ave_correlate_long.rst b/doc/src/fix_ave_correlate_long.rst index 003bdf897d..22fac89706 100644 --- a/doc/src/fix_ave_correlate_long.rst +++ b/doc/src/fix_ave_correlate_long.rst @@ -31,13 +31,14 @@ Syntax .. parsed-literal:: - *type* arg = *auto* or *upper* or *lower* or *auto/upper* or *auto/lower* or *full* + *type* arg = *auto* or *upper* or *lower* or *auto/upper* or *auto/lower* or *full* or *first* auto = correlate each value with itself upper = correlate each value with each succeeding value lower = correlate each value with each preceding value auto/upper = auto + upper auto/lower = auto + lower full = correlate each value with every other value, including itself = auto + upper + lower + first = correlate each value with the first value *start* args = Nstart Nstart = start accumulating correlations on this time step *file* arg = filename diff --git a/doc/src/fix_pimd.rst b/doc/src/fix_pimd.rst index 0c7f763ced..f29124e9aa 100644 --- a/doc/src/fix_pimd.rst +++ b/doc/src/fix_pimd.rst @@ -1,5 +1,7 @@ .. index:: fix pimd/langevin .. index:: fix pimd/nvt +.. index:: fix pimd/langevin/bosonic +.. index:: fix pimd/nvt/bosonic fix pimd/langevin command ========================= @@ -7,6 +9,12 @@ fix pimd/langevin command fix pimd/nvt command ==================== +:doc:`fix pimd/langevin/bosonic ` command +=========================================================== + +:doc:`fix pimd/nvt/bosonic ` command +====================================================== + Syntax """""" @@ -20,34 +28,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* - *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) @@ -62,6 +73,8 @@ Examples fix 1 all pimd/nvt method nmpimd fmass 1.0 sp 2.0 temp 300.0 nhc 4 fix 1 all pimd/langevin ensemble npt integrator obabo temp 113.15 thermostat PILE_L 1234 tau 1.0 iso 1.0 barostat BZP taup 1.0 +Example input files are provided in the examples/PACKAGES/pimd directory. + Description """"""""""" @@ -76,7 +89,14 @@ partition function for the original system to a classical partition function for a ring-polymer system is exploited, to efficiently sample configurations from the canonical ensemble :ref:`(Feynman) `. -The classical partition function and its components are given +.. versionadded:: 11Mar2025 + + Fix *pimd/langevin/bosonic* and *pimd/nvt/bosonic* were added. + +Fix *pimd/nvt* and fix *pimd/langevin* simulate *distinguishable* quantum particles. +Simulations of bosons, including exchange effects, are supported with the :doc:`fix pimd/langevin/bosonic ` and :doc:`fix pimd/nvt/bosonic ` commands. + +For distinguishable particles, the isomorphic classical partition function and its components are given by the following equations: .. math:: @@ -222,7 +242,7 @@ be *PILE_L* (path integral Langevin equation local thermostat, as described in : The keyword *tau* specifies the thermostat damping time parameter for fix style *pimd/langevin*. It is in time unit. It only works on the centroid mode. The keyword *scale* specifies a scaling parameter for the damping times of the non-centroid modes for fix style *pimd/langevin*. The default -damping time of the non-centroid mode :math:`i` is :math:`\frac{P}{\beta\hbar}\sqrt{\lambda_i\times\mathrm{fmass}}` (*fmmode* is *physical*) or :math:`\frac{P}{\beta\hbar}\sqrt{\mathrm{fmass}}` (*fmmode* is *normal*). The damping times of all non-centroid modes are the default values divided by *scale*. +damping time of the non-centroid mode :math:`i` is :math:`\frac{P}{\beta\hbar}\sqrt{\lambda_i\times\mathrm{fmass}}` (*fmmode* is *physical*) or :math:`\frac{P}{\beta\hbar}\sqrt{\mathrm{fmass}}` (*fmmode* is *normal*). The damping times of all non-centroid modes are the default values divided by *scale*. This keyword should be used only with *method*=*nmpimd*. The barostat parameters for fix style *pimd/langevin* with *npt* or *nph* ensemble is specified using one of *iso* and *aniso* keywords. A *pressure* value should be given with pressure unit. The keyword *iso* means couple all 3 diagonal components together when pressure is computed (hydrostatic pressure), and dilate/contract the dimensions together. The keyword *aniso* means x, y, and z dimensions are controlled independently using the Pxx, Pyy, and Pzz components of the stress tensor as the driving forces, and the specified scalar external pressure. @@ -334,8 +354,8 @@ it outputs multiple log files, and different log files contain information about different beads or modes (see detailed explanations below). If *ensemble* is *nve* or *nvt*, the vector has 10 values: - #. kinetic energy of the normal mode - #. spring elastic energy of the normal mode + #. kinetic energy of the bead (if *method*=*pimd*) or normal mode (if *method*=*nmpimd*) + #. spring elastic energy of the bead (if *method*=*pimd*) or normal mode (if *method*=*nmpimd*) #. potential energy of the bead #. total energy of all beads (conserved if *ensemble* is *nve*) #. primitive kinetic energy estimator @@ -412,12 +432,20 @@ variable, e.g. velocity all create 300.0 1234${ibead} rot yes dist gaussian +Related commands +"""""""""""""""" + +:doc:`pimd/nvt/bosonic `, :doc:`pimd/langevin/bosonic ` + Default """"""" The keyword defaults for fix *pimd/nvt* are method = pimd, fmass = 1.0, sp = 1.0, temp = 300.0, and nhc = 2. +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. + ---------- .. _Feynman: diff --git a/doc/src/fix_pimd_bosonic.rst b/doc/src/fix_pimd_bosonic.rst new file mode 100644 index 0000000000..67dd60ed39 --- /dev/null +++ b/doc/src/fix_pimd_bosonic.rst @@ -0,0 +1,237 @@ +.. index:: fix pimd/langevin/bosonic +.. index:: fix pimd/nvt/bosonic + +fix pimd/langevin/bosonic command +================================= + +fix pimd/nvt/bosonic command +============================ + +Syntax +"""""" + +.. code-block:: LAMMPS + + fix ID group-ID style keyword value ... + +* ID, group-ID are documented in :doc:`fix ` command +* style = *pimd/langevin/bosonic* or *pimd/nvt/bosonic* = style name of this fix command +* zero or more keyword/value pairs may be appended +* keywords for style *pimd/nvt/bosonic* + + .. parsed-literal:: + + *keywords* = *method* or *fmass* or *sp* or *temp* or *nhc* + *method* value = *pimd* or *nmpimd* + *fmass* value = scaling factor on mass + *sp* value = scaling factor on Planck constant + *temp* value = temperature (temperature units) + *nhc* value = Nc = number of chains in Nose-Hoover thermostat + +* keywords for style *pimd/langevin/bosonic* + + .. parsed-literal:: + + *keywords* = *integrator* or *ensemble* or *fmass* or *temp* or *thermostat* or *tau* or *fixcom* or *lj* or *esych* + *integrator* value = *obabo* or *baoab* + *ensemble* value = *nvt* or *nve* + *fmass* value = scaling factor on mass + *temp* value = temperature (temperature unit) + temperature = target temperature of the thermostat + *thermostat* values = style seed + style value = *PILE_L* + seed = random number generator seed + *tau* value = thermostat damping parameter (time unit) + *fixcom* value = *yes* or *no* + *lj* values = epsilon sigma mass planck mvv2e + epsilon = energy scale for reduced units (energy units) + sigma = length scale for reduced units (length units) + mass = mass scale for reduced units (mass units) + planck = Planck's constant for other unit style + mvv2e = mass * velocity^2 to energy conversion factor for other unit style + *esynch* value = *yes* or *no* + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all pimd/nvt/bosonic method pimd fmass 1.0 sp 1.0 temp 2.0 nhc 4 + fix 1 all pimd/langevin/bosonic integrator obabo temp 113.15 thermostat PILE_L 1234 tau 1.0 + +Example input files are provided in the examples/PACKAGES/pimd_bosonic directory. + +Description +""""""""""" + +These fix commands are based on the fixes :doc:`pimd/nvt and +pimd/langevin ` for performing quantum molecular dynamics +simulations based on the Feynman path-integral formalism. The key +difference is that fix *pimd/nvt* and fix *pimd/langevin* simulate +*distinguishable* particles, while fix *pimd/nvt/bosonic* and fix +*pimd/langevin/bosonic* perform simulations of bosons, including exchange +effects. The *bosonic* commands share syntax with the equivalent commands for distinguishable particles. +The user is referred to the documentation of :doc:`these commands ` +for a detailed syntax description and additional, general capabilities +of the commands. The major differences from fix *pimd/nvt* and fix *pimd/langevin* in terms of +capabilities are: + +* Fix *pimd/nvt/bosonic* only supports the "pimd" and "nmpimd" methods. Fix + *pimd/langevin/bosonic* only supports the "pimd" method, which is the default + in this fix. These restrictions are related to the use of normal + modes, which change in bosons. For similar reasons, *fmmode* of + *pimd/langevin* should not be used, and would raise an error if set to + a value other than *physical*. +* Fix *pimd/langevin/bosonic* currently does not support *ensemble* other than + *nve*, *nvt*. The barostat related keywords *iso*, *aniso*, + *barostat*, *taup* are not supported. +* Fix *pimd/langevin/bosonic* also has a keyword not available in fix + *pimd/langevin*: *esynch*, with default *yes*. If set to *no*, some + time consuming synchronization of spring energies and the primitive + kinetic energy estimator between processors is avoided. + +The isomorphism between the partition function of :math:`N` bosonic +quantum particles and that of a system of classical ring polymers at +inverse temperature :math:`\beta` is given by :ref:`(Tuckerman) +`: + +.. 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: + +.. math:: + + E^\sigma = \frac{mP}{2\beta^2 \hbar^2} \sum_{\ell=1}^N \sum_{j=1}^P \left(\mathbf{q}_\ell^j - \mathbf{q}_\ell^{j+1}\right)^2, + +where :math:`P` is the number of beads and :math:`\mathbf{q}_\ell^{P+1}=\mathbf{q}_{\sigma(\ell)}^1.` + +Hirshberg et. al. showed that the ring polymer potential +:math:`-\frac{1}{\beta}\textrm{ln}\left[ \frac{1}{N!} \sum_\sigma e ^ { +-\beta E^\sigma } \right]`, which scales exponentially with :math:`N`, +can be replaced by a potential :math:`V^{[1,N]}` defined through a +recurrence relation :ref:`(Hirshberg1) `: + +.. 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) `. + +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 `. + +Fix *pimd/nvt/bosonic* computes a global 4-vector, which can be accessed by +various :doc:`output commands `. The quantities in +the global vector are: + + #. the total spring energy of the quasi-beads, + #. the current temperature of the classical system of ring polymers, + #. the current value of the scalar virial estimator for the kinetic + energy of the quantum system :ref:`(Herman) ` (see the justification in the supporting information of :ref:`(Hirshberg2) `), + #. the current value of the scalar primitive estimator for the kinetic + energy of the quantum system :ref:`(Hirshberg1) `. + +The vector values calculated by fix *pimd/nvt/bosonic* are "extensive", except +for the temperature, which is "intensive". + +Fix *pimd/langevin/bosonic* computes a global 6-vector, which can be accessed +by various :doc:`output commands `. The quantities in the +global vector are: + + #. kinetic energy of the beads, + #. spring elastic energy of the beads, + #. potential energy of the bead, + #. total energy of all beads (conserved if *ensemble* is *nve*) if *esynch* is *yes* + #. primitive kinetic energy estimator :ref:`(Hirshberg1) ` + #. 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 +*pimd/langevin/bosonic* are "extensive". + +For both *pimd/nvt/bosonic* and *pimd/langevin/bosonic*, the contribution of the +exterior spring to the primitive estimator is printed to the first log +file. The contribution of the :math:`P-1` interior springs is printed +to the other :math:`P-1` log files. The contribution of the constant +:math:`\frac{PdN}{2 \beta}` (with :math:`d` being the dimensionality) is +equally divided over log files. + +Restrictions +"""""""""""" + +These fixes are part of the REPLICA package. They are only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +The restrictions of :doc:`fix pimd ` apply. + +Related commands +"""""""""""""""" + +:doc:`pimd/nvt `, :doc:`pimd/langevin ` + +Default +""""""" + +The keyword defaults for fix *pimd/nvt/bosonic* are method = pimd, fmass = 1.0, +sp = 1.0, temp = 300.0, and nhc = 2. + +The keyword defaults for fix *pimd/langevin/bosonic* are integrator = obabo, +method = pimd, ensemble = nvt, fmass = 1.0, temp = 298.15, thermostat = +PILE_L, tau = 1.0, fixcom = yes, esynch = yes, and lj = 1 for all its +arguments. + +---------- + +.. _book-Tuckerman: + +**(Tuckerman)** M. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford University Press, 2010) + +.. _Hirshberg: + +**(Hirshberg1)** B. Hirshberg, V. Rizzi, and M. Parrinello, "Path integral molecular dynamics for bosons," Proc. Natl. Acad. Sci. U. S. A. 116, 21445 (2019) + +.. _HirshbergInvernizzi: + +**(Hirshberg2)** B. Hirshberg, M. Invernizzi, and M. Parrinello, "Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality," J Chem Phys, 152, 171102 (2020) + +.. _Feldman: + +**(Feldman)** Y. M. Y. Feldman and B. Hirshberg, "Quadratic scaling bosonic path integral molecular dynamics," J. Chem. Phys. 159, 154107 (2023) + +.. _HigerFeldman: + +**(Higer)** J. Higer, Y. M. Y. Feldman, and B. Hirshberg, "Periodic Boundary Conditions for Bosonic Path Integral Molecular Dynamics," arXiv:2501.17618 (2025) + +.. _HermanBB: + +**(Herman)** M. F. Herman, E. J. Bruskin, B. J. Berne, J Chem Phys, 76, 5150 (1982). diff --git a/doc/src/fix_qeq.rst b/doc/src/fix_qeq.rst index 06a1f98375..ac0bca84f1 100644 --- a/doc/src/fix_qeq.rst +++ b/doc/src/fix_qeq.rst @@ -304,7 +304,7 @@ Chemistry, 95, 3358-3363 (1991). .. _CTIP1: **(CTIP)** G. Plummer, J. P. Tavenner, M. I. Mendelev, Z. Wu, J. W. Lawson, -in preparation +J Chemical Physics, 162, 054709 (2025) .. _vanDuin: diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index c8c1581c3c..c673561237 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -471,7 +471,7 @@ Phys, 110, 8254 (1999). .. _Plummer1: **(Plummer)** G. Plummer, J. P. Tavenner, M. I. Mendelev, Z. Wu, J. W. Lawson, -in preparation +J Chemical Physics, 162, 054709 (2025). .. _Jorgensen3: diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index a6548dbbcd..084f7c0ec1 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 @@ -1485,6 +1495,7 @@ hgrid hhmrr Hibbs Higdon +Higer hiID Hijazi Hilger @@ -2610,6 +2621,7 @@ nn nnodes npits npj +nmpimd nO Nocedal nocite @@ -2738,6 +2750,7 @@ nylo nz Nz nzlo +obabo ochre ocl octahedral @@ -2938,6 +2951,7 @@ Pieniazek Pieter pIm pimd +pimdb Piola pIp pipelining @@ -3289,6 +3303,7 @@ RiRj Risi Rix Riy +Rizzi rj Rj Rjinner @@ -3722,6 +3737,7 @@ Tanmoy Tartakovsky taskset taubi +taup Tavenner taylor tb diff --git a/examples/PACKAGES/pimd_bosonic/README b/examples/PACKAGES/pimd_bosonic/README new file mode 100644 index 0000000000..9e8e3cf527 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/README @@ -0,0 +1,11 @@ +This directory contains example input files for the pimd/nvt/bosonic and pimd/langevin/bosonic fix commands. +The demonstrated system is that of 3 non-interacting bosons in a harmonic trap. +The input files initiate a simulation of 100 steps using 4 beads. The expected output log files are included. + +The input files can be edited to provide a physically meaningful output, by extending the simulations to a few nanoseconds +and increasing the number of beads (32 beads should be enough). +The total energy can be computed by summing the potential energy and either the virial or primitive kinetic energy estimator, and averaging over time. +Notice that for pimd/langevin/bosonic the summed value of either kinetic estimator is printed across all log files. +It should be taken from a single log file, with all files providing the same output. +For pimd/nvt/bosonic, one has to sum over the output of all log files. For the potential energy, summing over log files is required for both fix commands. +The obtained energy should be compared to the analytical value of 0.00058031 Hartree. diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/in.lmp b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/in.lmp new file mode 100644 index 0000000000..f5d984eed4 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/in.lmp @@ -0,0 +1,57 @@ +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/langevin/bosonic ensemble nvt temp ${Temp} thermostat PILE_L 12345 tau 50 fixcom no + +# Outputs +variable prim_kinetic equal f_pimdb[5] +variable virial equal f_pimdb[6] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.0 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.0 new file mode 100644 index 0000000000..2e5f1b0225 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.0 @@ -0,0 +1,236 @@ +LAMMPS (4 Feb 2025) +Processor partition = 0 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188891 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/langevin/bosonic ensemble nvt temp ${Temp} thermostat PILE_L 12345 tau 50 fixcom no +fix pimdb all pimd/langevin/bosonic ensemble nvt temp 17.4 thermostat PILE_L 12345 tau 50 fixcom no + +# Outputs +variable prim_kinetic equal f_pimdb[5] +variable virial equal f_pimdb[6] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +Initializing PI Langevin equation thermostat... + Bead ID | omega | tau | c1 | c2 + 0 9.11206647e-03 5.00000000e+01 9.95012479e-01 9.97505201e-02 + 1 9.11206647e-03 5.00000000e+01 9.95012479e-01 9.97505201e-02 + 2 9.11206647e-03 5.00000000e+01 9.95012479e-01 9.97505201e-02 + 3 9.11206647e-03 5.00000000e+01 9.95012479e-01 9.97505201e-02 + PILE_L thermostat successfully initialized! + +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 2.801 | 2.801 | 2.801 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 1.3661449e-08 0.0009918329 + 1 1.3158398e-09 1.8742641e-09 0.00099182267 + 2 5.0818396e-09 7.1801177e-09 0.00099180053 + 3 1.2446727e-08 1.4372388e-08 0.00099176038 + 4 2.1754941e-08 3.1274072e-08 0.00099170266 + 5 3.2927703e-08 5.2470231e-08 0.00099162511 + 6 4.7293056e-08 7.2301291e-08 0.00099153 + 7 6.526455e-08 9.5440683e-08 0.00099141194 + 8 8.6974168e-08 1.1669162e-07 0.00099127168 + 9 1.1249595e-07 1.3551097e-07 0.00099109846 + 10 1.4190482e-07 1.6050328e-07 0.00099089759 + 11 1.7313555e-07 1.8645705e-07 0.0009906735 + 12 2.0672874e-07 2.1181418e-07 0.00099042271 + 13 2.4042615e-07 2.4218941e-07 0.00099012073 + 14 2.7524432e-07 2.7507982e-07 0.00098977876 + 15 3.1513375e-07 3.1771012e-07 0.00098940169 + 16 3.6161087e-07 3.6075963e-07 0.00098900826 + 17 4.0813191e-07 4.0059005e-07 0.00098859985 + 18 4.5194598e-07 4.3883585e-07 0.00098817535 + 19 4.9898075e-07 4.8036287e-07 0.00098770764 + 20 5.5458728e-07 5.2362451e-07 0.00098721266 + 21 6.1220768e-07 5.7705077e-07 0.00098670407 + 22 6.6573849e-07 6.3271105e-07 0.00098617081 + 23 7.2475089e-07 6.9382916e-07 0.00098557406 + 24 7.804837e-07 7.5979853e-07 0.0009849612 + 25 8.3194987e-07 8.232213e-07 0.00098434971 + 26 8.8723792e-07 8.9901705e-07 0.0009836645 + 27 9.3837241e-07 9.7033126e-07 0.0009830068 + 28 9.897077e-07 1.0444697e-06 0.00098229311 + 29 1.043431e-06 1.1208219e-06 0.00098150213 + 30 1.0900368e-06 1.1899315e-06 0.00098077686 + 31 1.1303314e-06 1.2560026e-06 0.00098003518 + 32 1.1725968e-06 1.3287355e-06 0.00097921099 + 33 1.2081344e-06 1.399856e-06 0.0009784037 + 34 1.2449387e-06 1.4769272e-06 0.00097754491 + 35 1.2852285e-06 1.5642459e-06 0.00097667479 + 36 1.3234189e-06 1.6571388e-06 0.00097572641 + 37 1.3666391e-06 1.7603492e-06 0.00097475881 + 38 1.4102743e-06 1.8666516e-06 0.00097373763 + 39 1.454172e-06 1.9696572e-06 0.0009727081 + 40 1.5031422e-06 2.0910539e-06 0.00097155807 + 41 1.5538491e-06 2.2198868e-06 0.00097030871 + 42 1.6079193e-06 2.3576614e-06 0.00096896762 + 43 1.6637051e-06 2.496284e-06 0.00096760456 + 44 1.7189714e-06 2.637657e-06 0.00096619083 + 45 1.7656688e-06 2.795845e-06 0.00096470693 + 46 1.8092153e-06 2.9727162e-06 0.00096303835 + 47 1.8415769e-06 3.1451673e-06 0.00096140406 + 48 1.8664597e-06 3.3241907e-06 0.00095971129 + 49 1.8940033e-06 3.5094834e-06 0.00095790964 + 50 1.9211844e-06 3.7049195e-06 0.00095603039 + 51 1.9533838e-06 3.9139587e-06 0.0009540918 + 52 1.9833776e-06 4.1289537e-06 0.0009520767 + 53 2.0106231e-06 4.3481079e-06 0.00095003676 + 54 2.0429292e-06 4.5702968e-06 0.00094802906 + 55 2.0747842e-06 4.7928135e-06 0.00094600351 + 56 2.1043981e-06 5.0253448e-06 0.0009438419 + 57 2.1474343e-06 5.2569331e-06 0.00094169256 + 58 2.2115003e-06 5.4997767e-06 0.00093948055 + 59 2.2773037e-06 5.7547729e-06 0.00093714003 + 60 2.3402972e-06 6.0074589e-06 0.00093480857 + 61 2.4107157e-06 6.2647899e-06 0.00093244785 + 62 2.4953095e-06 6.5384938e-06 0.00092998209 + 63 2.58923e-06 6.8124372e-06 0.00092751453 + 64 2.682021e-06 7.0913777e-06 0.00092505171 + 65 2.7538688e-06 7.3773956e-06 0.00092256944 + 66 2.8444661e-06 7.6811232e-06 0.00091984092 + 67 2.9532286e-06 7.9951107e-06 0.0009170037 + 68 3.0551339e-06 8.3096758e-06 0.00091414977 + 69 3.1625164e-06 8.6232285e-06 0.00091137219 + 70 3.2717759e-06 8.9388929e-06 0.00090857909 + 71 3.3924986e-06 9.2524846e-06 0.00090583473 + 72 3.5236727e-06 9.5688611e-06 0.00090308342 + 73 3.6440162e-06 9.8873586e-06 0.00090026584 + 74 3.7692169e-06 1.0203245e-05 0.00089755256 + 75 3.8946044e-06 1.0523755e-05 0.00089485804 + 76 4.0498738e-06 1.0848398e-05 0.00089217189 + 77 4.2099346e-06 1.1164286e-05 0.00088956967 + 78 4.3589564e-06 1.1466652e-05 0.00088700311 + 79 4.5069481e-06 1.177381e-05 0.00088439373 + 80 4.6635354e-06 1.2059606e-05 0.00088199635 + 81 4.8398269e-06 1.233837e-05 0.0008797438 + 82 5.0191539e-06 1.261531e-05 0.00087745162 + 83 5.205285e-06 1.2899972e-05 0.00087506863 + 84 5.3708116e-06 1.3160246e-05 0.00087288444 + 85 5.5534416e-06 1.3438712e-05 0.00087052617 + 86 5.7613e-06 1.372546e-05 0.00086818634 + 87 5.9645662e-06 1.4017631e-05 0.00086572407 + 88 6.1830289e-06 1.4328439e-05 0.000862992 + 89 6.423646e-06 1.4640964e-05 0.00086017364 + 90 6.688219e-06 1.4978644e-05 0.0008570595 + 91 6.9305603e-06 1.5307214e-05 0.00085411601 + 92 7.1857573e-06 1.5649328e-05 0.00085112325 + 93 7.4653538e-06 1.6015833e-05 0.0008479391 + 94 7.8132911e-06 1.6393131e-05 0.00084474237 + 95 8.1540941e-06 1.674573e-05 0.00084165639 + 96 8.4943231e-06 1.7094435e-05 0.00083858996 + 97 8.8504733e-06 1.7451741e-05 0.00083536142 + 98 9.2042324e-06 1.780703e-05 0.00083221292 + 99 9.5058078e-06 1.8141862e-05 0.00082913227 + 100 9.8087647e-06 1.8457846e-05 0.00082619877 +Loop time of 0.122922 on 1 procs for 100 steps with 3 atoms + +Performance: 35144210.362 fs/day, 0.000 hours/fs, 813.523 timesteps/s, 2.441 katom-step/s +70.1% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 1.1593e-05 | 1.1593e-05 | 1.1593e-05 | 0.0 | 0.01 +Comm | 9.2183e-05 | 9.2183e-05 | 9.2183e-05 | 0.0 | 0.07 +Output | 0.023243 | 0.023243 | 0.023243 | 0.0 | 18.91 +Modify | 0.099386 | 0.099386 | 0.099386 | 0.0 | 80.85 +Other | | 0.0001896 | | | 0.15 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 33 +Dangerous builds = 0 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.1 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.1 new file mode 100644 index 0000000000..a00ac90be8 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.1 @@ -0,0 +1,228 @@ +LAMMPS (4 Feb 2025) +Processor partition = 1 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188892 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/langevin/bosonic ensemble nvt temp ${Temp} thermostat PILE_L 12345 tau 50 fixcom no +fix pimdb all pimd/langevin/bosonic ensemble nvt temp 17.4 thermostat PILE_L 12345 tau 50 fixcom no + +# Outputs +variable prim_kinetic equal f_pimdb[5] +variable virial equal f_pimdb[6] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 2.801 | 2.801 | 2.801 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 1.3661449e-08 0.0009918329 + 1 7.3810761e-10 1.8742641e-09 0.00099182267 + 2 3.0654593e-09 7.1801177e-09 0.00099180053 + 3 6.4534507e-09 1.4372388e-08 0.00099176038 + 4 1.1681892e-08 3.1274072e-08 0.00099170266 + 5 1.8404818e-08 5.2470231e-08 0.00099162511 + 6 2.6171374e-08 7.2301291e-08 0.00099153 + 7 3.5990441e-08 9.5440683e-08 0.00099141194 + 8 4.9770229e-08 1.1669162e-07 0.00099127168 + 9 6.6159035e-08 1.3551097e-07 0.00099109846 + 10 8.7281867e-08 1.6050328e-07 0.00099089759 + 11 1.1505206e-07 1.8645705e-07 0.0009906735 + 12 1.4564491e-07 2.1181418e-07 0.00099042271 + 13 1.8674942e-07 2.4218941e-07 0.00099012073 + 14 2.3275077e-07 2.7507982e-07 0.00098977876 + 15 2.833552e-07 3.1771012e-07 0.00098940169 + 16 3.3078875e-07 3.6075963e-07 0.00098900826 + 17 3.7783787e-07 4.0059005e-07 0.00098859985 + 18 4.3007763e-07 4.3883585e-07 0.00098817535 + 19 4.8592613e-07 4.8036287e-07 0.00098770764 + 20 5.4429898e-07 5.2362451e-07 0.00098721266 + 21 6.0140068e-07 5.7705077e-07 0.00098670407 + 22 6.6036274e-07 6.3271105e-07 0.00098617081 + 23 7.3748601e-07 6.9382916e-07 0.00098557406 + 24 8.1086859e-07 7.5979853e-07 0.0009849612 + 25 8.8638234e-07 8.232213e-07 0.00098434971 + 26 9.8252709e-07 8.9901705e-07 0.0009836645 + 27 1.0768657e-06 9.7033126e-07 0.0009830068 + 28 1.1810654e-06 1.0444697e-06 0.00098229311 + 29 1.2964318e-06 1.1208219e-06 0.00098150213 + 30 1.4097624e-06 1.1899315e-06 0.00098077686 + 31 1.5368367e-06 1.2560026e-06 0.00098003518 + 32 1.6777082e-06 1.3287355e-06 0.00097921099 + 33 1.8186759e-06 1.399856e-06 0.0009784037 + 34 1.9730383e-06 1.4769272e-06 0.00097754491 + 35 2.1285065e-06 1.5642459e-06 0.00097667479 + 36 2.2809575e-06 1.6571388e-06 0.00097572641 + 37 2.4393024e-06 1.7603492e-06 0.00097475881 + 38 2.6212592e-06 1.8666516e-06 0.00097373763 + 39 2.8115312e-06 1.9696572e-06 0.0009727081 + 40 3.0115212e-06 2.0910539e-06 0.00097155807 + 41 3.2325069e-06 2.2198868e-06 0.00097030871 + 42 3.4582526e-06 2.3576614e-06 0.00096896762 + 43 3.6838101e-06 2.496284e-06 0.00096760456 + 44 3.9267475e-06 2.637657e-06 0.00096619083 + 45 4.1942569e-06 2.795845e-06 0.00096470693 + 46 4.4694976e-06 2.9727162e-06 0.00096303835 + 47 4.7570486e-06 3.1451673e-06 0.00096140406 + 48 5.0917213e-06 3.3241907e-06 0.00095971129 + 49 5.4581889e-06 3.5094834e-06 0.00095790964 + 50 5.8607381e-06 3.7049195e-06 0.00095603039 + 51 6.2673583e-06 3.9139587e-06 0.0009540918 + 52 6.6764632e-06 4.1289537e-06 0.0009520767 + 53 7.0828028e-06 4.3481079e-06 0.00095003676 + 54 7.4858144e-06 4.5702968e-06 0.00094802906 + 55 7.8711706e-06 4.7928135e-06 0.00094600351 + 56 8.2756726e-06 5.0253448e-06 0.0009438419 + 57 8.6667628e-06 5.2569331e-06 0.00094169256 + 58 9.0585083e-06 5.4997767e-06 0.00093948055 + 59 9.4634766e-06 5.7547729e-06 0.00093714003 + 60 9.8696653e-06 6.0074589e-06 0.00093480857 + 61 1.0308345e-05 6.2647899e-06 0.00093244785 + 62 1.0757919e-05 6.5384938e-06 0.00092998209 + 63 1.1183568e-05 6.8124372e-06 0.00092751453 + 64 1.1585073e-05 7.0913777e-06 0.00092505171 + 65 1.1993605e-05 7.3773956e-06 0.00092256944 + 66 1.2382612e-05 7.6811232e-06 0.00091984092 + 67 1.2768116e-05 7.9951107e-06 0.0009170037 + 68 1.3181941e-05 8.3096758e-06 0.00091414977 + 69 1.3605472e-05 8.6232285e-06 0.00091137219 + 70 1.4027012e-05 8.9388929e-06 0.00090857909 + 71 1.4471706e-05 9.2524846e-06 0.00090583473 + 72 1.4896833e-05 9.5688611e-06 0.00090308342 + 73 1.5345003e-05 9.8873586e-06 0.00090026584 + 74 1.5779905e-05 1.0203245e-05 0.00089755256 + 75 1.6189898e-05 1.0523755e-05 0.00089485804 + 76 1.6654778e-05 1.0848398e-05 0.00089217189 + 77 1.7131363e-05 1.1164286e-05 0.00088956967 + 78 1.7572176e-05 1.1466652e-05 0.00088700311 + 79 1.8000817e-05 1.177381e-05 0.00088439373 + 80 1.838117e-05 1.2059606e-05 0.00088199635 + 81 1.8762432e-05 1.233837e-05 0.0008797438 + 82 1.9099052e-05 1.261531e-05 0.00087745162 + 83 1.9507863e-05 1.2899972e-05 0.00087506863 + 84 1.9885648e-05 1.3160246e-05 0.00087288444 + 85 2.0266439e-05 1.3438712e-05 0.00087052617 + 86 2.0612316e-05 1.372546e-05 0.00086818634 + 87 2.0959061e-05 1.4017631e-05 0.00086572407 + 88 2.13065e-05 1.4328439e-05 0.000862992 + 89 2.1634898e-05 1.4640964e-05 0.00086017364 + 90 2.2023174e-05 1.4978644e-05 0.0008570595 + 91 2.2431573e-05 1.5307214e-05 0.00085411601 + 92 2.2872283e-05 1.5649328e-05 0.00085112325 + 93 2.338476e-05 1.6015833e-05 0.0008479391 + 94 2.3912199e-05 1.6393131e-05 0.00084474237 + 95 2.437725e-05 1.674573e-05 0.00084165639 + 96 2.4828999e-05 1.7094435e-05 0.00083858996 + 97 2.5286762e-05 1.7451741e-05 0.00083536142 + 98 2.5769956e-05 1.780703e-05 0.00083221292 + 99 2.624134e-05 1.8141862e-05 0.00082913227 + 100 2.6731735e-05 1.8457846e-05 0.00082619877 +Loop time of 0.122878 on 1 procs for 100 steps with 3 atoms + +Performance: 35156789.883 fs/day, 0.000 hours/fs, 813.815 timesteps/s, 2.441 katom-step/s +46.7% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 1.9787e-05 | 1.9787e-05 | 1.9787e-05 | 0.0 | 0.02 +Comm | 9.2033e-05 | 9.2033e-05 | 9.2033e-05 | 0.0 | 0.07 +Output | 0.0022584 | 0.0022584 | 0.0022584 | 0.0 | 1.84 +Modify | 0.12033 | 0.12033 | 0.12033 | 0.0 | 97.93 +Other | | 0.0001755 | | | 0.14 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 44 +Dangerous builds = 0 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.2 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.2 new file mode 100644 index 0000000000..9cc0033c8a --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.2 @@ -0,0 +1,228 @@ +LAMMPS (4 Feb 2025) +Processor partition = 2 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188893 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/langevin/bosonic ensemble nvt temp ${Temp} thermostat PILE_L 12345 tau 50 fixcom no +fix pimdb all pimd/langevin/bosonic ensemble nvt temp 17.4 thermostat PILE_L 12345 tau 50 fixcom no + +# Outputs +variable prim_kinetic equal f_pimdb[5] +variable virial equal f_pimdb[6] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 2.801 | 2.801 | 2.801 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 1.3661449e-08 0.0009918329 + 1 8.3480048e-10 1.8742641e-09 0.00099182267 + 2 2.8201389e-09 7.1801177e-09 0.00099180053 + 3 7.3605781e-09 1.4372388e-08 0.00099176038 + 4 1.4967393e-08 3.1274072e-08 0.00099170266 + 5 2.4978087e-08 5.2470231e-08 0.00099162511 + 6 3.5908849e-08 7.2301291e-08 0.00099153 + 7 4.8226614e-08 9.5440683e-08 0.00099141194 + 8 6.3621045e-08 1.1669162e-07 0.00099127168 + 9 8.0119736e-08 1.3551097e-07 0.00099109846 + 10 9.7965614e-08 1.6050328e-07 0.00099089759 + 11 1.1538858e-07 1.8645705e-07 0.0009906735 + 12 1.365443e-07 2.1181418e-07 0.00099042271 + 13 1.5920175e-07 2.4218941e-07 0.00099012073 + 14 1.8252589e-07 2.7507982e-07 0.00098977876 + 15 2.0858909e-07 3.1771012e-07 0.00098940169 + 16 2.3676046e-07 3.6075963e-07 0.00098900826 + 17 2.689757e-07 4.0059005e-07 0.00098859985 + 18 3.0022172e-07 4.3883585e-07 0.00098817535 + 19 3.3127637e-07 4.8036287e-07 0.00098770764 + 20 3.6044366e-07 5.2362451e-07 0.00098721266 + 21 3.8484646e-07 5.7705077e-07 0.00098670407 + 22 4.1435067e-07 6.3271105e-07 0.00098617081 + 23 4.4508994e-07 6.9382916e-07 0.00098557406 + 24 4.7879333e-07 7.5979853e-07 0.0009849612 + 25 5.0918223e-07 8.232213e-07 0.00098434971 + 26 5.4733661e-07 8.9901705e-07 0.0009836645 + 27 5.7932115e-07 9.7033126e-07 0.0009830068 + 28 6.1775401e-07 1.0444697e-06 0.00098229311 + 29 6.6883535e-07 1.1208219e-06 0.00098150213 + 30 7.0943445e-07 1.1899315e-06 0.00098077686 + 31 7.4477897e-07 1.2560026e-06 0.00098003518 + 32 7.8183156e-07 1.3287355e-06 0.00097921099 + 33 8.1987667e-07 1.399856e-06 0.0009784037 + 34 8.5514962e-07 1.4769272e-06 0.00097754491 + 35 8.9607784e-07 1.5642459e-06 0.00097667479 + 36 9.4986753e-07 1.6571388e-06 0.00097572641 + 37 1.0024729e-06 1.7603492e-06 0.00097475881 + 38 1.0526731e-06 1.8666516e-06 0.00097373763 + 39 1.103419e-06 1.9696572e-06 0.0009727081 + 40 1.1640352e-06 2.0910539e-06 0.00097155807 + 41 1.2292837e-06 2.2198868e-06 0.00097030871 + 42 1.3024454e-06 2.3576614e-06 0.00096896762 + 43 1.3788736e-06 2.496284e-06 0.00096760456 + 44 1.4503248e-06 2.637657e-06 0.00096619083 + 45 1.5334106e-06 2.795845e-06 0.00096470693 + 46 1.6246777e-06 2.9727162e-06 0.00096303835 + 47 1.715884e-06 3.1451673e-06 0.00096140406 + 48 1.8096683e-06 3.3241907e-06 0.00095971129 + 49 1.916714e-06 3.5094834e-06 0.00095790964 + 50 2.0374848e-06 3.7049195e-06 0.00095603039 + 51 2.1747488e-06 3.9139587e-06 0.0009540918 + 52 2.3236882e-06 4.1289537e-06 0.0009520767 + 53 2.4790586e-06 4.3481079e-06 0.00095003676 + 54 2.6478101e-06 4.5702968e-06 0.00094802906 + 55 2.8275135e-06 4.7928135e-06 0.00094600351 + 56 3.0338638e-06 5.0253448e-06 0.0009438419 + 57 3.2428318e-06 5.2569331e-06 0.00094169256 + 58 3.4414566e-06 5.4997767e-06 0.00093948055 + 59 3.6324968e-06 5.7547729e-06 0.00093714003 + 60 3.8269172e-06 6.0074589e-06 0.00093480857 + 61 4.0183315e-06 6.2647899e-06 0.00093244785 + 62 4.2263009e-06 6.5384938e-06 0.00092998209 + 63 4.4382735e-06 6.8124372e-06 0.00092751453 + 64 4.6488232e-06 7.0913777e-06 0.00092505171 + 65 4.8707403e-06 7.3773956e-06 0.00092256944 + 66 5.1178339e-06 7.6811232e-06 0.00091984092 + 67 5.377713e-06 7.9951107e-06 0.0009170037 + 68 5.63843e-06 8.3096758e-06 0.00091414977 + 69 5.8975743e-06 8.6232285e-06 0.00091137219 + 70 6.168921e-06 8.9388929e-06 0.00090857909 + 71 6.4488238e-06 9.2524846e-06 0.00090583473 + 72 6.7115203e-06 9.5688611e-06 0.00090308342 + 73 6.995134e-06 9.8873586e-06 0.00090026584 + 74 7.284199e-06 1.0203245e-05 0.00089755256 + 75 7.6029528e-06 1.0523755e-05 0.00089485804 + 76 7.9431391e-06 1.0848398e-05 0.00089217189 + 77 8.2768015e-06 1.1164286e-05 0.00088956967 + 78 8.6058917e-06 1.1466652e-05 0.00088700311 + 79 8.9683477e-06 1.177381e-05 0.00088439373 + 80 9.3324487e-06 1.2059606e-05 0.00088199635 + 81 9.6876036e-06 1.233837e-05 0.0008797438 + 82 1.0061754e-05 1.261531e-05 0.00087745162 + 83 1.0457019e-05 1.2899972e-05 0.00087506863 + 84 1.0866064e-05 1.3160246e-05 0.00087288444 + 85 1.1298832e-05 1.3438712e-05 0.00087052617 + 86 1.1762526e-05 1.372546e-05 0.00086818634 + 87 1.2264941e-05 1.4017631e-05 0.00086572407 + 88 1.2784125e-05 1.4328439e-05 0.000862992 + 89 1.332533e-05 1.4640964e-05 0.00086017364 + 90 1.3897966e-05 1.4978644e-05 0.0008570595 + 91 1.4473024e-05 1.5307214e-05 0.00085411601 + 92 1.5067033e-05 1.5649328e-05 0.00085112325 + 93 1.5677943e-05 1.6015833e-05 0.0008479391 + 94 1.6254294e-05 1.6393131e-05 0.00084474237 + 95 1.6791089e-05 1.674573e-05 0.00084165639 + 96 1.7362795e-05 1.7094435e-05 0.00083858996 + 97 1.7961498e-05 1.7451741e-05 0.00083536142 + 98 1.8568337e-05 1.780703e-05 0.00083221292 + 99 1.9188379e-05 1.8141862e-05 0.00082913227 + 100 1.9789011e-05 1.8457846e-05 0.00082619877 +Loop time of 0.112003 on 1 procs for 100 steps with 3 atoms + +Performance: 38570373.396 fs/day, 0.000 hours/fs, 892.833 timesteps/s, 2.678 katom-step/s +52.7% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 1.4356e-05 | 1.4356e-05 | 1.4356e-05 | 0.0 | 0.01 +Comm | 9.7936e-05 | 9.7936e-05 | 9.7936e-05 | 0.0 | 0.09 +Output | 0.0017373 | 0.0017373 | 0.0017373 | 0.0 | 1.55 +Modify | 0.10997 | 0.10997 | 0.10997 | 0.0 | 98.19 +Other | | 0.0001804 | | | 0.16 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 41 +Dangerous builds = 0 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.3 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.3 new file mode 100644 index 0000000000..941334c25c --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/log.lammps.3 @@ -0,0 +1,228 @@ +LAMMPS (4 Feb 2025) +Processor partition = 3 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188894 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/langevin/bosonic ensemble nvt temp ${Temp} thermostat PILE_L 12345 tau 50 fixcom no +fix pimdb all pimd/langevin/bosonic ensemble nvt temp 17.4 thermostat PILE_L 12345 tau 50 fixcom no + +# Outputs +variable prim_kinetic equal f_pimdb[5] +variable virial equal f_pimdb[6] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 2.801 | 2.801 | 2.801 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 1.3661449e-08 0.0009918329 + 1 8.448696e-10 1.8742641e-09 0.00099182267 + 2 3.5286528e-09 7.1801177e-09 0.00099180053 + 3 8.2936659e-09 1.4372388e-08 0.00099176038 + 4 1.4586045e-08 3.1274072e-08 0.00099170266 + 5 2.4161875e-08 5.2470231e-08 0.00099162511 + 6 3.8801085e-08 7.2301291e-08 0.00099153 + 7 5.9758384e-08 9.5440683e-08 0.00099141194 + 8 8.5838823e-08 1.1669162e-07 0.00099127168 + 9 1.1684437e-07 1.3551097e-07 0.00099109846 + 10 1.4954974e-07 1.6050328e-07 0.00099089759 + 11 1.8686687e-07 1.8645705e-07 0.0009906735 + 12 2.2747543e-07 2.1181418e-07 0.00099042271 + 13 2.807392e-07 2.4218941e-07 0.00099012073 + 14 3.3827867e-07 2.7507982e-07 0.00098977876 + 15 3.9918359e-07 3.1771012e-07 0.00098940169 + 16 4.6034561e-07 3.6075963e-07 0.00098900826 + 17 5.2315136e-07 4.0059005e-07 0.00098859985 + 18 5.8457073e-07 4.3883585e-07 0.00098817535 + 19 6.5273892e-07 4.8036287e-07 0.00098770764 + 20 7.2150647e-07 5.2362451e-07 0.00098721266 + 21 8.0596691e-07 5.7705077e-07 0.00098670407 + 22 8.946994e-07 6.3271105e-07 0.00098617081 + 23 9.9120142e-07 6.9382916e-07 0.00098557406 + 24 1.09722e-06 7.5979853e-07 0.0009849612 + 25 1.2053692e-06 8.232213e-07 0.00098434971 + 26 1.3196547e-06 8.9901705e-07 0.0009836645 + 27 1.4380033e-06 9.7033126e-07 0.0009830068 + 28 1.5700587e-06 1.0444697e-06 0.00098229311 + 29 1.7089133e-06 1.1208219e-06 0.00098150213 + 30 1.839944e-06 1.1899315e-06 0.00098077686 + 31 1.9718088e-06 1.2560026e-06 0.00098003518 + 32 2.1244847e-06 1.3287355e-06 0.00097921099 + 33 2.2823018e-06 1.399856e-06 0.0009784037 + 34 2.4607034e-06 1.4769272e-06 0.00097754491 + 35 2.6410061e-06 1.5642459e-06 0.00097667479 + 36 2.8395019e-06 1.6571388e-06 0.00097572641 + 37 3.0476467e-06 1.7603492e-06 0.00097475881 + 38 3.2661567e-06 1.8666516e-06 0.00097373763 + 39 3.4857766e-06 1.9696572e-06 0.0009727081 + 40 3.7310241e-06 2.0910539e-06 0.00097155807 + 41 3.9947356e-06 2.2198868e-06 0.00097030871 + 42 4.2722057e-06 2.3576614e-06 0.00096896762 + 43 4.5444932e-06 2.496284e-06 0.00096760456 + 44 4.8307956e-06 2.637657e-06 0.00096619083 + 45 5.1356773e-06 2.795845e-06 0.00096470693 + 46 5.4835294e-06 2.9727162e-06 0.00096303835 + 47 5.8235502e-06 3.1451673e-06 0.00096140406 + 48 6.1541132e-06 3.3241907e-06 0.00095971129 + 49 6.4796693e-06 3.5094834e-06 0.00095790964 + 50 6.8009213e-06 3.7049195e-06 0.00095603039 + 51 7.1253629e-06 3.9139587e-06 0.0009540918 + 52 7.4712081e-06 4.1289537e-06 0.0009520767 + 53 7.8322309e-06 4.3481079e-06 0.00095003676 + 54 8.1941694e-06 4.5702968e-06 0.00094802906 + 55 8.579188e-06 4.7928135e-06 0.00094600351 + 56 8.9656493e-06 5.0253448e-06 0.0009438419 + 57 9.351102e-06 5.2569331e-06 0.00094169256 + 58 9.7645407e-06 5.4997767e-06 0.00093948055 + 59 1.0230434e-05 5.7547729e-06 0.00093714003 + 60 1.0724041e-05 6.0074589e-06 0.00093480857 + 61 1.1207453e-05 6.2647899e-06 0.00093244785 + 62 1.171285e-05 6.5384938e-06 0.00092998209 + 63 1.221299e-05 6.8124372e-06 0.00092751453 + 64 1.2743686e-05 7.0913777e-06 0.00092505171 + 65 1.3295725e-05 7.3773956e-06 0.00092256944 + 66 1.3906023e-05 7.6811232e-06 0.00091984092 + 67 1.4511348e-05 7.9951107e-06 0.0009170037 + 68 1.5116776e-05 8.3096758e-06 0.00091414977 + 69 1.571005e-05 8.6232285e-06 0.00091137219 + 70 1.6297136e-05 8.9388929e-06 0.00090857909 + 71 1.6849099e-05 9.2524846e-06 0.00090583473 + 72 1.7417168e-05 9.5688611e-06 0.00090308342 + 73 1.7971613e-05 9.8873586e-06 0.00090026584 + 74 1.8509472e-05 1.0203245e-05 0.00089755256 + 75 1.9047825e-05 1.0523755e-05 0.00089485804 + 76 1.9510646e-05 1.0848398e-05 0.00089217189 + 77 1.9903304e-05 1.1164286e-05 0.00088956967 + 78 2.0300948e-05 1.1466652e-05 0.00088700311 + 79 2.0696905e-05 1.177381e-05 0.00088439373 + 80 2.1058218e-05 1.2059606e-05 0.00088199635 + 81 2.139726e-05 1.233837e-05 0.0008797438 + 82 2.1755155e-05 1.261531e-05 0.00087745162 + 83 2.2051018e-05 1.2899972e-05 0.00087506863 + 84 2.2294567e-05 1.3160246e-05 0.00087288444 + 85 2.2559076e-05 1.3438712e-05 0.00087052617 + 86 2.2815918e-05 1.372546e-05 0.00086818634 + 87 2.3062616e-05 1.4017631e-05 0.00086572407 + 88 2.3324206e-05 1.4328439e-05 0.000862992 + 89 2.3559323e-05 1.4640964e-05 0.00086017364 + 90 2.3792776e-05 1.4978644e-05 0.0008570595 + 91 2.3990667e-05 1.5307214e-05 0.00085411601 + 92 2.4198637e-05 1.5649328e-05 0.00085112325 + 93 2.4419398e-05 1.6015833e-05 0.0008479391 + 94 2.4625252e-05 1.6393131e-05 0.00084474237 + 95 2.4816363e-05 1.674573e-05 0.00084165639 + 96 2.4982897e-05 1.7094435e-05 0.00083858996 + 97 2.5153682e-05 1.7451741e-05 0.00083536142 + 98 2.5288512e-05 1.780703e-05 0.00083221292 + 99 2.5384836e-05 1.8141862e-05 0.00082913227 + 100 2.5401412e-05 1.8457846e-05 0.00082619877 +Loop time of 0.122921 on 1 procs for 100 steps with 3 atoms + +Performance: 35144393.915 fs/day, 0.000 hours/fs, 813.528 timesteps/s, 2.441 katom-step/s +88.8% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 1.5885e-05 | 1.5885e-05 | 1.5885e-05 | 0.0 | 0.01 +Comm | 9.4707e-05 | 9.4707e-05 | 9.4707e-05 | 0.0 | 0.08 +Output | 0.0027076 | 0.0027076 | 0.0027076 | 0.0 | 2.20 +Modify | 0.11993 | 0.11993 | 0.11993 | 0.0 | 97.57 +Other | | 0.0001738 | | | 0.14 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 42 +Dangerous builds = 0 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/run.sh b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/run.sh new file mode 100644 index 0000000000..c1329f646a --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_langevin/run.sh @@ -0,0 +1 @@ +mpirun -np $NBEADS $LMP -in in.lmp -partition ${NBEADS}x1 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/in.lmp b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/in.lmp new file mode 100644 index 0000000000..56d6d3932a --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/in.lmp @@ -0,0 +1,60 @@ +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Number of Nose-Hoover chains +variable nhc equal 4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/nvt/bosonic method pimd temp ${Temp} nhc ${nhc} + +# Outputs +variable virial equal f_pimdb[3] +variable prim_kinetic equal f_pimdb[4] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.0 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.0 new file mode 100644 index 0000000000..bf46982da6 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.0 @@ -0,0 +1,236 @@ +LAMMPS (4 Feb 2025) +Processor partition = 0 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Number of Nose-Hoover chains +variable nhc equal 4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188891 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/nvt/bosonic method pimd temp ${Temp} nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc 4 + +# Outputs +variable virial equal f_pimdb[3] +variable prim_kinetic equal f_pimdb[4] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No fixes with time integration, atoms won't move (../verlet.cpp:60) +Fix pimd/nvt -P/(beta^2 * hbar^2) = -2.2139311e-05 (kcal/mol/A^2) + +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 9.176 | 9.176 | 9.176 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 0 0.00024794798 + 1 8.5849317e-10 4.5664581e-10 0.00024794215 + 2 3.4338094e-09 1.3424585e-09 0.00024793325 + 3 7.7255501e-09 2.6573351e-09 0.00024792129 + 4 1.3733112e-08 4.4011211e-09 0.00024790626 + 5 2.1455686e-08 6.5736113e-09 0.00024788818 + 6 3.0892259e-08 9.1745487e-09 0.00024786704 + 7 4.2041614e-08 1.2203626e-08 0.00024784285 + 8 5.4902328e-08 1.5660483e-08 0.0002478156 + 9 6.9472776e-08 1.9544712e-08 0.00024778532 + 10 8.5751128e-08 2.385585e-08 0.00024775199 + 11 1.0373535e-07 2.8593386e-08 0.00024771564 + 12 1.234232e-07 3.3756758e-08 0.00024767626 + 13 1.4481225e-07 3.9345352e-08 0.00024763386 + 14 1.6789985e-07 4.5358504e-08 0.00024758844 + 15 1.9268316e-07 5.1795501e-08 0.00024754003 + 16 2.1915914e-07 5.8655576e-08 0.00024748862 + 17 2.4732452e-07 6.5937915e-08 0.00024743423 + 18 2.7717588e-07 7.3641654e-08 0.00024737686 + 19 3.0870956e-07 8.1765876e-08 0.00024731653 + 20 3.4192172e-07 9.0309618e-08 0.00024725324 + 21 3.7680831e-07 9.9271863e-08 0.00024718701 + 22 4.1336509e-07 1.0865155e-07 0.00024711785 + 23 4.5158762e-07 1.1844756e-07 0.00024704578 + 24 4.9147126e-07 1.2865874e-07 0.0002469708 + 25 5.3301119e-07 1.3928387e-07 0.00024689293 + 26 5.7620236e-07 1.5032169e-07 0.00024681218 + 27 6.2103957e-07 1.6177089e-07 0.00024672857 + 28 6.675174e-07 1.7363012e-07 0.00024664212 + 29 7.1563024e-07 1.8589798e-07 0.00024655283 + 30 7.6537229e-07 1.9857299e-07 0.00024646073 + 31 8.1673756e-07 2.1165368e-07 0.00024636584 + 32 8.6971988e-07 2.2513849e-07 0.00024626816 + 33 9.2431288e-07 2.3902582e-07 0.00024616772 + 34 9.8051001e-07 2.5331402e-07 0.00024606453 + 35 1.0383045e-06 2.6800143e-07 0.00024595862 + 36 1.0976895e-06 2.8308629e-07 0.00024585001 + 37 1.1586578e-06 2.9856682e-07 0.0002457387 + 38 1.2212021e-06 3.1444121e-07 0.00024562473 + 39 1.2853151e-06 3.3070757e-07 0.00024550812 + 40 1.350989e-06 3.4736399e-07 0.00024538888 + 41 1.4182159e-06 3.644085e-07 0.00024526703 + 42 1.4869879e-06 3.8183911e-07 0.00024514261 + 43 1.5572969e-06 3.9965375e-07 0.00024501562 + 44 1.6291343e-06 4.1785034e-07 0.0002448861 + 45 1.7024917e-06 4.3642673e-07 0.00024475406 + 46 1.7773604e-06 4.5538075e-07 0.00024461954 + 47 1.8537315e-06 4.7471016e-07 0.00024448255 + 48 1.9315959e-06 4.9441271e-07 0.00024434312 + 49 2.0109445e-06 5.1448608e-07 0.00024420127 + 50 2.0917677e-06 5.3492792e-07 0.00024405703 + 51 2.1740562e-06 5.5573584e-07 0.00024391043 + 52 2.2578001e-06 5.769074e-07 0.00024376149 + 53 2.3429895e-06 5.9844014e-07 0.00024361023 + 54 2.4296146e-06 6.2033154e-07 0.00024345669 + 55 2.5176649e-06 6.4257904e-07 0.00024330088 + 56 2.6071302e-06 6.6518005e-07 0.00024314285 + 57 2.6980001e-06 6.8813194e-07 0.00024298261 + 58 2.7902637e-06 7.1143203e-07 0.0002428202 + 59 2.8839103e-06 7.3507763e-07 0.00024265563 + 60 2.9789289e-06 7.5906598e-07 0.00024248895 + 61 3.0753084e-06 7.833943e-07 0.00024232018 + 62 3.1730376e-06 8.0805977e-07 0.00024214934 + 63 3.2721049e-06 8.3305953e-07 0.00024197647 + 64 3.372499e-06 8.5839069e-07 0.00024180161 + 65 3.474208e-06 8.8405032e-07 0.00024162477 + 66 3.5772201e-06 9.1003545e-07 0.00024144598 + 67 3.6815235e-06 9.3634309e-07 0.00024126529 + 68 3.7871058e-06 9.6297021e-07 0.00024108272 + 69 3.893955e-06 9.8991374e-07 0.0002408983 + 70 4.0020586e-06 1.0171706e-06 0.00024071206 + 71 4.1114041e-06 1.0447376e-06 0.00024052403 + 72 4.2219789e-06 1.0726116e-06 0.00024033425 + 73 4.3337702e-06 1.1007895e-06 0.00024014274 + 74 4.4467651e-06 1.1292679e-06 0.00023994955 + 75 4.5609507e-06 1.1580437e-06 0.0002397547 + 76 4.6763137e-06 1.1871136e-06 0.00023955822 + 77 4.792841e-06 1.2164741e-06 0.00023936014 + 78 4.9105191e-06 1.2461221e-06 0.00023916051 + 79 5.0293346e-06 1.276054e-06 0.00023895935 + 80 5.149274e-06 1.3062666e-06 0.00023875669 + 81 5.2703234e-06 1.3367563e-06 0.00023855257 + 82 5.392469e-06 1.3675198e-06 0.00023834703 + 83 5.515697e-06 1.3985535e-06 0.00023814009 + 84 5.6399934e-06 1.4298539e-06 0.00023793178 + 85 5.7653439e-06 1.4614174e-06 0.00023772215 + 86 5.8917344e-06 1.4932406e-06 0.00023751123 + 87 6.0191505e-06 1.5253199e-06 0.00023729904 + 88 6.1475779e-06 1.5576515e-06 0.00023708563 + 89 6.277002e-06 1.590232e-06 0.00023687102 + 90 6.4074081e-06 1.6230576e-06 0.00023665526 + 91 6.5387817e-06 1.6561246e-06 0.00023643837 + 92 6.6711079e-06 1.6894294e-06 0.00023622039 + 93 6.8043718e-06 1.7229682e-06 0.00023600135 + 94 6.9385585e-06 1.7567372e-06 0.00023578129 + 95 7.0736529e-06 1.7907328e-06 0.00023556024 + 96 7.2096401e-06 1.8249512e-06 0.00023533824 + 97 7.3465046e-06 1.8593885e-06 0.00023511531 + 98 7.4842314e-06 1.8940408e-06 0.0002348915 + 99 7.622805e-06 1.9289045e-06 0.00023466684 + 100 7.76221e-06 1.9639756e-06 0.00023444136 +Loop time of 0.00940749 on 1 procs for 100 steps with 3 atoms + +Performance: 459208566.791 fs/day, 0.000 hours/fs, 10629.828 timesteps/s, 31.889 katom-step/s +90.3% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 8.466e-06 | 8.466e-06 | 8.466e-06 | 0.0 | 0.09 +Comm | 7.8365e-05 | 7.8365e-05 | 7.8365e-05 | 0.0 | 0.83 +Output | 0.0012482 | 0.0012482 | 0.0012482 | 0.0 | 13.27 +Modify | 0.0079193 | 0.0079193 | 0.0079193 | 0.0 | 84.18 +Other | | 0.0001532 | | | 1.63 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 24 +Dangerous builds = 0 +Total wall time: 0:00:00 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.1 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.1 new file mode 100644 index 0000000000..f129f10c09 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.1 @@ -0,0 +1,234 @@ +LAMMPS (4 Feb 2025) +Processor partition = 1 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Number of Nose-Hoover chains +variable nhc equal 4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188892 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/nvt/bosonic method pimd temp ${Temp} nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc 4 + +# Outputs +variable virial equal f_pimdb[3] +variable prim_kinetic equal f_pimdb[4] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No fixes with time integration, atoms won't move (../verlet.cpp:60) +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 9.176 | 9.176 | 9.176 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 0 0.00024796164 + 1 8.5850843e-10 5.313724e-10 0.00024796083 + 2 3.4339027e-09 1.4919466e-09 0.00024795841 + 3 7.7257901e-09 2.8816048e-09 0.00024795438 + 4 1.3733516e-08 4.7001635e-09 0.00024794873 + 5 2.1456164e-08 6.9473741e-09 0.00024794148 + 6 3.0892555e-08 9.6229225e-09 0.00024793261 + 7 4.2041252e-08 1.2726429e-08 0.00024792214 + 8 5.4900553e-08 1.6257449e-08 0.00024791006 + 9 6.9468499e-08 2.0215473e-08 0.00024789638 + 10 8.5742867e-08 2.4599926e-08 0.0002478811 + 11 1.0372118e-07 2.9410168e-08 0.00024786422 + 12 1.2340069e-07 3.4645495e-08 0.00024784576 + 13 1.447784e-07 4.0305136e-08 0.0002478257 + 14 1.6785106e-07 4.6388259e-08 0.00024780406 + 15 1.9261515e-07 5.2893965e-08 0.00024778084 + 16 2.190669e-07 5.9821291e-08 0.00024775605 + 17 2.4720228e-07 6.7169211e-08 0.00024772969 + 18 2.7701701e-07 7.4936635e-08 0.00024770176 + 19 3.0850655e-07 8.312241e-08 0.00024767228 + 20 3.416661e-07 9.1725317e-08 0.00024764125 + 21 3.7649064e-07 1.0074408e-07 0.00024760868 + 22 4.1297485e-07 1.1017735e-07 0.00024757457 + 23 4.511132e-07 1.2002372e-07 0.00024753893 + 24 4.9089989e-07 1.3028173e-07 0.00024750177 + 25 5.3232887e-07 1.4094984e-07 0.0002474631 + 26 5.7539386e-07 1.5202647e-07 0.00024742293 + 27 6.2008832e-07 1.6350996e-07 0.00024738126 + 28 6.6640547e-07 1.7539859e-07 0.00024733811 + 29 7.1433828e-07 1.876906e-07 0.00024729348 + 30 7.6387949e-07 2.0038414e-07 0.00024724738 + 31 8.150216e-07 2.1347732e-07 0.00024719983 + 32 8.6775685e-07 2.2696818e-07 0.00024715084 + 33 9.2207727e-07 2.4085472e-07 0.0002471004 + 34 9.7797464e-07 2.5513485e-07 0.00024704855 + 35 1.0354405e-06 2.6980645e-07 0.00024699528 + 36 1.0944662e-06 2.8486732e-07 0.00024694061 + 37 1.1550427e-06 3.0031522e-07 0.00024688456 + 38 1.2171611e-06 3.1614785e-07 0.00024682713 + 39 1.2808118e-06 3.3236283e-07 0.00024676833 + 40 1.3459853e-06 3.4895776e-07 0.00024670818 + 41 1.4126717e-06 3.6593016e-07 0.0002466467 + 42 1.4808611e-06 3.832775e-07 0.0002465839 + 43 1.5505431e-06 4.009972e-07 0.00024651978 + 44 1.6217074e-06 4.1908663e-07 0.00024645437 + 45 1.6943431e-06 4.3754308e-07 0.00024638767 + 46 1.7684394e-06 4.5636381e-07 0.00024631971 + 47 1.8439852e-06 4.7554603e-07 0.0002462505 + 48 1.9209691e-06 4.9508689e-07 0.00024618005 + 49 1.9993796e-06 5.1498348e-07 0.00024610839 + 50 2.079205e-06 5.3523286e-07 0.00024603551 + 51 2.1604333e-06 5.5583202e-07 0.00024596145 + 52 2.2430524e-06 5.7677791e-07 0.00024588621 + 53 2.32705e-06 5.9806742e-07 0.00024580981 + 54 2.4124135e-06 6.1969741e-07 0.00024573228 + 55 2.4991302e-06 6.4166468e-07 0.00024565362 + 56 2.5871872e-06 6.6396597e-07 0.00024557386 + 57 2.6765714e-06 6.8659801e-07 0.00024549301 + 58 2.7672695e-06 7.0955745e-07 0.00024541108 + 59 2.859268e-06 7.328409e-07 0.00024532811 + 60 2.9525533e-06 7.5644494e-07 0.00024524409 + 61 3.0471115e-06 7.8036609e-07 0.00024515907 + 62 3.1429287e-06 8.0460083e-07 0.00024507304 + 63 3.2399907e-06 8.2914562e-07 0.00024498603 + 64 3.3382831e-06 8.5399683e-07 0.00024489807 + 65 3.4377915e-06 8.7915084e-07 0.00024480916 + 66 3.5385011e-06 9.0460395e-07 0.00024471933 + 67 3.6403972e-06 9.3035245e-07 0.0002446286 + 68 3.7434648e-06 9.5639257e-07 0.00024453698 + 69 3.8476886e-06 9.8272051e-07 0.0002444445 + 70 3.9530535e-06 1.0093324e-06 0.00024435118 + 71 4.0595439e-06 1.0362244e-06 0.00024425703 + 72 4.1671443e-06 1.0633926e-06 0.00024416208 + 73 4.275839e-06 1.0908331e-06 0.00024406635 + 74 4.3856121e-06 1.1185418e-06 0.00024396985 + 75 4.4964475e-06 1.1465147e-06 0.00024387262 + 76 4.6083292e-06 1.1747478e-06 0.00024377466 + 77 4.7212408e-06 1.2032371e-06 0.00024367601 + 78 4.835166e-06 1.2319783e-06 0.00024357668 + 79 4.9500883e-06 1.2609674e-06 0.00024347669 + 80 5.0659909e-06 1.2902001e-06 0.00024337607 + 81 5.1828572e-06 1.3196723e-06 0.00024327483 + 82 5.3006702e-06 1.3493797e-06 0.000243173 + 83 5.419413e-06 1.379318e-06 0.0002430706 + 84 5.5390685e-06 1.4094831e-06 0.00024296765 + 85 5.6596194e-06 1.4398705e-06 0.00024286417 + 86 5.7810485e-06 1.4704759e-06 0.00024276019 + 87 5.9033384e-06 1.501295e-06 0.00024265573 + 88 6.0264715e-06 1.5323234e-06 0.00024255081 + 89 6.1504303e-06 1.5635565e-06 0.00024244545 + 90 6.275197e-06 1.5949902e-06 0.00024233967 + 91 6.400754e-06 1.6266197e-06 0.0002422335 + 92 6.5270834e-06 1.6584408e-06 0.00024212696 + 93 6.6541672e-06 1.6904488e-06 0.00024202007 + 94 6.7819875e-06 1.7226394e-06 0.00024191285 + 95 6.9105262e-06 1.7550078e-06 0.00024180534 + 96 7.0397651e-06 1.7875497e-06 0.00024169754 + 97 7.169686e-06 1.8202604e-06 0.00024158948 + 98 7.3002707e-06 1.8531354e-06 0.00024148118 + 99 7.4315008e-06 1.88617e-06 0.00024137268 + 100 7.563358e-06 1.9193596e-06 0.00024126398 +Loop time of 0.00941353 on 1 procs for 100 steps with 3 atoms + +Performance: 458913876.206 fs/day, 0.000 hours/fs, 10623.006 timesteps/s, 31.869 katom-step/s +50.9% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 8.215e-06 | 8.215e-06 | 8.215e-06 | 0.0 | 0.09 +Comm | 7.7692e-05 | 7.7692e-05 | 7.7692e-05 | 0.0 | 0.83 +Output | 0.0047662 | 0.0047662 | 0.0047662 | 0.0 | 50.63 +Modify | 0.004407 | 0.004407 | 0.004407 | 0.0 | 46.82 +Other | | 0.0001545 | | | 1.64 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 23 +Dangerous builds = 0 +Total wall time: 0:00:00 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.2 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.2 new file mode 100644 index 0000000000..074acc9d50 --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.2 @@ -0,0 +1,234 @@ +LAMMPS (4 Feb 2025) +Processor partition = 2 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Number of Nose-Hoover chains +variable nhc equal 4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188893 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/nvt/bosonic method pimd temp ${Temp} nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc 4 + +# Outputs +variable virial equal f_pimdb[3] +variable prim_kinetic equal f_pimdb[4] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No fixes with time integration, atoms won't move (../verlet.cpp:60) +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 9.176 | 9.176 | 9.176 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 0 0.00024796164 + 1 8.5850845e-10 -3.6831179e-10 0.0002479594 + 2 3.4338783e-09 -3.0736303e-10 0.00024795267 + 3 7.7256433e-09 1.8277487e-10 0.00024794145 + 4 1.3733026e-08 1.1019528e-09 0.00024792576 + 5 2.145494e-08 2.449944e-09 0.00024790558 + 6 3.0889985e-08 4.226444e-09 0.00024788093 + 7 4.2036455e-08 6.4310711e-09 0.00024785181 + 8 5.4892332e-08 9.0633658e-09 0.00024781823 + 9 6.9455288e-08 1.2122792e-08 0.00024778019 + 10 8.572269e-08 1.5608734e-08 0.0002477377 + 11 1.0369159e-07 1.9520504e-08 0.00024769077 + 12 1.2335875e-07 2.3857331e-08 0.00024763941 + 13 1.4472059e-07 2.8618373e-08 0.00024758363 + 14 1.6777327e-07 3.3802708e-08 0.00024752344 + 15 1.9251262e-07 3.940934e-08 0.00024745886 + 16 2.1893415e-07 4.5437194e-08 0.00024738989 + 17 2.4703311e-07 5.1885122e-08 0.00024731656 + 18 2.768044e-07 5.8751899e-08 0.00024723887 + 19 3.0824267e-07 6.6036227e-08 0.00024715685 + 20 3.4134222e-07 7.3736729e-08 0.00024707051 + 21 3.7609709e-07 8.1851958e-08 0.00024697986 + 22 4.1250102e-07 9.0380389e-08 0.00024688494 + 23 4.5054742e-07 9.9320426e-08 0.00024678575 + 24 4.9022945e-07 1.086704e-07 0.00024668231 + 25 5.3153996e-07 1.1842856e-07 0.00024657466 + 26 5.744715e-07 1.285931e-07 0.00024646281 + 27 6.1901635e-07 1.3916212e-07 0.00024634678 + 28 6.6516648e-07 1.5013367e-07 0.0002462266 + 29 7.1291361e-07 1.6150571e-07 0.00024610229 + 30 7.6224914e-07 1.7327614e-07 0.00024597388 + 31 8.1316422e-07 1.8544279e-07 0.0002458414 + 32 8.6564969e-07 1.9800342e-07 0.00024570487 + 33 9.1969614e-07 2.109557e-07 0.00024556432 + 34 9.7529388e-07 2.2429727e-07 0.00024541978 + 35 1.0324329e-06 2.3802567e-07 0.00024527128 + 36 1.0911031e-06 2.5213838e-07 0.00024511886 + 37 1.1512938e-06 2.6663282e-07 0.00024496253 + 38 1.2129943e-06 2.8150634e-07 0.00024480234 + 39 1.2761936e-06 2.9675622e-07 0.00024463831 + 40 1.3408803e-06 3.1237968e-07 0.00024447049 + 41 1.407043e-06 3.2837388e-07 0.0002442989 + 42 1.4746698e-06 3.447359e-07 0.00024412359 + 43 1.5437487e-06 3.6146277e-07 0.00024394458 + 44 1.6142674e-06 3.7855146e-07 0.00024376192 + 45 1.6862132e-06 3.9599886e-07 0.00024357564 + 46 1.7595735e-06 4.1380182e-07 0.00024338577 + 47 1.8343351e-06 4.3195711e-07 0.00024319237 + 48 1.9104848e-06 4.5046147e-07 0.00024299547 + 49 1.9880091e-06 4.6931154e-07 0.00024279511 + 50 2.0668943e-06 4.8850394e-07 0.00024259132 + 51 2.1471263e-06 5.0803521e-07 0.00024238416 + 52 2.228691e-06 5.2790184e-07 0.00024217367 + 53 2.311574e-06 5.4810027e-07 0.00024195988 + 54 2.3957606e-06 5.6862687e-07 0.00024174285 + 55 2.481236e-06 5.8947797e-07 0.00024152261 + 56 2.5679851e-06 6.1064984e-07 0.00024129922 + 57 2.6559926e-06 6.3213871e-07 0.00024107271 + 58 2.7452431e-06 6.5394073e-07 0.00024084313 + 59 2.8357209e-06 6.7605203e-07 0.00024061054 + 60 2.9274101e-06 6.9846868e-07 0.00024037498 + 61 3.0202947e-06 7.2118669e-07 0.00024013649 + 62 3.1143583e-06 7.4420204e-07 0.00023989513 + 63 3.2095846e-06 7.6751065e-07 0.00023965094 + 64 3.3059569e-06 7.911084e-07 0.00023940398 + 65 3.4034585e-06 8.1499113e-07 0.00023915429 + 66 3.5020722e-06 8.3915461e-07 0.00023890193 + 67 3.601781e-06 8.6359461e-07 0.00023864694 + 68 3.7025676e-06 8.8830681e-07 0.00023838939 + 69 3.8044144e-06 9.132869e-07 0.00023812931 + 70 3.9073038e-06 9.3853048e-07 0.00023786677 + 71 4.0112181e-06 9.6403314e-07 0.00023760182 + 72 4.1161392e-06 9.8979043e-07 0.00023733451 + 73 4.2220491e-06 1.0157978e-06 0.0002370649 + 74 4.3289295e-06 1.0420509e-06 0.00023679303 + 75 4.4367621e-06 1.068545e-06 0.00023651897 + 76 4.5455283e-06 1.0952755e-06 0.00023624277 + 77 4.6552095e-06 1.1222378e-06 0.00023596449 + 78 4.7657869e-06 1.1494273e-06 0.00023568418 + 79 4.8772416e-06 1.1768393e-06 0.0002354019 + 80 4.9895545e-06 1.2044689e-06 0.00023511771 + 81 5.1027067e-06 1.2323116e-06 0.00023483166 + 82 5.2166788e-06 1.2603625e-06 0.00023454381 + 83 5.3314515e-06 1.2886168e-06 0.00023425422 + 84 5.4470054e-06 1.3170697e-06 0.00023396295 + 85 5.5633208e-06 1.3457162e-06 0.00023367006 + 86 5.6803783e-06 1.3745516e-06 0.0002333756 + 87 5.798158e-06 1.403571e-06 0.00023307963 + 88 5.9166402e-06 1.4327693e-06 0.00023278222 + 89 6.0358049e-06 1.4621418e-06 0.00023248341 + 90 6.1556323e-06 1.4916833e-06 0.00023218329 + 91 6.2761023e-06 1.5213889e-06 0.00023188189 + 92 6.3971948e-06 1.5512537e-06 0.00023157929 + 93 6.5188897e-06 1.5812727e-06 0.00023127554 + 94 6.6411667e-06 1.6114407e-06 0.0002309707 + 95 6.7640056e-06 1.6417527e-06 0.00023066484 + 96 6.8873862e-06 1.6722038e-06 0.00023035801 + 97 7.011288e-06 1.7027888e-06 0.00023005028 + 98 7.1356907e-06 1.7335027e-06 0.0002297417 + 99 7.2605738e-06 1.7643404e-06 0.00022943234 + 100 7.3859169e-06 1.7952968e-06 0.00022912226 +Loop time of 0.00941372 on 1 procs for 100 steps with 3 atoms + +Performance: 458904516.311 fs/day, 0.000 hours/fs, 10622.790 timesteps/s, 31.868 katom-step/s +24.7% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 8.785e-06 | 8.785e-06 | 8.785e-06 | 0.0 | 0.09 +Comm | 7.9921e-05 | 7.9921e-05 | 7.9921e-05 | 0.0 | 0.85 +Output | 0.0071119 | 0.0071119 | 0.0071119 | 0.0 | 75.55 +Modify | 0.0020558 | 0.0020558 | 0.0020558 | 0.0 | 21.84 +Other | | 0.0001572 | | | 1.67 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 24 +Dangerous builds = 0 +Total wall time: 0:00:00 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.3 b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.3 new file mode 100644 index 0000000000..edebd7096c --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/log.lammps.3 @@ -0,0 +1,234 @@ +LAMMPS (4 Feb 2025) +Processor partition = 3 +# Units and dimensions +units electron +dimension 3 +boundary p p p +atom_style atomic +atom_modify map yes +pair_style none + +# Time step (in femtoseconds) +timestep 0.5 + +# Temperature (in Kelvin) +variable Temp equal 17.4 + +# Number of Nose-Hoover chains +variable nhc equal 4 + +# Force constant (in Hartree energies per Bohr radius squared) +variable k equal 1.2154614120000001e-08 + +# Number of beads +variable Nbeads equal 4 +variable ibead uloop ${Nbeads} pad +variable ibead uloop 4 pad +variable seed equal 18889 + +# Create box and atoms. All distances are in Bohr +region box block -1500 1500 -1500 1500 -1500 1500 +create_box 1 box +Created orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + 1 by 1 by 1 MPI processor grid +variable a loop 3 +label loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14.5 ${y1} ${z1} +create_atoms 1 single -14.5 9.5 ${z1} +create_atoms 1 single -14.5 9.5 1.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -14 ${y1} ${z1} +create_atoms 1 single -14 9 ${z1} +create_atoms 1 single -14 9 2 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop +variable x1 equal -15.0+0.5*v_a +variable y1 equal 10.0-0.5*v_a +variable z1 equal 1+0.5*v_a +create_atoms 1 single ${x1} ${y1} ${z1} +create_atoms 1 single -13.5 ${y1} ${z1} +create_atoms 1 single -13.5 8.5 ${z1} +create_atoms 1 single -13.5 8.5 2.5 +Created 1 atoms + using lattice units in orthogonal box = (-1500 -1500 -1500) to (1500 1500 1500) + create_atoms CPU = 0.000 seconds +next a +jump SELF loop + +# Electron mass (in amu) +mass 1 0.00054858 + +# Initialize velocities +velocity all create ${Temp} ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 ${seed}${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 18889${ibead} mom yes rot yes dist gaussian +velocity all create 17.4 188894 mom yes rot yes dist gaussian + +# Add harmonic external force +fix harm all spring/self ${k} +fix harm all spring/self 1.215461412e-08 + +# Add harmonic potential energy to total energy and potential energy +fix_modify harm energy yes + +# PIMD command +fix pimdb all pimd/nvt/bosonic method pimd temp ${Temp} nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc ${nhc} +fix pimdb all pimd/nvt/bosonic method pimd temp 17.4 nhc 4 + +# Outputs +variable virial equal f_pimdb[3] +variable prim_kinetic equal f_pimdb[4] + +thermo_style custom step pe v_virial v_prim_kinetic +thermo 1 + +run 100 +WARNING: No fixes with time integration, atoms won't move (../verlet.cpp:60) +WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (../atom.cpp:2444) +WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (../comm_brick.cpp:212) +Per MPI rank memory allocation (min/avg/max) = 9.176 | 9.176 | 9.176 Mbytes + Step PotEng v_virial v_prim_kinetic + 0 0 0 0.00024796164 + 1 8.5847879e-10 4.4463719e-10 0.00024796024 + 2 3.4336684e-09 1.3183829e-09 0.00024795602 + 3 7.7250055e-09 2.6210625e-09 0.000247949 + 4 1.3731669e-08 4.3524371e-09 0.00024793917 + 5 2.1452581e-08 6.5122035e-09 0.00024792653 + 6 3.0886407e-08 9.0999942e-09 0.0002479111 + 7 4.2031554e-08 1.2115377e-08 0.00024789287 + 8 5.4886174e-08 1.5557858e-08 0.00024787185 + 9 6.9448162e-08 1.9426875e-08 0.00024784804 + 10 8.5715159e-08 2.3721806e-08 0.00024782146 + 11 1.0368455e-07 2.8441963e-08 0.0002477921 + 12 1.2335346e-07 3.3586595e-08 0.00024775998 + 13 1.4471878e-07 3.9154888e-08 0.0002477251 + 14 1.6777712e-07 4.5145965e-08 0.00024768748 + 15 1.9252485e-07 5.1558884e-08 0.00024764712 + 16 2.1895809e-07 5.8392642e-08 0.00024760403 + 17 2.4707272e-07 6.5646173e-08 0.00024755823 + 18 2.7686434e-07 7.3318348e-08 0.00024750973 + 19 3.0832834e-07 8.1407976e-08 0.00024745854 + 20 3.4145982e-07 8.9913804e-08 0.00024740467 + 21 3.7625367e-07 9.8834518e-08 0.00024734813 + 22 4.127045e-07 1.0816874e-07 0.00024728895 + 23 4.5080671e-07 1.1791504e-07 0.00024722714 + 24 4.9055443e-07 1.2807191e-07 0.00024716271 + 25 5.3194155e-07 1.3863779e-07 0.00024709568 + 26 5.7496174e-07 1.4961107e-07 0.00024702607 + 27 6.1960841e-07 1.6099007e-07 0.00024695389 + 28 6.6587473e-07 1.7277305e-07 0.00024687917 + 29 7.1375365e-07 1.8495821e-07 0.00024680191 + 30 7.6323786e-07 1.975437e-07 0.00024672216 + 31 8.1431984e-07 2.1052759e-07 0.00024663991 + 32 8.6699183e-07 2.2390793e-07 0.0002465552 + 33 9.2124582e-07 2.3768267e-07 0.00024646805 + 34 9.7707361e-07 2.5184974e-07 0.00024637847 + 35 1.0344667e-06 2.6640698e-07 0.0002462865 + 36 1.0934165e-06 2.813522e-07 0.00024619215 + 37 1.1539141e-06 2.9668313e-07 0.00024609545 + 38 1.2159503e-06 3.1239747e-07 0.00024599643 + 39 1.2795157e-06 3.2849286e-07 0.00024589511 + 40 1.344601e-06 3.4496686e-07 0.00024579151 + 41 1.4111961e-06 3.61817e-07 0.00024568566 + 42 1.4792913e-06 3.7904076e-07 0.0002455776 + 43 1.5488762e-06 3.9663555e-07 0.00024546734 + 44 1.6199404e-06 4.1459874e-07 0.00024535492 + 45 1.6924734e-06 4.3292764e-07 0.00024524035 + 46 1.7664643e-06 4.5161952e-07 0.00024512369 + 47 1.841902e-06 4.7067159e-07 0.00024500494 + 48 1.9187753e-06 4.90081e-07 0.00024488415 + 49 1.9970729e-06 5.0984488e-07 0.00024476135 + 50 2.076783e-06 5.2996028e-07 0.00024463656 + 51 2.1578939e-06 5.5042422e-07 0.00024450981 + 52 2.2403935e-06 5.7123366e-07 0.00024438115 + 53 2.3242696e-06 5.9238553e-07 0.0002442506 + 54 2.4095099e-06 6.138767e-07 0.00024411819 + 55 2.4961017e-06 6.3570399e-07 0.00024398396 + 56 2.5840323e-06 6.5786419e-07 0.00024384795 + 57 2.6732887e-06 6.8035404e-07 0.00024371018 + 58 2.7638579e-06 7.0317021e-07 0.0002435707 + 59 2.8557265e-06 7.2630938e-07 0.00024342953 + 60 2.9488811e-06 7.4976814e-07 0.00024328672 + 61 3.0433081e-06 7.7354305e-07 0.0002431423 + 62 3.1389936e-06 7.9763065e-07 0.00024299631 + 63 3.2359238e-06 8.2202741e-07 0.00024284878 + 64 3.3340844e-06 8.4672977e-07 0.00024269976 + 65 3.4334612e-06 8.7173415e-07 0.00024254927 + 66 3.5340397e-06 8.970369e-07 0.00024239736 + 67 3.6358055e-06 9.2263435e-07 0.00024224406 + 68 3.7387437e-06 9.4852279e-07 0.00024208942 + 69 3.8428394e-06 9.7469849e-07 0.00024193348 + 70 3.9480777e-06 1.0011576e-06 0.00024177626 + 71 4.0544433e-06 1.0278964e-06 0.00024161782 + 72 4.1619211e-06 1.054911e-06 0.0002414582 + 73 4.2704954e-06 1.0821976e-06 0.00024129742 + 74 4.3801508e-06 1.1097521e-06 0.00024113554 + 75 4.4908716e-06 1.1375706e-06 0.0002409726 + 76 4.6026419e-06 1.1656492e-06 0.00024080863 + 77 4.7154459e-06 1.1939839e-06 0.00024064367 + 78 4.8292673e-06 1.2225706e-06 0.00024047778 + 79 4.9440901e-06 1.2514052e-06 0.00024031099 + 80 5.059898e-06 1.2804837e-06 0.00024014333 + 81 5.1766745e-06 1.309802e-06 0.00023997487 + 82 5.2944032e-06 1.3393558e-06 0.00023980562 + 83 5.4130674e-06 1.369141e-06 0.00023963565 + 84 5.5326505e-06 1.3991534e-06 0.00023946499 + 85 5.6531355e-06 1.4293887e-06 0.00023929369 + 86 5.7745057e-06 1.4598429e-06 0.00023912178 + 87 5.8967439e-06 1.4905114e-06 0.00023894931 + 88 6.0198332e-06 1.5213902e-06 0.00023877633 + 89 6.1437564e-06 1.5524748e-06 0.00023860287 + 90 6.2684962e-06 1.5837609e-06 0.00023842898 + 91 6.3940353e-06 1.6152441e-06 0.00023825471 + 92 6.5203562e-06 1.6469202e-06 0.0002380801 + 93 6.6474416e-06 1.6787847e-06 0.00023790518 + 94 6.7752739e-06 1.7108332e-06 0.00023773001 + 95 6.9038354e-06 1.7430613e-06 0.00023755463 + 96 7.0331085e-06 1.7754645e-06 0.00023737908 + 97 7.1630755e-06 1.8080385e-06 0.00023720341 + 98 7.2937187e-06 1.8407787e-06 0.00023702765 + 99 7.4250201e-06 1.8736806e-06 0.00023685186 + 100 7.5569619e-06 1.9067398e-06 0.00023667607 +Loop time of 0.00939597 on 1 procs for 100 steps with 3 atoms + +Performance: 459771778.655 fs/day, 0.000 hours/fs, 10642.865 timesteps/s, 31.929 katom-step/s +25.2% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0 | 0 | 0 | 0.0 | 0.00 +Neigh | 8.404e-06 | 8.404e-06 | 8.404e-06 | 0.0 | 0.09 +Comm | 8.6872e-05 | 8.6872e-05 | 8.6872e-05 | 0.0 | 0.92 +Output | 0.0071309 | 0.0071309 | 0.0071309 | 0.0 | 75.89 +Modify | 0.0020085 | 0.0020085 | 0.0020085 | 0.0 | 21.38 +Other | | 0.0001612 | | | 1.72 + +Nlocal: 3 ave 3 max 3 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 0 +Ave neighs/atom = 0 +Neighbor list builds = 24 +Dangerous builds = 0 +Total wall time: 0:00:00 diff --git a/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/run.sh b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/run.sh new file mode 100644 index 0000000000..c1329f646a --- /dev/null +++ b/examples/PACKAGES/pimd_bosonic/harmonic_trap_nvt/run.sh @@ -0,0 +1 @@ +mpirun -np $NBEADS $LMP -in in.lmp -partition ${NBEADS}x1 diff --git a/src/EXTRA-FIX/fix_ave_correlate_long.cpp b/src/EXTRA-FIX/fix_ave_correlate_long.cpp index aaae2624f7..4b4ba462b1 100644 --- a/src/EXTRA-FIX/fix_ave_correlate_long.cpp +++ b/src/EXTRA-FIX/fix_ave_correlate_long.cpp @@ -42,7 +42,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using MathSpecial::powint; -enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL }; +enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL, FIRST }; static const char cite_fix_ave_correlate_long[] = "fix ave/correlate/long command: doi:10.1063/1.3491098\n\n" @@ -142,6 +142,8 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS *lmp, int narg, char **arg) : type = AUTOLOWER; else if (strcmp(arg[iarg + 1], "full") == 0) type = FULL; + else if (strcmp(arg[iarg + 1], "first") == 0) + type = FIRST; else error->all(FLERR, iarg_orig + 1, "Unknown fix ave/correlate/long type: {}"); iarg += 2; @@ -262,7 +264,7 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS *lmp, int narg, char **arg) : // npair = # of correlation pairs to calculate - if (type == AUTO) npair = nvalues; + if (type == AUTO || type == FIRST) npair = nvalues; if (type == UPPER || type == LOWER) npair = nvalues * (nvalues - 1) / 2; if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues * (nvalues + 1) / 2; if (type == FULL) npair = nvalues * nvalues; @@ -299,6 +301,9 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < nvalues; i++) for (int j = 0; j < nvalues; j++) fprintf(fp," %s*%s",earg[i],earg[j]); + else if (type == FIRST) + for (int i = 0; i < nvalues; i++) + fprintf(fp," %s*%s",earg[0],earg[i]); fprintf(fp,"\n"); } if (ferror(fp)) @@ -612,6 +617,8 @@ void FixAveCorrelateLong::accumulate() if (i == j) add(ipair++,cvalues[i]); else add(ipair++,cvalues[i],cvalues[j]); } + } else if (type == FIRST) { + for (i=0; i < nvalues; i++) add(i,cvalues[0],cvalues[i]); } last_accumulated_step = update->ntimestep; } diff --git a/src/REPLICA/bosonic_exchange.cpp b/src/REPLICA/bosonic_exchange.cpp new file mode 100644 index 0000000000..6ffbeb8eb7 --- /dev/null +++ b/src/REPLICA/bosonic_exchange.cpp @@ -0,0 +1,431 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Package BosonicExchange + + Purpose Handle Bosonic Exchange in Path Integral Molecular Dynamics + according to J. Chem. Phys. 159, 154107 (2023) (doi:10.1063/5.0173749) + Copyright Hirshberg lab @ Tel Aviv University + Authors Ofir Blumer, Jacob Higer, Yotam Feldman (yotam.feldman at gmail.com), Barak Hirshberg (hirshb at tauex.tau.ac.il) + + Updated Jan-06-2025 + Version 1.0 +------------------------------------------------------------------------- */ + +#include "bosonic_exchange.h" + +#include "domain.h" +#include "error.h" +#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), 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; + + 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); +} + +/* ---------------------------------------------------------------------- */ + +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; +} + +/* ---------------------------------------------------------------------- */ + +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; + + 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]; +} + +/* ---------------------------------------------------------------------- */ + +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; + + 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)); + } + } +} + +/* ---------------------------------------------------------------------- */ + +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)); + } + } + + V_backwards[0] = 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_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::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; + + 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); + + 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; + } +} + +/* ---------------------------------------------------------------------- */ + +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]; + + diff_two_beads(x_first_bead, l, x_last_bead, prev_l, diff_prev); + + 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; + } +} + +/* ---------------------------------------------------------------------- */ + +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_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; + } +} + +/* ---------------------------------------------------------------------- */ + +double BosonicExchange::prim_estimator() +{ + // 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. 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 boolean variable to false amounts to adopting the reduced-beta convention. + double convention_correction = (physical_beta_convention ? 1 : 1.0 / np); + + // 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()); + } + + // For the first bead, add the contribution of the exterior bead (See Equations 4-5 in the + // Supporting Information of Hirshberg et. al., doi.org/10.1073/pnas.1913365116) + 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 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]); } + + 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)); + } + + double sig_denom_m = m * exp(-beta * (V[m] - Elongest)); + + temp_nbosons_array[m] = sig / sig_denom_m; + } + + 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 new file mode 100644 index 0000000000..be36e2910e --- /dev/null +++ b/src/REPLICA/bosonic_exchange.h @@ -0,0 +1,71 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef BOSONIC_EXCHANGE_H +#define BOSONIC_EXCHANGE_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +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); + + double get_potential() const; + double get_bead_spring_energy() const; + + void spring_force(double **f) const; + + 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(); + + const int nbosons; + const int np; + const int bead_num; + const bool apply_minimum_image; + const bool physical_beta_convention; + + 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; +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index 271725e238..037af71990 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -53,17 +53,12 @@ using MathConst::MY_SQRT2; using MathConst::THIRD; using MathSpecial::powint; -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 }; - -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"}}; /* ---------------------------------------------------------------------- */ @@ -248,7 +243,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : removecomflag = 1; else if (strcmp(arg[i + 1], "no") == 0) removecomflag = 0; - } else { + } else if (strcmp(arg[i], "") != 0) { error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix {}", arg[i], style)); } } @@ -445,7 +440,7 @@ void FixPIMDLangevin::init() } planck *= sp; hbar = planck / (MY_2PI); - double beta = 1.0 / (force->boltz * temp); + beta = 1.0 / (force->boltz * temp); double _fbond = 1.0 * np * np / (beta * beta * hbar * hbar); omega_np = np / (hbar * beta) * sqrt(force->mvv2e); @@ -522,7 +517,7 @@ void FixPIMDLangevin::setup(int vflag) else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]); } else if (method == PIMD) { - inter_replica_comm(x); + prepare_coordinates(); if (cmode == SINGLE_PROC) spring_force(); else if (cmode == MULTI_PROC) @@ -682,6 +677,13 @@ void FixPIMDLangevin::final_integrate() /* ---------------------------------------------------------------------- */ +void FixPIMDLangevin::prepare_coordinates() +{ + inter_replica_comm(atom->x); +} + +/* ---------------------------------------------------------------------- */ + void FixPIMDLangevin::post_force(int /*flag*/) { int nlocal = atom->nlocal; @@ -715,7 +717,7 @@ void FixPIMDLangevin::post_force(int /*flag*/) if (mapflag) { for (int i = 0; i < nlocal; i++) { domain->unmap(x[i], image[i]); } } - inter_replica_comm(x); + prepare_coordinates(); spring_force(); compute_spring_energy(); compute_t_prim(); diff --git a/src/REPLICA/fix_pimd_langevin.h b/src/REPLICA/fix_pimd_langevin.h index 869281243f..4f1bac2532 100644 --- a/src/REPLICA/fix_pimd_langevin.h +++ b/src/REPLICA/fix_pimd_langevin.h @@ -29,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; @@ -94,8 +103,9 @@ class FixPIMDLangevin : public Fix { int *counts, *displacements; void comm_init(); + virtual void prepare_coordinates(); void inter_replica_comm(double **ptr); - void spring_force(); + void virtual spring_force(); /* normal-mode operations */ @@ -165,12 +175,12 @@ class FixPIMDLangevin : public Fix { class Compute *c_pe; class Compute *c_press; - void compute_totke(); // 1: kinetic energy - void compute_spring_energy(); // 2: spring elastic energy - void compute_pote(); // 3: potential energy - void compute_tote(); // 4: total energy: 1+2+3 for all the beads + void compute_totke(); // 1: kinetic energy + virtual void compute_spring_energy(); // 2: spring elastic energy + void compute_pote(); // 3: potential energy + void compute_tote(); // 4: total energy: 1+2+3 for all the beads void compute_stress_tensor(); - void compute_t_prim(); + virtual void compute_t_prim(); void compute_t_vir(); void compute_t_cv(); void compute_p_prim(); diff --git a/src/REPLICA/fix_pimd_langevin_bosonic.cpp b/src/REPLICA/fix_pimd_langevin_bosonic.cpp new file mode 100644 index 0000000000..140dd1afc5 --- /dev/null +++ b/src/REPLICA/fix_pimd_langevin_bosonic.cpp @@ -0,0 +1,193 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Package FixPIMDBLangevin + + Purpose Path Integral Molecular Dynamics of Bosons with Langevin Thermostat + Copyright Hirshberg lab @ Tel Aviv University + Authors Ofir Blumer, Jacob Higer, Yotam Feldman (yotam.feldman at gmail.com), Barak Hirshberg (hirshb at tauex.tau.ac.il) + + Updated Jan-06-2025 + Version 1.0 +------------------------------------------------------------------------- */ + +#include "fix_pimd_langevin_bosonic.h" + +#include "bosonic_exchange.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "universe.h" +#include "update.h" + +#include +#include +#include + +using namespace LAMMPS_NS; +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 = 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."); + } + + method = PIMD; + + size_vector = 6; + + memory->create(f_tag_order, nbosons, 3, "FixPIMDBLangevin:f_tag_order"); +} + +/* ---------------------------------------------------------------------- */ + +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) +{ + 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; +} + +/* ---------------------------------------------------------------------- */ + +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]; + + bosonic_exchange->prepare_with_coordinates(me_bead_positions, last_bead_positions, + next_bead_positions, beta_np, ff); +} + +/* ---------------------------------------------------------------------- */ + +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); + + 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(); + + 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(); + } +} + +/* ---------------------------------------------------------------------- */ + +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; +} diff --git a/src/REPLICA/fix_pimd_langevin_bosonic.h b/src/REPLICA/fix_pimd_langevin_bosonic.h new file mode 100644 index 0000000000..fb1f7e7af9 --- /dev/null +++ b/src/REPLICA/fix_pimd_langevin_bosonic.h @@ -0,0 +1,55 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(pimd/langevin/bosonic, FixPIMDBLangevin); +// clang-format on +#else + +#ifndef FIX_PIMDB_LANGEVIN_H +#define FIX_PIMDB_LANGEVIN_H + +#include "fix_pimd_langevin.h" + +namespace LAMMPS_NS { + +class FixPIMDBLangevin : public FixPIMDLangevin { + public: + FixPIMDBLangevin(class LAMMPS *, int, char **); + ~FixPIMDBLangevin(); + + double compute_vector(int) override; + void compute_spring_energy() override; + void compute_t_prim() override; + + char **filtered_args; + int filtered_narg; + + protected: + void prepare_coordinates() override; + void spring_force() override; + + 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_pimd_nvt.cpp b/src/REPLICA/fix_pimd_nvt.cpp index e899c75a27..09c2258c17 100644 --- a/src/REPLICA/fix_pimd_nvt.cpp +++ b/src/REPLICA/fix_pimd_nvt.cpp @@ -43,8 +43,6 @@ using namespace MathConst; using MathSpecial::powint; -enum { PIMD, NMPIMD, CMD }; - /* ---------------------------------------------------------------------- */ FixPIMDNVT::FixPIMDNVT(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) @@ -220,7 +218,7 @@ void FixPIMDNVT::init() const double Plank = force->hplanck; double hbar = Plank / (2.0 * MY_PI) * sp; - double beta = 1.0 / (Boltzmann * nhc_temp); + beta = 1.0 / (Boltzmann * nhc_temp); double _fbond = 1.0 * np / (beta * beta * hbar * hbar); omega_np = sqrt((double) np) / (hbar * beta) * sqrt(force->mvv2e); @@ -271,12 +269,20 @@ void FixPIMDNVT::final_integrate() /* ---------------------------------------------------------------------- */ +void FixPIMDNVT::prepare_coordinates() +{ + comm_exec(atom->x); +} + +/* ---------------------------------------------------------------------- */ + void FixPIMDNVT::post_force(int /*flag*/) { for (int i = 0; i < atom->nlocal; i++) for (int j = 0; j < 3; j++) atom->f[i][j] /= np; - comm_exec(atom->x); + prepare_coordinates(); + pre_spring_force_estimators(); spring_force(); if (method == CMD || method == NMPIMD) { @@ -533,6 +539,24 @@ void FixPIMDNVT::nmpimd_transform(double **src, double **des, double *vector) /* ---------------------------------------------------------------------- */ +void FixPIMDNVT::pre_spring_force_estimators(){ + vir_estimator(); +} + +/* ---------------------------------------------------------------------- */ + +void FixPIMDNVT::vir_estimator() { + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + + virial = 0.0; + for (int i = 0; i < nlocal; i++) { + virial += -0.5 * (x[i][0] * f[i][0] + x[i][1] * f[i][1] + x[i][2] * f[i][2]); + } +} +/* ---------------------------------------------------------------------- */ + void FixPIMDNVT::spring_force() { spring_energy = 0.0; @@ -546,8 +570,6 @@ void FixPIMDNVT::spring_force() double *xlast = buf_beads[x_last]; double *xnext = buf_beads[x_next]; - virial = 0.0; - for (int i = 0; i < nlocal; i++) { double delx1 = xlast[0] - x[i][0]; double dely1 = xlast[1] - x[i][1]; @@ -567,8 +589,6 @@ void FixPIMDNVT::spring_force() double dy = dely1 + dely2; double dz = delz1 + delz2; - virial += -0.5 * (x[i][0] * f[i][0] + x[i][1] * f[i][1] + x[i][2] * f[i][2]); - f[i][0] -= (dx) *ff; f[i][1] -= (dy) *ff; f[i][2] -= (dz) *ff; diff --git a/src/REPLICA/fix_pimd_nvt.h b/src/REPLICA/fix_pimd_nvt.h index 2cd1397562..a0bff062ab 100644 --- a/src/REPLICA/fix_pimd_nvt.h +++ b/src/REPLICA/fix_pimd_nvt.h @@ -30,6 +30,8 @@ class FixPIMDNVT : public Fix { FixPIMDNVT(class LAMMPS *, int, char **); ~FixPIMDNVT() override; + enum { PIMD, NMPIMD, CMD }; + int setmask() override; void init() override; @@ -53,7 +55,6 @@ class FixPIMDNVT : public Fix { void unpack_forward_comm(int, int, double *) override; protected: - int method; int np; double inverse_np; @@ -62,8 +63,10 @@ class FixPIMDNVT : public Fix { double omega_np, fbond, spring_energy, sp, virial; int x_last, x_next; - - void spring_force(); + virtual void prepare_coordinates(); + virtual void pre_spring_force_estimators(); + virtual void spring_force(); + void vir_estimator(); /* fictitious mass */ @@ -101,6 +104,7 @@ class FixPIMDNVT : public Fix { int nhc_nchain; bool nhc_ready; double nhc_temp, dtv, dtf, t_sys; + double beta; double **nhc_eta; /* coordinates of NH chains for ring-polymer beads */ double **nhc_eta_dot; /* velocities of NH chains */ diff --git a/src/REPLICA/fix_pimd_nvt_bosonic.cpp b/src/REPLICA/fix_pimd_nvt_bosonic.cpp new file mode 100644 index 0000000000..1458a4b72f --- /dev/null +++ b/src/REPLICA/fix_pimd_nvt_bosonic.cpp @@ -0,0 +1,99 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Package FixPIMDBNVT + Purpose Path Integral Molecular Dynamics of Bosons with Nose-Hoover Thermostat + Copyright Hirshberg lab @ Tel Aviv University + Authors Ofir Blumer, Jacob Higer, Yotam Feldman (yotam.feldman at gmail.com), Barak Hirshberg (hirshb at tauex.tau.ac.il) + + Updated Jan-06-2025 + Version 1.0 +------------------------------------------------------------------------- */ + +#include "fix_pimd_nvt_bosonic.h" + +#include "bosonic_exchange.h" + +#include "atom.h" +#include "error.h" +#include "force.h" +#include "universe.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +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.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"); + } +} + +/* ---------------------------------------------------------------------- */ + +FixPIMDBNVT::~FixPIMDBNVT() +{ + delete bosonic_exchange; +} + +/* ---------------------------------------------------------------------- */ + +void FixPIMDBNVT::pre_spring_force_estimators() +{ + FixPIMDNVT::pre_spring_force_estimators(); + spring_energy = bosonic_exchange->get_bead_spring_energy(); + prim = bosonic_exchange->prim_estimator(); +} + +/* ---------------------------------------------------------------------- */ + +void FixPIMDBNVT::prepare_coordinates() +{ + comm_exec(atom->x); + double **x = atom->x; + 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); +} + +void FixPIMDBNVT::spring_force() +{ + double **f = atom->f; + + bosonic_exchange->spring_force(f); +} + +/* ---------------------------------------------------------------------- */ + +double FixPIMDBNVT::compute_vector(int 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!"); + } + + return 0.0; +} diff --git a/src/REPLICA/fix_pimd_nvt_bosonic.h b/src/REPLICA/fix_pimd_nvt_bosonic.h new file mode 100644 index 0000000000..3d4adb2a7d --- /dev/null +++ b/src/REPLICA/fix_pimd_nvt_bosonic.h @@ -0,0 +1,46 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(pimd/nvt/bosonic,FixPIMDBNVT); +// clang-format on +#else + +#ifndef FIX_PIMDB_NVT_H +#define FIX_PIMDB_NVT_H + +#include "fix_pimd_nvt.h" + +namespace LAMMPS_NS { + +class FixPIMDBNVT : public FixPIMDNVT { + public: + 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; + + private: + class BosonicExchange *bosonic_exchange; + double prim; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/fix_ave_correlate.cpp b/src/fix_ave_correlate.cpp index 10280682df..a48788b777 100644 --- a/src/fix_ave_correlate.cpp +++ b/src/fix_ave_correlate.cpp @@ -36,7 +36,7 @@ using namespace LAMMPS_NS; using namespace FixConst; enum { ONE, RUNNING }; -enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL }; +enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL, FIRST }; /* ---------------------------------------------------------------------- */ @@ -108,6 +108,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+1],"auto/upper") == 0) type = AUTOUPPER; else if (strcmp(arg[iarg+1],"auto/lower") == 0) type = AUTOLOWER; else if (strcmp(arg[iarg+1],"full") == 0) type = FULL; + else if (strcmp(arg[iarg+1], "first") == 0) type = FIRST; else error->all(FLERR, iarg+1, "Unknown fix ave/correlate type: {}"); iarg += 2; } else if (strcmp(arg[iarg],"ave") == 0) { @@ -203,7 +204,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) : // npair = # of correlation pairs to calculate - if (type == AUTO) npair = nvalues; + if (type == AUTO || type == FIRST) npair = nvalues; if (type == UPPER || type == LOWER) npair = nvalues*(nvalues-1)/2; if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues*(nvalues+1)/2; if (type == FULL) npair = nvalues*nvalues; @@ -242,6 +243,9 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < nvalues; i++) for (int j = 0; j < nvalues; j++) fprintf(fp," %s*%s",earg[i],earg[j]); + else if (type == FIRST) + for (int i = 0; i < nvalues; i++) + fprintf(fp," %s*%s",earg[0],earg[i]); fprintf(fp,"\n"); } if (ferror(fp)) @@ -579,6 +583,16 @@ void FixAveCorrelate::accumulate() m--; if (m < 0) m = nrepeat-1; } + } else if (type == FIRST) { + m = n = lastindex; + for (k = 0; k < nsample; k++) { + ipair = 0; + for (i = 0; i < nvalues; i++) { + corr[k][ipair++] += cvalues[m][0]*cvalues[n][i]; + } + m--; + if (m < 0) m = nrepeat-1; + } } }