Crashing optimization

[susilehtola]

See https://github.com/psi4/psi4/issues/2213

This happens with optking 16066d0daef18f6c33138fb6950994b689e37a51

[susilehtola]

Error is

[INFO]: Creating a Compute Wrapper
[INFO]: psi4 has been initialized
[INFO]: Running optimize(o_molsys,computer)
[INFO]:     Augmenting connectivity matrix to join fragments.
[CRITICAL]: 1
[WARNING]: This method should be adjusted after dimer fragments
[ERROR]: AlgError: Exception created.

[ERROR]: AlgError: Exception created.

[CRITICAL]:     A non-optimization-specific error has occurred.

[CRITICAL]:     Resetting all optimization options for potential queued jobs.

[ERROR]: Error Type:  <class 'optking.exceptions.AlgError'>
Traceback (most recent call last):
  File "/tmp/optking/optking/tors.py", line 155, in q
    tau = v3d.tors(geom[self.A], geom[self.B], geom[self.C], geom[self.D])
  File "/tmp/optking/optking/v3d.py", line 227, in tors
    raise AlgError(
optking.exceptions.AlgError: Interior angle of   0.1 or  93.5  can't work in good torsion.

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

Traceback (most recent call last):
  File "/tmp/optking/optking/optimize.py", line 56, in optimize
    logger.debug(str(o_molsys))
  File "/tmp/optking/optking/molsys.py", line 46, in __str__
    s += F.__str__()
  File "/tmp/optking/optking/frag.py", line 49, in __str__
    s += "\t%-18s=%17.6f%19.6f\n" % (x, x.q(self._geom), x.q_show(self._geom))
  File "/tmp/optking/optking/tors.py", line 157, in q
    raise AlgError("Tors.q: unable to compute torsion value") from error
optking.exceptions.AlgError: Tors.q: unable to compute torsion value
[ERROR]: Error caught:Tors.q: unable to compute torsion value
Traceback (most recent call last):
  File "/tmp/optking/optking/tors.py", line 155, in q
    tau = v3d.tors(geom[self.A], geom[self.B], geom[self.C], geom[self.D])
  File "/tmp/optking/optking/v3d.py", line 227, in tors
    raise AlgError(
optking.exceptions.AlgError: Interior angle of   0.1 or  93.5  can't work in good torsion.

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

Traceback (most recent call last):
  File "/tmp/optking/optking/optimize.py", line 56, in optimize
    logger.debug(str(o_molsys))
  File "/tmp/optking/optking/molsys.py", line 46, in __str__
    s += F.__str__()
  File "/tmp/optking/optking/frag.py", line 49, in __str__
    s += "\t%-18s=%17.6f%19.6f\n" % (x, x.q(self._geom), x.q_show(self._geom))
  File "/tmp/optking/optking/tors.py", line 157, in q
    raise AlgError("Tors.q: unable to compute torsion value") from error
optking.exceptions.AlgError: Tors.q: unable to compute torsion value

[psi-rking]

https://github.com/psi4/psi4/issues/2213#issuecomment-1583231316

[psi-rking]

Today I looked at this ~D3h Mn(NO)(CO)4 complex. It is pathologically interesting.

1. There are no reasonable torsions defined among the bonded atoms in this molecule, as each of the 5 spokes out from the metal involve 2 ~collinear atoms.
2. The reason that optking was trying to use torsions is that the input geometry is so crowded that the covalent radii * 1.3 formula was bonding the Mn to the atom once removed (beta?). That is, in the linear Mn-N-O, both Mn-N and Mn-O bonds were being created. So I tried covalent_connect=1.2, which results in only linear bends being defined and no torsions. This raises the question of whether we should use something beyond a linear scaling, or count on the user to provide a reasonable structure. (The clean tool in Spartan lengthened the Mn-N bond by 25%.). For the moment, I am not changing.
3. Using these coordinates (all stretches and bends), the first step blew up to 100 Angstroms Cartesian step. This was not due to forces, which were large but not that ridiculous.
4. The numerical test of the B matrix passed. (a non-trivial accomplishment for this thing) I figured out that the pseudoinverse of (B B^T) was blowing up, due to inversion of very small values. It is customary to invert only those eigenvalues whose magnitude is >1e-10 or so. A numerical problem was introduced when my custom linear algebra function was replaced by numpy.linalg.pinv which uses a default of 1e-14. I'm actually surprised if this threshold isn't causing problems more commonly. Anyway, I increased this to 1e-12 and added a keyword to increase it further. More testing is needed to see if we can make the default larger. This particular case needs 1e-8 and then the optimization iterates successfully downhill.The geometry is very far from equilibrium and may be changing electronic states, so I'm not pursuing its STO-3G to a minimum whatever that might be.

So I am going to commit some related things but I'm not going to attend to this particular case further.