Computing peratom contributions from different individual pair styles

LAMMPS LAMMPS Mailing List Mirror
#1

Hello all,

I have a group of atoms of interest experiencing a sum of LJ energy+forces, Coulomb energy+forces, and “test” EAM potential energy.

The test comes from a duplicate of the pair_eam files. I have named the files and necessary terms to correspond to “pair_eamtest”, commented out the force additions, and successfully compiled them with the source code.

I was wondering if you have any pointers on to how to create peratom vectors that can individually tally the total LJ, total Coulomb, and total test EAM energies of each atom in the specified group. I looked through the compute tally and compute pair pages, but could not find something that quite matched this.

Any help would be greatly appreciated.

Thank you,
Anne

#2

There is no existing source code in LAMMPS that can do this directly.

If you need to do this data collection “online”, i.e. during the simulation, you might be able to do this by hacking compute pe/tally to collect and output 3 columns per atom instead of 2, but you would have to come up with a way to tell the EAM style contributions apart from evdwl and ecoul contributions when the tally callback is called. E.g. by setting both atom indices i and j negative and then flip their sign when tallying their (evdwl) data into the 3rd column. This would not be a modification that we would want to integrate into the LAMMPS distribution unless you can figure out a less hackish way to do this.

if you can do the analysis “offline”, an option would be to record a trajectory and then use the rerun command for analysis where you can collect evdwl and ecoul in the first round but do not define your custom EAM style and then unset the LJ and coulomb potentials and collect the EAM data. that should be possible without further code modifications.

axel.

#3

Hi Axel, thank you very much for the advice. I tried something similar to what what was suggested, flagging this particular pair computation and modifying the tally callback.

As a less hackish approach, would it be possible to replace the call to ev_tally with a callback to a custom tallying function in a new fix? The plan was for the fix to then access the vector of stored eam_test values for further operations after the force computations, and for eam_test to have no influence on the remaining energies/forces in the simulation.

Thanks,
Anne

#4

You have the source code, so you can change it any way you want as you already do. But we won’t do something like that in the release version of LAMMPS.

#5

Understood; are there any prior examples of a fix’s function being called to during force computation, saving the data, and then working on this data after the force computation? I tried searching online, but could not find something similar.

#6

Not that I know of

#7

But you could also just allocate and store the data inside the pair style on added arrays and then access it later by retrieving pointers to the internal data via the extract() function. That is used in several places.

Axel

#8

one more comment.

perhaps it would be more helpful, if you first explain what you ultimately want to achieve with your modifications instead of asking for specific solutions to problems that you encounter while trying to implement what you have devised.

from many years or answering on mailing lists, i’ve noticed that people sometimes get lost in solving problems with solutions for other problems that were solutions for some issue with their original approach and so on and then forget that there perhaps are other approaches to the original problem, which may be more straightforward to implement.

axel.

#9

Hi Axel,

Thank you for your feedback. To give more context, I want to calculate the energy that “would be” experienced by each individual atom of one type (Type 9 in an earlier email) if it was part of a greater lattice (Type 11) at its given position near the lattice. I only want to accumulate the influence of Type 11 atoms on Type 9, without: Type 9 affecting Type 11 back, Type 9s affecting each other, or this calculation being combined with other LJ/Coulomb forces experienced by Type 9. I took your advice and used “fix property/atom” to create a custom vector of doubles named “eTest”. I then successfully restructured the custom eam_test files to accumulate the energy accordingly:

int tmp_eTest = -1;
double *eTest = atom->dvector[atom->find_custom(“eTest”,tmp_eTest)];
[ modified EAM loops]

eTest[i] += …
[ modified EAM loops]

However, I have the following questions on accessing eTest in fixes separate from the force evaluation:

  • What exactly would be the difference between using find_custom() and extract() ? It appears both provide access to the custom array.

  • Is there an example for what syntax should be used to run extract() on a vector create by fix property/atom?

Thank you,

Anne Brant

Graduate Student

Micro and Nano Manufacturing Laboratory

Mechanical Engineering

University of Cincinnati

#10

I don’t fully understand your explanation, but if you want to run dynamics but be calculating
other quantities (from pair styles) that don’t contribute forces to the dynamics, e.g.
your test EAM energy … Then you could look at using the rerun command to
process dump snapshots from a previous simulation (w/out test EAM), and define
whatever pair styles you want in the rerun script. E.g. turn on EAM but only between
types you want, and have no potential defined (pair_style none) for other types.

Steve

#11

Thank you, Steve. I understand my application is a bit weird. However, I got it to work during force computation in the run, which is where it was needed.

The part where I am stuck is how to access the values of an array defined by fix property/atom in a custom fix after the force computation. I was wondering how extract() differs from find_custom(), and what syntax should be used for extract().

Thank you,
Anne Brant

#12

Hi Axel,

Thank you for your feedback. To give more context, I want to calculate the energy that “would be” experienced by each individual atom of one type (Type 9 in an earlier email) if it was part of a greater lattice (Type 11) at its given position near the lattice. I only want to accumulate the influence of Type 11 atoms on Type 9, without: Type 9 affecting Type 11 back, Type 9s affecting each other, or this calculation being combined with other LJ/Coulomb forces experienced by Type 9. I took your advice and used “fix property/atom” to create a custom vector of doubles named “eTest”. I then successfully restructured the custom eam_test files to accumulate the energy accordingly:

this was not my advice. i suggested something much simpler, i.e. to declare and allocate an array within the pair style and then make a pointer to it accessible via its extract() function. using fix property/atom is a much more convoluted and complex process that requires additional work in the input and passing this information along or creating the fix instance from within the pair style or something even more dirty and hackish.

int tmp_eTest = -1;
double *eTest = atom->dvector[atom->find_custom(“eTest”,tmp_eTest)];
[ modified EAM loops]

eTest[i] += …
[ modified EAM loops]

However, I have the following questions on accessing eTest in fixes separate from the force evaluation:

  • What exactly would be the difference between using find_custom() and extract() ? It appears both provide access to the custom array.

those are two different things with different purposes. the major difference is which class instance is allocating and maintaining the data. it would be very, very bad idea to re-export a pointer allocated though fix property/atom with an extract function. besides, you already have a clean mechanism to access this data, which you have already been using.

  • Is there an example for what syntax should be used to run extract() on a vector create by fix property/atom?

no. as i already pointed out, it is a bad idea to re-export this.

in comparison, i think this is also a worse solution than what we had discussed previously using a modification of compute pe/tally.

axel.

#13

Thank you Axel, I had misinterpreted your advice earlier. I have now instead declared eTest as a public vector in pair_eam_test.h

Using the extract function of atom.cpp as an example, I added this to the extract function in pair_eam_test.cpp:

if(strcmp(name,“eTest”)==0) return (void *) eTest;

In the fix, I believe I should then access the eTest vector as:

(?) ->eTest

However, I am stuck on what goes in the (?) section.

Also, if a post-force fix such as fix setforce uses atom style arrays like atom->x, atom->f, atom->type, etc, why would there be complications with a custom fix/property vector? I had thought that the communications of custom vectors are done in the same way the existing atom style vectors are.

#14

Thank you Axel, I had misinterpreted your advice earlier. I have now instead declared eTest as a public vector in pair_eam_test.h

it does not have to be public. more importantly, you have to allocate/resize/free it as necessary using memory->create() or memory->destroy() (see the array PairEam::rho for example)

Using the extract function of atom.cpp as an example, I added this to the extract function in pair_eam_test.cpp:

if(strcmp(name,“eTest”)==0) return (void *) eTest;

that is not sufficient, you also need to set dim as necessary to indicate the dimensionality. the existing example is for a 2d array, but you have only 1d data.

In the fix, I believe I should then access the eTest vector as:

(?) ->eTest

no. you have to check on your pair style and call force->pair->extract() with the corresponding settings to obtain the pointer. check whether it is not NULL and whether the dimensionality matches. force->pair->extract() is used in several places.

However, I am stuck on what goes in the (?) section.

considering the kinds of dirty hacks that you are doing, i am surprised of your limited ability to read and search through source code. if you keep doing what you are doing, you have to acquire and practice this kind of skill or you will be stuck all the time and there may not be people around that are willing and able to help you.
with grep you should be able quickly identify cases where the Pair::extract() method is used and learn from it.

Also, if a post-force fix such as fix setforce uses atom style arrays like atom->x, atom->f, atom->type, etc, why would there be complications with a custom fix/property vector? I had thought that the communications of custom vectors are done in the same way the existing atom style vectors are.

they are. i just said that it was rather complex and convoluted to have a (hacked) pair style that is dependent on you defining a specific property with fix property/atom in your input and then consume it with yet another fix. at the very least. this is just an invitation to get crashes. perhaps i am thinking too much in terms of what is a meaningful and maintainable modification in LAMMPS and less in terms of how to just hack your way to some result no matter how dirty and hackish.

you didn’t say that you want to access this data from fix setforce. in that case, you fix property/atom would have the advantage to already being exposed to the LAMMPS script interface and (atom style) variable expressions. but then again, i don’t understand why the compute pe/tally version wasn’t sufficient, that would be the cleanest approach in my opinion and also produce a per-atom property that is accessible in atom style variables. using Pair::extract() on the other hand would be most useful when trying to access it directly from a inside a fix, which is what you had specified you would be doing.

besides there is no communication necessary here. all data is only stored and accessed only for local atoms on the respective MPI subdomain.

at any rate, this is starting to get too hackish and too dirty for me. you are moving too far into uncharted territory with an approach that has no chance to ever been included into the LAMMPS distribution and thus useful for others but you.

feel free to continue on your path, but understand that i rather spend my time on modifications and improvements to LAMMPS that are of use for more than just one person. we are a very small team of people and have to budget our efforts carefully. up to this point, i had the hope that there would be something in here that could be later cleaned up and refactored to be usable in a more general way. but that is now gone and consequently so is my desire to spend time on it.

axel.

#15

Hi Axel, I understand, I really appreciate the help you have given on this.

Thanks,

Anne