Fix bi-rad singlet composite jobs silently running a restricted reference

arc/job/adapters/common.py

@@ -660,21 +660,28 @@ def is_species_restricted(obj: JobAdapter,
         bool: Whether to run as restricted (``True``) or not (``False``).
     """
 
-    if obj.level.method_type in REFERENCE_AGNOSTIC_METHOD_TYPES:
-        return True
-
     multiplicity = obj.multiplicity if species is None else species.multiplicity
     species_obj = obj.species[0] if species is None else species
     number_of_radicals = species_obj.number_of_radicals
     species_label = species_obj.label
-    if multiplicity > 1 or (number_of_radicals is not None and number_of_radicals > 1):
-        # run an unrestricted electronic structure calculation if the spin multiplicity is greater than one,
-        # or if it is one but the number of radicals is greater than one (e.g., bi-rad singlet)
-        # don't run unrestricted for composite methods such as CBS-QB3, it'll be done automatically if the
-        # multiplicity is greater than one, but do specify uCBS-QB3 for example for bi-rad singlets.
-        if number_of_radicals is not None and number_of_radicals > 1:
-            logger.info(f'Using an unrestricted method for species {species_label} which has '
-                        f'{number_of_radicals} radicals and multiplicity {multiplicity}.')
+
+    if number_of_radicals is not None and number_of_radicals > 1:
+        # A declared number_of_radicals > 1 with multiplicity 1 (e.g., a bi-rad singlet) always
+        # needs an unrestricted reference, *even for a reference-agnostic method type* such as a
+        # composite method (e.g., CBS-QB3, G4): for those, an unrestricted reference is applied
+        # automatically when the multiplicity is greater than one, but a bi-rad singlet's
+        # multiplicity of one gives the method no such signal, so ARC must still specify it
+        # explicitly (e.g., uCBS-QB3, uG4). This check must therefore run before the
+        # REFERENCE_AGNOSTIC_METHOD_TYPES early return below, not after it.
+        logger.info(f'Using an unrestricted method for species {species_label} which has '
+                    f'{number_of_radicals} radicals and multiplicity {multiplicity}.')
+        return False
+
+    if obj.level.method_type in REFERENCE_AGNOSTIC_METHOD_TYPES:
+        return True
+
+    if multiplicity > 1:
+        # run an unrestricted electronic structure calculation if the spin multiplicity is greater than one
         return False
     if adopted_reference_is_unrestricted(species_obj):
         if not level_admits_a_broken_symmetry_reference(obj.level):

arc/job/adapters/common_test.py

@@ -5,6 +5,8 @@
 This module contains unit tests of the arc.job.adapters.common module
 """
 
+import glob
+import os
 import shutil
 import tempfile
 import unittest
@@ -185,6 +187,50 @@ def test_a_composite_level_ignores_a_derived_verdict(self):
         self.assertEqual(job.level.method_type, 'composite')
         self.assertTrue(common.is_species_restricted(job))
 
+    def test_a_bi_rad_singlet_composite_job_is_unrestricted(self):
+        """Test that a declared bi-rad singlet (multiplicity 1, number_of_radicals 2) is run
+        unrestricted even for a composite method, whose reference-agnostic early return would
+        otherwise make it restricted (the u prefix, e.g. uCBS-QB3/uG4, must still be specified)"""
+        for method in ['cbs-qb3', 'g4']:
+            job = self._singlet_job(method=method, number_of_radicals=2)
+            self.assertEqual(job.level.method_type, 'composite')
+            self.assertFalse(common.is_species_restricted(job),
+                             msg=f'a bi-rad singlet at {method} was not run unrestricted')
+
+    def test_a_closed_shell_singlet_composite_job_stays_restricted(self):
+        """Test that an ordinary closed-shell singlet composite job is unaffected by the fix"""
+        job = self._singlet_job(method='cbs-qb3', number_of_radicals=None)
+        self.assertEqual(job.level.method_type, 'composite')
+        self.assertTrue(common.is_species_restricted(job))
+
+    def test_a_composite_job_above_singlet_multiplicity_is_unchanged(self):
+        """Test that a composite job at multiplicity > 1 keeps returning restricted, as it did
+        before the fix: the method itself is expected to go unrestricted automatically"""
+        job = self._singlet_job(method='cbs-qb3', multiplicity=3, number_of_radicals=None)
+        self.assertEqual(job.level.method_type, 'composite')
+        self.assertTrue(common.is_species_restricted(job))
+
+    def test_a_bi_rad_singlet_composite_job_writes_the_u_prefix_in_the_route_section(self):
+        """End-to-end: a bi-rad singlet composite job's generated Gaussian input carries the u prefix"""
+        species = ARCSpecies(label='O2_a1Dg', xyz=['O 0 0 0', 'O 0 0 1.2'], multiplicity=1, number_of_radicals=2)
+        project_directory = tempfile.mkdtemp(prefix='arc_test_common_')
+        self.addCleanup(shutil.rmtree, project_directory, ignore_errors=True)
+        job = GaussianAdapter(execution_type='incore',
+                              job_type='composite',
+                              level=Level(method='g4'),
+                              project='test',
+                              project_directory=project_directory,
+                              species=[species],
+                              testing=True,
+                              )
+        job.write_input_file()
+        input_files = glob.glob(os.path.join(project_directory, '**', 'input.gjf'), recursive=True)
+        self.assertEqual(len(input_files), 1)
+        with open(input_files[0], 'r') as f:
+            content = f.read()
+        self.assertIn('ug4', content)
+        self.assertNotIn(' g4', content)
+
     def test_a_correlated_single_point_is_not_flipped_by_an_adopted_verdict(self):
         """Test that an adopted verdict decides no reference for a correlated wavefunction level"""
         for method in ['dlpno-ccsd(t)', 'ccsd(t)-f12', 'ccsd(t)', 'mp2']:

arc/species/species.py

@@ -1008,6 +1008,17 @@ def from_dict(self, species_dict):
                 self.mol = rmg_mol_from_inchi(inchi)
             elif smiles is not None:
                 self.mol = Molecule(smiles=smiles)
+        # Read the multiplicity and charge implied by the structure the user actually GAVE (a `mol`
+        # in the dict, or an adjacency list / InChI / SMILES) BEFORE mol_from_xyz() replaces self.mol
+        # with one perceived from the coordinates. Perception routinely lands on a different
+        # electronic state for an open-shell species -- `C[C]C#N` (multiplicity 3) perceives as
+        # `C[C][C][N]` (multiplicity 1) -- and reading the multiplicity off the perceived molecule
+        # afterwards silently runs the entire job on that other state, with no warning anywhere.
+        # ARCSpecies.__init__ already reads both here, before its own (guarded) mol_from_xyz call;
+        # from_dict differed from it by statement order alone, so the same input gave the two
+        # constructors different answers.
+        structure_multiplicity = self.mol.multiplicity if self.mol is not None else None
+        structure_charge = self.mol.get_net_charge() if self.mol is not None else None
         # Perceive molecule from xyz coordinates. This also populates the .mol attribute of the Species.
         # It overrides self.mol generated from adjlist or smiles so xyz and mol will have the same atom order.
         if self.final_xyz or self.initial_xyz or self.most_stable_conformer or self.conformers or self.ts_guesses:
@@ -1023,9 +1034,10 @@ def from_dict(self, species_dict):
                     logger.debug(f'TS species {self.label}: using xyz-based multiplicity '
                                  f'{self.multiplicity} (ignored mol.multiplicity)')
                 else:
-                    self.multiplicity = self.mol.multiplicity
+                    self.multiplicity = structure_multiplicity if structure_multiplicity is not None \
+                        else self.mol.multiplicity
             if self.charge is None:
-                self.charge = self.mol.get_net_charge()
+                self.charge = structure_charge if structure_charge is not None else self.mol.get_net_charge()
         if 'conformers' in species_dict:
             self.conformers = [str_to_xyz(conf) for conf in species_dict['conformers']]
             self.conformer_energies = species_dict['conformer_energies'] if 'conformer_energies' in species_dict \

arc/species/species_test.py

@@ -966,6 +966,34 @@ def test_stability_sequencing_state_round_trip(self):
         self.assertEqual(restored.stability_pending_opt_job, 'opt_a7')
         self.assertTrue(restored.stability_reoptimized)
 
+    def test_from_dict_reads_multiplicity_from_the_given_structure(self):
+        """Test that from_dict() and the keyword constructor agree on multiplicity and charge.
+
+        Both constructors replace .mol with a molecule perceived from the coordinates so that atom
+        orders match, and neither may read the multiplicity off that perceived molecule: perception
+        picks an electronic state, and for an open-shell species it routinely picks a different one
+        from the structure that was given. 'C[C]C#N' is a triplet, and these same atoms perceive as
+        the closed-shell 'C[C][C][N]'. Reading the multiplicity after the override therefore ran a
+        triplet as a singlet with no warning anywhere - a wrong answer rather than a failure.
+        """
+        xyz = """C       1.83783800   -0.00006500   -0.00024000
+                 C       0.39175700    0.00019900    0.00069900
+                 C      -0.89358800    0.00004200    0.00009600
+                 N      -2.10264600   -0.00008000   -0.00026000
+                 H       2.23488200    0.14322800    1.01178500
+                 H       2.23365500   -0.94854900   -0.38246700
+                 H       2.23394500    0.80482700   -0.63082500"""
+        kwargs_spc = ARCSpecies(label='CCCN_triplet', smiles='C[C]C#N', xyz=xyz)
+        dict_spc = ARCSpecies(species_dict={'label': 'CCCN_triplet', 'smiles': 'C[C]C#N', 'xyz': xyz})
+        self.assertEqual(kwargs_spc.multiplicity, 3)
+        self.assertEqual(dict_spc.multiplicity, kwargs_spc.multiplicity)
+        self.assertEqual(dict_spc.charge, kwargs_spc.charge)
+
+        # An explicitly declared multiplicity still wins over the structure.
+        explicit = ARCSpecies(species_dict={'label': 'CCCN_triplet', 'smiles': 'C[C]C#N',
+                                            'xyz': xyz, 'multiplicity': 3})
+        self.assertEqual(explicit.multiplicity, 3)
+
     def test_from_dict(self):
         """Test Species.from_dict()"""
         species_dict = self.spc2.as_dict()