I have a small system which is equilibrated at 5 GPa pressure and 1K, and now I make it a hybrid model with a mixture of rigid bodies and non-rigid particles. To do so, I define the top part (rigid1) and the bottom part (rigid2) of the model as rigid bodies, and the middle part as non-rigid (mobile), like a sandwich. Then, the following fixes are applied:
fix 1 mobile npt temp 1 1 100 x $xval xval 2000 y {ypressVal} {ypressVal} 2000 z {zpressVal} ${zpressVal} 2000 xy 0.0 0.0 2000 xz 0.0 0.0 2000 yz 0.0 0.0 2000 dilate all
fix 2 boundary rigid/nve group 2 rigid1 rigid2
fix 3 rigid1 setforce 0.0 0.0 NULL
fix 4 rigid2 setforce 0.0 0.0 NULL
run 10000
Since the system is well equilibrated before, I do not expect the system to oscillate too much. However, the model separates into several parts during the simulation (as you can see in the video), the pressure and temperature vary a lot as well.
Another question is, although I set the x and y force of the rigid bodies as 0, they still can move along these two directions.
The version of LAMMPS is 29Oct 2020. The log file and the video are attached below.
If you could provide any ideas, that’s appreciated.
After going though more posts in the mailing list, I tried following combinations to make sure the problem doesn’t come from my system:
fix 1 mobile nve
fix 2 boundary rigid/nve group 2 rigid1 rigid2
fix 3 rigid1 setforce 0.0 0.0 NULL
fix 4 rigid2 setforce 0.0 0.0 NULL
run 10000
and
fix 1 mobile nvt temp 1 1 1000
fix 2 boundary rigid/nve group 2 rigid1 rigid2
fix 3 rigid1 setforce 0.0 0.0 NULL
fix 4 rigid2 setforce 0.0 0.0 NULL
run 10000
either combination works.
However I want to control the pressure of the system, and once I use ‘fix npt’ for the mobile part, the temperature becomes so high that the simulation stops. Since my system is triclinic, I cannot use ‘fix rigid/npt’ to control the pressure.
Note, someone mentioned ‘fix setforce’ doesn’t work with ‘fix rigid’ in the mailing list, it seems not a problem for me.
There is a lot of artifice in your model so it is somewhat questionable if all these complexities are needed to get results that are significantly improved over just using two wall fixes.
Using fix rigid/npt, even if it would work for triclinic cells, would still be a major problem, since you are doing a variable cell simulation with two(!!) particles (the two rigid bodies). Try setting up a system with two argon atoms and do fix npt on it and see if you can achieve proper pressure control. You’ll very likely run into trouble with that, too.
Furthermore, you are treating your system as periodic, but conceptually, you have a non-periodic setup, you want to use the two rigid bodies as pistons that are adjusted to achieve a target pressure in the system in between them.
Since you have equilibrated your system without keeping those boundaries separate, there is likely a possibility that molecules have migrated and mixed and thus are now unphysically interlinked (difficult to tell from remote).
Next issue is that by keeping your boundary objects rigid you are now suddenly freezing them to 0K and thus also creating some unphysical situation.
If you consider this an approximate approximation, you probably should look into making your system a) non-periodic in the direction of the boundaries and b) try using fix aveforce (after setting velocities to zero of those atoms) as that would be the way to implement a piston-like setup with a constant force/pressure applied.
close to 0 as well. Thus the sudden freezing should be ok.
you are missing my point.
Following your suggestion and some posts in the mailing list, I make my system non-periodic and use ‘fix aveforce’ to apply the pressure on the rigid body with following fixes:
velocity boundary set 0.0 0.0 0.0
fix 2 rigid1 setforce 0.0 0.0 NULL
fix 3 rigid2 setforce 0.0 0.0 0.0
fix 8 rigid1 aveforce 0.0 0.0 v_zpressVal
fix 5 mobile npt temp 1 1 1000 x {xpressVal} {xpressVal} 5000 y {ypressVal} {ypressVal} 5000 xy 0.0 0.0 5000 xz 0.0 0.0 5000 yz 0.0 0.0 5000
fix 4 boundary nve
However, the problem is still there, the simulation can go well once I apply a non-zero fz using ‘fix aveforce’. I’ve attached the input.script, log.lammps and the .dat file.
you have not paid attention to what I suggested. this is wrong. please re-read my previous point.
when you keep your boundary immobile then you must not use fix npt. in fact, I have yet to see a good argument for you why what you are trying to do is worth doing and what you can learn from that that is meaningful and not tainted by your bad choices.
I have one more question about ‘fix aveforce’. Since I have many different types of atoms with different masses, is it reasonable to use ‘fix aveforce’ to apply the same force to each atom? Can ‘fix aveforce’ take care of this automatically?
check the documentation or do the math/physics allows for it, or simply set up a small test system and observe.
Sorry, I wasn’t capable of capturing the highlights in your previous email and misunderstood them.
What I’m trying to do is to calculate the energy of stacking faults in my system at a specific pressure and the current model represents a perfect model without stacking faults.
According to your statement, once there are pistons in a system, barostat is not allowed. This is in agreement with my tests and I appreciate that you pointed it out to me again, patiently. Then I have a question about how a piston is defined in LAMMPS. If a boundary is immobile along one or more directions, that is a piston, am I right? What if I use ‘fix addforce’ to move a boundary which may or may not be rigid, is that boundary a piston and is NPT still not allowed for the rest of the system? If I understand this correctly, ‘fix right/npt’ works for (small) rigid bodies which are not part of the boundary of the system.
Is it necessary to keep a system with pistons non-periodic? I saw that the example system in example/in.flow.pois is periodic with pistons.
And now I figured out how to use ‘fix aveforce’ for molecules.
Is it necessary to keep a system with pistons non-periodic? I saw that the example system in example/in.flow.pois is periodic with pistons.
No, it is not. It is periodic in the direction to the flow, but it uses shrinkwrap orthogonal to it, i.e. in the direction of the added force.
You really have to train yourself to pay better attention to detail.
I should say once the boundary is a piston in a system, barostat is not allowed.
it is not so much a question of what is “allowed” but rather a matter of what makes sense.
once you have a rigid boundary object, it doesn’t make sense to apply a variable cell integrator, as it would have to do unphysical transformations.
in a corresponding macroscopic experiment the rigid boundary would block any deformation orthogonal to it.
Technically speaking, In the pressure direction you may have periodicity, but you need to have vacuum to decouple the upper boundary object from the lower boundary object, and then it makes more sense to not use a periodic boundary at all. with a periodic boundary, it would only make sense to have a single rigid boundary object.