Interpolation in internal coordinates

[LeonardoBarneschi]

Dear Mr,

May I ask how I would go about a linear interpolation in internal coordinates given 2 geometries in chemcoord? I guess it is possible to generate the 2 corresponding z-matrix and interpolate externally, but I was wondering whether is there any built-in tool to do that.

Also, on the same note, I would like to ask if given two conformers of the same molecule we can enforce the 2 z-matrix to be generated using the same reference/order.

Thank you in advance, Leonardo

[ledragna]

Dear Leonardo,

I'm just an enthusiastic user of chemcoord, and probably someone is going to give you a better reply. However, I used chemcoord in the past to do exactly what you are describing and I don't know if there is now a builtin function for it.

As starting point, you can get inspiration from the tutorial jupyter notebooks (i.e. Transformation.ipynb).

For example you can get the same zmat for two structures like this

    file1xyz = cc.Cartesian.read_xyz(file1, start_index=1)
    table = file1xyz.get_construction_table()
    file1zmat = file1xyz.get_zmat(table)
    file2xyz = cc.Cartesian.read_xyz(file2, start_index=1)
    file2zmat = file2xyz.get_zmat(table)

Now, you can compute the difference between the two structures and with some unsafe_loc get a smooth transition between the two structures.

    dist = file1zmat.copy()
    dist.unsafe_loc[:, ['bond', 'angle', 'dihedral']] = 0
    dist.unsafe_loc[:, ['bond', 'angle', 'dihedral']] = file2zmat.loc[:, ['bond', 'angle', 'dihedral']] - file1zmat.loc[:, ['bond', 'angle', 'dihedral']]
    dist.unsafe_iloc[0,[2,4,6]] = 0
    dist.unsafe_iloc[1,[4,6]] = 0
    dist.unsafe_iloc[2,[6]] = 0

Maybe you can check the difference between the dihedrals and force the smaller rotation angle or force the rotation direction.

There is probably a more clean way to get the same result, but this can be a starting point.

Marco

[mcocdawc]

Dear Leonardo thank you for the question and dear Marco thank you for already giving a good answer,

There are other possibilities too though: The binary operators are also overloaded for Z-matrices. This means you can in principle write something like

STEPS = 10
[zm1 + i * (zm2 - zm1) / (STEPS - 1) for i in range(STEPS)]

to obtain a list of interpolated zmatrices between zm1 and zm2. There is a big caveat though. Usually the difference between two Z-matrices (zm2 - zm1) is not a valid¹ Z-matrix anymore.² Also the whole expression i * (zm2 - zm1) / (STEPS - 1) is usually not valid, only after adding this difference to another Z-matrix zm1 we assume and require the result to be valid again. Similar to safe_loc I added a test to every binary operations that asserts that the result is a valid Z-matrix. Hence you usually cannot write the above list comprehension.

There is an explicit way around it though. With the context manager TestOperators you can switch the test on and off. Hence the cleanest and recommended way to write such an interpolation is:

STEPS = 10
with TestOperators(False):
    Deltas = [i * (zm2 - zm1) / (STEPS - 1) for i in range(STEPS)]
# We want the final conversion to be tested for validity
[zm1 + delta for delta in Deltas]

¹ Valid defined as: "can be converted to cartesian coordinates. ² This is in contrast to cartesian coordinates where (m2 - m1) could be visualized and gives you the movement of atoms directly.

[LeonardoBarneschi]

Thank you very much Oscar and Marco, works like a charm!