Hello everyone,
I’m trying to simulate flow of water molecules through a slit between two sheets of an organic substance. I place a reservoir of water molecules on both ends of the slit and at the extreme end I have sheets of graphite on both sides. I want to essentially apply a pressure differential on the graphite slabs and the impact of pressure differential on the flow of water through the slit. I treat the graphite slabs as rigid molecules and use addforce on it to create a pressure differential. The problem I’m facing is that even with addforce there seems to be no change in the position of the slabs. The slabs remain stationary and do not move and there seems to be no force acting on the water molecule. Could you please let me know what I am doing wrong?
Another question I have is when I want to treat the graphite slabs as just walls, should I be describing a force field for graphite? (If I don’t that the slabs will pass right through the water molecules…)
graphite 1 is the slab on the left and graphite 2 is slab on the right
The wall is essentially confining the molecules to a region so that water molecules don’t escape.
reset_timestep 0
thermo 1000
timestep 0.001
fix wall all wall/region graphitezone lj126 1.0 1.0 2.5
fix 4 graphite12 rigid group 1 graphite12
fix 14 graphite1 addforce 50.0 0.0 0.0
fix 15 graphite2 addforce -50.0 0.0 0.0
fix move water12 nve
run 200000
using fix rigid + fix addforce may not be the most elegant way of implementing the behavior you are looking for (it does not allow for good parallel scaling).
you should be able to get the intended behavior by using fix aveforce instead and including those atoms into the group to which fix nve is applied.
I also suggest you first test this without the water.
axel.
Thank you, Axel! That was helpful.
I used aveforce and nve on the rigid group. I tried testing this with just the graphite slabs as you suggested. When I apply the aveforce in the x direction, the rigid slabs are not moving in the x direction, but at an angle to the z axis. So I tried making the Y and Z force components in aveforce NULL, but even now the slabs seem to rotate instead? Would you have an idea what could be going on?
Also I used the fix rigid/nve, using fix nve on the slabs seems to cause lost atoms issue.
timestep 0.011
fix 4 graphite12 rigid/nve group 1 graphite12
fix 5 graphite1 aveforce 10.0 NULL NULL
fix 6 graphite2 aveforce -10.0 NULL NULL
Thank you
Thank you, Axel! That was helpful.
I used aveforce and nve on the rigid group. I tried testing this with just the graphite slabs as you suggested. When I apply the aveforce in the x direction, the rigid slabs are not moving in the x direction, but at an angle to the z axis. So I tried making the Y and Z force components in aveforce NULL, but even now the slabs seem to rotate instead? Would you have an idea what could be going on?
you probably have initialized the atoms in the slabs with some velocities other than 0.0
please remind me, why do you want to apply a force and not just have those slabs move at a prescribed fixed velocity?
Also I used the fix rigid/nve, using fix nve on the slabs seems to cause lost atoms issue.
you must not apply two time integrating fixes to the same atoms. LAMMPS will print a warning about this, which you should take seriously.
axel.
Thank you for the suggestion, Axel.
I tried moving the slabs with constant velocity as you suggested and it works pretty well. However, what I want is the pressure differential applied between the slabs. That’s why I was aiming to control the forces applied on the slabs (which still doesn’t move the slabs linearly, and I haven’t figured out why…).
With the velocity method, the issue I’m facing is with the calculation of pressure. I used the compute stress/atom and reduce it to get the total stress. But since this requires calculation of volume of the atoms as well, I tried the voronoi/atom and reduce that to get the total volume of the atoms in the slab. This gives me an insanely high pressure for very low velocities which I’m questioning…
I suppose I will try computing the voronoi volume and stress on smaller systems of which I know the property, and debug it.
Thanks again for all your help!