[lammps-users] rigid/small quaternions - read data and restart

Dear Lammps-Users,

I have some trouble understanding how quaternions are calculated/set using rigid/small bodies.
I simulate rigid cubes in an explicit lennard jones solvent. After orienting the cubes randomly and writing the coordinates in a data file (as molecules), I add the cubes to my lj-fluid simulation by reading this data file and deleting overlapping lj particles.
The question is now: after reading the data file, how are the quaternions of the cubes defined? I would have guessed they are normalized to 0 for every cube, but they arent. Here is a snippet of the quaternions (dump.0.quat) after reading the data file at t=0

ITEM: TIMESTEP
0
ITEM: NUMBER OF ENTRIES
96
ITEM: BOX BOUNDS pp pp ff
0.0000000000000000e+00 8.0000000000000000e+01
0.0000000000000000e+00 8.0000000000000000e+01
0.0000000000000000e+00 1.6000000000000000e+02
ITEM: ENTRIES c_quat[1] c_quat[2] c_quat[3] c_quat[4] c_quat[5]
33 0 0.707107 0.707107 0
69 0 0.707107 0.707107 0
7 0.5 -0.5 -0.5 -0.5
61 0.707107 0.707107 -0 0
85 0.707107 0.707107 -0 0
73 0 0.707107 0.707107 0
92 1 0 0 0
45 0.707107 0.707107 -0 0
66 0.707107 0 -0.707107 0
51 0.996283 0 -0.0861436 0
55 0.921232 -0.381587 -0.0289496 0.0698904
22 0.704753 0.704753 0.0576413 -0.0576413

As you can see some of them are set to 0, some to 0.5, some to sqrt(2)/2 and some to something different…
A minimal input script for reading the data file:

units lj
atom_style full
atom_modify map array first cubes
dimension 3
boundary p p f
pair_style lj/cut/opt 2.5

read_restart restart.dat.10000000 # equilibrated lj fluid
group solvent type 1
read_data data.cubes group cubes add append # read data file
neighbor 0.4 bin
neigh_modify exclude molecule/intra cubes

pair_modify shift yes
pair_coeff 1 1 1.0 1.0 #solvent
pair_coeff 2 2 1.0 1.0 1.12 #np
pair_coeff 3 3 0.0 0.0 0.0 #dummy
pair_coeff 2 3 0.0 0.0 0.0 #np dummy
pair_coeff 1 2 ${sol} 1.0 #np solvent
pair_coeff 1 3 0.0 0.0 2.5 #dummy solvent

delete_atoms overlap 1.0 solvent cubes #delete overlapping solvent
pair_coeff 1 3 0.0 0.0 0.0 #dummy solvent=
comm_modify mode multi cutoff/multi 1 2.5 cutoff/multi 2 2.5 cutoff/multi 3 6.95
reset_timestep 0
fix 1 cubes rigid/nve/small molecule
fix 2 solvent nve
compute quat cubes rigid/local 1 mol quatw quati quatj quatk
dump quat cubes local 100000 dump..quat c_quat[]

thermo 1000
thermo_style custom step temp epair evdwl etotal press pe ke atoms
timestep 0.005
timer full
run 10000000

Would it be possible to use the set command or other to define just the quaternion values, without rotating the molecule?

Also after running and reading the restart from the simulation above i have a similar problem: the quaternions of the last dump before the restart (t=10000000) are not the same as the quaternions of the first dump (t=0) after the restart.
e.g.
before restart (dump.10000000.quat):
7 0.543479 0.112023 0.28107 0.782995
76 -0.0397375 -0.373768 -0.30478 -0.875116
32 0.411682 -0.406197 0.538993 0.61238
after restart (dump.0.quat)
7 0 0.5226 0.852578 0
76 0.770062 0.562899 0.177185 -0.242394
32 0.9928 0 -0.119787 0

Where does this discrepancy come from?
Of course I’ll galdly provide more scripts or data if needed.

Thanks in advance!

Best regards
Konstantin

A couple other things which may be helpful.

The set quat commands are for aspherical individual particles, like ellipsoids.
Nothing to do with rigid bodies.

I’m not clear why a restart would cause the quats to change. Nothing about the rigid
bodies is written into the restart file. You have to re-specify your fix rigid/small in the new
simulation, which then repeats what I mentioned in the previous email. I’d be curious
to know if that same behavior occurs using the infile keyword with fix rigid/small.
Which does write an auxiliary file as output with a restart. And it can be read back in
for the new restarted simulation. But it has nothing about quats in it, just space frame
quantities.

If your rigid bodies are degenerate in that a different quat could essentially represent
the same body (within round-off error) then LAMMPS may be flipping quats in that manner on a restart.
But I don’t think that should affect dynamics.

Steve

When you use the fix rigid (or rigid/small) command it identifies each rigid body and sets
its initial quaternion as well as angular momentum from the particle coords and velocities.

This is done by diagonalizing a small matrix to get the directions of the principal axes for the body
in the space frame (simulation box). Those 3 axes are converted to
an initial quaternion by the method exyz_to_q() in math_extra.h.
So no reason that quat would be zero.

Not sure if that is answering your Q.

Steve