Question about reax/c/species

Hello all,

I am working with ReaxFF to simulate metal/water. I have two questions about “fix reax/c/species”:

1- I use “fix reax/c/species” command to get the information about different species for a zone. As I need to have species in different zones I need to use “fix reax/c/species” command more than one time in my input script, but because of “Reuse of compute ID” error I cannot do that. Right now I have to run my simulation over and over with one “fix reax/c/species” command to get species for different zones. Could you please let me know if there is a better way to get information of species for different zones? (I attached my input script in case you wanted to take a look at it)

2- When I simulate iron/water if I need species just with H and O (like H2O, H3O, O, H and etc.) for the whole domain I could do it by defining a “water” group including H and O and use it in “fix reax/c/species” command (you could see in the attached input script). But I am not sure how I can do it for a special zone. For example in the second attachment I would like to have species with just O and H in the black rectangle, and not O, H and Fe. Do you have any idea about the way I can do it.

Thank you in advance for your time and consideration,

Regards,
Hossein,

in.Ferun (2.18 KB)

Iron-Water_simulation.docx (232 KB)

Hi Hossein,

Replies below, thanks.

Ray

Hi Ray,

Thank you for your reply.

I tried to use “position” keyword in “reax/c/species” command but I get a blank file. I could not find the reason. Could you please let me know what the reason could be. It should be noted that I get the species output file but “HOFe.pos” is blank. The input file has been attached.

Thanks
Hossein,

in.Ferun (2.21 KB)

There is no indication that you used the position keyword anywhere in your input script.

Ray

There is actually “position” keyword at the end of “fix reax/c/species” command as follows in the input script I attached (based on LAMMPS manual):

fix 4 all reax/c/species 1 1 1000 species1.out element H O Fe position 1000 HOFe.pos

Isn’t it the correct way to use it?

Thanks,

Hossein

Okay, I was looking at the first input script you attached. If I add the position keyword to the LAMMPS/examples/reax/in.reaxc.tatb script (fix 3 all reax/c/species 1 5 5 species.tatb position 25 pos.out), it works as advertised on the doc page.

Which LAMMPS version are you using? Does adding the position keyword to the reaxc/in.reaxc.tatb example work for you? When the filepos output was blank, did you start your run with the structure file or the restart file?

Ray

Thank you Ray.

I realized that the “position” part of “reax/c/species” command does not work with “minimize” command. I checked it without minimize command and it worked.

I am working with the LAMMPS version of 25 Sep 2015, as I changed the code for my research. So my question is: does “position” work with minimize command in the newer version? and if not how could I use it with minimize command? because I need the minimization, too.

Another question that I would like to know is that: what is the reason that the “species” part works with minimize command but the position part does not work?

Thank you again,

Hossein,

1- I use “fix reax/c/species” command to get the information about different species for a zone. As I need to have species in different zones I need to use “fix reax/c/species” command more than one time in my input script, but because of “Reuse of compute ID” error I cannot do that. Right now I have to run my simulation over and over with one “fix reax/c/species” command to get species for different zones. Could you please let me know if there is a better way to get information of species for different zones? (I attached my input script in case you wanted to take a look at it)

2- When I simulate iron/water if I need species just with H and O (like H2O, H3O, O, H and etc.) for the whole domain I could do it by defining a “water” group including H and O and use it in “fix reax/c/species” command (you could see in the attached input script). But I am not sure how I can do it for a special zone. For example in the second attachment I would like to have species with just O and H in the black rectangle, and not O, H and Fe. Do you have any idea about the way I can do it.

— Just bit different approach, but how about bond connection table? I’m not sure what is your goal of your analysis, but if your goal is track the number of particular type of product, than connection table could help. If you want to limit the analysis for specific region in your system, you could write some script to use connection table and xyz trejectory together to get what you want. This way, you could solve question 1 and 2 all together.

Hi Jejoon,

Please start a new thread for a new set of questions – though yours are similar to Hossein’s.

Thanks,

Ray

Hi Hossein,

The entire fix reax/c/species is not invoked during minimization, not just the position keyword. It does not in the previous and it does not in the current version.

25Sep2015 is a very outdated version – run the examples/reax/in.reaxc.tatb example and see if you get postion output. Also, what part of LAMMPS did you change? If you changed your source code and got a different behavior, then it is very difficult for me to assist you.

Thanks,

Ray

It is not a question, it is my suggestion for Hossein. With the connection table and xyz trajectory, one could track and count the specific type of product from simulation, in the specific region of the simulation box. But this depends on what is the goal of the analysis, and what is the target molecule / species.

Thank you Ray.

I run the “examples/reax/in.reaxc.tatb” example and I could get the position output, using 25Sep2015 version which I am working with. I changed “fix_qeq_reax.cpp” and “reaxc_nonbonded.cpp”. But I will compile and work with new version soon.

Now the question is: if I need to use minimization, how could I get the positions? One method that I just checked is: I run my simulation to minimize it for a couple of thousands steps and then I used restart file without minimize command and I could get positions using restart file and reax/c/species command. Is it a correct way to do? or do you have any other recommendations?

Thank you again,

Hossein,

Thank you JeJoon.

This is exactly my goal you just mentioned: “track and count the specific type of product from simulation, in the specific region of the simulation box.”

Could you please let me know in a bit more details on the method you said (using connection table and xyz trajectory) as I haven’t done it with “bond” and “trajectory” files of ReaxFF yet. I also could not find much information about these two files in the manual.

In more details: I have modeled metal/water using ReaxFF. My goal is to count number of H, O, OH, H2O and H3O species near the interface (let’s say about 2 Angstrom) in different times.

Thank you again and regards,

Hossein,

For the bond connection table:
http://lammps.sandia.gov/doc/fix_reax_bonds.html

And you could print out xyz trajectory or xyz coordinate by dump command:
http://lammps.sandia.gov/doc/dump.html

It doesn’t matter which dump style to use (custom style dump or xyz style dump or else), but you can print out xyz coordinate of all atoms.
I usually use xyz format because it is the simplest one for coordinate data.
Oh, and don’t forget to set the same print out frequency for both xyz coordinates and connection table.

After that, you can read both files, count the number of all specific species.
For example, if you want to count the number of water, you can count the number of O atom, which have only have two H atoms as neighbor.
If you want to count H3O, count the number of O atom, whose neighbor number is 3, and all of neighbors are H.
Larger products are little bit tricky to use this method, but it is not impossible. I used this method to count the number of n-decane and even bigger molecules.

If you want to set the specific region, you can limit the O number counting only in specific region, by setting the limit of x, y, and z coordinates.
You could choose the other way instead of x, y, and z coordinate, such as distance from specific atom type or etc…

I use Fortran to do this but you can use Matlab, Python, TCL, or any other simple script.

Thank you JeJoon, it was a big help. Just some minor questions on “bond” and “trajectory” files:

1- In the “bond” output file; what are “abo” and “nlp” parameters?
2- In the “.trj” output file (please see attached image); after atom id we could see x,y and z coordinate of the atom. What are another 3 quantities after that (before the last number which is the charge of the atom)?

Best,

Hossein,

trj_file.jpg

Thank you JeJoon, it was a big help. Just some minor questions on “bond” and “trajectory” files:

1- In the “bond” output file; what are “abo” and “nlp” parameters?

Averaged bond order and number of lone pairs. You basically don’t need those but just the bo between the main atom and its neighbors. You can choose a bond order cutoff for determining a bond.

2- In the “.trj” output file (please see attached image); after atom id we could see x,y and z coordinate of the atom. What are another 3 quantities after that (before the last number which is the charge of the atom)?

These are atomic forces or velocities, depending on your control file settings:

atom_info 1 ! 0: no atom info, 1: print basic atom info in the trajectory file

atom_forces 1 ! 0: basic atom format, 1: print force on each atom in the trajectory file

atom_velocities 0 ! 0: basic atom format, 1: print the velocity of each atom in the trajectory file

For example, the above settings they would be atomic forces.

Ray

Thank you Ray and JeJoon, big help.

Ray,

I hope it is my last question in this thread (I asked it before but I think it was missed):

I run the “examples/reax/in.reaxc.tatb” example and I could get species and position outputs, using 25Sep2015 version which I am working with. I changed “fix_qeq_reax.cpp” and “reaxc_nonbonded.cpp”. But I will compile and work with new version soon.

Now the question is: If I need to use minimization, how could I get the species and positions? One method that I just checked is: I run my simulation to minimize it for a couple of thousands steps and then I used restart file without minimize command and I could get positions using restart file and reax/c/species command. Is it a correct way to do? or do you have any other recommendations?

Thanks again and regards,

Hossein,

Yes, that seems like a valid approach though I would automate it with a shell script.

One recommendation is that: it doesn’t seem to make a lot of sense to me to identify molecular species during minimization, not to mention that structure that needs to be minimized “for a couple of thousand steps” is already suspicious. What you want to do is to compare the initial state and the final, minimized state, not the intermediate states along the minimization paths. The whole purpose of “fix reax/c/species” is to identify molecular species and chemical reactions for real thermodynamic states (different T, P, time, etc), not along hypothetical, minimization paths. Just a recommendation more on the philosophical side…

Ray