compute pair reaxff for group of atoms

Dear all,

I wanted to compute hydrogen bond energy of a group of atoms using the ReaxFF force field.

I have selected the atoms using the group command as follows:

group water1 id 13351:13353

and calculated the hydrogen bond energy for the group using the compute pair reax/c command, as follows:

compute reax1 water1 pair reax/c

and used variable command to save the output, as follows:

variable ehb1 equal c_reax1[8]

However the calculated result was the total hydrogen bond energy of the system not the group of atoms I selected.

I have been using lammps version 31Mar17

How do I calculate properly the hydrogen bond energy of a group atoms?

Thank you.

Kind regards,
Iwan

Though it is not obvious from the doc page, it does state that “This pair style tallies a breakdown of the total ReaxFF potential energy into sub-categories, which can be accessed via the compute pair command as a vector of values of length 14.” This feature is applied to group all.

Use re-run and delete all other atoms to obtain the hydrogen bond energy of a group of atoms.

Ray

FYI, the compute pair docs explicitly say that the group id is ignored, but it is not very visible. i’ll submit a patch to make it more prominent.

I would be a bit concerned with using re-run, as the removed atoms will not contribute to charge equilibration and other dynamically computed properties.
I think the only way to get exactly the desired output would be to modify the sources. however, that is very non-trivial, as the reax/c code is not “native” LAMMPS code but an independent code that was (for the most part) adapted with some wrapper functions to be run from inside of LAMMPS.

axel.

Dear Ray,

Thank you for your reply.

I need to calculate the hydrogen bond energy of a group of atoms considering the effect of the surrounding atoms, therefore I cannot delete all the other atoms.
It seems, thus I need to write my own compute to do that.

Kind regards,
Iwan

Dear Axel,

Thank you for your reply. I did not read the compute pair docs properly thus I did not notice that the group id is ignored.

Kind regards,
Iwan

FYI, the compute pair docs explicitly say that the group id is ignored, but it is not very visible. i’ll submit a patch to make it more prominent.

I would be a bit concerned with using re-run, as the removed atoms will not contribute to charge equilibration and other dynamically computed properties.
I think the only way to get exactly the desired output would be to modify the sources. however, that is very non-trivial, as the reax/c code is not “native” LAMMPS code but an independent code that was (for the most part) adapted with some wrapper functions to be run from inside of LAMMPS.

axel.

I think it is possible without changing the source code, but it is a little bit tricky – it requires rerunning the simulation with a modified potential file. Basically, you have to get rid of all energy contributions except hydrogen bonds. In order to do it, modify the potential file in the following manner: