# How to thermostat a piston wall with given external pressure

Dear all,

I’ve got a problem. I am doing liquid bubble nucleation on solid wall and I want to keep my system at constant pressure. Therefore, I used a solid wall at the top working as a piston and applied constant force using “fix aveforce” to mimic constant external pressure which is shown in the figure. And to make the piston atoms not blown up, I used the “fix rigid” command to make the piston wall move as a whole.(I am not sure if it is correct to use “fix rigid”, but without this command the piston wall will explode.)

I want to maintain the whole system temperature at 300K during the equilibration period so I used a langevin thermostat for both the solid and liquid atoms as shown below. Besides, to keep the solid atoms under thermal vibration, I used the “fix spring” to let the solid atoms vibrating around their initial positions. However, I monitored the temperature of the piston wall during equilibration. The temperature was not 300K but a pretty large value. It doesn’t make sense and I don’t know where is wrong.

The commands that I used and the temperature output information were pasted below as well. As you can see, the thermostat worked well for liquid and the bottom solid wall atoms but the piston temperature rise up to more than 10000K. It would be much appreaciated if anyone could point out my problem and show me the solution.

Best regards,
Huaqiang

fix spring solid spring/self \${springcons}
fix pistonrigid piston rigid single

variable pressure equal 1.01325 # atmosphere pressure, 1.01325 bar
variable force equal -{pressure}*lx*ly*6.2415e-7 #1 eV/A3 = 1.60218e6 bar, 1 bar = 6.2415e-7 eV/A3 fix conspress piston aveforce NULL NULL {force}

fix fNVE gthermo nve
fix fNVT gthermo langevin 300 300 \${Tdamp} 54687 tally yes

Step PotEng c_liquidtemp c_ltempcom c_bottomsolidtemp c_pistontemp
0 -6072.8409 300.67666 300.60388 296.35187 314.35183
1000 -5950.8854 306.04658 306.03024 298.59288 11979.62
2000 -5955.0416 299.24039 299.17544 305.31793 11687.528
3000 -5960.9367 298.44887 298.39681 302.05843 11613.145
4000 -5966.0016 298.28059 298.22337 307.98911 11766.455
5000 -5963.1173 299.58155 299.57656 300.30442 11690.392
6000 -5961.362 295.74873 295.68325 304.9265 11692.116
7000 -5964.2282 305.99951 305.90512 303.0286 11884.135
8000 -5964.8743 305.23263 305.11932 297.68174 11579.761
9000 -5972.4619 299.64192 299.51238 303.11088 11585.04
10000 -5969.8647 297.36952 297.30997 306.85762 11481.543
11000 -5966.6676 294.67152 294.65507 306.8742 11589.772
12000 -5970.8856 301.88411 301.75024 299.63099 12114.777
13000 -5971.2732 303.98168 303.92532 307.06429 11689.798
14000 -5975.3776 299.39211 299.35435 303.44743 11871.216
15000 -5970.0912 296.15279 296.14327 305.45128 12044.773
16000 -5970.2396 304.1463 303.94733 302.33427 11676.827
17000 -5972.2242 303.73806 303.72924 298.91539 11611.167
18000 -5974.9747 298.75664 298.75383 299.46743 12047.719
19000 -5979.7322 300.94428 300.93173 298.99499 12202.711
20000 -5977.3123 300.64568 300.63493 297.20767 11834.675
21000 -5980.0247 297.63133 297.58785 289.16902 11480.162
22000 -5979.0115 304.83932 304.83261 302.70581 11815.479
23000 -5983.2138 303.22544 303.18585 298.15083 11645.94
24000 -5975.351 301.54933 301.49702 307.79553 11692.013
25000 -5977.7105 299.99199 299.9294 298.85417 12082.813
26000 -5977.6303 299.29086 299.25458 301.28117 11631.49
27000 -5978.8198 301.38791 301.34431 303.52968 12140.948
28000 -5984.5888 304.80038 304.75096 285.31906 12342.003
29000 -5983.2467 300.56778 300.5384 304.09003 11900.625
30000 -5980.0742 299.52052 299.49614 297.81838 12122.636

The information you provide is not sufficient to make any specific recommendations.
If your setup would be correct, then fix aveforce should be sufficient to maintain the structure of the piston wall.
please see example/flow/in.flow.pois for a 2D example. It should work in 3d if adjusted properly and fix rigid would not be needed.
Since applying fix rigid will result in a system with very limited degrees of freedom, you cannot trust the temperature compute for the atoms of the piston. in fact, there would be only one degree of freedom left, if done properly, but LAMMPS cannot easily determine this internally.

Axel.

There is a simple example of this in examples/flow/in.flow.pois using fix aveforce.
You can visualize what it does by uncommenting the dump image commands.

I wouldn’t make the top red atoms rigid. Fix aveforce effectively does that for you,
b/c it adds a variable force to each red atom so that all experience the same force.

Steve

Dear Steve and Axel,

Thanky you so much for your advices. I have checked the example script carefully and tried different methods which are described below:

• If only fix aveforce is used without thermostat which is similar to the example script, it works fine with the following command. In this way, the piston was like “frozen” atoms with displacement only in z direction.

fix conspress piston aveforce NULL NULL \${force}
fix pistonkeep piston setforce 0.0 0.0 NULL
fix fNVE piston nve

• If the piston temperature is controlled using lagevin thermostat using the following command, it reports error “Domain too large for neighbor bins”. I think it is because the system collapsed.

fix conspress piston aveforce NULL NULL \${force}

fix pistonkeep piston setforce 0.0 0.0 NULL

fix fNVE piston nve

fix fNVTlan piston langevin {Ts} {Ts} \${Tdamp} 54687 tally yes

• Then, I delete the constrain in x and y direction (the command " fix pistonkeep piston setforce 0.0 0.0 NULL "). The piston is still in langevin thermostat. The piston soon exploded as shown in the below image.

• If Nose-Hoover thermostat (i.e. fix nvt) is used instead of langevin thermostat, the piston slowly converges to around 200K. The reason I think is that I constrain the piston atoms to xy plane using “fix spring/self xy” command so that the degree of freedom in z direction is gone. That seems fine as long as the correct degrees of freedom is used to calculate the temperature.

fix springpis piston spring/self 2.9332 xy
fix conspress piston aveforce NULL NULL {force} fix fNVTlan piston nvt temp {Ts} {Ts} {Tdamp}

In summary. it seems that my option is either to control the piston as “frozen” atoms with no thermostat, or using Nose-Hoover thermostat and calculate the temperature myself. But I still have some questions:

1. Why Nose-Hoover thermostat works while langevin thermostat doesn’t? Is it because the langevin thermostat includes a random force that makes the stoms couldn’t move as a whole?
2. I calculate the forces applied on the piston atoms using the methods below but I don’t know if the force should be divided by the atom number since it is average force on each atom (I saw someone did in literature). I have read the manual. I think the divion marked in red is not needed but I am not sure.

what seems to be missing is that you are also setting the z velocity of the atoms in the piston to zero. please find below a simple demonstration based on the melt example for a setup where both bottom and top layer are force to keep the overall crystal structure with fix spring/self, but are allow some movement and then it uses fix aveforce in z direction to apply a pressure on the liquid part.

axel.

units lj
atom_style atomic
boundary p p m

lattice fcc 0.8442
region box block 0 10 0 10 0 10
create_box 3 box
create_atoms 1 box
mass * 1.0

region bottom block -1 11 -1 11 -0.1 1.5
region top block -1 11 -1 11 8.5 11
group bottom region bottom
group top region top

set group bottom type 2
set group top type 3

velocity all create 3.0 87287

pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0 2.5

neighbor 0.3 bin
neigh_modify every 2 delay 2

fix 1 all nvt temp 1.5 1.5 0.1
#fix 1 all nve
#fix 0 all langevin 1.5 1.5 0.1 67123456

velocity top set NULL NULL 0.0

fix 2 bottom spring/self 50.0 xyz
fix 3 top spring/self 50.0 xy
fix 4 top aveforce NULL NULL -3.0

dump 3 all movie 10 movie.mpg type type &
axes yes 0.8 0.02 view 60 -30 size 1024 1024