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.0*nper/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