Maintaing temperature of a group of particle using fix temp/rescale

Dear all,
I am simulating in isenthalpic ensemble and trying to maintain the temperature of a group of particles (part of a rigid molecule) to some fixed T. I first used temp/rescale to rescale velocity to T, for some run lengths. I know temp/rescale will ramp the temperature of the group from T_start to T. I am using variable for T_start. The T_group goes nicely to T, but now I want to maintain that T for that group using temp/rescale only, as I am using fix nph. After the ramping, I again ran the simulation with same fixes thinking that now temperature is T = T_start and T_last =T, so the temperature should be maintained at T. But, for my surprise I am not able to maintain the temperature to T, but again T_start took a lower value and was starting the ramping to T.
So, in summary I am asking how to maintain the T_group using fix temp/rescale.

From your message, I am guessing that your code looks like that :

fix 3 flow temp/rescale 100 ${Tstart} ${T} 0.02 0.5
run 1000
run 1000

But what you want is :

fix 3 flow temp/rescale 100 ${Tstart} ${T} 0.02 0.5 # perform a temperature ramp from Tstart to T
run 1000
unfix 3 # optional
fix 3 flow temp/rescale 100 ${T} ${T} 0.02 0.5 # maintain the temperature at T
run 1000


Thanks, Simon, but I did not know that I could take T_stop as a variable too. I was trying with
fix 3 flow temp/rescale 100 2.0 2.0 0.02 0.5
but in this case, I was getting exactly 2.0 for the complete run, with no fluctuations.

No fluctuations in temperature at all is probably because you have no fix updating the positions and velocities of your atoms (like fix nve).

I am using fix nph.
For a region near the edges, the rescaling is running fine with temperature fluctuations. The problem is with a hemispherical shell of a Janus particle. I am trying to fix the temperature for this hemisphere at T=2. I started with uniform T for my complete system, then I rescaled box edges to 0.75 and the hemisphere to 2. The ramping was all fine. The commands are pasted here:

**variable Tedges equal temp**
**variable Themishell equal temp**
**fix 3 g_edges temp/rescale 100 ${Tedges} 0.75 0.02 1.0**
**fix 4 g_hemishell temp/rescale 100 ${Themishell} 2.0 0.02 1.0**
**thermo_style custom step temp c_3_temp c_4_temp epair emol press vol**
**fix nph1 all nph iso 0.0 0.01 $(1000*dt)**
**run 10000**

The issue is when I order to keep the T_hemisphere to be about 2. And that I am doing as pasted below:

fix 3 g_edges temp/rescale 100 ${Tedges} ${Tedges} 0.02 1.0
fix 4 g_hemishell temp/rescale 100 ${Themishell} ${Themishell} 0.02 1.0
thermo_style custom step temp c_3_temp c_4_temp epair emol press vol
fix nph1 all nph iso 0.0 0.01 $(1000*dt)
run 10000

And the result that I get after the last run is pasted below:

Step Temp c_3_temp **c_4_temp** E_pair E_mol Press Volume 
   15004   0.76908167   0.75536175    **1.6751742**   -5.0305912   0.74976659  0.011457015    146607.96 
   15100   0.76516441   0.75439046   **0.76908167**   -5.0267566   0.74870991  0.015455284    146649.03 
   15200   0.76280784   0.75298959   0.76908167   -5.0232352   0.74670651 0.0041983502    146765.09 
   15300   0.76244456   0.75581584   0.76908167   -5.0228891   0.74563291 -0.015837387    146898.55 
   15400    0.7698248   0.76908167   0.76908167    -5.018928   0.74431351 -0.012700605    146979.82 
   15500   0.76749191   0.75757706   0.76908167   -5.0156567   0.74345156 -0.0080016577    147039.49

The problem seems the value that is being calculated for the hemisphere. After ramping it to around 2, the next calculated value by the last fix 4 is not around 2 but is 0.76, and that too remains constant throughout.

The atoms of the Janus particle are bounded by fene potential.

Helping you is very difficult.

First you say that you have no fluctuation in temperature, but you obviously have some.

Then you copy a small part of your code with too many variable ${} (which makes it very difficult to read) and not enough context.

Hi @syed

This set of commands if weird compared to what you say you want to achieve:

You are taking the current temperature of the whole system (the thermo keyword temp) in your variables and use it in both your fix commands. I don’t think this is what you want. If your system ends the first simulation with a global temperature of roughly 0.75, this is the temperature that you will use in both your following commands.

To be clear, I think that putting these lines:

is not the good way to proceed. (Plus, as @simongravelle said, these info are not enough to understand clearly what is happening so I am only making a guess here. The log output with variables evaluation would be more useful.)

I suggest you take a closer look at your log.lammps output file and see which values are passed to your fix commands to understand what is wrong compared to what you expected. As far as I can tell, LAMMPS does what you told it to do.

Okay sorry for the confusion. I will try to put my query clearly.
So, I have a Janus nanoparticle, and I am heating the hemisphere of the Janus by rescaling its T_hemisphere.
I ramped the T_hemisphere to 2 using fix temp/rescale while updating the positions, velocities, etc via fix nph.
After ramping, I make a second run, where I need the T_hemisphere to be fixed at 2. I again used fix temp/rescale with equal style variable temp, to maintain T_hemisohere at 2. But, the T_hemisphere quickly drops to the surrounding temperature of 0.76. And, then ramps again to 2.
But, If I use fix temp/rescale 2.0 2.0, then the T_hemishpgere remains constant (no fluctuation in T_hemisphere) for the whole run.

Dear @Germain
You are correct, indeed both the fixes are taking global temperature. It occurred to me that using group ID in equal style variable will make fix temp/rescale calculate group temperature. So, what to use instead of equal style variable to make fix 3 and 4 to compute group temperature and not the global one?

No you are not. You are just using fix nph, which is not the same. The fix selects the integrator algorithm. An ensemble is the result of all fixes that manipulate a simulation.

This is a very problematic project. Rigid body fixes maintain internal data structures for center of mass velocity, and rotational velocity and atomic velocities are extracted from them after updating them. So using a thermostat collides with that.

Fix temp/rescale is a very bad choice for any production calculation. It does not sample any meaningful statistical mechanical ensemble, and has many times been proven to yield unphysical results since it gives your system a systematic “kick” when it is applied.

For the reason described above, it may look as if it would manage your temperatures, but it will just change velocities after the fact because the rescale fix will be always applied after the rigid body time integration, yet the rigid fix will ignore those updated velocities in the next rigid body velocity update.

In general, what you are trying to do is a flawed model. The temperature of a part of a rigid object is an ill-defined property. What is it that you want to achieve by this? What is the physics that motivates your model?

As others have commented. Unless you provide a complete input deck (best of a smaller, simplified system that more clearly demonstrates its purpose and does not have to be meaningful for production), everybody has to speculate about the missing parts and thus every advice given is tainted by the guesses people have to make. Good advice requires well laid out problems and questions. None of that has happened in this discussion.

1 Like

Dear @akohlmey , I am trying to study the dynamics of a self-propelled Janus particle in a fluid. This work is motivated by the study done in " Dynamics of Self-Propelled Janus Particles in Viscoelastic Fluids". I have also read papers that have computationally(MD) investigated such a system.

So in an attempt to achieve diffusion through thermophoresis, I start by building a nanoparticle and filling the simulation box with the simple fluid. I used FENE to bind atoms in the Janus nano-particle. Then
equilibrated temperature and pressure of the system using fix npt. I did not treat the Janus separately as a rigid body.
Then to mimic the heating of the hemisphere of the Janus, I was trying rescaling the velocities of that group to 2. And the edges of the simulation box have to be maintained at 0.75.
I am pasting here the minimal script that I am using here to achieve such a system.

# Box and units  (use LJ units and periodic boundaries)

units lj                 # use lennard-jones (i.e. dimensionless) units
atom_style bond        #  atom bounded by bonds 
bond_style fene

boundary p p p           # all boundaries are periodic

# Pair interactions require lists of neighbors to be calculated
neighbor 1.9 bin
neigh_modify every 1 delay 1 check yes
comm_modify mode single cutoff 5.0 vel yes

#Generating initial configuration

lattice fcc 1.465 spacing 0.72 0.72 0.72 
region bx1 block 0.0 51.0 0.0 51.0 0.0 51.0
region rg1 sphere 25.5 25.5 25.5 5 side in
region rg11 sphere 25.5 25.5 25.5 6 side in
region rg2 sphere 25.5 25.5 25.5 3.5 side out
region bx2 block 0.0 51.0 0.0 51.0 25.5 51.0
region rg3 intersect 3 rg1 rg2 bx2
region bxrg1 intersect 2 bx1 rg11 side out
region bx3 block 5.0 46.0 5.0 46.0 5.0 46.0 side out
create_box 4 bx1 bond/types 1 extra/bond/per/atom 15 
create_atoms 1 region rg1 
lattice sc 0.8
create_atoms 1 region bxrg1 basis 1 1 

#/ grouping
group g_sphere region rg1
group g_hemishell region rg3
group g_argon region bxrg1
group g_edges region bx3

set group g_sphere type 2
set group g_hemishell type 3

mass * 1.0 

# Set up interaction potentials
pair_style lj/cut 2.5
pair_coeff      * * 1.0 1.0 2.5

#/create bonds
special_bonds fene
bond_coeff 1 30.0 1.5 1.0 1.0
create_bonds many g_sphere g_sphere 1 0.8 1.2

###########do energy minimization for a better starting point #####################
min_style       sd
minimize        1.0e-12 0.5e4 1000000 1000000
reset_timestep 0

#/setting NPT
timestep	0.01
fix npt1 all npt temp 0.75 0.75 $(100*dt) iso 0.0 0.01 $(1000*dt)

# set timestep of integrator
timestep 0.01 

#### Run the simulation for some time steps 
run 5000
unfix npt1

#/ramping T for edge and hemisphere 
variable Tedges equal temp
variable Themishell equal temp
fix 3 g_edges temp/rescale 100 ${Tedges} 0.75 0.02 1.0
fix 4 g_hemishell temp/rescale 100 ${Themishell} 2.0 0.02 1.0
thermo_style custom step temp c_3_temp c_4_temp epair emol press vol 
fix nph1 all nph iso 0.0 0.01 $(1000*dt)
run 1000
unfix 3
unfix 4
unfix nph1

#/maintaining T for edge and hemisphere 
variable Tedges equal temp
variable Themishell equal temp
fix 3 g_edges temp/rescale 1 ${Tedges} ${Tedges} 0.02 1.0
fix 4 g_hemishell temp/rescale 1 ${Themishell} ${Themishell} 0.02 1.0
thermo_style custom step temp c_3_temp c_4_temp epair emol press vol
fix nph1 all nph iso 0.0 0.01 $(1000*dt)
run 10000

I do not want to control the temperature of the rest of the box. So I want my system to evolve under these conditions while maintaining the pressure. So I use fix nph as an integrator for all of the atoms in the box. **

It is to say that Janus is on its own.**

@akohlmey @simongravelle

What is the physics behind this? I know about particles moving differently due to a external temperature gradient, but how would one physically (i.e. experimentatlly) implement different temperature areas inside a particle, considering that the concept of temperature as such is not so well defined on the nanoscale area in the first place.

To my knowledge this is usually achieved by adding a force. The easiest way to do this in LAMMPS is by just using extended individual particles and then fix propel/self. But for nanoparticles built from explicit atoms, this is doable as well: either pick two reference atoms or split the sphere in half and give each half of the atoms a different atom type, then the orientation is defined by the vector of the two center of masses. That can be normalized and multiplied by the desired force. It should be possible to do that with LAMMPS input file scripting, but it would be complex. I would consider implementing a custom time integration fix based on fix nve that would then also determine the orientation of each particle (grouped by molecule ID, a chunk of code that does that can be found in fix rigid) and then apply the force (or velocity if so desired).

It is not clear at this point why fix temp/rescale would be a good choice for this, since it is a completely unphysical way to adjust kinetic energy of a group of particles. It can be useful during system preparation and equilibration, but most MD codes only have it implemented because it is the easiest temperature control method to implement and thus people tend to do it first and then just don’t remove it (even though they should since it sets such a bad example). If you need to move kinetic energy in and out of a group of particles, you should look for a dissipative thermostat like one of the multiple langevin style thermostats.

This is in contrast to what you wrote in your original post:

which underlines the importance of providing complete and consistent information and best full examples to demonstrate your point.

This also makes very little sense. The self-propelling will add energy to your system and thus heat it up. Unless there is a way to remove that excess energy, you will have a non-equilibrium system and trying to maintain a constant pressure for that would seem strange (to put it mildly).

1 Like

Here are a few lines of LAMMPS input code that turn your spherical particle into a self-propelled one based on my outline, if added to the input example your provided (and after ripping out all the temp/rescale stuff):

# get particle position and direction
variable jpx equal xcm(g_sphere,x)
variable jpy equal xcm(g_sphere,y)
variable jpz equal xcm(g_sphere,z)
variable jdx equal v_jpx-xcm(g_hemishell,x)
variable jdy equal v_jpy-xcm(g_hemishell,y)
variable jdz equal v_jpz-xcm(g_hemishell,z)

# compute per-atom force in direction
variable jpf equal 0.001*count(g_sphere)/sqrt(v_jdx*v_jdx+v_jdy*v_jdy+v_jdz*v_jdz)
variable jfx equal v_jpf*v_jdx
variable jfy equal v_jpf*v_jdy
variable jfz equal v_jpf*v_jdz

# apply force to particle
fix propel g_sphere addforce v_jfx v_jfy v_jfz

There is no such thing as a “box edge” in periodic boundaries. You only see those due to visualization, but since the system is translation invariant, what is “at the edge” can just as well be considered to be “in the middle”.
More importantly, since you define those edge atoms with a (static) group, they will not remain in their “edge” region, but start diffusing around and mixing with the atoms that are not coupled to a thermostat.

I’ll answer to the question you re-asked to me first but cannot help to add something on @akohlmey insights afterwards:

What you are looking for here is the use of compute temp command combined with fix_modify. I suggest you read the documentation of both commands, as well as the documentation of thermalisation fixes (such as nvt). However, please also consider the following.

Thermophoresis is an effect due to temperature and mass gradients relationship. As @akohlmey correctly pointed out, the gradients are usually external to your particles (using HEX algorithms or some variations). This is not what you are trying to achieve since you are trying to make an active swimmer which is different (a light gradient induces a local heating around only half of the particle). I think I understand what you are trying to do but I do not see the point of some of your choices. The physics of rigid objects is indeed quite specific and I do not see how using one relates to your problem. Also since the heating is dependent on a light source (and its direction) this problem seems already tricky to translate to a well-defined MD simulation of a particles with several atoms and its surrounding fluid using only integrators and velocity scaling. And finally, as @akohlmey pointed out, the surrounding real fluid also acts as a heat bath, so the use of the nph integrator seems also weird to me.

But this scopes out of LAMMPS and I might not be of very good help with that. Good luck.

1 Like

I was aware of the fact that the group (by the edge of the box I meant the region near the box boundaries ) that I defined at the start will not be the same throughout the simulation. But this issue I would have tackled separately.

I am trying for self-phoresis of the double-faced Janus particle. I am not trying to create a temperature gradient inside the nano-particle but using the strategy to atomistically realize the double-faced structure of the experimentally employed Janus particles by
giving the atoms on the two hemispheres different thermal resistances to the solvent. The different thermal resistance will be employed by parameterizing the two faces differently. For the current script, I am not modeling the two hemispheres differently. But, I am just trying to make the face of the hemisphere heated.

The heating of the face is done using a laser. Once I get an equilibrated phase, I push my Janus+lj fluid system into a heated phase by switching off the global thermostat and rescaling the face of one of Janus’s hemispheres and maintaining the boundaries at T=0.75. The situation is completely far from equilibrium. Then I leave my system to evolve under NPH (constantly scaling the hemisphere face to T=2 and boundaries to T=0.75) and once the steady-state will reach, I will collect the trajectories for post-processing.

By heating the face of the hemisphere, I just induce a temperature gradient in the solvent around the Janus particle. The face of the heated hemisphere is defined as the outermost shell of width 1.5 sigma.

This again spells a different story than what one could infer from what you wrote before due to lack of detail. What you describe as experiment, I can follow, but how you derive your model from it is a mystery to me. It makes even less sense now than anything you wrote before.
You talk about reaching an equilibrium but then acknowledge that you are far from an equilibrium at the same time. You say that you want to let your system evolve, but then talk about adding a barostat and keep rescaling velocities.
At any rate, I am tired of trying to discuss science when the details of the story keep changing in every iteration. As @Germain stated, this is getting off-topic anyway.

Summarizing your inquiries and misconceptions with respect to thermostats only:

  • if you apply a thermostat to a group of atoms, it will thermalize this group of atoms only. the fix will create its own compute temp instance for that purpose and you can access that. details are in the LAMMPS manual.
  • if you define a temperature ramp, it will be repeated for every “run” command. this is also documented behavior.
  • if you provide arguments as ${name} variables, they are replaced when pre-processing the input and thus they are not recognized by any fix as a variable. so the corresponding section in the manual must be disregarded unless you pass the variable as a variable reference with v_name (when and where supported by the fix as stated by the documentation.).
  • fix temp/rescale is a bad choice for any application in a production simulation

Well, It just happened. It never came to me that I had to write down my whole project. But in the flow of the discussion, I wrote my actual plan.

I may not be very clear about this in the initial messages. I say that I start with a starting point where I use NPT for the initial run. After this equilibration, I start my heating phase (not adding any heat), where both the velocities in the boundaries of the box and the face of the Janus hemisphere are rescaled to match a prescribed Temperature. This way I generate a temperature gradient around the Janus particle. My story will also include the scenarios where the faces of the Janus have different thermal capacities or have different wetting(in case I take binary fluid mixture). These scenarios will have their own physics to investigate.
This heating phase is my only concern and here the global thermostat is off, maintaining the face of Janus at T=2 and the boundaries at T=0.75. For investigating the science part I wait till the heating phase reaches a steady-state (sorry for the “evolve” term)

I beg your pardon for such a statement, but I am trying to recreate a situation of a heated Janus particle in a liquid and I run my simulation under constant pressure with no enthalpic changes. And this is a non-equilibrium MD study.