I am trying to write a plugin in LAMMPS that would need to compute the energy between two atoms as a function of the first atoms coordination number. From the LAMMPS GitHub repository, I found that there is a compute command that can do the calculation of the coordination number. This brings me to the following question: Can I call that from a pair style (let’s say gaussian)?
The specifics of the problem are as follows:
I need some way to bias the energies obtained from using pair_style gauss/cut to be dependent on the current coordination number of the central atom. This means that I need to modify the hgauss term in the gaussian on the fly.
Yes, there is no precedent, but I don’t see why this should not be possible.
There are two ways to deal with this.
- You (re-)create the compute directly inside the pair style constructor or the init_style() method of the pair style and store the created compute ID.
- You create the compute explicitly in your input and then pass the compute ID as an argument to the pair style and store it.
Inside the compute() method (or single(), if available) you then can use:
auto icompute = modify->get_compute_by_id(id_coord);
to get a pointer to the compute and then use icompute->compute_peratom(); to trigger the computation. Please keep in mind that this will only compute information for local atoms. If you need the coordination number also on ghost atoms, you will have to implement a forward communication.
you’ll likely also need to add:
modify->addstep_compute(update->ntimestep + 1);
at the end of the respective functions to avoid an error.
Since the coordination number is rather simple to compute, you could also do that directly inside the pair style. You would just loop the neighbor lists twice (for optimization purposes, you can also request a second neighbor list from the pair style that has a shorter cutoff). Collect the per-atom data directly and then use a reverse communication to have this information available (possibly with a forward communication to send the information also to ghost atoms as mentioned above). You can look at pair_eam.cpp as an example where first the property “rho” is computed followed by a reverse communication to collect contributions from ghost atoms and from the derivative of the embedding energy “fp” which is then subject to a forward communication.
Thank you, Prof. Kohlmeyer, for your response. I was considering writing the section corresponding to the calculation of the coordination number in the computation of the pair_gauss_cut.cpp file directly. I do not have ghost atoms in my simulation, and I was trying to reuse the pairlist information that is available in the atom object. To be precise, I was modifying the section that computes the geometric distance between the atoms, because that is the information that I also require to compute the coordination number.
Once you run in parallel, there will be ghost atoms.
This makes no sense. There is no “pairlist” stored in “atom”.
The neighbor list according to the neighbor list request is stored in the pair style during the call to the method init_list(). Please also note that compute coord/atom requests a full neighbor list while pair style gauss/cut uses a half neighbor list (same as eam, BTW).
There is a detailed discussion of the various computational steps in a pair style in the LAMMPS manual in the programmer’s guide section at: 4.8.1. Writing new pair styles — LAMMPS documentation
As I mentioned before, you must follow the EAM pair style for a direct implementation, since with a half neighbor list, you also have information stored on ghost atoms and that has to be collected back to their local counterparts with a reverse communication and then you will have to use a forward communication. Also, you cannot have the complete information about neighbors without walking the neighbor list twice.
If you bypass the steps I have outlined you are likely to have inconsistent information and then bogus results and ultimately have wasted your time.
Thank you Prof. Kohlmeyer. I will follow the protocol for EAM pair style.
I cannot stress enough how important it is to take the time to fully understand what you are doing and to exploit the available information to the max by reading the manual and studying the source code.
1 Like