[lammps-users] compute_modify extra/dof command for a system in an external field

lammps-users · May 20, 2021, 8:15pm

Dear LAMMPS developers,

We are studying a system made of a few particles into a framework that is kept frozen in our simulation. We are aware that LAMMPS is imposing the conservation of the total linear momentum[1], assuming the system under study does not have any external force applied. This is performed using the variable “extra_dof” that, by default, is equal to the dimension of the system, and the degree of freedom is 3N-3.


extra_dof = domain->dimension

dof = domain->dimension * natoms_temp

dof -= extra_dof + fix_dof

For very few particles under an external field, the “extra_dof” is changing the dynamics that we want to study, leading to lower kinetic energy. We have tried “compute_modify extra/dof 0” and “compute modify extra/dof 0 dynamic/dof yes” from the manual[2]. Those modifications help in retrieving the correct value for the degree of freedom (3N) but we still get the same lower kinetic energy and a lower temperature than what we expect from canonical distribution. Results with/without “compute_modify extra/dof” are attached below, and a brief input file is also attached.

We are wondering if the command “compute_modify extra/dof” should be applied differently, or some other commands can help us to get the correct canonical distribution for our case. We will be happy to give more details on our tests if it is needed. We would be happy to modify the source code if this is the only way to get what we need.

We really appreciate your help and the time.

Best greetings,

Henglu

[1] https://lammps.sandia.gov/doc/fix_nh.html

[2] https://lammps.sandia.gov/doc/compute_modify.html#compute-modify-command

Canonical distribution using “lammps-stable_3Mar2020”

Canonical distribution

temperature 298.00 K

kinetic energy 3.72 kJ/mol

LAMMPS SETTINGS

compute allTemp all temp

compute CH4Temp CH4 temp

compute framTemp fram temp

compute CH4ke CH4 ke

fix 1 CH4 nvt temp 298 298 100*dt

Average properties from log.lammps

temp 14.00 K

alltemp 14.00 K

CH4temp 297.84 K

framtemp 0.00 K

ke 3.48 kJ/mol

CH4ke 3.48 kJ/mol

Average Ek from CH4 velocities

ke 3.47 kJ/mol

Ek distribution (shown in the fig)

Canonical distribution using “lammps-stable_3Mar2020” with compute_modify extra/dof

Canonical distribution

temperature 298.00 K

kinetic energy 3.72 kJ/mol

LAMMPS SETTINGS

compute allTemp all temp

compute CH4Temp CH4 temp

compute framTemp fram temp

compute CH4ke CH4 ke

compute_modify CH4Temp extra/dof 0 dynamic/dof yes

fix 1 CH4 nvt temp 298 298 100*dt

Average properties from log.lammps

temp 14.00 K

alltemp 14.00 K

CH4temp 279.22 K

framtemp 0.00 K

ke 3.48 kJ/mol

CH4ke 3.48 kJ/mol

Average Ek from CH4 velocities

ke 3.47 kJ/mol

Ek distribution (shown in the fig)

akohlmey · May 20, 2021, 8:59pm

I am copying Aidan, who is much more an expert in statistical thermodynamics than me.

The removal of the 3 global translational degrees of freedom is based on the system being invariant to translation.
When you apply an external field or some other external force or do not have periodic boundary conditions, you don’t have this invariance anymore. You also don’t have a canonical ensemble (despite using fix nvt).

There is not enough information given to make any more specific statements, at least at my level of understanding.

Axel.

lammps-users · June 8, 2021, 2:13pm

Dear Axel,

Thanks for your reply.

We are studying the diffusion of particles inside a framework that is kept fixed during the simulation. By default, LAMMPS reduces the total degrees of freedom by 3 for linear momentum conservation. Since we are interested in the diffusion of particles inside the fixed framework, we don’t want to subtract the center-of-mass motions.

As far as we understand from the manual[1], this could be achieved by “compute modify extra/dof 0” or “compute modify extra/dof 0 dynamic/dof yes”. We have added some print statements to the code in compute.cpp and compute_temp.cpp[2] to check if we get the correct degree of freedom. We noticed that the values of extra_dof are not consistent in the first iteration. If instead, we remove the “extra_dof” variable from the equation for computing the degrees of freedom [2] (enforcing a value of 0 all the time including the first iteration), we obtain the results we expect. Is there a way to use “compute modify extra/dof” in the inputs to achieve what we are looking for?

You find the lammps input file we used below. Any hints/suggestions would be greatly appreciated.

Best wishes,
Henglu

[1] https://lammps.sandia.gov/doc/compute_modify.html#compute-modify-command
[2] https://github.com/lammps/lammps/compare/master…HengluXu:master_print

Possibly related:
https://lammps.sandia.gov/threads/msg58151.html
https://lammps.sandia.gov/threads/msg65987.html

in.MFI_CH4 (1.97 KB)

data.MFI_CH4 (20.2 KB)

akohlmey · June 8, 2021, 2:56pm

Henglu,

After looking at your inputs I believe the problem is a different one.

When you define your nvt fix, it will create its own temperature compute which is used to determine the difference to the target temperature.
Now you change the extra DOFs on your own temperature compute, that will adjust the results for that compute but not for the thermostat.
You either need to apply the compute_modify to the temperature compute for the fix (compute_modify eqCH4_temp extra/dof 0) or use
the custom temperature compute that you have defined (fix_modify eqCH4 temp CH4Temp).

Since the thermostat will have no impact on the initial step, you don’t see a difference there.

Axel.

p.s.: there doesn’t seem to be a need for using dynamic/dof yes.

akohlmey · June 9, 2021, 2:29pm

the key piece of information is in the fix nvt/npt documentation (or any other thermostat) where it states that the fixes will create their own temperature compute and that you can replace it with a compute of your own (e.g. to apply a temperature bias). compute_modify applies to all computes and that includes the ones created by other fixes. those computes just have a name that includes the fix-id. they will show up in the list of computes with “info compute” just the same.

please note that both of my suggestions have the same effect, only the route of where the compute is created is different.

axel.

lammps-users · June 9, 2021, 2:00pm

Dear Axel,

Thanks a lot for your help.

compute_modify fix-ID_temp extra/dof 0 solved my problem, giving the correct degree of freedom and particle kinetics. And I use thermo_style custom ... c_fix-ID_temp ... to print the temperature of the thermostat.

Besides, I didn’t find the compute_modify fix-ID_temp keyword in the LAMMPS manual. I got a bit confused from the manual and thought compute_modify only works for compute-ID, the same for fix_modify and fix-ID.[1][2]

Best,
Henglu

[1] https://docs.lammps.org/compute_modify.html
[2] https://docs.lammps.org/fix_modify.html

lammps-users · June 9, 2021, 8:17pm

Dear Axel,

Thank you very much for the discussion.

In fix_nh.cpp, t_current is used in the velocity scaling and to calculate kecurrent. In my problem, to retrieve the correct particle kinetics, t_current needs to be calculated using 3N degree of freedom, so compute_modify extra/dof needs to be applied on the exact compute-ID_temp of t_current because t_current is up to date.

In the two approaches you suggested, which are using different computes, compute-ID_temp of t_current is eqCH4_temp in the first case and CH4Temp in the second case, these two approaches give equivalent results.


## APPROACH 1
fix eqCH4 CH4 nvt temp ...
#BY DEFAULT# compute eqCH4_temp CH4 temp
compute_modify eqCH4_temp extra/dof 0
## APPROACH 2
fix eqCH4 CH4 nvt temp ...
compute CH4Temp CH4 temp
compute_modify CH4Temp extra/dof 0
fix_modify eqCH4 temp CH4Temp ## `eqCH4_temp` vanishes

Best,
Henglu