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}
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.
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.
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.
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.
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.
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:
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?
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.
