[lammps-users] compute temperature of lubricant constrained by fix shake

Hello, everyone, this is rufer. I’m modeling the temperature of lubricant which is under shearing in the x-direction and pressing in the z-direction. the lubricant which I used is Hexadecane, and the forcefield is OPLS-AA force field. The C-H bond is constrained by fix shake. when I compute the entire temperature of the whole lubricant by compute temp/profile, the temperature is 330K. the thermostat is 300K, so 330K is right when lubricant is under friction. But when I use fix ave/chunk to compute the distribution of temperature in the z-direction,the temperature is only 250K which is not correct. I also use bias keyword in fix ave/chunk to correct the computation of temperature distribution, It dose not works. So can anyone give me some advice?
Best regards.

From the documentation of compute temp/chunk:

Note that currently the global and per-chunk temperatures calculated by this compute only include translational degrees of freedom for each atom. No rotational degrees of freedom are included for finite-size particles. Also no degrees of freedom are subtracted for any velocity bias or constraints that are applied, such as compute temp/partial, or fix shake or fix rigid. This is because those degrees of freedom (e.g. a constrained bond) could apply to sets of atoms that are both included and excluded from a specific chunk, and hence the concept is somewhat ill-defined. In some cases, you can use the adof and cdof keywords to adjust the calculated degrees of freedom appropriately, as explained below.

Thanks for your advice. Actually, I saw the description as follows:
The adof and cdof keywords define the values used in the degree of freedom (DOF) formulas used for the global or per-chunk temperature, as described above. They can be used to calculate a more appropriate temperature for some kinds of chunks. Here are 3 examples:

If spatially binned chunks contain some number of water molecules and fix shake is used to make each molecule rigid, then you could calculate a temperature with 6 degrees of freedom (DOF) (3 translational, 3 rotational) per molecule by setting adof to 2.0.

If compute temp/partial is used with the bias keyword to only allow the x component of velocity to contribute to the temperature, then adof = 1.0 would be appropriate.

If each chunk consists of a large molecule, with some number of its bonds constrained by fix shake or the entire molecule by fix rigid/small, adof = 0.0 and cdof could be set to the remaining degrees of freedom for the entire molecule (entire chunk in this case), e.g. 6 for 3d, or 3 for 2d, for a rigid molecule.

but I still not know how to determine the adof and cdof of Hexadecane, should I change the adof to 2 for atoms in Hexadecane? But H atom is bond to C atom, it is only constrained once, for CH3, the C atom is constrained three times. So I am very confused about the definement of adof? Any advice will be appreciated.

Best regards.

Axel Kohlmeyer <[email protected]> 于2021年6月16日周三 下午2:33写道:

There is no simple answer to that. The main issue is that the temperature compute has no information about which atoms are part of constrained bonds and thus have to have half a DOF subtracted to give the total number of DOFs per chunk to compute the temperature.

Given that temperature is not a well defined property at this scale (only kinetic energy is, we compute temperature based on the assumptions that all degrees of freedom contribute equally and that averaging over time is equivalent to averaging over space) and that you have a homogeneous system you may be able to get a suitable approximation by determining what the average number of degrees of freedom per atom with constraints applied is and thus how many DOFs must be subtracted per atom. “adof” need not be an integer.


Thank you for the clarification.

Axel Kohlmeyer <[email protected]> 于2021年6月16日周三 下午8:57写道: