I'm trying to run a simulation with a mobile wall in the x-direction, which should move in response to its total force, composed of (a constant external force) + (the force exerted by all atoms on the wall). This seems possible using the "fix wall/region" command together with a dynamic region definition. However, I cannot seem to find a way of accessing the total force exerted on the wall, which is stored in ewall[1:3] in the fix wall/region routine. Is there a simple way of extracting this variable for use in my input script (through, for example, a compute command), or should I explicitly make this variable available by modifying the code? Thank you!
I'm trying to run a simulation with a mobile wall in the x-direction,
which should move in response to its total force, composed of (a
constant external force) + (the force exerted by all atoms on the wall).
This seems possible using the "fix wall/region" command together with a
dynamic region definition. However, I cannot seem to find a way of
accessing the total force exerted on the wall, which is stored in
ewall[1:3] in the fix wall/region routine. Is there a simple way of
extracting this variable for use in my input script (through, for
example, a compute command), or should I explicitly make this variable
available by modifying the code? Thank you!
fix wall/region stores these data in a global vector that can be accessed,
e.g. with an equal style variable in a way like this.
Thank you for this information, Axel! I can now indeed access the wall force, but it seems this instead leads to a difficult problem when it comes to the order in which I evaluate things… Here’s my specific example, which is supposed to move the wall region in response to the force exerted by the particles on the wall:
variable fwx equal f_w1[1]
region mobilewall block -10.0 0.0 0.0 100.0 -0.25 0.25 side out move v_fwx NULL NULL
fix w1 all wall/region mobilewall lj126 1.0 1.0 1.122462
When I print out v_fwx in the above setting I get zero, while if I remove the mobility of the region (i.e. change v_fwx --> NULL in the region command), I get a sensible value of the wall force. I assume this is because, when I call the variable command (v_fwx) inside the region command, fwx is evaluated before the fix wall/region command is called, and thus fwx is still 0. On the other hand, I cannot put the region definition after the fix wall/region command, as the latter relies on knowledge of the former. So it seems I’m stuck! Do you (or someone else) have a neat idea about how to solve this issue with as little hacking as possible? Thank you!
Thank you for this information, Axel! I can now indeed access the wall force, but it seems this instead leads to a difficult problem when it comes to the order in which I evaluate things… Here’s my specific example, which is supposed to move the wall region in response to the force exerted by the particles on the wall:
variable fwx equal f_w1[1]
region mobilewall block -10.0 0.0 0.0 100.0 -0.25 0.25 side out move v_fwx NULL NULL
fix w1 all wall/region mobilewall lj126 1.0 1.0 1.122462
When I print out v_fwx in the above setting I get zero, while if I remove the mobility of the region (i.e. change v_fwx → NULL in the region command), I get a sensible value of the wall force. I assume this is because, when I call the variable command (v_fwx) inside the region command, fwx is evaluated before the fix wall/region command is called, and thus fwx is still 0. On the other hand, I cannot put the region definition after the fix wall/region command, as the latter relies on knowledge of the former. So it seems I’m stuck! Do you (or someone else) have a neat idea about how to solve this issue with as little hacking as possible? Thank you!
Please check the mailing lists archive for a discussion on storing properties from a previous step. We had this discussion multiple times, but I don’t recall the exact solution presently.
Thank you for this information, Axel! I can now indeed access the wall
force, but it seems this instead leads to a difficult problem when it comes
to the order in which I evaluate things... Here's my specific example,
which is supposed to move the wall region in response to the force exerted
by the particles on the wall:
variable fwx equal f_w1[1]
region mobilewall block -10.0 0.0 0.0 100.0 -0.25 0.25 side
out move v_fwx NULL NULL
fix w1 all wall/region mobilewall lj126 1.0 1.0 1.122462
When I print out v_fwx in the above setting I get zero, while if I remove
the mobility of the region (i.e. change v_fwx --> NULL in the region
command), I get a sensible value of the wall force. I assume this is
because, when I call the variable command (v_fwx) inside the region
command, fwx is evaluated *before* the fix wall/region command is called,
and thus fwx is still 0. On the other hand, I cannot put the region
definition *after* the fix wall/region command, as the latter relies on
knowledge of the former. So it seems I'm stuck! Do you (or someone else)
have a neat idea about how to solve this issue with as little hacking as
possible? Thank you!
nothing like being stuck on the subway without internet to jog your memory:
you can use fix ave/time to get access to the data before it is destroyed.
try this:
variable fwx equal f_w2[1]
region mobilewall block -10.0 0.0 0.0 100.0 -0.25 0.25 side
out move v_fwx NULL NULL
fix w1 all wall/region mobilewall lj126 1.0 1.0 1.122462
fix w2 all ave/time 1 1 1 f_w1[1]
fix ave/time will compute the average over 1 step (i.e. make a copy) of the
data from fix wall/region and store it in a global vector, so that it is
accessible in the next iteration.
Thanks for the suggestion! It seems neat, but still seems to give me exactly the same problem - i.e., that fwx (and w2) will evaluate to zero as soon as I put the variable into the region command. Are you positive that fix ave/time will actually make a difference? This command still relies on f_w1, which in turn relies on the region command, which is where v_fwx is needed, so we’re back at square one…
Thanks a lot for the input! I know that the “move” option was not designed for this particular implementation, but I was hoping that I might still work. After all, the amount of force on the wall is, to me, indeed a function that can give me a displacement as a function of time - albeit a fairly complex one. As you say, I will of course need to define an equation of motion for the wall, including an effective mobility which, together with an assumption of overdamped dynamics, can directly translate the total force on the wall to a displacement (after multiplying with the timestep). It seems to me that this idea is in principle working, apart from the small caveat that the force variable is not fully accessible to me.
Regarding this issue, it seems moving the initialisation line doesn’t make any difference, neither does Axel’s suggestion with fix ave/time do the trick. I was considering if making a “dummy” copy of the wall/lj126 fix (under a different name) which only calculates the force on the wall without updating the particle forces might work. I would then use this fix as input to the region command, which will then be unlinked from the fix wall/region command. Is this a silly idea? It seems that this specific variable-access problem is simple enough that I should be able to solve it without a Python script, but it might turn out that this would anyway be the best solution in the end!
I now realised what you mean - I slightly misunderstood the documentation of the region command. I thought that the “displacement” meant the displacement per timestep rather than the total displacement relative the position at the start of the run. Might still be possible to write this as a variable that’s continuously updated in the input script though, but significantly more tricky! I might have a look at the python interfacing as you suggested to see if that’s an easier option. Thanks for your help and patience!
I now realised what you mean - I slightly misunderstood the documentation
of the region command. I thought that the "displacement" meant the
displacement per *timestep* rather than the total displacement relative the
position at the start of the run. Might still be possible to write this as
a variable that's continuously updated in the input script though, but
significantly more tricky! I might have a look at the python interfacing as
you suggested to see if that's an easier option. Thanks for your help and
patience!
steve's suggestion of using a python style variable for your purpose is
very good and one of the best use cases that i have seen for using python
style variables.
I have had a read through this thread, but I may be barking up the wrong tree.
It seems to me that you want to move the wall in response to the force on that wall. This is easily achieved if you use the velocity aspect of the ‘fix move’ command, rather than the displacement. This is because the velocity does not directly move the wall, so there will be no bogus behaviour at the next step.
I have built a ‘2 term PI controller’ which uses the force between 2 groups, compares that to the force I would like to have between those groups, and changes the velocity of one of the groups (which are otherwise not time integrated, and so act as a fixed body) to bring the force towards what I am asking for. This does not need a python variable.
Thanks for this tip! However, it seems to me that the “fix move” command will apply a prescribed displacement to a group of atoms rather than to a wall, which is what I want to achieve. Or would your solution imply building up a wall from actual particles? Furthermore, it seems that fix/move should not be combined with explicit time integration, which is exactly what I want to do (i.e., to have an equation of motion for the wall). I’m not sure if your routine would help me for the same reason (i.e., that it seems to apply to atoms rather than walls), but I might have misunderstood it! What do you think?
It’s a simple use of a couple of variables, and a modified form of fix ave/time so that it give the sum up to that point. You may want to make the mod – we (I say we; Duncan Harris gets the credit/blame) did it by copying the ave/time and forcing the ‘normalising divider’ to always be 1.
I’ve now tried to implement a basic python script controlling the wall movement, which indeed seems like the way to go. I do have some thoughts on the basic implementation though, which I was hoping that you might be able to give some feedback on. Here is a simplified version of the script, which so far only pushes the wall at a constant speed (the actual integration should be easy to fix once I get a good base solution):
from lammps import lammps
infile = ‘in.dipole’
nstep = 10000
lmp = lammps()
lmp.file(infile)
wallpos = 0.
for istep in range (0,nstep):
wallpos += 0.001
cmd = “fix w1 all wall/lj126 xlo " + str(wallpos) + " 1.0 1.0 1.122462 units box”
lmp.command(cmd)
lmp.command(“run 1 pre no post no”)
This works! The main issue is, however, that the “run 1” inside the loop makes LAMMPS output huge amounts of information (in spite of the “pre no post no” options), which I assume might also affect performance. As I need to update the wall position at every timestep I can’t see immediately to avoid this, neither in this setting, nor by using python-style variables inside the input script instead. The “run every” option seemed like a viable option, but that will not allow me to integrate the wall position in the script, as the command inside the run command will have to be updated on each timestep. Any ideas on how to circumvent this?
I've now tried to implement a basic python script controlling the wall
movement, which indeed seems like the way to go. I do have some thoughts on
the basic implementation though, which I was hoping that you might be able
to give some feedback on. Here is a simplified version of the script, which
so far only pushes the wall at a constant speed (the actual integration
should be easy to fix once I get a good base solution):
from lammps import lammps
infile = 'in.dipole'
nstep = 10000
lmp = lammps()
lmp.file(infile)
wallpos = 0.
for istep in range (0,nstep):
wallpos += 0.001
cmd = "fix w1 all wall/lj126 xlo " + str(wallpos) + " 1.0 1.0 1.122462
units box"
lmp.command(cmd)
lmp.command("run 1 pre no post no")
This works! The main issue is, however, that the "run 1" inside the loop
makes LAMMPS output huge amounts of information (in spite of the "pre no
post no" options), which I assume might also affect performance. As I need
to update the wall position at every timestep I can't see immediately to
avoid this, neither in this setting, nor by using python-style variables
inside the input script instead. The "run every" option seemed like a
viable option, but that will not allow me to integrate the wall position in
the script, as the command inside the run command will have to be updated
on each timestep. Any ideas on how to circumvent this?
you misunderstood steve's suggestion. you are using the LAMMPS library
wrapper and the python module built on top of it, but what steve and i
were suggesting was to use python style variables and the python *package*
that allows to embed python functions into LAMMPS. that bypasses your
problem entirely. so to do what you are doing above, you can use the
following LAMMPS code (which can be easily expanded to do time integration
or a fictitious wall position object):
# declare variables for later use
variable step equal step
variable wallpos python movewall
# code to initizlize and reset python globals
python initwall input 1 -5.0 format f here """
oldstep = 0
pos = 0.0
def initwall(initpos):
global oldstep
global pos
pos = initpos
oldstep = -1
"""
# initialize python data
python initwall invoke
# python code to determine location of the wall
python movewall input 2 0.001 v_step return v_wallpos format fff here """
def movewall(delta,step):
global oldstep
global pos
if (oldstep != int(step)):
pos += delta*(step-oldstep)
oldstep = step
return pos
"
fix w1 all wall/lj126 xlo v_wallpos 1.0 1.0 1.122462 units box
thermo_style custom step v_wallpos
thermo
10
0
run
1000
Thanks for this! I actually didn’t misunderstand Steve’s suggestion - I did give a thorough thought, but didn’t see how it’d solve my initial problem with accessing the wall force from fix wall in the previous timestep. I managed this with the Python wrapper, albeit in a non-elegant way. Nevertheless, I did give your suggestion a try, which reads something like this:
python movewall input 4 v_mobility v_pres f_w1[1] v_timestep return v_wallpos format fffff here “”"
def movewall(mob,pres,force,timestep):
global pos
pos += timestepmob(pres*100.0 + force)
return pos “”"
fix w1 all wall/lj126 xlo v_wallpos 1.0 1.0 1.122462 units box
fix w2 all wall/lj126 xhi EDGE 1.0 1.0 1.122462 units box
This gives me an error message “ERROR: Expected floating point parameter in input script or data file (…/python.cpp:344)”, but when I change f_w1[1] for a numerical value it works fine; thus, I’m back to the problem I tried to solve initially, and I still don’t see a way around this without integrating one step at a time, saving the force to a variable, using this to integrate the wall position and then reiterate. I assume I can do this with a “run 1” construction inside the input script rather than using the wrapper, but it’d still leave me with the problem I emailed you about yesterday. Or am I being stupid now?
Thanks for this! I actually didn't misunderstand Steve's suggestion - I
did give a thorough thought, but didn't see how it'd solve my initial
problem with accessing the wall force from fix wall in the previous
timestep. I managed this with the Python wrapper, albeit in a non-elegant
way. Nevertheless, I did give your suggestion a try, which reads something
like this:
python movewall input 4 v_mobility v_pres f_w1[1] v_timestep return
v_wallpos format fffff here """
def movewall(mob,pres,force,timestep):
global pos
pos += timestep*mob*(pres*100.0 + force)
return pos """
fix w1 all wall/lj126 xlo v_wallpos 1.0 1.0 1.122462 units box
fix w2 all wall/lj126 xhi EDGE 1.0 1.0 1.122462 units box
This gives me an error message "ERROR: Expected floating point parameter
in input script or data file (../python.cpp:344)", but when I change
f_w1[1] for a numerical value it works fine; thus, I'm back to the problem
I tried to solve initially, and I still don't see a way around this without
integrating one step at a time, saving the force to a variable, using this
to integrate the wall position and then reiterate. I assume I can do this
with a "run 1" construction inside the input script rather than using the
wrapper, but it'd still leave me with the problem I emailed you about
yesterday. Or am I being stupid now?
python input (currently) only accepts numbers or variables. if you want to
access data from a fix to have to access it via an equal style variable.