I am simulating a system of atoms in which I have two reservoirs connected by a cylindrical pore. I am using hybrid/overlay to use both the Yukawa potential and lj/cut as pair styles describing the interactions between the atoms. I am also using fix wall/table and fix langevin to apply additional forces to the atoms.
I want to scale the forces from the pair styles (specifically the Yukawa potential) to the atoms that are only within the pore. For example, if atoms 27,214, and 319 are within the pore, I want to compute the pair-wise forces and scale it by X10, and addforce it to the respective atoms. So, for atom 27, I would sum the forces fz(27,214) and fz(27,319), scale them by X10, and add this force to atom 27, and so for the other two atoms.
I was searching for a command such as compute group/group but something like compute group/atom but I am not aware of such a command. My thought process is to define a dynamic group containing the atoms within the pore. Then, I use compute pair/local to calculate the pair forces between atoms strictly within this dynamic group. The output of this compute is an array of pair forces (and/or distances, energies), but I couldn’t find a way to map those pair-wise interactions to the pairs of atoms—there are no atom identifiers that I am aware of. Is there a logic behind how the order of the pair interactions appear such that I can deduce the corresponding atoms? Is it for example ordered by the atom ID? Or perhaps the order in the neighbor list? I compared the trajectory of my system with the output from the compute pair/local but I couldn’t find a trend of how the ordering of the interactions works.
I appreciate any suggestions to the above, and perhaps if there is a neater way to achieving force scaling through a different pipeline.
The order is given by the order of the pairs in the neighbor list.
How to determine the identity of the atoms in the pair is explicitly given in the compute pair/local
documentation: you use compute property/local and output the atom ids alongside the force info.
Looks like you need to spend more time reading the LAMMPS manual.
Doing this kind of thing through LAMMPS scripting is always going to be cumbersome and slow. The more effective approach would be to write a custom fix, request and walk the neigbor list for affected atoms only and compute the adjustment required for the full pair force, accumulate and apply it.
Indeed, it slipped from me although I read it a couple of times. Thanks for pointing it out.
I used the property/local to extract the atom IDs corresponding to pair interactions calculated in pair/local. Is there a way to natively aggregate forces from pair/local and property/compute computes into atom-style variables within the input script? While such variables are evaluated per-atom, they cannot loop over all pairwise interactions, but I was wondering if there is leeway around this. Perhaps casting the vector output from those computes into per-atom variables. Other options include:
(i) cumbersome looping within the input script using criterion like id==c_pairatomid[i],
(ii) using fix external + Python scripting to do the whole operation, or
(iii)
While I am not expecting to using this selective force scaling frequently in my simulations, I am tending to walk down the options from (i) to (iii) to avoid over-engineering my problem unless needed. Any suggestions? Thanks.
P.S. I am assuming this is relevant to this thread, but I could open a new topic if required.
Any kind of post-processing of pairwise interactions is going to be costly and also has limitations for applications to certain pair styles and atom styles.
The lowest overhead is to create a compute like the ones in the TALLY package.
The next option is to create a new pair style by creating a derived class from the original and modifying the compute() and single() (if available) functions.
Then I would consider writing a fix (the fix can create the computes and access their data internally), but please note that this will recompute pairwise forces in a rather inefficient way.
I don’t think the other two options you mentioned are viable, particularly not using LAMMPS scripting. That will “unparallelize” the process and is restricted by the limited capabilities of the LAMMPS input file script syntax.
In your situation I would write a fix that (1) detects whether particles have entered or left a given region at the end of step and (2) changes their particle type from “outside” to “inside” and vice versa as appropriate.
You can then set the potential parameters differently for the “inside” and “outside” types in the simulation setup, and then let the usual LAMMPS structures do all the work for you.
Here’s an example of what your eventual input script could look like:
pair_style lj/cut 10.0
pair_coeff 1 1 1.0 3.0
pair_coeff 1 2 1.0 3.0
pair_coeff 2 2 0.1 3.0 # scale down interactions to 0.1 within cylinder
...
fix MUTATE all mutate/region 1 cylinder type 1 in 2 out 1
# 1: every 1 steps
# cylinder: check if particles have entered or left region "cylinder"
# type: change particles' types upon entering or leaving -- you may want other keywords such as "charge"
# 1: we'll only check one type
# in: starts list of types -- there must be only 1 entry based on the previous number
# 2: a type 1 particle (see later) going inside gets changed to type 2
# out: starts list of types -- as with "in", there must only be 1 entry
# 1: a type 2 particle going outside gets changed to type 1
While this is still a non-trivial exercise, it is still relatively much simpler than directly modifying forces, energies, and so on.