Dear Lammps Users:
I am trying to develop a python script which creates a LAMMPS instance with OPLS force field.
OPLS is a non reactive model, so it requires bonding information at the beginning. I would like to mimic chemical reactions with the help of call back function method.
During the MD simulation “ fix all callback myfunction” will call “myfunction” right before the verlet step.
In myfunction I would like to check distances between particular atoms and if they are so close or certain conditions met I would like to create/break bonds at the simulation.
I have 2 questions at this point.
1- Is it possible to access inter atomic distances which is already calculated by LAMMPS at that point?
2- How to create/break or basically update bond topology of lammps instance.
As far as I know there is no functions for these two at the python interface and probably I need to write one. Can someone guide me how to do this?
Thanks in advance.
Dear Lammps Users:
I am trying to develop a python script which creates a LAMMPS instance with OPLS force field.
OPLS is a non reactive model, so it requires bonding information at the beginning. I would like to mimic chemical reactions with the help of call back function method.
During the MD simulation “ fix all callback myfunction” will call “myfunction” right before the verlet step.
In myfunction I would like to check distances between particular atoms and if they are so close or certain conditions met I would like to create/break bonds at the simulation.
I have 2 questions at this point.
1- Is it possible to access inter atomic distances which is already calculated by LAMMPS at that point?
2- How to create/break or basically update bond topology of lammps instance.
As far as I know there is no functions for these two at the python interface and probably I need to write one. Can someone guide me how to do this?
Thanks in advance.
Rather than do this using “fix external”, is there a reason you don’t want to use “fix bond/react” ?
Here is the paper describing their method:
Modeling chemical reactions in classical molecular dynamics simulationsJacob R.Gissinger, Benjamin D.Jensen, Kristopher E.Wise, Polymer, 2017
From the description you provided, this fix should do everything you want.
Andrew
P.S. We are also working on a LAMMPS fix to achieve these same goals, but it’s not quite ready to release yet. When writing your own fix, being able to access the atom positions and bonds is not the difficult part. If you are writing a custom fix, the atom ID-numbers, positions, and bonds owned by each processor are located in, atom->tag, atom->x, neighbor->bondlist, neighbor->nbondlist, respectively. (See “fix_bond_break.cpp” and “fix_bond_create.cpp” for details.) If you are writing an external call-back function, then it would have to access these variables. But this isn’t the difficult part. It would also have to inform each processor whenever an bond is broken, or an atom type has changed, and also tell LAMMPS to pack/unpack this information so that it can be shared with neighboring processors. I don’t know if that can be done using an external callback function. It seems better to write a fix to do this. (That is, assuming you don’t want to use “fix bond/react”.) My own experience writing these fixes has taught me that it is easy to make a custom fix code work on one processor. But it can be tricky to make them work robustly in parallel.