mcocdawc / chemcoord

A python module for manipulating cartesian and internal coordinates.
GNU Lesser General Public License v3.0
72 stars 19 forks source link

Allow specifying first three atom positions for give_cartesian() #14

Closed sgill2 closed 7 years ago

sgill2 commented 7 years ago

When using give_cartesian() (previously to_xyz()) the function takes the convention of setting the first atom in the build list to be centered at the origin. In the cases where you are converting back and forth between xyz and zmat representations though, it might be helpful to keep the xyz coordinates more consistent with the original coordinates, which may not be centered at the origin. You can effectively do this by specifying the coordinates of the first three atoms.

I'd propose to add an additional optional argument to give_carrtesian() that takes up to a 3x3 numpy array that corresponds to the x,y,z coordinates of the first three build list atoms. If this argument is specified then the add_first/second/third_atom() would draw from these coordinates instead. (Of course, if these atom positions don't satisfy the zmat bonds/angles then the resulting representation will be messed up, so checking the bonds and angles given by the first three atoms and throwing an warning or error if they are not consistent would probably be good too.)

I've written up some code that mostly does this already (in the old chemcoord framework at least), and would be happy to share it if you would find it helpful.

mcocdawc commented 7 years ago

Thank you for your interest! And for writing in order to avoid double work. It seems as if we stumbled over the same problem.

Short answer: I am implementing it at the moment.

Long answer: I don't think that a manual insertion of first three atoms into cartesian space is the best option. Chemically speaking translational and rotational degrees were ommitted before. (in total six) If you specify now three locations, you fix nine degrees of freedom. Somehow the nine degrees of freedom have to be coupled to reduce to only six and I think this will usually bite you. Because even if they are consistently inserted and don't mess up the resulting representation, they could still throw an exception because of floating point noise. And it would become quite involved to consistently check for this, because the accumulation of errors due to float noise would be distance dependent...

I would prefer two options, which are:

  1. Take implicitly the canonical unitvectors as reference for the first three atoms. Fill the "upper right triangle" of the zmatrix with bond, angle and dihedral in reference to the canonical unitvectors. (Thats also how they did it in this paper: http://www.tandfonline.com/doi/abs/10.1080/08927020600900337)

  2. Specify the location of the barycenter of the molecule in the first row using spherical coordinates. Use the remaining three angles as euler angles to specify how the axis of inertia are rotated. This would allow it, to rotate a molecule along principal axes of inertia by just assigning an angle.

The problem with option two is, that the rotation of the molecule is still undefined, if you have degenerate inertia moments. And this is not rare. At least as far as I thought about it. For this reason I am implementing at the moment option 1.

sgill2 commented 7 years ago

Yeah, that definitely sounds like a better approach. Thanks for the response and the reference!

mcocdawc commented 7 years ago

@sgill2 You can install the new master branch and try the new functionality. I am happy about feedback