diff --git a/arc/job/adapters/common.py b/arc/job/adapters/common.py index b3472c251b..7a03888d29 100644 --- a/arc/job/adapters/common.py +++ b/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): diff --git a/arc/job/adapters/common_test.py b/arc/job/adapters/common_test.py index e4eb970056..0507c3ddc4 100644 --- a/arc/job/adapters/common_test.py +++ b/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']: diff --git a/arc/species/species.py b/arc/species/species.py index 3dae4e4051..b3cdf3123f 100644 --- a/arc/species/species.py +++ b/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 \ diff --git a/arc/species/species_test.py b/arc/species/species_test.py index 815c4f426b..dd9cebb4ac 100644 --- a/arc/species/species_test.py +++ b/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()