From ab074273391675b89581e12c4b5b2bbf64fcfd62 Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 09:52:05 +1000 Subject: [PATCH 01/11] Corrected default extra_dof --- src/compute.h | 2 +- src/compute_temp_profile.cpp | 6 ++++++ src/compute_temp_profile.h | 1 + 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/compute.h b/src/compute.h index 71c07737d4..810a1464d3 100644 --- a/src/compute.h +++ b/src/compute.h @@ -104,7 +104,7 @@ class Compute : protected Pointers { Compute(class LAMMPS *, int, char **); virtual ~Compute(); void modify_params(int, char **); - void reset_extra_dof(); + virtual void reset_extra_dof(); virtual void init() = 0; virtual void init_list(int, class NeighList *) {} diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 6938560359..5f95ef29e5 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -575,6 +575,12 @@ void ComputeTempProfile::bin_assign() /* ---------------------------------------------------------------------- */ +void ComputeTempProfile::reset_extra_dof() { + extra_dof = 0.0; +} + +/* ---------------------------------------------------------------------- */ + double ComputeTempProfile::memory_usage() { double bytes = (double)maxatom * sizeof(int); diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h index f0c07bbd48..92e76f9095 100644 --- a/src/compute_temp_profile.h +++ b/src/compute_temp_profile.h @@ -34,6 +34,7 @@ class ComputeTempProfile : public Compute { void compute_vector(); void compute_array(); + void reset_extra_dof(); void remove_bias(int, double *); void remove_bias_thr(int, double *, double *); void remove_bias_all(); From c643389ec4db55b77a789afa38a3a6f4d95a2d2c Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 10:01:45 +1000 Subject: [PATCH 02/11] Treat extra_dof as system-wide in compute_array for consistency, and include fix_dof --- src/compute_temp_profile.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 5f95ef29e5..5bb415f5ed 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -333,11 +333,15 @@ void ComputeTempProfile::compute_array() MPI_Allreduce(tbin,tbinall,nbins,MPI_DOUBLE,MPI_SUM,world); - int nper = domain->dimension; + double totcount = 0.0; for (i = 0; i < nbins; i++) { array[i][0] = binave[i][ncount-1]; + totcount += array[i][0]; + } + double nper = domain->dimension - (extra_dof + fix_dof)/totcount; + for (i = 0; i < nbins; i++) { if (array[i][0] > 0.0) { - dof = nper*array[i][0] - extra_dof; + dof = nper*array[i][0] - domain->dimension; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; array[i][1] = tfactor*tbinall[i]; From 5e4dd5321c423fa0abf80e18a1a3a3930c24dc10 Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 11:43:25 +1000 Subject: [PATCH 03/11] Using local dof and tfactor for compute_array to prevent overwrite --- src/compute_temp_profile.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 5bb415f5ed..52f2f7de6f 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -339,12 +339,13 @@ void ComputeTempProfile::compute_array() totcount += array[i][0]; } double nper = domain->dimension - (extra_dof + fix_dof)/totcount; + double dofbin, tfactorbin; for (i = 0; i < nbins; i++) { if (array[i][0] > 0.0) { - dof = nper*array[i][0] - domain->dimension; - if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); - else tfactor = 0.0; - array[i][1] = tfactor*tbinall[i]; + dofbin = nper*array[i][0] - domain->dimension; + if (dofbin > 0) tfactorbin = force->mvv2e / (dofbin * force->boltz); + else tfactorbin = 0.0; + array[i][1] = tfactorbin*tbinall[i]; } else array[i][1] = 0.0; } } From b1b7f7a248c8c6c13917f98ae750b6e9971d0180 Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 11:47:32 +1000 Subject: [PATCH 04/11] Fixed treatment of DoF when streaming velocity not subtracted in some dimensions --- src/compute_temp_profile.cpp | 13 ++++++++++--- src/compute_temp_profile.h | 1 + 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 52f2f7de6f..9790faf9b5 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -118,6 +118,9 @@ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : nbins = nbinx*nbiny*nbinz; if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command"); + nconstraints = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1); + nconstraints *= nbins; + memory->create(vbin,nbins,ncount,"temp/profile:vbin"); memory->create(binave,nbins,ncount,"temp/profile:binave"); @@ -196,9 +199,10 @@ void ComputeTempProfile::dof_compute() natoms_temp = group->count(igroup); dof = domain->dimension * natoms_temp; - // subtract additional d*Nbins DOF, as in Evans and Morriss paper + // subtract additional Nbins DOF for each adjusted direction, + // as in Evans and Morriss paper - dof -= extra_dof + fix_dof + domain->dimension*nbins; + dof -= extra_dof + fix_dof + nconstraints; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -342,7 +346,7 @@ void ComputeTempProfile::compute_array() double dofbin, tfactorbin; for (i = 0; i < nbins; i++) { if (array[i][0] > 0.0) { - dofbin = nper*array[i][0] - domain->dimension; + dofbin = nper*array[i][0] - nconstraints/nbins; if (dofbin > 0) tfactorbin = force->mvv2e / (dofbin * force->boltz); else tfactorbin = 0.0; array[i][1] = tfactorbin*tbinall[i]; @@ -582,6 +586,9 @@ void ComputeTempProfile::bin_assign() void ComputeTempProfile::reset_extra_dof() { extra_dof = 0.0; + if (xflag == 0) extra_dof++; + if (yflag == 0) extra_dof++; + if (zflag == 0 && domain->dimension == 3) extra_dof++; } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h index 92e76f9095..5e62d923b7 100644 --- a/src/compute_temp_profile.h +++ b/src/compute_temp_profile.h @@ -48,6 +48,7 @@ class ComputeTempProfile : public Compute { int nbinx,nbiny,nbinz,nbins; int ivx,ivy,ivz; double tfactor; + double nconstraints; int box_change,triclinic; int *periodicity; From aad0a9a0f303f62bee0fdd41c6afbc3ac48831aa Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 14:14:33 +1000 Subject: [PATCH 05/11] Updated documentation to reflect changes in compute_temp_profile --- doc/src/compute_temp_profile.rst | 32 +++++++++++++++++++++----------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/doc/src/compute_temp_profile.rst b/doc/src/compute_temp_profile.rst index 0fa8ce5807..3f045f2bda 100644 --- a/doc/src/compute_temp_profile.rst +++ b/doc/src/compute_temp_profile.rst @@ -76,13 +76,20 @@ velocity for each atom. Note that if there is only one atom in the bin, its thermal velocity will thus be 0.0. After the spatially-averaged velocity field has been subtracted from -each atom, the temperature is calculated by the formula KE = (dim\*N -- dim\*Nx\*Ny\*Nz) k T/2, where KE = total kinetic energy of the group of -atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the -simulation, N = number of atoms in the group, k = Boltzmann constant, -and T = temperature. The dim\*Nx\*Ny\*Nz term are degrees of freedom -subtracted to adjust for the removal of the center-of-mass velocity in -each of Nx\*Ny\*Nz bins, as discussed in the :ref:`(Evans) ` paper. +each atom, the temperature is calculated by the formula +KE = (dim\*N - stream\*Nx\*Ny\*Nz - extra ) k T/2, where KE = total +kinetic energy of the group of atoms (sum of 1/2 m v\^2), dim = 2 +or 3 = dimensionality of the simulation, stream = 0, 1, 2 or 3 for +streaming velocity subtracted in 0, 1, 2 or 3 dimensions, extra = extra +degrees-of-freedom, N = number of atoms in the group, k = Boltzmann +constant, and T = temperature. The stream\*Nx\*Ny\*Nz term is degrees +of freedom subtracted to adjust for the removal of the center-of-mass +velocity in each direction of the Nx\*Ny\*Nz bins, as discussed in the +:ref:`(Evans) ` paper. The extra term defaults to (dim - stream) +and accounts for overall conservation of center-of-mass velocity across +the group in directions where streaming velocity is not subtracted. This +can be altered using the *extra* option of the +:doc:`compute_modify ` command. If the *out* keyword is used with a *tensor* value, which is the default, a kinetic energy tensor, stored as a 6-element vector, is @@ -123,10 +130,13 @@ needed, the subtracted degrees-of-freedom can be altered using the .. note:: When using the *out* keyword with a value of *bin*\ , the - calculated temperature for each bin does not include the - degrees-of-freedom adjustment described in the preceding paragraph, - for fixes that constrain molecular motion. It does include the - adjustment due to the *extra* option, which is applied to each bin. + calculated temperature for each bin includes the degrees-of-freedom + adjustment described in the preceding paragraph for fixes that + constrain molecular motion, as well as the adjustment due to + the *extra* option (which defaults to dim - stream as described above), + by fractionally applying them based on the fraction of atoms in each + bin, so that the degrees-of-freedom summed over all bins is consistent + with the degrees-of-freedom in the scalar temperature calculation. See the :doc:`Howto thermostat ` doc page for a discussion of different ways to compute temperature and perform From 471cfa8ac316f2afe3a6fe826c62186947ea7bc4 Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 14:31:17 +1000 Subject: [PATCH 06/11] Fixed inconsistent default extra_dof value --- src/compute_temp_profile.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index 9790faf9b5..a2dfb741e4 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -120,6 +120,7 @@ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : nconstraints = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1); nconstraints *= nbins; + reset_extra_dof(); memory->create(vbin,nbins,ncount,"temp/profile:vbin"); memory->create(binave,nbins,ncount,"temp/profile:binave"); From 97f90f114697595e221bd725aa77fd60f5c28aa6 Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 15:56:48 +1000 Subject: [PATCH 07/11] Cleaned up math for clarity --- src/compute_temp_profile.cpp | 7 +++---- src/compute_temp_profile.h | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index a2dfb741e4..a64b0e0321 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -118,8 +118,7 @@ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : nbins = nbinx*nbiny*nbinz; if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command"); - nconstraints = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1); - nconstraints *= nbins; + nstreaming = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1); reset_extra_dof(); memory->create(vbin,nbins,ncount,"temp/profile:vbin"); @@ -203,7 +202,7 @@ void ComputeTempProfile::dof_compute() // subtract additional Nbins DOF for each adjusted direction, // as in Evans and Morriss paper - dof -= extra_dof + fix_dof + nconstraints; + dof -= extra_dof + fix_dof + nstreaming*nbins; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -347,7 +346,7 @@ void ComputeTempProfile::compute_array() double dofbin, tfactorbin; for (i = 0; i < nbins; i++) { if (array[i][0] > 0.0) { - dofbin = nper*array[i][0] - nconstraints/nbins; + dofbin = nper*array[i][0] - nstreaming; if (dofbin > 0) tfactorbin = force->mvv2e / (dofbin * force->boltz); else tfactorbin = 0.0; array[i][1] = tfactorbin*tbinall[i]; diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h index 5e62d923b7..c6d2243148 100644 --- a/src/compute_temp_profile.h +++ b/src/compute_temp_profile.h @@ -48,7 +48,7 @@ class ComputeTempProfile : public Compute { int nbinx,nbiny,nbinz,nbins; int ivx,ivy,ivz; double tfactor; - double nconstraints; + double nstreaming; int box_change,triclinic; int *periodicity; From a3e204a99df75dcafb23be283f8b49cf00cf95ea Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Wed, 17 Mar 2021 16:23:04 +1000 Subject: [PATCH 08/11] Simplified reset_extra_dof --- src/compute_temp_profile.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index a64b0e0321..da0363c09e 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -585,10 +585,7 @@ void ComputeTempProfile::bin_assign() /* ---------------------------------------------------------------------- */ void ComputeTempProfile::reset_extra_dof() { - extra_dof = 0.0; - if (xflag == 0) extra_dof++; - if (yflag == 0) extra_dof++; - if (zflag == 0 && domain->dimension == 3) extra_dof++; + extra_dof = domain->dimension - nstreaming; } /* ---------------------------------------------------------------------- */ From 907bcd16a84c17d40a30e85915611eab2beb1d7a Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 22 Apr 2022 14:21:57 -0600 Subject: [PATCH 09/11] Update compute_temp_profile.rst --- doc/src/compute_temp_profile.rst | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/doc/src/compute_temp_profile.rst b/doc/src/compute_temp_profile.rst index 3f045f2bda..7ccf403eb9 100644 --- a/doc/src/compute_temp_profile.rst +++ b/doc/src/compute_temp_profile.rst @@ -77,17 +77,17 @@ bin, its thermal velocity will thus be 0.0. After the spatially-averaged velocity field has been subtracted from each atom, the temperature is calculated by the formula -KE = (dim\*N - stream\*Nx\*Ny\*Nz - extra ) k T/2, where KE = total -kinetic energy of the group of atoms (sum of 1/2 m v\^2), dim = 2 -or 3 = dimensionality of the simulation, stream = 0, 1, 2 or 3 for -streaming velocity subtracted in 0, 1, 2 or 3 dimensions, extra = extra -degrees-of-freedom, N = number of atoms in the group, k = Boltzmann -constant, and T = temperature. The stream\*Nx\*Ny\*Nz term is degrees +*KE* = (*dim\*N* - *Ns\*Nx\*Ny\*Nz* - *extra* ) *k* *T*/2, where *KE* = total +kinetic energy of the group of atoms (sum of 1/2 *m* *v*\^2), *dim* = 2 +or 3 = dimensionality of the simulation, *Ns* = 0, 1, 2 or 3 for +streaming velocity subtracted in 0, 1, 2 or 3 dimensions, *extra* = extra +degrees-of-freedom, *N* = number of atoms in the group, *k* = Boltzmann +constant, and *T* = temperature. The *Ns\*Nx\*Ny\*Nz* term is degrees of freedom subtracted to adjust for the removal of the center-of-mass -velocity in each direction of the Nx\*Ny\*Nz bins, as discussed in the -:ref:`(Evans) ` paper. The extra term defaults to (dim - stream) +velocity in each direction of the *Nx\*Ny\*Nz* bins, as discussed in the +:ref:`(Evans) ` paper. The extra term defaults to (*dim* - *Ns*) and accounts for overall conservation of center-of-mass velocity across -the group in directions where streaming velocity is not subtracted. This +the group in directions where streaming velocity is *not* subtracted. This can be altered using the *extra* option of the :doc:`compute_modify ` command. @@ -95,9 +95,9 @@ If the *out* keyword is used with a *tensor* value, which is the default, a kinetic energy tensor, stored as a 6-element vector, is also calculated by this compute for use in the computation of a pressure tensor. The formula for the components of the tensor is the -same as the above formula, except that v\^2 is replaced by vx\*vy for -the xy component, etc. The 6 components of the vector are ordered xx, -yy, zz, xy, xz, yz. +same as the above formula, except that *v*\^2 is replaced by *vx\*vy* for +the xy component, etc. The 6 components of the vector are ordered *xx, +yy, zz, xy, xz, yz.* If the *out* keyword is used with a *bin* value, the count of atoms and computed temperature for each bin are stored for output, as an @@ -133,10 +133,12 @@ needed, the subtracted degrees-of-freedom can be altered using the calculated temperature for each bin includes the degrees-of-freedom adjustment described in the preceding paragraph for fixes that constrain molecular motion, as well as the adjustment due to - the *extra* option (which defaults to dim - stream as described above), + the *extra* option (which defaults to *dim* - *Ns* as described above), by fractionally applying them based on the fraction of atoms in each - bin, so that the degrees-of-freedom summed over all bins is consistent - with the degrees-of-freedom in the scalar temperature calculation. + bin. As a result, the bin degrees-of-freedom summed over all bins exactly + equals the degrees-of-freedom used in the scalar temperature calculation, + :math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature + is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T` See the :doc:`Howto thermostat ` doc page for a discussion of different ways to compute temperature and perform From 16dd89c641dc804f6d39e5c87bd9eacd20fae9eb Mon Sep 17 00:00:00 2001 From: Stephen Sanderson Date: Sat, 23 Apr 2022 10:30:44 +1000 Subject: [PATCH 10/11] Added qualifier to assertions + note about inhomogeneous rigid systems --- doc/src/compute_temp_profile.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/doc/src/compute_temp_profile.rst b/doc/src/compute_temp_profile.rst index 7ccf403eb9..1be82efbdd 100644 --- a/doc/src/compute_temp_profile.rst +++ b/doc/src/compute_temp_profile.rst @@ -138,7 +138,12 @@ needed, the subtracted degrees-of-freedom can be altered using the bin. As a result, the bin degrees-of-freedom summed over all bins exactly equals the degrees-of-freedom used in the scalar temperature calculation, :math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature - is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T` + is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`, so long as + all bins contain more degrees of freedom than the number of constraints per bin. + Note that as degrees-of-freedom subtracted due to fixes that contstrain + molecular motion are distributed equally between bins, the reported + temperature within a bin may not be accurate for systems containing + inhomogeneously distributed rigid molecules. See the :doc:`Howto thermostat ` doc page for a discussion of different ways to compute temperature and perform From 5aebd151b69c6a3b534d9963b8d27ce758d9fa66 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Fri, 22 Apr 2022 18:47:13 -0600 Subject: [PATCH 11/11] Update compute_temp_profile.rst --- doc/src/compute_temp_profile.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/src/compute_temp_profile.rst b/doc/src/compute_temp_profile.rst index 1be82efbdd..b773ee67ee 100644 --- a/doc/src/compute_temp_profile.rst +++ b/doc/src/compute_temp_profile.rst @@ -138,12 +138,12 @@ needed, the subtracted degrees-of-freedom can be altered using the bin. As a result, the bin degrees-of-freedom summed over all bins exactly equals the degrees-of-freedom used in the scalar temperature calculation, :math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature - is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`, so long as - all bins contain more degrees of freedom than the number of constraints per bin. - Note that as degrees-of-freedom subtracted due to fixes that contstrain - molecular motion are distributed equally between bins, the reported - temperature within a bin may not be accurate for systems containing - inhomogeneously distributed rigid molecules. + is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`. + These relations will breakdown in cases where the adjustment + exceeds the actual number of degrees-of-freedom in a bin. This could happen + if a bin is empty or in situations where rigid molecules + are non-uniformly distributed, in which case the reported + temperature within a bin may not be accurate. See the :doc:`Howto thermostat ` doc page for a discussion of different ways to compute temperature and perform