I am trying to write a compute function to find the truncated radius of a granule. By that I mean, if one granule is considered, depending on how that granule is compressed, it will expand from its initial nominal radius. The new radius is given by a closed-form solution which is dependent on the center to center distance between touching/neighboring granules. If you look at the figure below, essentially I’d like to loop over neighboring granules, find the little h, and use the closed-form equation to update and find the new large R.
There are a few caveats though. One is that this compute involves history-dependent values (i.e. the value of little h influences terms in the closed-form solution). How can I keep track of that within the compute? For example, if the instantaneous little h is smaller than its previous value, that contact has less influence on the closed-form expression.
However, I believe this can be obtained by writing a compute. I have several questions I hope you don’t mind me asking.
When a compute is invoked, is it calculated at the beginning of the timestep?
How can I keep track of history-dependent values within compute? i.e h and the new truncated radius R
Is there a similar compute function that I can use as a template?
Potentially, what files would I have to modify when writing this compute?
Thanks in advance! Appreciate any help as I’m still learning how to modify LAMMPS.
Computes are invoked when they are consumed. If they are referenced a second time during the same timestep, the computation is not repeated, but the unchanged results are presented.
You need to create an internal fix or custom property that can store your data and migrate it with the atoms as they move between processors. If you need that data also be associated with ghost atoms (and not just the local atom), then you probably need to use fix property/atom with the ghost yes flag.
I cannot think of anything very similar. The computes msd, vacf, and displace/atom store data with fix STORE.
Ideally, you would only need to write your compute class with its header and implementation files. When configuring LAMMPS for recompilation, it should pick it up automatically if you put it into the src/ folder. Longer term, you may want to move it to the GRANULAR folder. But that can be done later.
Within the compute function of your pair style, you would loop over the list of neighbors and identify those that are in contact (that is R_i + R_j is larger than the distance between the two, d_ij). If I am not mistaken, then h_ij would be 0.5*(R_i + R_j - d_ij).
This can then be stored in a per-atom array in the pair style itself, and if you provide the suitable pack/unpack functions it would also migrate with the atoms. You can also have a h_old that would be the value from the previous step. And you can have multiple such h values per particle, for convenience you would just set a maximum supported number of contacts and then allocate and manage another per-atom property that is the count of contacts (and thus h_ij values).
One caveat would be that you must not update the per-particle radii until after the force computation is complete, or else the results would depend on the order in which you process them. Another is that updating the radius would also require to update the per-particle mass or the per particle density.
Thanks for both your detailed replies, and for thinking deeply about this. I agree with using a pair style to do this as it stores contact level properties similar to what i need. What you mentioned about h_ij is correct but the problem is R_i and R_j dynamically changes based on the state of compression of the granule at a given timestep. So my thought process was to use a compute to find the new R for all the granules and then calculate h_ij.
Eventually after the calculation of the new truncated radius (R), I use the new R to calculate the force. For example if the truncated radius increases, there will be more overlap using the equation you provided and it will change the force.
My concern is the following which I think you alluded in this comment:
If I were to implement your suggestion, how will the compute function within the pair style know the R of both R_i, and R_j? During the i for loop for example in pair gran hertz history, I can use similar syntax to calculate R_i, but the force calculation would require the use of R_j as well.
I guess, I can have two sets of i and j for loops. The first i and j for loop is to compute the force for all the granules, and the second is for the computation of the new radii? All of this within the compute function of the pair style. Then I can probably store contact level, history-dependent properties like h_old, similar to the shear vector in gran hertz history?
Yes, you can look at pair style EAM for an example. It first computes a “density” property rho from all pairs and stores it with the individual atoms. That needs to be “collected” with newton_pair on from the ghost atoms with a reverse communication. And from that density it computes an embedding energy term fp which is then used in a second loop over all pairs to compute the force. To have the fp data available to all atoms, it is distributed to the ghost atoms with a forward communication.
If you don’t need to store the “shear” properties that the hookian or hertzian styles with history use, you could “hijack” that exact same mechanism (implemented via a fix) and store up to 3 floating point properties per pair of atoms with it. Rolling your own or storing more properties would be tricky since those properties need to be “managed” during neighbor list updates and thus would require substantial changes to a part of LAMMPS that is not so easy to change.
Thanks so much Dr Axel. I leave you with two last questions, do I need to be concerned about ghost atoms if I only intend on using the serial version of LAMMPS? Or would it mess up the way the rest of LAMMPS is written if I don’t include it?
Do you also mind sharing the name of the .cpp file which has pair style EAM you referenced? Is it the same as that found in the MEAM folder in /src?
Yes, always. For the same reason as for granular pair styles you need to use comm_modify vel yes with granular pair styles, since they need updated atom velocities on ghost atoms.
You can learn more about ghost atoms here: 4.4. Parallel algorithms — LAMMPS documentation
I started working on this pair style and have run into a problem trying to understand the pack communications. Until now, I did 2 things.
The first is to create a new atom style variable to store the continuously updated truncated radius (defined as trunc_rad). I did that by modifying atom.cpp and atom_vec_sphere.cpp by looking at how existing variables are created.
I also allocated larger array sizes to store shear history info for pair_gran styles by editing fix_shear_history. So instead of storing 3 floating point properties, I changed it to accommodate 11 by inferring how the previous shear vector values are stored and changing appropriately. I thought this might be better than hijacking the 3 variables as I need these 3 original variables.
If I am to do:
Must I define the following methods like in pair_eam.cpp (as shown below but instead of rho, I use trunc_rad?). I’m struggling to see what buf, m, and last are.
int PairEAM::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) buf[m++] = rho[i];
return 1;
}
/* ---------------------------------------------------------------------- */
void PairEAM::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
rho[j] += buf[m++];
}
}
You need the reverse communication only, if you have Newton pair set to on and need to compute a property that is the sum of contributions from neighbors (like the sum of density contributions for EAM.
The meaning of “n” is the number of relevant ghost atoms for this sub-domain part, “buf” is a pointer to a communication buffer, first is the local index of the first local ghost atom and then last the last+1 index.
I’ve made some progress writing this however I have one bottle neck. Is there a way to retrieve the original radius I specify in the input file as I dynamically keep updating the radius? I get NaNs after a few 1000 timesteps and I speculate it to be a result of this issue.
To elaborate, I’d like to update my radius=truncated radius to help build future neighbor lists but my problem is that the truncated radius itself is a function of the original radius. Is there a way I can call the initial radius and still be able to set radius=truncated radius,
You have to make a copy. E.g. by creating an instance of fix STORE internally and assigning data to it and then retrieve that. Or use fix store/state externally and then passing the fix ID to the pair style and then look it up and access the data.
Thanks Dr Axel. I like the second option more. So basically I would have to use ‘fix store’ prior to a run command in my input script, and then access the data from that fix?
Is there a template I can follow that can retrieve data from a fix? I tried searching for this in the archive and didn’t see anything of importance showing up.
All examples I know use fix STORE and none of them is a pair style.
The principle is straightforward:
you get a pointer to the fix by using modify->get_fix_by_id() and assign it to a local Fix pointer variable, say Fix *fix_store;
You then query the array dimensions in the compute via looking at fix_store->size_peratom_cols and fix_store->size_peratom_rows. “cols” should be 0 in your case.
You can then access the stored data via fix_store->vector_atom
I have written the ::compute() method for the pair style, and now I want to write the ::single() method to access the pair-wise forces so that it can be written into a dump file. My question pertains to when is the single() method invoked? Is it activated whenever ::compute() is?
Asking this because I have defined variables in compute() that I would like to access in single()… Can I simply pass these variables/pointers as arguments to my single() method if so?
Only in special cases. E.g. when you use pair_write, compute group/group ,or compute pair/local in some use cases of fix gcmc and related and so on. You can easily find most of those locations by running the command:
grep -- 'pair->single(' src/*/*.cpp src/*.cpp
in the lammps tree.
That is a very bad idea™. The single function has to stand on its own. If you want to share data, you can make those class members, but the single() may be used out of context or at any time, so it must not depend on anything done in the compute() function previously. There is one exception for that in pair style eam, and it is creating a lot of headaches and inconsistent behavior.
Of course not. The interface must remain unchanged or it won’t be recognized as overload of the (virtual) base class function and thus not be called.
Thanks for the advice. After looking at it some more, I would like to use compute pair/local.
Since my potential is very similar to the eam pairstyle, I guess I can borrow from that? Are there any examples of using class members to share data that I can look at?
Is it also possible for me to store desired output in the shear vector and then just convert it to the svector within single()?
