DifferentPrimitivesException in some COS calculations

[RaphaelRobidas]

Describe the bug
A small, but constant portion of growing string calculations I launch end up crashing with this exception:

Traceback (most recent call last):
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/intcoords/helpers.py", line 36, in get_tangent
    tangent = prims2 - prims1
              ~~~~~~~^~~~~~~~
ValueError: operands could not be broadcast together with shapes (341,) (340,) 

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/robidasr/goldenv_3.11/bin/pysis", line 8, in <module>
    sys.exit(run())
             ^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/run.py", line 2054, in run
    run_result = run_from_dict(run_dict, **run_kwargs)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/run.py", line 1993, in run_from_dict
    run_result = main(run_dict, restart, cwd, scheduler)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/run.py", line 1516, in main
    opt_result = run_opt(
                 ^^^^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/drivers/opt.py", line 151, in run_opt
    opt.run()
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/optimizers/Optimizer.py", line 928, in run
    reparametrized = self.geometry.reparametrize()
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/cos/GrowingString.py", line 474, in reparametrize
    self.reparam_dlc(desired_param_density, thresh=self.reparam_tol)
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/cos/GrowingString.py", line 299, in reparam_dlc
    cur_param_density = self.get_cur_param_density()
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/cos/GrowingString.py", line 88, in get_cur_param_density
    diffs = [
            ^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/cos/GrowingString.py", line 89, in <listcomp>
    image - self.images[max(i - 1, 0)] for i, image in enumerate(self.images)
    ~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/Geometry.py", line 360, in __sub__
    diff = -get_tangent(
            ^^^^^^^^^^^^
  File "/home/robidasr/goldenv_3.11/lib/python3.11/site-packages/pysisyphus/intcoords/helpers.py", line 38, in get_tangent
    raise DifferentPrimitivesException
pysisyphus.intcoords.exceptions.DifferentPrimitivesException

This is very similar to a bug I previously reported (#218)

To Reproduce
The crash occurs repeatably using these input files:
input_files.zip

Expected behavior
I would expect either the calculation to work or Pysisyphus to handle this exception and stop gracefully.

OS and Python:

• OS: Debian trixie
• Python version: 3.11.6

Pysisyphus version
This kind of crash has been occurring randomly for at least two months (although most calculations worked fine) and I confirmed that the behavior is the same with the latest commit on master (5950976)

[eljost]

Thanks for the report. This path looks wild. The issue is that at one point one internal coordinate gets invalid ((<PrimTypes.PROPER_DIHEDRAL: 8>, 15, 13, 12, 7) or if you open cycle_023.trj in jmol you can issue measure 16 14 13 8 to see the coordinate). The reason is the flipping methyl group and one hydrogen that was formerly a terminal atom in a dihedral now moves inbetween the oxygen and the carbon, leading to a collinear arrangement. Following this, pysisyphus rebuilds the internal coordinates, but does it in an incosistent way. The set of internal coordinates used to rebuild them on all images is probably not valid for all images, leading to different number of coordinates.

Locally, I just tried a patch that makes it work, but this is not a comprehensive solution and I'll have to take a closer look at this before I can push the patch.

diff --git a/pysisyphus/cos/GrowingString.py b/pysisyphus/cos/GrowingString.py
index 41829acf..ead07439 100644
--- a/pysisyphus/cos/GrowingString.py
+++ b/pysisyphus/cos/GrowingString.py
@@ -296,7 +296,11 @@ class GrowingString(GrowingChainOfStates):
         # is above or below the desired param_density on the normalized arc.
         #
         # The reparametrization is done in micro cycles, until it is converged.
-        cur_param_density = self.get_cur_param_density()
+        try:
+            cur_param_density = self.get_cur_param_density()
+        except DifferentPrimitivesException:
+            self.reset_geometries(self.images[0])
+            cur_param_density = self.get_cur_param_density()

For now this seems to converge reasonably using cart coordinates, at least on the current dev branch:

geom:
 type: cart
 fn: [start.xyz, end.xyz]
calc:
 type: xtb
 pal: 1
 charge: 1
 mult: 1
 alpb: ch2cl2
 acc: 0.1
preopt:
  max_cycles: 5
cos:
 type: gs
 max_nodes: 21
 climb: true
 reparam_every: 2
 reparam_every_full: 1
opt:
 type: string
 stop_in_when_full: 50
 align: true  # note that aling is enabled for 'cart' coords
 coord_diff_thresh: 0.005

The Cartesian GSM stops after 50 cycles after growing fully

       -----------------------------------------------------------------------------------
         60    -0.005100*    0.005560     0.000474*    0.022240*    0.002228*        4.38 
        String= Full HEI=11/23 (E_hei-E_0)= 344.1 kJ/mol norm(forces_true,hei)=0.009951 E_h/a_0
       Converged!

The patched GSM in DLC also stops after 50 cycles after being fully grown:

      -----------------------------------------------------------------------------------          
         60    -0.003524*    0.022589     0.002031     0.054787*    0.005184*        4.69           
        String= Full HEI=10/23 (E_hei-E_0)= 318.0 kJ/mol norm(forces_true,hei)=0.071273 E_h/a_0     
       Converged!  

[RaphaelRobidas]

Thanks for your very quick reply and patch! Indeed, the path is pretty odd and I didn't necessarily expect Pysisyphus to nicely converge to a transition state, but the exception seemed unexpected. The patch indeed solves the issue, thank you very much!

[eljost]

Yeah, it solves this issue but the patch relies on the fact that self.images[0] contains a set of internal coordinates that are valid for all images, which may not be the case.
I'll try to provide a more general patch in the near future.

[runtian97]

Hi, is this bug fixed and the new code is pushed?