Geometry optimization fails due to changing comp. point group -- follow up

[Geoff]

Similar to this issue from a number of years ago: Geometry optimization fails due to changing comp. point group

I’m generating some coordinates for planar conjugated molecules using rdkit and using the resultant cartesians for the optimizations, so the geoms look something like this:

0 1
C -2.1962244741105725 -0.46169909698923206 -0.20891270685180444
C -0.8758009187983459 0.20980442180480696 -0.08232387781514809
N -0.7315608609335387 1.5462856680004695 -0.0635822140001667
C 0.4377411643912381 2.179525760223909 0.04869127100579021
N 1.5951030640963157 1.4780825215574294 0.15465206209955645
C 1.5897841087193778 0.13417294785168185 0.14892273545322615
C 2.553443832943142 -0.8558486769060368 0.2355541801052097
N 1.921333468436344 -2.031727904283795 0.17160855449182943
C 0.5833492989447032 -1.8720755405092093 0.046589129399722345
C 0.29513849732858344 -0.5271363216806334 0.024769718248766878
H -2.3800922299955465 -1.19430775651866 0.5986956005595973
H -2.19592941119346 -0.9975468624678573 -1.1639239326072979
H -2.9820300302821745 0.3128091427674146 -0.13924187648261496
H -1.5874317691678135 2.160902681628291 -0.1415789955399663
H 0.5125207891983367 3.2689098054760373 0.05992780257781991
H 3.633431674364055 -0.7075268537503937 0.3375740409229852
H -0.17277620394060184 -2.6426239362042336 -0.027421491567487515

There’s probably something to do with the fact that it’s C1 according to the coordinates, but the optimizaer fails with:

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psikit/psikit.py:72, in Psikit.optimize(self,     basis_sets, return_wfn, name, multiplicity, maxiter)
     70 self.psi4.set_options({'GEOM_MAXITER':maxiter})
 71 try:
---> 72     scf_energy, wfn = self.psi4.optimize(basis_sets, return_wfn=return_wfn)
     73     self.wfn = wfn
     74 except self.psi4.OptimizationConvergenceError as cError:

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/driver.py:1341, in optimize(name, **kwargs)
   1339 # Compute the gradient - preserve opt data despite core.clean calls in gradient
   1340 core.IOManager.shared_object().set_specific_retention(1, True)
-> 1341 G, wfn = gradient(lowername, return_wfn=True, molecule=moleculeclone, **kwargs)
   1342 thisenergy = core.variable('CURRENT ENERGY')
   1344 # above, used to be getting energy as last of energy list from gradient()
   1345 # thisenergy below should ultimately be testing on wfn.energy()
   1346 
   1347 # Record optimization steps
   1348 # Add wavefunctions later

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/driver.py:700, in gradient(name, **kwargs)
    697 lowername = name.lower()
    698 if dertype == 1:
    699     # Bounce to CBS in pure-gradient mode if "method/basis" name and all parts have analytic grad. avail.
--> 700     return driver_cbs._cbs_gufunc(gradient, name, ptype='gradient', **kwargs)
    701 else:
    702     # Set method-dependent scf convergence criteria (test on procedures['energy'] since that's guaranteed)
    703     optstash = driver_util._set_convergence_criterion('energy', cbs_methods[0], 8, 10, 8, 10, 8)

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/driver_cbs.py:1965, in     _cbs_gufunc(func, total_method_name, **kwargs)
   1963 optstash = p4util.OptionsState(['BASIS'])
   1964 core.set_global_option('BASIS', basis)
    -> 1965 ptype_value, wfn = func(method_name, return_wfn=True, molecule=molecule, **kwargs)
   1966 if core.get_option("SCF", "DF_INTS_IO") != "SAVE":
   1967     core.clean()

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/driver.py:741, in gradient(name, **kwargs)
    738     core.print_out("""gradient() will perform analytic gradient computation.\n""")
    740     # Perform the gradient calculation
--> 741     wfn = procedures['gradient'][lowername](lowername, molecule=molecule, **kwargs)
    743 else:
744     core.print_out("""gradient() will perform gradient computation by finite difference of analytic energies.\n""")

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/procrouting/proc.py:2485, in run_scf_gradient(name, **kwargs)
   2483 ref_wfn = kwargs.get('ref_wfn', None)
   2484 if ref_wfn is None:
-> 2485     ref_wfn = run_scf(name, **kwargs)
   2487 if core.get_option('SCF', 'REFERENCE') in ['ROHF', 'CUHF']:
   2488     ref_wfn.semicanonicalize()

File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/procrouting/proc.py:2390, in run_scf(name, **kwargs)
   2387 if not core.has_global_option_changed('SCF_TYPE'):
   2388     core.set_global_option('SCF_TYPE', 'DF')
-> 2390 scf_wfn = scf_helper(name, post_scf=False, **kwargs)
   2391 returnvalue = scf_wfn.energy()
   2393 ssuper = scf_wfn.functional()

    File ~/opt/anaconda3/envs/psi4/lib/python3.9/site-packages/psi4/driver/procrouting/proc.py:1526, in scf_helper(name, post_scf, **kwargs)
   1523 Cb_occ = old_wfn.Cb_subset("SO", "OCC")
   1525 if old_wfn.molecule().schoenflies_symbol() != scf_molecule.schoenflies_symbol():
    -> 1526     raise ValidationError("Cannot compute projection of different symmetries.")
   1528 if old_wfn.basisset().name() == scf_wfn.basisset().name():
   1529     core.print_out(f"  Reading orbitals from file {read_filename}, no projection.\n\n")

ValidationError: Cannot compute projection of different symmetries.

Is there someway to turn symmetry off or is this still the same bug as before? I’m using:

psi4 1.5+e9f4d6d py39ha809fef_0 psi4

[jmisiewicz]

Add symmetry c1 to the molecule specification.

[Geoff]

I think a great feature request would be having the code compute the Abelian group at the start and parse that to the optimizer most other QM codes do this.

[jmisiewicz]

I don’t understand what you’re requesting. “parse that to” means some kind of datatype conversion, but converting an abelian group into a geometry optimizer is not what you mean. Maybe you mean “pass that to,” but then I’m not sure what the geometry optimizer is supposed to do with that information.

Are you suggesting that the code compute the (maximal finite) Abelian group at the start and then keeps the molecule in that point group throughout optimization?

[Geoff]

" Are you suggesting that the code compute the (maximal finite) Abelian group at the start and then keeps the molecule in that point group throughout optimization?"

Yes! That would be great. If you start in C1 why would you want to convert to Cs and fall over? Gaussian, Q-Chem and others do this at the start of an optimization and keeps that symmetry throughout

[jmisiewicz]

The maximum finite Abelian group is already computed at the start, so the sticking point is “keep the molecule in its original point group by default.”

This is a useful default for some users (like yourself), but not for others. Whether a weakly bound complex was in a Cs or a C2v point group was extremely important for a PhD project of mine, and I would have wanted to know “the point group changed, look this over before doing any more computations” at the earliest possible opportunity.

There’s absolutely no way to satisfy all users, and in cases where there is a legitimate use for the existing default and a way to get the non-default behavior (specifying “symmetry c1”), I’m inclined to keep the existing default.

If you want to push the issue further, the best thing you can do is to file an issue on our GitHub page. That’ll be more likely to attract the notice of other developers, who may have different opinions than I do.

[crawdad]

I’m a bit confused by this topic, because, if a molecule starts in a certain point group, a geometry optimization should not be able to break that symmetry unless the gradients are computed with inadequate numerical accuracy. It is possible to find a stationary point with higher symmetry than the initial geometry, but the symmetry should not drop if one is using analytical gradients. Am I misunderstanding the issue here?

[jmisiewicz]

Yes. In this particular example, we start with a molecule in C1 symmetry, and optimization (correctly) finds a Cs solution.

Also in this particular example, the problem occurs during the orbital read step, but I have a simple example where the point group change is deliberately trapped by the driver:

molecule {
O 0 0 0 
H 0 1 0 
H 0 0.0001 1
}

set g_convergence gau_verytight

optimize('scf/cc-pvdz')

yields

!----------------------------------------------------------------------------------!
!                                                                                  !
!  Point group changed! (c2v <-- cs) You should restart using the last geometry in !
!     the output, after carefully making sure all symmetry-dependent input, such   !
!     as DOCC, is correct.                                                         !
!                                                                                  !
!----------------------------------------------------------------------------------!

Why the reporter observes a different error message isn’t obvious to me.

[Geoff]

I’m not sure if I’m following this either. Your water example changed symmetry but the optimization converged? If so, does optimization output strictly change to C2v or does it stay as Cs, in your example?

[jmisiewicz]

No, the optimization didn’t converge. I posted the error message I got.

[Geoff]

OK thanks! It wasn’t obvious that was error or warning.