Feature request: reorder atoms

[Jellby]

Maybe this is already possible, but today I needed to re-create my RasOrb file because i had to change the order of the atoms in the molecule (so that the geometry can be interpolated with another conformer), and I thought that would be a useful feature in penny.

(Also, in molcas it should be in Tools, if anywhere, not in sbin, don't you think?)

[steabert]

Thanks for the suggestion, this is indeed something that could be useful and will be easy with penny. The idea is to be able to easily reorder orbitals any way you like. Right now I'm refactoring the entire code to make everything more structured.

The release in Molcas is a single large script that won't be further updated, I decided to shift to a separate development process. So, in the future, any new additions will only appear here. When I think the github version is ready to replace it, I will remove the script from the Molcas source.

[steabert]

When needing to re-order atoms, it would only work with HDF5, because the RasOrb file itself doesn't carry basis set information necessary for the reordering. Would that be enough?

[Jellby]

Can't you just reorder the MO coefficients in the RasOrb file?

[steabert]

No, how do I know which coefficients match which center?

[Jellby]

Doh! Very right... you could take the information from the molden file... But well, I guess the HDF5 file is enough, that should be what we use in the future.

[steabert]

How would one specify the new ordering? I'm thinking about an option to read the new ordering from an xyz file?

[Jellby]

An xyz file may not be enough, since xyz files may have only the atom symbols, and no indication of which carbon is which, for instance. I'd say the easiest is to provide an index list, like "1 3 2 4" (meaning exchange the 2nd and 3rd atoms). For bonus points there could be an option to make the reverse reordering (treat "1 4 2 3" as "1 3 4 2"), or use two xyz files as input and reorder atoms matching whatever reordering was done in the xyz files.

Now, something possibly more complicated would be modifying the orbitals to account for rotations... 😉

[steabert]

It can be done by sorting the coordinates, that way you get a mapping of identical positions between the provided xyz file and the coordinates stored in HDF5. This is similar to your last suggestion. So no big problems there, this is the solution I would go for. An index list can be added too, but that is more work I think, I mean, for the person using it ;)

What do you mean with accounting for rotations? If you just rotate the molecule, then I would think the orbitals keep the same composition. Or ...?

[Jellby]

Except for the angular components. If you rotate 90 degrees in the x-y plane, the coefficient for px should now be for py...

[steabert]

Oops, yes of course :) On the other hand, couldn't Molcas reorient the molecule?

[Jellby]

What do you mean? Reorient the coordinates to match the orbitals (or HDF5 file)? Sometimes you don't want that. My use case is the following.

1. Do a calculation of your problem molecule to get the orbitals
2. Get a similar structure from elsewhere (a different calculation, a published structure, a QM/MM system).
3. Try to use the orbitals in (1) for the structure in (2).
4. Oops! The atom order is different, the orientation of the molecule is different (even though the internal structure is almost the same).

Two possible fixes without re-computing initial orbitals: reorder and reorient the structure in (2) to match the orbitals in (1); or reorder and reorient the orbitals in (1) to match the structure in (2).

The first solution is usually easy, unless the structure is part of a larger system and a larger calculation (like QM/MM dynamics, where you want to communicate between different programs). The second could in theory be done with molpy, which is why this is a feature request :)

[steabert]

It seems that for either 1 or 2, one would need to first write a code that figures out molecular overlap. If you have to do that anyway, then it is much easier to rotate/reorder the xyz file than it is to transform the orbitals...

[Jellby]

Yes, the code for overlapping two structures is not very complicated. Rotating the whole system in (2) is sometimes such a hassle that one prefers re-computing the orbitals (been there, done that), for instance, if the system is periodic, it is not rotation-invariant, so the periodic conditions need to be rotated as well. I don't know how easy it would be to rotate the orbitals, it seems to me it should be almost as easy once you have the rotation matrix, if you know how the different basis functions relate to each other.

If you decide to implement this, I could give you a nice reference for the structure overlap. Or maybe I'll implement it some day.