Dear Ray,
hi. I think there is a problem with the new fix reax/c/species command. Maybe you already know about it. It does not take into account periodic boundary conditions when using “position” keyword. I’m attaching the example to this message. In our own composition analisys tool we solved the same problem by checking if the distance between neighbors is higher than maximum possible bond distance: then they are connected through PBC, so atomic positions should be corrected for this molecule. I think the same idea may be used in the fix.
I also have another questions,
- am I right that you calculate CoM there as geometric center of the molecule?
- will the fix work if nfreq != neighbor->every, but ((nfreq % neighbor->every) == 0)? Do we have to reset neighboring settings in this case?
Sincerely yours,
Oleg Sergeev,
junior researcher, VNIIA.
data2.reax (238 Bytes)
ffield.reax.rdx (12.3 KB)
in.example (457 Bytes)
A little addition: LAMMPS version is that of 14 May.
Oleg.
15.05.2013, 16:33, "Oleg Sergeev" <[email protected]...>:
Oleg,
This is not a bug, rather it depends on the treatment of center-of-mass (COM) positions with periodic boundary conditions (PBC).
Your test case is a C2 molecule across the PBC, and the fix reax/c/species has no problem identifying that. But how would you treat the COM? Shift the COM to either one of the atoms does not make sense. Other LAMMPS COM functions estimate COM in the same fashion.
For example, adding this “variable xcom equal xcm(all,x)” to your input and use thermo_style to output this variable, you will see the COM of the C2 molecule is also at the center - same as reax/c/species.
Thanks for testing this out.
Ray
Ray,
15.05.2013, 19:24, "Ray Shan" <[email protected]...>:
Oleg,
This is not a bug, rather it depends on the treatment of center-of-mass (COM) positions with periodic boundary conditions (PBC).
Your test case is a C2 molecule across the PBC, and the fix reax/c/species has no problem identifying that. But how would you treat the COM? Shift the COM to either one of the atoms does not make sense. Other LAMMPS COM functions estimate COM in the same fashion.
Why shouldn't the COM be at 0.0 (or 9.5 at some timesteps) in my example? I think it would be no problem if it sometimes jump from one boundary to another, same way as atom coordinates under PBC. The problem with present implementation is that if you have a system of homogeneously distributed large molecules that span across PBC, you will not obtain homogeneous distribution of their concentration with fix reax/c/species, but dense region near the cell center. On the other hand, if some of the molecules would "be" near x_lo and some near x_hi, you'd observe the same concentration across all simulation cell.
For example, adding this "variable xcom equal xcm(all,x)" to your input and use thermo_style to output this variable, you will see the COM of the C2 molecule is also at the center - same as reax/c/species.
May I ask you to briefly comment two other questions?
Regards,
Oleg.
Ray,
15.05.2013, 19:24, "Ray Shan" <[email protected]...>:
Oleg,
This is not a bug, rather it depends on the treatment of center-of-mass (COM) positions with periodic boundary conditions (PBC).
Your test case is a C2 molecule across the PBC, and the fix reax/c/species has no problem identifying that. But how would you treat the COM? Shift the COM to either one of the atoms does not make sense. Other LAMMPS COM functions estimate COM in the same fashion.
Why shouldn't the COM be at 0.0 (or 9.5 at some timesteps) in my example? I think it would be no problem if it sometimes jump from one boundary to another, same way as atom coordinates under PBC. The problem with present implementation is that if you have a system of homogeneously distributed large molecules that span across PBC, you will not obtain homogeneous distribution of their concentration with fix reax/c/species, but dense region near the cell center. On the other hand, if some of the molecules would "be" near x_lo and some near x_hi, you'd observe the same concentration across all simulation cell.
oleg,
the problem in your example is that you don't have the proper image
flags. if you had *real* molecules with a bond in the topology, lammps
would complain that your bond is far too long. the group->xcm() method
considers image flags, so if a molecule is very long (up to the point
of crossing the PBC multiple times) its COM will be computed
correctly.
axel.
I suppose you can; adding something along the lines of using keyword “shift_pos yes”. But currently it is treated in the same way as group->xcm(). Please feel free to modify it and contribute if you would like - I’d be happy to acknowledge you.
Replies to your Qs:
1). Yes.
2). The fix resets neighboring criteria automatically and writes output when “timestep % nfreq == 0”
Ray
Does it work for reaxff too? It does not have predefined topologies. I'll check it out, thanks. Though I don't yet understand why this C2 is not a real molecule.
Oleg.
15.05.13, 19:57, "Axel Kohlmeyer" <[email protected]>":
In atomic style potentials, such as ReaxFF, there are only atoms and no bonds at all. LAMMPS does not consider it to be bonded hence not a molecule - no bonds, no true LAMMPS molecules. Only ReaxFF thinks it is a molecule because there is a strong bond order between the two.
Ray
Does it work for reaxff too?
image flags work everywhere. since you don't use a bond topology,
lammps cannot check for inconsistencies, that is all. even for
conventional force fields that require a bond topology, you can ignore
the warning about inconsistent image flags and lammps will do most of
the time the right thing. only when you use a feature that *requires*
consistent image flags (e.g. when using replicate or computing the
com) lammps will produce the wrong answer.
axel.
15.05.13, 20:23, "Ray Shan" <[email protected]...>":
In atomic style potentials, such as ReaxFF, there are only atoms and no bonds at all. LAMMPS does not consider it to be bonded hence not a molecule - no bonds, no true LAMMPS molecules. Only ReaxFF thinks it is a molecule because there is a strong bond order between the two.
Ray
Thanks, I thought so.
I'll see if I can modify the fix. Regarding neighbor resetting, can't it be a problem if species output is done rarely, say, each 10000th step? Building neighbor lists this rare is bad for reactive systems. That was the reason of my second question.
Oleg.
It is not a good idea to do time-averaging on reactive systems with timescale that span this long, so it was designed to, and thus only appropriate, to be used frequently, e.g. “1 100 100”.
Ray
Sorry for many questions, but would "1 100 10000" be incorrect? Nfreq is 10000 in this case, but time averaging is only done for 100 steps of a 10000-step period.
Oleg.
15.05.13, 21:12, "Ray Shan" <[email protected]...>":
This setting is correct, but may not work well with the implementation.
I see. The system is intended to be reactive, so predefined bonding cannot be used.
Oleg.
15.05.13, 20:39, "Axel Kohlmeyer" <[email protected]>":
Thanks for the answers 
Oleg.
15.05.13, 21:20, "Ray Shan" <[email protected]...>":
I see. The system is intended to be reactive, so predefined bonding cannot be used.
sure, but when your initial configuration has inconsistent image
flags, they will *stay* inconsistent and the COM will *always* be
wrong. if you have a reactive system and reactions do happen and stuff
diffuses in or away, then it may take quite a while until things
*perhaps* become inconsistent. in other words, there are no
guarantees, but you are already starting on the wrong foot.
axel.