Cs symmetry is not detected

[TiborGY]

Psi4 seems to be unable to recognize the symmetry of molecules that should be Cs. Minimal input:

molecule {
pubchem:fluoroethane               
}
set basis cc-pVDZ
energy('scf')

Output is attached. Cs_test.log I have tried overriding this by manually specifying a Cs symmetry, but that results in an error:

!----------------------------------------------------------------------------------!
!                                                                                  !
! Fatal Error: User specified point group (Cs(X)) is not a subgroup of the highest !
!     detected point group (C1). If this is because the symmetry increased, try to !
!     start the calculation again from the last geometry, after checking any       !
!     symmetry-dependent input, such as DOCC.                                      !
! Error occurred in file: /scratch/psilocaluser/conda-                             !
!     builds/psi4-multiout_1645476997931/work/psi4/src/psi4/libmints/molecule.cc   !
!     on line: 2011                                                                !
! The most recent 5 function calls were:                                           !
! psi::Molecule::find_point_group(double) const                                    !
! psi::Molecule::update_geometry()                                                 !
! from_dict(pybind11::dict)                                                        !
!                                                                                  !
!----------------------------------------------------------------------------------!

[susilehtola]

The problem is that the pubchem structure does not appear to have Cs symmetry.

The pubchem geometry gives

  ==> Geometry <==

    Molecular point group: c1
    Full point group: C1

    Geometry (in Angstrom), charge = 0, multiplicity = 1:

       Center              X                  Y                   Z               Mass       
    ------------   -----------------  -----------------  -----------------  -----------------
         F            1.054979085535     0.279491800552     0.000000000000    18.998403162730
         C           -0.000320914465    -0.577108199448     0.000000000000    12.000000000000
         C           -1.279120914465     0.226491800552     0.000000000000    12.000000000000
         H            0.072379085535    -1.206508199448    -0.891000000000     1.007825032230
         H            0.072279085535    -1.206508199448     0.891000000000     1.007825032230
         H           -2.153920914465    -0.429408199448    -0.000100000000     1.007825032230
         H           -1.322020914465     0.874191800552     0.881000000000     1.007825032230
         H           -1.321920914465     0.874291800552    -0.880900000000     1.007825032230

  Running in c1 symmetry.

  Rotational constants: A =      1.21162  B =      0.32319  C =      0.28216 [cm^-1]
  Rotational constants: A =  36323.59189  B =   9689.02155  C =   8459.03982 [MHz]
  Nuclear repulsion =   80.019990391667150

Reading this into IQmol and running "Symmetrize geometry", I get

  ==> Geometry <==

    Molecular point group: cs
    Full point group: Cs

    Geometry (in Angstrom), charge = 0, multiplicity = 1:

       Center              X                  Y                   Z               Mass       
    ------------   -----------------  -----------------  -----------------  -----------------
         F           -0.264271602373     1.058895456489     0.000000000000    18.998403162730
         C            0.577048397627    -0.008634543511     0.000000000000    12.000000000000
         C           -0.244891602373    -1.275724543511     0.000000000000    12.000000000000
         H            1.207428397627     0.054945456489    -0.891000000000     1.007825032230
         H            1.207428397627     0.054945456489     0.891000000000     1.007825032230
         H            0.398348397627    -2.159884543511     0.000000000000     1.007825032230
         H           -0.893191602373    -1.309244543511     0.880950000000     1.007825032230
         H           -0.893191602373    -1.309244543511    -0.880950000000     1.007825032230

  Running in cs symmetry.

  Rotational constants: A =      1.21161  B =      0.32319  C =      0.28216 [cm^-1]
  Rotational constants: A =  36323.26201  B =   9689.02376  C =   8459.02361 [MHz]
  Nuclear repulsion =   80.019847421417836

However, the HF/cc-pVDZ energies appear almost identical: -178.08964068212876 for the pubchem structure vs -178.08964078908636 for the symmetrized structure. (The symmetrized structure is 1e-7 Eh lower.)

[susilehtola]

Maybe there's a threshold that should be loosened?

[TiborGY]

According to the docs psi::Molecule::find_point_group is the function that does the detection, and it can take a threshold argument, but I cannot find a keyword that would change that. Have not looked at the source yet, but I think the threshold may be hardcoded at the moment.

[hokru]

The symmetry detection does not try to symmetrize the geometry to my knowledge. But there is a symmetrize function than can be explicitly called like this:

molecule {
pubchem:fluoroethane
}
set basis cc-pVDZ
psi4.get_active_molecule().symmetrize(1e-2)
# or <molecule_name>.symmetrize(1e-2) if a name was set
energy('scf')

edit: Otherwise finding the point group with a lower tolerance can be done like this, I think. (Doesn't work in this case)

molecule {
pubchem:fluoroethane
}
set basis cc-pVDZ

mol = psi4.get_active_molecule()
tol = 1e-2
PG = mol.find_point_group(tol)
print(PG.full_name())
mol.set_point_group(PG)
energy('scf')

[JonathonMisiewicz]

The action item here is to have some way to adjust the symmetrization threshold.

[TiborGY]

The action item here is to have some way to adjust the symmetrization threshold.

I think updating the docs would also be a good idea.

The symmetry detection does not try to symmetrize the geometry to my knowledge. But there is a symmetrize function than can be explicitly called

It has been a while, but I cannot recall seeing this anywhere when I was reading the docs.