Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BUG] cursor.rotate should rotate all vector properties #103

Closed
lohedges opened this issue Sep 7, 2023 · 8 comments · Fixed by #120
Closed

[BUG] cursor.rotate should rotate all vector properties #103

lohedges opened this issue Sep 7, 2023 · 8 comments · Fixed by #120
Assignees
Labels
bug Something isn't working
Milestone

Comments

@lohedges
Copy link
Contributor

lohedges commented Sep 7, 2023

Currently the cursor.rotate function only rotates the coordinates property, when it should really rotate all vector properties, i.e. velocity too. At present this can be achieved by calling the function multiple times, using the map option to remap coordinates to velocities, e.g.:

# Here m is a rotation matrix.

# First the coordinates.
cursor = cursor.rotate(matrix=m)

# Now the velocities.
cursor = cursor.rotate(matrix=m, map={"coordinates": "velocity"})
@lohedges lohedges added the bug Something isn't working label Sep 7, 2023
@lohedges
Copy link
Contributor Author

lohedges commented Sep 7, 2023

Actually, the above approach doesn't work because, unless the center keyword argument is used, the coordinates property is used to compute the center-of-mass to use for the rotation, e.g. here.

@lohedges
Copy link
Contributor Author

lohedges commented Sep 7, 2023

This is relevant to this feature branch that I'm working on, which naively uses the workaround in the original post. For now I can just compute the center of mass before calling the function and pass it in.

@lohedges
Copy link
Contributor Author

lohedges commented Sep 7, 2023

Which also makes me realise that I'll need to be careful to choose the correct center when applying the rotation, since this won't be (necessarily) the same as what is used for the box vectors themselves which would be rotated about (0, 0, 0).

@chryswoods chryswoods added this to the 2023.4.0 milestone Sep 11, 2023
@lohedges
Copy link
Contributor Author

I think I've now got this working (using a workaround) on my feature_reduce branch. I've also added a unit test to check that the rotation is performed correctly and that the coordinates property is updated.

The two things that I'm not sure of are:

  • For the vector properties, I aim performing the rotation about the box origin, which can be specified by the user (defaulting to (0,0,0)). It would be nice to be able to determine this automatically.
  • Following the rotation of the box, the z component of the third lattice vector needs to be flipped if it is negative. I'm not sure how this should impact the properties.

@chryswoods chryswoods modified the milestones: 2023.4.0, 2023.5.0 Oct 7, 2023
@chryswoods
Copy link
Contributor

I've been thinking about this - am I right that the centre of rotation for each velocity vector would be the origin (0,0,0)? This is because velocities are all vectors that represent change in coordinates - their origin represents the "no change - atom not moving", and so they all have an origin at (0,0,0) regardless of where the atom is located or where the centre of the box sits (effectively a velocity vector is always in a frame of reference centered on its atom).

So, if an atom's velocity is (1,0,0), then rotating the atom by 90 degrees about z would change the velocity to (0,1,0), i.e. rotated the velocity vector about its origin. And this is the same regardless of the position of the atom or the choice of centre of box.

@chryswoods chryswoods self-assigned this Oct 25, 2023
@chryswoods
Copy link
Contributor

I'm also not sure how the above fix worked, because velocity should be of type AtomVelocities, which cannot be rotated by that function (it only works on AtomCoords).

For example;

>>> import sire as sr
>>> # Creating some velocities by running MD
>>> mols = sr.load(sr.expand(sr.tutorial_url, "ala.top", "ala.traj"), silent=True)
>>> mols = mols.minimisation().run().commit()
>>> d = mols.dynamics(timestep="4fs", temperature="25oC")
>>> d.run("10ps")
>>> mols = d.commit()
>>> mols[0].property("velocity")
SireMol::AtomVelocities( size=22
0: ( -14.35 Å ps-1, -13.7221 Å ps-1, -3.16664 Å ps-1 )
1: ( 3.83649 Å ps-1, 0.782261 Å ps-1, -1.09576 Å ps-1 )
2: ( -1.63303 Å ps-1, 30.2243 Å ps-1, 19.0436 Å ps-1 )
3: ( -6.70623 Å ps-1, -15.6312 Å ps-1, -14.4375 Å ps-1 )
4: ( -1.49417 Å ps-1, 5.25292 Å ps-1, 0.825426 Å ps-1 )
...
17: ( 11.5151 Å ps-1, 26.6377 Å ps-1, -27.9718 Å ps-1 )
18: ( 0.576416 Å ps-1, 1.62035 Å ps-1, 2.71101 Å ps-1 )
19: ( -3.38965 Å ps-1, 7.1897 Å ps-1, 8.87909 Å ps-1 )
20: ( -8.95596 Å ps-1, -22.0357 Å ps-1, -5.77348 Å ps-1 )
21: ( 0.297291 Å ps-1, 0.673374 Å ps-1, 19.2852 Å ps-1 )
)
>>> c = mols[0].cursor()
>>> c.rotate(5*sr.units.degrees, map={"coordinates":"velocity"})
TypeError: SireError::invalid_cast: Cannot cast from an object of class "SireMol::AtomVelocities" to an object of class "SireMol::AtomCoords". (call 
sire.error.get_last_error_details() for more info)

@lohedges
Copy link
Contributor Author

Yes, I think (0, 0, 0) is correct and is what my current workaround uses.

@chryswoods
Copy link
Contributor

This commit is the beginning of adding in the option of rotating velocities.

@chryswoods chryswoods mentioned this issue Oct 25, 2023
@chryswoods chryswoods linked a pull request Oct 25, 2023 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants