I use LAMMPS with reax/c pair-style to simulate thermal decomposition of
various materials, and currently trying to characterize the products in my
system.
Particularly, I would like to calculate the radius of gyration of each
identified molecule so to have an estimate of the size of the molecules and
agglomerates at each time step.
Doing it with 'compute gyration' isn't easy since there are no molecules in
reaxff, so i don't have the group_id's to use with that fix.
compute gyration or compute gyration/molecule? Because compute
gyration does not require molecule IDs but group IDs. You can define
the group IDs how ever you want to, and it certainly works with
reax/c. It just does not output the gyration for each identified
molecule with fix reax/c/species. "compute gyration/molecule" require
molecule IDs so that it does not work with pair_style reax/c.
So, I think that leaves me with the option of editing the source code where
molecule identification takes place (fix species / reaxc_bonds) and maybe
add a subroutine
To modify the source code, you just add the formula defined on the
compute gyration doc page to somewhere after the FindMolecule()
function in fix_reaxc_species.cpp by looping through Nmole (number of
molecules and each molecule has a sequenced molecule ID).
Besides modifying the source code, you can also dump the per-atom
molecule ID from fix reax/c/species to a per-atom dump. Eg,
fix 1 all reax/c/species 1 100 100 species.out
dump 1 all custom 100 *.dump id type x y z f_1
This way you have molecule ID in the dump file. Then write a simple
converter to convert the dump file to a atom_style molecule data file
to be read in using read_data. Then you can use compute
gyration/molecule and perhaps "rerun" or "run 0" to obtain the
gyration of each molecules.
Dear Ray,
Thank you for your kind response and guidance.
I tried to dump the molecule ID of that fix, however it seems to me that it
doesn't support outputting per-atom quantities.
I got the error:
ERROR: Dump custom fix does not compute per-atom info
(../dump_custom.cpp:1193)
And reading the fix reax/c/species command manual it says:
"No global or per-atom quantities are stored by this fix for access by
various output commands"
Does it mean that the source code of this fix should be modified so that it
stores the molecule id's ?
I picked what seemed as the easiest method you offered, which is to dump the
per-atom molecule ID from fix/reaxc/species, but
i didn't get the molecule ID in the dump file, but rather a column of
floating point numbers (positive and negative), and
I'm not sure what they mean. They obviously not the mol id's.
I used the following 2 commands:
fix bonds all reax/c/species 1 100 500 species.out cutoff 1 4 0.3 cutoff 1 1
0.55 cutoff 1 3 0.65 cutoff 1 2 0.4 cutoff 3 3 0.65 cutoff 4 3 0.4 cutoff 3
2 0.4 cutoff 2 2 0.55 cutoff 2 4 0.55 cutoff 4 4 0.55 element C H O N
position 1000 pos.out
dump 1 all custom 5000 dump.cook_${name} id f_bonds
I used 21Nov13 version and did not see what you saw. Which version
are you using? If you are using a new version, you can attached a
simple input deck that reproduces the problems.
Fix reax/c/species only determines molecules and species at timesteps
that are multiples of Nfreq, which is 50 in your script. Fix
reax/c/species has no information on on-the-fly molecules IDs before
50 MD steps. I will have to check the source code to see what it is
outputting to the global quantity before 50 steps, but it is expected
that you have to wait after the first Nfreq steps to access this
information. As you can you in your dump file, you have on-the-fly
molecule IDs after 50 steps.
Please also note that these molecule IDs are calculated on-the-fly,
due to the nature of an atomistic reactive potential. So please keep
in mind that a molecule with a certain ID can have another ID at the
next Nfreq, and a same ID number at different Nfreqs may very well
represent a totally different molecule.