Identifying molecules and changing atom types

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

I am running a script where I have one molecule type that can react with all other molecule types in the same way. The configuration of the first molecule type looks like:

X-X (two particles of type X)

The configuration of all other molecules looks like:

A1-B1-A1 (particles of type A1 and B1)

A2-B2-A2

A3-B3-A3

…

The X-particles can only react with A-particles and not with B-particles. I would like to implement this by using the fix bond/create command, but it appears that this will not work properly since I would have to use this command multiple times on the same atom type (type X). Since I’ve read that subsequent fix bond/create commands don’t communicate within the same time step, this will probably not give the result that I’m looking for. Therefore I’m thinking about doing it differently. At the start of the simulation I would assign all particles with atom types A1, A2, A3, etc to the same atom type, that is atom type A. This would allow me to use a single fix bond/create command to form all bonds between A and X. After all bonds are formed I would then change the atom types back to their original atom type. This should still be possible, since I didn’t change the atom type of the B-particles, so in principle I should be able to identify each A-particle. The only question now is how to do this last step of identifying the B-particles and then assigning the respective A-particles to their original atom type?

Kind regards,

Tim

2

Therefore I'm thinking about doing it differently. At the start of the simulation
I would assign all particles with atom types A1, A2, A3, etc to the same atom
type, that is atom type A. This would allow me to use a single fix bond/create
command to form all bonds between A and X. After all bonds are formed
I would then change the atom types back to their original atom type.
This should still be possible, since I didn't change the atom type of the
B-particles, so in principle I should be able to identify each A-particle.
The only question now is how to do this last step of identifying the B-particles
and then assigning the respective A-particles to their original atom type?

I wrote a fix that enables atom type changes in response to bonds
breaking and forming, (So did Tim Sirk, and Jacob Gissinger, see
below). Alas apparently this new fix does not yet work reliably in
parallel. I've been using it anyway (running slow simulations on 1
CPU). Getting this to work in parallel has simply not been a priority
for my employer so far. (That is about to change.) Anyway, I can
share with you what I have.

but it appears that this will not work properly since I would have to use this
command (fix bond/create) multiple times on the same atom type (type X).
Since I've read that subsequent fix bond/create commands don't communicate
within the same time step, this will probably not give the result that I'm looking for.

I think a simple work-around would be to just use multiple fix
bond/create commands, and hack the fix_bond_create.cpp source code to
add an extra "delay" argument:

fix f1 all bond/create 100 0 X A1 Rmin
fix f2 all bond/create 100 1 X A2 Rmin
fix f3 all bond/create 100 2 X A3 Rmin

(...where "X", "A1", "Rmin",... are numbers)

The "0", "1", "2" in the example above are the extra arguments, and
they mean the following:

f1: when t = 0, 100, 200, 300, ... check to see if we need to create
bonds between X and A1

f2: when t = 1, 101, 201, 301, ... check to see if we need to create
bonds between X and A2

f3: when t = 2, 102, 202, 302, ... check to see if we need to create
bonds between X and A3

In your case, the order of the fixes ("f1", "f2", "f3") do not matter.
Some people might consider this approach an ugly kluge. The important
thing is that the system is in a consistent state after every time
step. By the time we need to decide whether to for a bond between "X"
and "A2", all of the bonds between "X" and "A1" have already been
formed, and all of the processors know about these new bonds.

You can implement this in the code by editing the LAMMPS code. Here's how

(If anyone cares, I'm happy to put in a pull request with this feature.)

1) Edit "fix_bond_create.h", and add the following private member to
the "FixBondCreate" class:

   int nevery_delay;

2) Near the beginning of the "FixBondCreate::FixBondCreate()"
function, I would add this code:

   nevery_delay = 0;

(...or add "nevery_delay(0)" to the end of the class initializer list,
if you prefer that)

3) Near the end of the "FixBondCreate::FixBondCreate()" function, I
would add this code:

    } else if (strcmp(arg[iarg],"delay") == 0) {
      if (iarg+1 >= narg) error->all(FLERR,"Illegal fix bond/change command");
      nevery_delay = force->inumeric(FLERR,arg[iarg+1]);
      iarg += 2;
    }

4) At the beginning of "FixBondCreate::post_integrate()", I would
replace this line:

  if (update->ntimestep % nevery) return;

... with this line:

  if ((update->ntimestep - nevery_delay) % nevery) return;

-- fix bond/react --

I suggest you definitely take a look at fix bond/react. There's a web
site with a great animated gif showing what this fix can do but I
forgot the link. Here are some other links:

https://www.sciencedirect.com/science/article/pii/S0032386117309114

It's cool. It can probably do what you need, but it might be overkill.

I hope this helps.

Andrew

3

Therefore I'm thinking about doing it differently. At the start of the simulation
I would assign all particles with atom types A1, A2, A3, etc to the same atom
type, that is atom type A. This would allow me to use a single fix bond/create
command to form all bonds between A and X. After all bonds are formed
I would then change the atom types back to their original atom type.
This should still be possible, since I didn't change the atom type of the
B-particles, so in principle I should be able to identify each A-particle.
The only question now is how to do this last step of identifying the B-particles
and then assigning the respective A-particles to their original atom type?

I think a simple work-around would be to just use multiple fix
bond/create commands, and hack the fix_bond_create.cpp source code to
add an extra "delay" argument:

fix f1 all bond/create 100 0 X A1 Rmin
fix f2 all bond/create 100 1 X A2 Rmin
fix f3 all bond/create 100 2 X A3 Rmin

My apologies, that should have been:

  fix f1 all bond/create 100 delay 0 X A1 Rmin
  fix f2 all bond/create 100 delay 1 X A2 Rmin
  fix f3 all bond/create 100 delay 2 X A3 Rmin

(Hopefully replying to my own post does not cause the rest of it to be
hidden from view. If so, please check the previous message on this
thread. Again, if anyone wants, I'm happy to put in a pull request to
add the "delay" argument to the "fix bond/create" and "fix bond/break"
commands.)

Cheers
Andrew

4

Therefore I'm thinking about doing it differently. At the start of the simulation
I would assign all particles with atom types A1, A2, A3, etc to the same atom
type, that is atom type A. This would allow me to use a single fix bond/create
command to form all bonds between A and X. After all bonds are formed
I would then change the atom types back to their original atom type.
This should still be possible, since I didn't change the atom type of the
B-particles, so in principle I should be able to identify each A-particle.
The only question now is how to do this last step of identifying the B-particles
and then assigning the respective A-particles to their original atom type?

I think a simple work-around would be to just use multiple fix
bond/create commands, and hack the fix_bond_create.cpp source code to
add an extra "delay" argument:

My apologies, that should have been:

  fix f1 all bond/create 100 delay 0 X A1 Rmin
  fix f2 all bond/create 100 delay 1 X A2 Rmin
  fix f3 all bond/create 100 delay 2 X A3 Rmin

(Incredible.) My apologies -again-. That should have been:

   fix f1 all bond/create 100 X A1 Rmin Btype delay 0
   fix f2 all bond/create 100 X A2 Rmin Btype delay 1
   fix f3 all bond/create 100 X A3 Rmin Btype delay 2

(Sorry, my last two messages were typed on a cell phone. This time
I'm on a computer and I checked the input files for my own version of
fix bond/create. This time commands I just posted above use the
syntax which is compatible with the code I posted in my original
reply. This works. Again, hopefully replying to my own post does not
cause the rest of it to be hidden from view. If so, please check the
previous message on this thread. Again, if anyone wants, I'm happy to
put in a pull request to add the "delay" argument to the "fix
bond/create" and "fix bond/break" and "fix atom/swap" commands.)

Hi

5

Dear Andrew,

Thank you very much for your reply (well, actually replies). This helps a lot.

It might indeed be useful to add the delay to this command in the future.

Kind regards,Tim

6

Does this “delay” option force the code to reneighbor on each of steps 100,101,102, etc?

Steve

7

Hi Steve
   For reference, in Tim's example where there are 3 fix bond/create commands:

   fix f1 all bond/create 100 X A1 Rmin Btype delay 0
   fix f2 all bond/create 100 X A2 Rmin Btype delay 1
   fix f3 all bond/create 100 X A3 Rmin Btype delay 2

For fix-ID f1, post_integrate()'s contents are run at timesteps 100,
200, 300, ...
For fix-ID f2, post_integrate()'s contents are run at timesteps 101,
201, 301, ...
For fix-ID f3, post_integrate()'s contents are run at timesteps 102,
202, 302, ...
Any of them can set "next_reneighbor = update->ntimestep;"

The change in "fix_bond_create.cpp" to make this work (attached) is
minimal. Here is the diff output.

44,45c44
< probability(NULL), created(NULL), copy(NULL), random(NULL), list(NULL),
< nevery_delay(0)

src_fix_bond_create_with_delay.tar.gz (41.4 KB)