Storing old atomic positions

Hello All,

I am writing a compute that requires the instantaneous displacement of each atom during the run as a part of its operation. That is, I want to subtract the position of the atom at the previous time step from its current position after integration. For that, I need to add an array to the atom class to store the old atomic positions, because as far as I know (I am using version 10-Aug-2015), there is no such array anywhere. I believe the PERI package has something like that but I can’t just wrap my head around it.

so I have a couple of questions here :

  1. Where should I allocate this array ? I assume in the AtomVec classes, is this correct ?
  2. When packing/unpacking the atomic data for interprocess communication, the old positions need to be packed/unpacked, however, I believe I need to modify pack_exchange and unpack_exchange for that, the problem is I am not sure where the size of the buffer passed to these functions is set, if I simply pack the extra data, will the communication message size adjust automatically since the pack_exchange function returns the index of the last packed data item ?
  3. Are there any other locations where I need to worry about communications ? I will store the old atomic positions in the Verlet class before the first call to any integrations but I just want to make sure that the atomic positions are not modified somewhere else.

I understand this is a somewhat involved question, I would be very grateful if you can give me rough directions.

Thanks

Hello All,

I am writing a compute that requires the instantaneous displacement of
each atom during the run as a part of its operation. That is, I want to
subtract the position of the atom at the previous time step from its
current position after integration. For that, I need to add an array to the
atom class to store the old atomic positions, because as far as I know (I
am using version 10-Aug-2015), there is no such array anywhere. I believe
the PERI package has something like that but I can't just wrap my head
around it.

so I have a couple of questions here :
1. Where should I allocate this array ? I assume in the AtomVec classes,
is this correct ?

​no.​

2. When packing/unpacking the atomic data for interprocess communication,
the old positions need to be packed/unpacked, however, I believe I need to
modify pack_exchange and unpack_exchange for that, the problem is I am not
sure where the size of the buffer passed to these functions is set, if I
simply pack the extra data, will the communication message size adjust
automatically since the pack_exchange function returns the index of the
last packed data item ?

​you don't need to worry about any of that.​

3. Are there any other locations where I need to worry about
communications ? I will store the old atomic positions in the Verlet class
before the first call to any integrations but I just want to make sure that
the atomic positions are not modified somewhere else.

​you should not need to make any code modifications outside of your
compute. just look at what compute msd does and create an instance of fix
store/state with suitable settings instead in your compute.

if you cannot make that work for you in the way you want, they you should
instead create an instance of fix property/atom with 3 double precision
entries (one for x, y, and z each) and then update and retrieve that data
instead.

adding any additional storage to any AtomVec classes for the purpose of a
compute or fix should be avoided at all cost. the fix mentioned fixes will
automatically take care of migrating the per-atom data across processors.

axel.

Thanks Axel

So how can I make LAMMPS call the compute every time step even if I don’t output anything (dump or thermo) every time step ? I am writing an integrator and I need to keep the cumulative sum of the integration kernel.

Thanks Axel

So how can I make LAMMPS call the compute every time step even if I don't
output anything (dump or thermo) every time step ? I am writing an
integrator and I need to keep the cumulative sum of the integration kernel.

compute_instance_variable->compute_peratom()

axel.

one question remains, though: why do you need a compute, when you can instantiate fix store/state directly from your integrator just as well?

But be sure not to instantiate another fix from the constructor of your fix.

The logic in modify.cpp won’t support that. You’ll have to do it later, like in init().

Steve

I am bit confused here, initially, I didn’t write an integrator, I only wrote a compute that keeps track of the old position and uses it when needed. This is working well now.

However, to keep track of old position in every time step, I need to call the compute every time step, but to my knowledge, this computes are only called when evaluating variables, dumping output or outputting thermodynamic states. I don’t want a dump file or a thermo line every step.

So what I did is creating a new fix, called it fix compute (just a name for now) that takes in any number of computes as arguments and invokes them every n time steps where n is another argument passed to the fix. I guess at this point I have “an integrator” which basically calls the computes I need to update. The fix does not store any data or pass anything back to anything else, it just invokes the computes in the end_of_step() function when necessary.

Does this sound like a good way to approach this problem ? Did I end up rewriting fix store/state ?

Two things are unclear about what you say:

I only wrote a compute that keeps track of the old position and uses it when needed.

What is “when needed”. I.e. when do you need them and what other command
in LAMMPS will use them? Also, what

does “old position” mean? Do you mean only the most recent timestep, or N timesteps

ago?

However, to keep track of old position in every time step, I need to call the compute every time step

Again, what LAMMPS command will call the compute every time step?

Did you write that command as well?

I don’t think your first email explained clearly what it is you actually want to do. Forget

about computes or fixes. If you just explain it at a high level you’ll probably

get some suggestions about how to do it in LAMMPS.

Steve

I am bit confused here, initially, I didn't write an integrator, I only
wrote a compute that keeps track of the old position and uses it when
needed. This is working well now.

​the confusion is in part due to not providing much context and asking for
specific solutions to specific problems, that - as more context becomes
visible - ​don't seem very suitable to the overall approach. even more so,
that all input we have is your "say so". you are basically asking us to
debug/develop your code without even seeing it.

However, to keep track of old position in every time step, I need to call
the compute every time step, but to my knowledge, this computes are only
called when evaluating variables, dumping output or outputting
thermodynamic states. I don't want a dump file or a thermo line every step.

​there are ways to signal on which timesteps a compute is needed to be
executed. but i am not going to share any more recommendations without
knowing more context.

So what I did is creating a new fix, called it fix compute (just a name
for now) that takes in any number of computes as arguments and invokes them
every n time steps where n is another argument passed to the fix. I guess
at this point I have "an integrator" which basically calls the computes I
need to update. The fix does not store any data or pass anything back to
anything else, it just invokes the computes in the end_of_step() function
when necessary.

Does this sound like a good way to approach this problem ? Did I end up
rewriting fix store/state ?

how are you going to handle atoms migrating to subdomains?
fix store state already does what you were asking for *with* proper support
of atoms migrating between subdomains. so yes, it currently sounds a lot,
that you've re-implemented available functionality, but in a much less
efficient way and possibly even incorrectly when running in parallel.

axel.

Ok, here is the full story.

I want to calculate the mechanical work done by an atom, the integral of the force multiplied by the displacement increment. This is a per atom compute.

For that, I need the position in the last time step to get the displacement increment, and I need to increment the cumulative sum of F.dx every time step whether I call the compute the does that or not every time step.

So I wrote the compute, tested it, it works fine, but in the test I dumped the data every time step. I called it compute power because it calculates the mechanical power as well as the mechanical work. Now to make things simpler, I don’t really store the old positions anymore, I calculate F.v and Sum_{time step 0}^{current time step} of F.v , this last sum is the mechanical work divided by the time step.

Now I need to call compute_peratom() every time step whether I dump the data or not. The immediate solution I thought of is to write a fix (basically, taking fix ave/atom and stripping it down to something that calls the computes list without storing anything). I called it fix compute (bad naming, I know). Now in the LAMMPS input script, I call compute power … then fix compute … , one of the arguments of fix compute is the power compute.

So I was a bit imprecise in my last post, I don’t explicitly keep track of old positions, I infer them from, or more precisely I replaced them in the calculation by, the current velocities.

To answer Axel’s question, no subdomain migration or periodicity issues can arise here, the velocities take care of it all.

In my last email I only wanted to ask how to invoke a compute every time step without outputting anything.

I would just make the whole thing a fix, e.g. fix_work.cpp

You can trigger it to be called every step, whenever you want

during the timestep. If you need to store some persistent

per-atom quantities (from timestep to timestep), you can

do that in fix work itself, or you can have it use fix store/state to

do it for you.

Fixes can output per-atom values (e.g. to a dump file or a variable), just

like computes can.

What you describe is really not what a compute is designed for.

All the code you wrote for your compute should basically drop

into a method in the fix.

Steve

I was just thinking about that, I noticed that the fix does everything and the compute is nothing more than a sum and a product operation.

Now speaking of fixes, I noticed that when I run with the fix, it gets called every N time steps as required, but also every dump operation where the compute is defined. I understand that it is better to move everything to the fix, but in general, if I choose to leave it a compute and a fix not just a fix, how to call the compute once per time step. I do

modify->clearstep_compute()

and I follow it by

if(!compute->invoked_flag |= INVOKED_PERATOM)

before I do

compute->compute_peratom()

and I checked the dump, it does the same thing. It seems to me that one doesn’t need a modify->clearstep_compute() before the if condition because it resets everything, however, I see it before almost every fix that uses computes, the logic here confuses me.

I apologize for the so many questions, I am a new LAMMPS user and it is the first time I extend LAMMPS and the complexity of the program is a bit confusing to me.

All you are really doing is a time integration of mechanical power F.v for one atom. LAMMPS already contains a lot of functionality related to this. For example, the compute heat/flux contains a very similar expression. Also, since you are primarily interested in the time integral of the power (work), you can use fix ave/time to do the integration. A good example of this is examples/VISCOSITY/in.einstein.2d:

fix avstress all ave/time $s $p $d v_pxy ave one

which captures the time integral of the momentum flux = f_avstress*$d*dt

I will add a note documenting this in fix_eve_time.html.

Aidan

I checked compute heat/flux, it is very similar to what I want, the main issue however is that it is not a per atom compute, hence, you don’t need to handle data exchange, an MPI_Allreduce call at the end takes care of summing up the vector J.

In my case, I need to compute the mechanical work per atom, I have this ready and working well and everything, but only on one processor, I implemented pack_exchange and unpack_exchange but they only get called every now and then (don’t know why, but there is a condition int he verlet integrator that calls them selectively).

So my problem now is, should I implement pack_comm_forward so that the cumulative sums (the running integration) get sent along with the coordinates ? Is there another solution ?

Thanks

I was just thinking about that, I noticed that the fix does everything and
the compute is nothing more than a sum and a product operation.

Now speaking of fixes, I noticed that when I run with the fix, it gets
called every N time steps as required, but also every dump operation where
the compute is defined. I understand that it is better to move everything
to the fix, but in general, if I choose to leave it a compute and a fix not
just a fix, how to call the compute once per time step. I do

modify->clearstep_compute()

and I follow it by

if(!compute->invoked_flag |= INVOKED_PERATOM)

before I do

compute->compute_peratom()

and I checked the dump, it does the same thing. It seems to me that one
doesn't need a modify->clearstep_compute() before the if condition because
it resets everything, however, I see it before almost every fix that uses
computes, the logic here confuses me.

I apologize for the so many questions, I am a new LAMMPS user and it is
the first time I extend LAMMPS and the complexity of the program is a bit
confusing to me.

​then you should even more follow the advice being given instead of
insisting to do something else.

axel.​

Axel,

I cannot describe how delighted I become every time I hear your kind, heart warming words. Your masterfully crafted advice never fail to amaze me. The amount of consideration loaded in every statement you make cannot reflect anything but deep, pure and unparalleled desire to offer the best of help to me in particular and the scientific community in general. We are yet to see a soldier with more generosity and patience than yourself in the battlefield of scientific endeavor.

Having said that, and as grateful as I can get, I would like to remind you that you don’t have to respond to people that you deem stupid, ignorant and stubborn. I understand that you are a very busy individual with a lot more important things to take care of than to respond to questions from people who are hopelessly rejecting your sharp advice and trying to think on their own, or worse, those who post questions without complete information.

In this regard I like to bring a fact that you know well back under the spotlight of your attention, this community is largely made up of scientists, who most of them happen to be curious. Asking vague questions, or thinking of a less than optimal solution to a particular problem, should be highly tolerated, even more, encouraged. This kind of intellectual diversity is a precursor for the natural selection process spontaneously taking place among scholars, which in turn constitutes a cornerstone of scientific and technological advancements. I am confident that you are more than aware of the history of human struggle against those who shun intellectual diversity.

So, speaking of the specific matter in hand, and speaking for the future matters that might fall in the same category of “questions asked by ignorant stubborns”, I would like to kindly ask you to be more patient with those who post to this mailing list. Some of them don’t know LAMMPS inside out, some of them like to experiment with different solutions, some of them are just taken aback by the complexity of the code, and some of them have a combination of all of these. And you know what is the most exciting thing of all, maybe some of them are here just to share their ideas with the broader community and engage in fruitful discussions rather than the “our way or the highway” approach. People are here primarily to learn, not to get the job done and go home.

I wholeheartedly like to emphasize that I did benefit from this very discussion, for real, but the vibe I am getting from this exchange makes me concerned about the extent of “freedom of expression” allowed on this mailing list. So please, if there is absolutely no way that you can tolerate the question being asked or the person asking it, you can clear the stage for other less experienced community members to give the asker a lower quality advice.

Hopefully, this might help bringing humanity one step closer to world peace.

Thank You

Ahmed

i am trying to not drag this into a discussion about philosophy of science. i’ll try to reduce what could be a long lecture to a few points.

  • the scientific method is at its very core based on precise observation and careful and complete description of those observations. vagueness has no room in here (not to be confused with creativity). i think it is not too much to expect people that identify themselves as scientists to expect them to be able to apply the scientific method.

  • i think i am very well justified to criticize somebody that asks for advice, identifies himself as a LAMMPS beginner, has been pointed to a viable solution by the expert, i.e. steve, and then turns around and says, “this is all fine, but i choose to do things differently, please now explain how to do this as well”. considering that people here are volunteering their time, what entitles you to get even more help, after you have been given advice that you yourself have identified as satisfactory? at the very least, you can show some respect (to steve in this case) by giving his suggestion a try.

  • if you feel challenged by my way of response, then you show exactly the intended reaction. if you would observe more carefully, i am quite patient with people asking for advice the first time. however, as the exchange develops, and people are changing their stories, ignoring advice, not improving their approach or leaving out vital information, i am cranking up the volume as well and become more challenging and critical.

  • you wrote that “People are here primarily to learn, not to get the job done and go home”. that is exactly my intention. i am quite successfully training people in different areas of scientific computing for over 15 years now, and even if it is sometimes not pleasant and it may make you feel uncomfortable at times (you are not alone there), the consensus of the large majority of people that i have worked with is, that they learned more this way than with any other (less challenging) approach. the fact that they are not 100% happy in the process is a hard to avoid side effect. no pain, no gain.

axel.

Axel,

I agree, while I would be very interested in a philosophy of science discussion, this is not the appropriate place for it. Maybe next time I am in Philly, we can grab a cup of coffee together and discuss philosophy at length.

A couple of short notes I want to make to take this discussion back to the technical nature it enjoyed 24 hours ago.

I did not make any vague scientific observations, conclusions or explanations. I did not make any observations or conclusions at all. The possible vagueness is in the questions asked by people who are less experienced with LAMMPS, and the vagueness in the question does not stem from a lack of desire to be less vague, it it only because beginners, while they know a lot about the problem in hand and the solutions they have tried, they don’t usually have a “list of most relevant details” ready in their head when asking questions and discussing technical issues. What you deem an important detail, and is an important detail, might be less important because I don’t have such an experience you have.

Whether I am a LAMMPS beginner or not, whether you have been teaching LAMMPS for three weeks or since the fall of the Roman Empire, you are free to criticize me, or any one for that matter, without justification.

Nothing really entitles me to get more help, nothing entitles me to get any help actually. You are are volunteering your time and I am very aware of and grateful for that. This is why I said if you cannot tolerate the question or the person asking it, you don’t have to respond. I am nothing but a peasant who went to the town center and called out loud asking for advice on how to milk his cow. You can be a considerate advisor or a silent passerby.

And by the way, I did give Steve’s solution a try, in fact, Steve’s suggested solution is what I currently have. I am not sure if you received all the communications in this discussion, because for some reason, my emails stopped appearing on the mailing list (the one on sourceforge), but anyway, I implemented what Steve suggested, and I also looked up an example suggested by Aidan. In my last email I asked a couple of questions related to that particular solution, because I really want to know how exchanging fix information between processors works, and I really wanted to know the logic behind clearstep_compute(), and a couple of other things. I apologize if you missed these emails, I will make sure to copy you on all emails.

I have to admit that my story changed over the course of the discussion. I believe it is better to have started a new thread or kept this as a private communication away from the mailing list. But you know how persistent this is in scientific discussions, where you start somewhere and the discussion drifts as one of the parties learns more about the subject.

Finally, I am not asking you to change your ways, and I am not claiming that your challenging approach is this or that. I am just saying that your approach might turn someone truly interested in learning off. Since you are not obligated to respond, and since you strongly believe in your approach, I am just suggesting that, in the interest of benefiting that part of the community who would shy away from asking further questions after observing the attitude, that you don’t bother communicating with the stubborn few that prefer a larger community discussion over a “we did it before and we know how it works” approach.

At the end, I am really thankful to you, Steve and Aidan for all the time you guys wasted on this. I truly benefited from every piece of advice, as for someone as beginner as myself, every piece of information is new information.

Best,

Ahmed

ok, I suggest we all just drop this thread and just move on.

As a general rule, when things become personal, I suggest folks take their

email discussion off-line. Ahmed has gotten some ideas

for what he wants to do, and hopefully can make progress.

Thanks,

Steve

Agreed. Thank you all.

I believe I am almost there, I will try to push things forward and get back to you once I start tripping over.