I am trying to implement a simulation where a molecule consists of two bonded particles — a protein and a patch. These exist in two possible states: Active and Inactive, defined as follows:
group activeprotein type 5
group activepatch type 6
group inactiveprotein type 3
group inactivepatch type 4
The state switching is currently implemented as:
# type 5/6 switches to type 3/4
set group activeprotein type/fraction 3 ${f_fracswitching} ${switchseed}
set group activepatch type/fraction 4 ${f_fracswitching} ${switchseed}
However, this approach causes independent switching of the protein (type 5) and patch (type 6).
What I would like instead is coupled switching — that is, when one molecule is chosen to switch, both bonded particles (the protein and its patch partner) should simultaneously change their types (5→3 and 6→4).
Could you please suggest how this can be implemented in LAMMPS? For example, is there a way to perform type switching based on bonded pairs rather than independent particles, or to synchronise the switching of bonded partners?
I cannot think of an easy and reliable way to do this kind of manipulation from the LAMMPS script interface, especially not during a simulation (if that is required).
It would be possible to process this from within C++ code. For example, you can access the list of atoms connected to any given atom in the Atom::special list. This is, for example, used in compute fragment/atom to determine clusters of atoms connected by bonds.
If you need this as a one-time operation, this can be implemented as a command style, otherwise it would be a fix style (compare to create_bonds and fix bond/create).
Suppose you know for sure that each bonded pair of particles has a unique molecule number. (You may be able to enforce this using reset_atoms mol ...)
Then you can perform the desired manipulation using good old chunks:
group inactiveprotein type 3
group activeprotein type 4
group inactivepatch type 5
group activepatch type 6
# assuming all proteins and patches are initially active:
set group activeprotein type/fraction 3 ...
# set up molecular chunks
compute molchunk chunk/atom all molecule
# calculate minimum of type for all molecules ...
variable atype atom type
compute mintype reduce/chunk all molchunk min v_atype
# ... and spread back to all atoms, including patches
compute spreadtype chunk/spread/atom all molchunk c_mintype
# this variable will be 5 for molecules with inactiveproteins
# and 6 for molecules with activeproteins
variable patchtype atom c_spreadtype+2
# which makes this possible:
set group activepatch type v_patchtype
Note that you need to reorder your types to make this work (if inactivepatch’s type is numerically smaller than activeprotein’s type, this won’t work). Also the usual caveats apply (please make sure you understand what this does, you may need a run 0 to properly populate variables and computes, etc etc).
In particular, please be very aware that LAMMPS groups are not dynamic by default. That is, the four groups you’ve defined will not change membership just because particles’ types have changed. You either need to actively make these groups dynamic, or to redefine the groups based on particle type after particle types have changed (if this is just once at the start of the simulation rather than happening periodically).
An alternative could be to use compute coord/atom to identify type 5 atoms that have type 4 atoms as close neighbor. This would assume some geometry constrains in that there would be only one pair of relevant atoms within the cutoff radius.
One should point out in this context that the group command is incremental, i.e. it always adds to a group. To completely redefine a group, it must first be deleted.
Thanks for your insights. Following your ideas, I implemented it a bit differently (Thanks to LAMMPS documentation module!).
# All proteins and patches are initially active
variable a loop ${switches}
label swloop # start of loop
#switching a fraction of type 5 to type 3 (Inactivating)
set group activeprotein type/fraction 3 ${f_fracswitching} ${switchseed}
# Grouping all type 3 (new/old) to a temporary group
group switch type 3
# Extending group to include the bonded partners (all bonded pairs have separate mol id)
group switch include molecule
# adding freshly inactivated patch-partners
group inactivepatch intersect activepatch switch
# setting up inactive patch type
set group inactivepatch type 4
# regrouping
group switch delete
group activeprotein delete
group activepatch delete
group inactiveprotein delete
group inactivepatch delete
group activeprotein type 5
group activepatch type 6
group inactiveprotein type 3
group inactivepatch type 4
# --- Run steps for this switching cycle ---
run ${switchsteps}
next a
jump SELF swloop
Please let me know your thoughts. Here is a log-file snippet-- how it performs:
set group activeprotein type/fraction 3 ${f_fracswitching} ${switchseed}
set group activeprotein type/fraction 3 0.01 ${switchseed}
set group activeprotein type/fraction 3 0.01 60214
Setting atom values ...
1 settings made for type/fraction
group switch type 3
333 atoms in group switch
group switch include molecule
666 atoms in group switch
group inactivepatch intersect activepatch switch
333 atoms in group inactivepatch
set group inactivepatch type 4
Setting atom values ...
333 settings made for type
As I wrote in my first response, to implement what you were asking for reliably requires C++ code since you need to access data structures that are not otherwise accessible.
Everything else would require heuristics based on specific properties of your system. You have not shared any such information, so it is impossible to comment in a meaningful way on your approach. In the end, it is up to you anyway to be certain that your choices will work as intended, just like it is for everybody else on their simulations.