Td Molecule Breaks Symmetry-Alignment with Ref. Axes During Optimization

[JonathonMisiewicz]

Transferred from the forums and possibly related to #1177.

In the following optimization, Psi quits after reporting a tetrahedral molecule falls from C2v to Cs molecular point group while maintaining Td full point group after the first step.

import psi4
mol_string = """ 
     C
     H 1 %f
     H 1 %f 2 TDA
     H 1 %f 2 TDA 3 120
     H 1 %f 2 TDA 4 120
"""
mol = lambda x: mol_string%(x,x,x,x)
psi4.set_options({"BASIS": "def2-qzvppd", "scf_type": "df"})
geom = psi4.geometry(mol(1.106))
e,wfn = psi4.optimize("PBE",molecule=geom, return_wfn="yes")

The first step appears to break the alignment of the molecular symmetry axes with the reference axes. For comparison, here are the starting coordinates (Bohr) and the coordinates after step 1:

[[ 0.          0.          0.        ]
 [ 0.          1.70650814 -1.20668348]
 [ 1.70650814  0.          1.20668348]
 [-1.70650814  0.          1.20668348]
 [ 0.         -1.70650814 -1.20668348]]

[[-6.64087882e-10  2.94128033e-09  0.00000000e+00]
 [-6.64087882e-10  2.06989822e+00  0.00000000e+00]
 [-9.75759370e-01 -6.89966085e-01  1.69006481e+00]
 [-9.75759370e-01 -6.89966085e-01 -1.69006481e+00]
 [ 1.95151875e+00 -6.89966086e-01  0.00000000e+00]]

The original reporter said this problem was observed with this basis set but not others. To me, this suggests numerical noise determines whether the bug occurs.

[susilehtola]

This sounds similar to the issue I ran into with the first version of the SAD guess rewrite: CH4 switched symmetry and the program crashed due to memory errors.

[dgasmith]

Can you print the gradient at the first step, it should not be possible to break that symmetry as we symmetrize the gradient to cleanup small noise. Looking at the molecule you have a completely different orientation than the starting one, where a YZ plan switch and shift appears to have happened.

If you track the point group does it change per iteration?

[JonathonMisiewicz]

Gradient and geometry are as follows:

	 C           0.0000000000        0.0000000000        0.0000000000
	 H           0.0000000000        1.7065081419       -1.2066834793
	 H           1.7065081419        0.0000000000        1.2066834793
	 H          -1.7065081419        0.0000000000        1.2066834793
	 H           0.0000000000       -1.7065081419       -1.2066834793
	            -0.0000000000        0.0000000000       -0.0000000011
	            -0.0000000000        0.0054395963       -0.0038463766
	             0.0054395971        0.0000000000        0.0038463771
	            -0.0054395971       -0.0000000000        0.0038463771
	             0.0000000000       -0.0054395963       -0.0038463766

The YZ plane shift surprised me at first, but then I remembered orientation is discontinuous in Psi. The shift still surprises me.

I should have been clearer: it is the very first step that reduces the computational point group. So yes, the point group does change per geometry iteration.

[JonathonMisiewicz]

I decided after the previous post to repeat the gradient check for def2-tzvp, which the original reporter said doesn't display this problem. The results are:

	---Fragment 1 Geometry and Gradient---
	 C           0.0000000000        0.0000000000        0.0000000000
	 H           0.0000000000        1.7065081419       -1.2066834793
	 H           1.7065081419        0.0000000000        1.2066834793
	 H          -1.7065081419        0.0000000000        1.2066834793
	 H           0.0000000000       -1.7065081419       -1.2066834793
	            -0.0000000000        0.0000000000       -0.0000000001
	            -0.0000000000        0.0048395198       -0.0034220573
	             0.0048395199        0.0000000000        0.0034220573
	            -0.0048395199       -0.0000000000        0.0034220573
	             0.0000000000       -0.0048395198       -0.0034220573

In the bugged example, the -z component of the gradients disagree in the last printed decimal place, which they don't in the non-bugged example. So am I correct in assuming something has gone wrong with the gradient?

[dgasmith]

Hmm, dollars to donuts we are only symmetrizing to the closest Abelian symmetry, not the full point group symmetry. I assume Psi find the molecule as Td, but guessing only Cs molecular symmetry?

In matrix.cc:958 I believe we want molecule->find_highest_point_group not molecule->point_group I think? @loriab @andysim

[loriab]

Been a while, but I think we don't have character table for non-Abelian point groups. https://github.com/psi4/psi4/blob/master/psi4/src/psi4/libmints/maketab.cc

Symmetry detection happens in molecule.cc, so we do have non-Abelian there, but only for reporting.

[dgasmith]

Should be closed by #1483 correct?

[loriab]

Should, but worth another check.

[JonathonMisiewicz]

It's fixed now, but how? #1483 was just a DFMP2 fix. PBE isn't an MP2 hybrid, and it certainly isn't MP2.

[JonathonMisiewicz]

Testing indicates 1483 was not the fix. Whatever the fix was, it was in-between dev359 84ee495 and dev377 6b93fa5. At least it's fixed?