Cartesian.get_zmat() fails

pjknowles commented 1 month ago

Code Sample, a copy-pastable example if possible

import chemcoord
print(chemcoord.__version__)
chemcoord.show_versions()
xyzfile='chemcoord-test.xyz'
with open(xyzfile,'w') as f:
    f.write("""
    4

    C 0 0 0
    O 0 0 -1
    H 0 1 1
    H 0 -1 1
    """)
xyz=chemcoord.Cartesian.read_xyz(xyzfile)
xyz.get_zmat()

Problem description

It crashes while trying to take the square root of a numpy array:

2.1.2

INSTALLED VERSIONS
------------------
python: 3.12.4.final.0
python-bits: 64
OS: Linux
OS-release: 6.2.0-39-generic
machine: x86_64
processor: x86_64
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: en_GB.UTF-8

chemcoord: 2.1.2
numpy: 1.26.4
pandas: 2.2.2
numba: 0.60.0
sortedcontainers: 2.4.0
sympy: 1.12.1
pytest: None
pip: 24.0
setuptools: 70.1.1
IPython: 8.26.0
sphinx: None
AttributeError: 'float' object has no attribute 'sqrt'

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/peterk/Google Drive/Research/Optg/chemcoord-test.py", line 15, in <module>
    xyz.get_zmat()
  File "/home/peterk/miniconda3/envs/pymolpro/lib/python3.12/site-packages/chemcoord/cartesian_coordinates/_cartesian_class_get_zmat.py", line 608, in get_zmat
    c_table = self.get_construction_table(use_lookup=use_lookup)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peterk/miniconda3/envs/pymolpro/lib/python3.12/site-packages/chemcoord/cartesian_coordinates/_cartesian_class_get_zmat.py", line 263, in get_construction_table
    full_table = fragment._get_frag_constr_table(use_lookup=use_lookup)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peterk/miniconda3/envs/pymolpro/lib/python3.12/site-packages/chemcoord/cartesian_coordinates/_cartesian_class_get_zmat.py", line 95, in _get_frag_constr_table
    molecule = self.get_distance_to(self.get_centroid())
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peterk/miniconda3/envs/pymolpro/lib/python3.12/site-packages/chemcoord/cartesian_coordinates/_cartesian_class_core.py", line 1127, in get_distance_to
    new['distance'] = norm((new - origin).loc[:, ['x', 'y', 'z']],
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/peterk/miniconda3/envs/pymolpro/lib/python3.12/site-packages/numpy/linalg/linalg.py", line 2583, in norm
    return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: loop of ufunc does not support argument 0 of type float which has no callable sqrt method

Expected Output

Output of `cc.show_versions()`

# Paste the output of cc.show_versions() here INSTALLED VERSIONS ------------------ python: 3.12.4.final.0 python-bits: 64 OS: Darwin OS-release: 23.5.0 machine: arm64 processor: arm LC_ALL: None LANG: None LOCALE: en_GB.UTF-8 chemcoord: 2.1.2 numpy: 1.26.4 pandas: 2.2.2 numba: 0.60.0 sortedcontainers: 2.4.0 sympy: 1.12.1 pytest: None pip: 24.0 setuptools: 70.1.1 IPython: 8.26.0 sphinx: None

pjknowles commented 1 month ago

Further investigation shows that the problem goes away if all of the coordinates in the xyz file have an explicit . rather than being integers. I guess it is debatable whether or not xyz files that do not conform to this should be supported.

mcocdawc commented 1 month ago

Thank you for the bug report and the already found solution.

I am pretty sure that this would not have been a problem some time ago. Overall numpy and pandas got much more strict about types over time.

I guess it is debatable whether or not xyz files that do not conform to this should be supported.

I always write my own xyz files with explicit floats, but the library should be either permissive or give better error messages. I think I will add an explicit cast to float.

PS:

all of the coordinates in the xyz file have an explicit .

Probably it is already sufficient if one entry per column has the explicit dot, isn't it?

pjknowles commented 1 month ago

Yes. It's pretty unlikely that you will get an xyz with no '.' (except hand written like mine) but it could happen if written from python str().

mcocdawc / chemcoord