problem with compute temp/partial?

Dear all,

I'm getting strange results from the compute temp/partial command,
where the temperature computed along the three directions
(temp/partial 1 1 1) is not the average of the temperature computed
along each of the three directions (temp/partial 1 0 0, temp/partial 0
1 0, temp/partial 0 0 1). It seems that their ratio remains constant
though (in my case, 1.01015...).

My system is a graphene sheet (200 atoms) in a triclinic box, modeled
with the tersoff potential, see enclosed archive with input and output
files. I'm using the current svn version of lammps that I just
downloaded and compiled.

Could there be a problem with compute temp/partial or am I missing something?

Best regards,
Laurent

temp_partial.tgz (16.3 KB)

Dear all,

I'm getting strange results from the compute temp/partial command,
where the temperature computed along the three directions
(temp/partial 1 1 1) is not the average of the temperature computed
along each of the three directions (temp/partial 1 0 0, temp/partial 0
1 0, temp/partial 0 0 1). It seems that their ratio remains constant
though (in my case, 1.01015...).

​since you have 200 atoms and 1.01015 is 203/200, this is indeed suspicious
and hints at the global translation DOFs not being removed consistently.
the following change should correct that:​

diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp
index 845b10e..9ae7184 100644
--- a/src/compute_temp_partial.cpp
+++ b/src/compute_temp_partial.cpp
@@ -77,7 +77,7 @@ void ComputeTempPartial::dof_compute()
   natoms_temp = group->count(igroup);
   int nper = xflag+yflag+zflag;
   dof = nper * natoms_temp;
- dof -= (1.0*nper/domain->dimension)*fix_dof + extra_dof;
+ dof -= (1.0*nper/domain->dimension)*fix_dof + (extra_dof -
domain->dimension) + nper;
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
}

an alternate approach would be:

diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp
index 845b10e..9ae7184 100644
--- a/src/compute_temp_partial.cpp
+++ b/src/compute_temp_partial.cpp
@@ -77,7 +77,7 @@ void ComputeTempPartial::dof_compute()
   natoms_temp = group->count(igroup);
   int nper = xflag+yflag+zflag;
   dof = nper * natoms_temp;
- dof -= (1.0*nper/domain->dimension)*fix_dof + extra_dof;
+ dof -= (1.0*nper/domain->dimension)*fix_dof +
extra_dof/domain->dimension*nper;
   if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
   else tfactor = 0.0;
}

the problem here is, that there is no unambiguous way to know how the DOFs
in extra_dof should be divided across the directions. i personally would
prefer the first option, as that means that the explicitly added/removed
DOFs that are given in fix_modify will be handled by this compute as given
in the input, whereas in the second case, they would be evenly distributed
across the dimensions (but consistent with automatically removed DOFs from
fixes like fix rigid).

axel.

With the current code, a user like Laurent who is paying attention can get the behavior they want by doing:

compute_modify TempPartialX extra 1
compute_modify TempPartialAll extra 3

and most users won’t care at all, so I am not sure it worth changing. If the code was changed, I would recommend Axel’s second choice, but written as:

dof -= (1.0nper/domain->dimension)(fix_dof + extra_dof);

In that case, the following sentence should be added to the compute_modify doc

“The default is 2 or 3 for 2d or 3d systems which is a correction factor for an ensemble of velocities with zero total linear momentum. For compute temp/partial, if one or more velocity components are excluded, the default value is reduced accordingly.”

Aidan