Setting atom quaternions via Python

LAMMPS LAMMPS General Discussion
1

I am attempting to rapidly compute the energy of interaction between pairs of ellipsoidal particles in specific positions & orientations using the lammps python module (e.g. up to 10k pairs of particles and up to 1k unique sets parameters).

I initially hoped to manually set quaternions like one would set atom positions or properties but it does not appear that python can access the pointers for the quaternions based on the conversation here: Cannot extract quaternions with python - #2 by akohlmey . I did attempt to compute the quaternion and then set the compute pointers to the desired values, but it didn’t update the quaternions (unsurprisingly).
Does anyone know - Is there another backdoor way to directly set an ellipsoid’s quaternion?

As a 2nd (slightly slower) successful approach, I created two ellipsoids via in the desired positions & orientations, calculated their energy, deleted the ellipsoids, and create two new ellipsoids in the desired positions & orientations using read_data, and repeated.
Here is an example single script that shows how to do it for future reference. You can decompose the script and run the commands from python.

#2 GB Particles

dimension	3
units 		real	
boundary 	f f f

atom_style	ellipsoid
pair_style 	gayberne 1.0 1.0 1.0 20.0 #gamma, upsilon, mu, cutoff

read_data "GB_template1.data"

pair_coeff 1 1 1.0 1.0 2.0 1.0 3.0 2.0 1.0 3.0
#          i j eps sig e1a e1b e1c e2a e2b e2c
compute gbpair all pair gayberne epair
compute quat all property/atom quati quatj quatk quatw
group GB id 1:2 #necessary for the delete_atoms

timestep 1

#Run
fix 1 all nve
run 0
unfix 1

delete_atoms group GB

read_data "GB_template2.data" add append

pair_coeff 1 1 1.0 1.0 2.0 1.0 3.0 2.0 1.0 3.0
#          i j eps sig e1a e1b e1c e2a e2b e2c

timestep 1

#Run
fix 1 all nve
run 0
unfix 1

and the GB_template#.data files are:

Lammps Description

2 atoms
1 atom types
2 ellipsoids

0 40.0 xlo xhi
0 40.0 ylo yhi
0 40.0 zlo zhi

Atoms # atom-ID atom-type ellipsoidflag density x y z

1 1 1 1 0.0 0.0 0.0 
2 1 1 1 5.0 0.0 0.0

Ellipsoids # atom-ID shapex shapey shapez quatw quati quatj quatk **Diameters**

1 2.0 1.5 1.0 1.0 0.0 0.0 0.0
2 2.0 1.5 1.0 1.0 0.0 0.0 0.0
2

Changing the data in the compute is pointless, since that is a copy. The quaternion parameters cannot be changed that way, because they are not stored as per-atom parameters. The per-atom data related to this are two integers: one that is a flag indicating whether a particle has the ellipsoid data or not and a second containing an index into a separate data structure which contains the shape information (as a “bonus” struct). Thus any attempt to manipulate this data directly is doomed.

Since you have only two particles, parallel performance should not be an issue. Thus the straightforward way to set the ellipsoid shape would be through using the set command. With set atom you can set parameters for each individual particle.

Why that? Why use python at all? Or why not use the lammps.file() method?

3

Ah, thank you for the set command. It simplifies the overall script.
It would be nice if set could directly set a particle’s quaternion (not the a, b, c, theta values of the rotation, requiring getting the initial quat), but it’s a small price to pay.

Big picture, I’m doing force field optimization of a rigid body via force matching, so python is managing the optimization (shape, GayBerne parameters) and comparing the Lammps-computed energies to the DFT-calculated (BSSE-corrected) energy.

4

That would require adding a new function to the LAMMPS C-library interface and then wrap it with Python. The problem here is to implement this cleanly in C since internally.
I have to think about how to implement this, especially in Python.