Atom Type Change Based on Neighbor in a Polymer

LAMMPS LAMMPS General Discussion
,
1

I am working with a polymer system in LAMMPS where each polymer consists of beads of three types: 1, 2, and 3.

I would like to implement the following logic during the simulation:

  1. Randomly pick an atom (say, index ‘i’)
  2. If that atom is of type 2, check whether the ‘(i+1)th atom’ along the polymer is of type 1
  3. If so, then with 50% probability, change the type of atom ‘i’ to 3

This would ideally run during the simulation.

I explored the ‘set’ command, but it seems insufficient for:

  1. Accessing neighboring atoms along the polymer (like ‘i+1’)
  2. Doing probabilistic checks
  3. Dynamically changing type based on conditions

Is there a way to implement this using built-in LAMMPS commands?
Or would this require a custom fix?

Any suggestions, hacks, or references would be greatly appreciated!

I really appreciate any help you can provide.

2

Directly there is probably nothing. It could be possible to (ab-)use fix bond/react for this, but I don’t know enough about it to say for certain. @jrgissing can comment. He wrote that fix.

Writing a fix in C++ would be the straightforward approach. There are some challenges there, too, like identifying what constitutes a molecule and that they may be distributed over multiple MPI ranks.

Alternatives would be to call python from LAMMPS and use the LAMMPS python module to get access to the data of the simulation, but that is complicated and is subject to the same problems as implementing it in C++.

Another option would be to write a custom program that can transform a LAMMPS data file in the way you want and then run LAMMPS in small chunks, write out a data file, use the shell command to transform it, replace all atoms with those from the modified data file and continue for another chunk.

3

Is your polymer chain directional? Or, can the (i+1)th atom be in either direction?

Both should be possible with ‘fix bond/react’. Depending on other details of your system (e.g., whether or not it contains angles, dihedrals, etc.) it may be easier (require fewer reaction templates) to use a simplified version of ‘fix bond/react’ that turns off various sanity checks. Please let me know if interested.