Resonance structure generation is creating radical catalyst surface sites

[rwest]

Bug Description

See cfgoldsmith#70 for further history, consequences, and debugging efforts, but the tl;dr:
rmgpy.molecule.resonance.generate_allyl_delocalization_resonance_structures is creating surface site atoms (element X) with radical electrons. These then can't be found in the adsorption thermo database, because surface sites (of a metal) shouldn't have unpaired electrons.

I'm not yet sure how best to address this.

• perhaps let it generate the resonance structure, then reset the radical count to 0 on the metal?
• perhaps prevent it from generating that resonance structure?
  • perhaps do this by having update_atomtypes raise an AtomTypeError if there's a radical on an X?
• perhaps prevent entire resonance forms from occurring in adsorbates?

Any suggestions?

How To Reproduce

Put an assert not self.is_surface_site() in Atom.increment_radical. (eg. rwest@b6b8751 )
then build a model with heterogeneous catalysis.

For me the resonance structure that causes the problem is *N=O

Error: multiplicity 2
1 O u0 p2 c0 {2,D}
2 N u0 p1 c0 {1,D} {3,S}
3 X u1 p0 c0 {2,S}

though I'm not sure what the form was before this structure was generated.
Presumably it was something like...

Error: multiplicity 2
1 O u1 p2 c0 {2,S}
2 N u0 p1 c0 {1,S} {3,D}
3 X u0 p0 c0 {2,D}

??

Expected Behavior

It wouldn't make radical metals and, more importantly, wouldn't crash.

Installation Information

Describe your installation method and system information.

• OS: Mac OS and CentOS
• Installation method: from source in anaconda environment
• RMG version information:
  • RMG-Py: 2.4.1-508-gb6b8751ef
  • RMG-database: 2.4.0-37-gcc0c77f9

Additional Context

Another observation, perhaps related, is that for many (perhaps all?) recent runs with catalysis, a large fraction of the overall CPU time is spent on resonance structure generation (see cfgoldsmith#69 (comment) ). This may be just that we're making molecules with many resonance structures (as discussed in that thread), but it struck me that perhaps the resonance structure algorithms perform poorly for adsorbates. This may be a red herring though.

[rwest]

The previous species added to the core was

*N[O]
multiplicity 2
1 O u1 p2 c0 {2,S}
2 N u0 p1 c0 {1,S} {3,S} {4,S}
3 H u0 p0 c0 {2,S}
4 X u0 p0 c0 {2,S}

so it's quite likely that this lost an H forming

Error: multiplicity 2
1 O u1 p2 c0 {2,S}
2 N u0 p1 c0 {1,S} {3,D}
3 X u0 p0 c0 {2,D}

as hypothesized above

[rwest]

I tried just making it so that the action items increment_radical, decrement_radical, increment_lone_pairs, and decrement_lone_pairs do nothing to surface sites. That way a resonance structure generation that tries to put a radical on the surface will just have no effect on it, but it'll make the structure anyway (like the option "let it generate the resonance structure, then reset the radical count to 0 on the metal" above).
rwest@b727616

That seemed to solve things and the job ran for further, but about 20 species later, it crashed with

For reaction generation 4 processes are used.
Exception in thread Thread-255:
Traceback (most recent call last):
  File "/Users/rwest/anaconda/envs/rmg3/lib/python3.7/threading.py", line 926, in _bootstrap_inner
    self.run()
  File "/Users/rwest/anaconda/envs/rmg3/lib/python3.7/threading.py", line 870, in run
    self._target(*self._args, **self._kwargs)
  File "/Users/rwest/anaconda/envs/rmg3/lib/python3.7/multiprocessing/pool.py", line 470, in _handle_results
    task = get()
  File "/Users/rwest/anaconda/envs/rmg3/lib/python3.7/multiprocessing/connection.py", line 251, in recv
    return _ForkingPickler.loads(buf.getbuffer())
  File "rmgpy/species.py", line 130, in rmgpy.species.Species.__init__
rmgpy.exceptions.SpeciesError: Multiplicities of molecules in species  do not match.

I guess there's now a species where some resonance structures have radicals (on the adsorbate) and some do not (because the unpaired electron was transferred to the metal, then disappeared). There's a check when making a new Species that the multiplicities of all its constituent molecules (resonance structures) match, which would usually make sense (so a species is either a radical or is not). This is only enforced when calling __init__, which is why we got away with it for a while, but that happens sometimes when passing things back and forth through the multiprocessing (I was generating reactions on 4 cores in parallel.

So either...

• making resonance structures the way I did was a bad idea, or
• we need to relax the "all resonance structures have the same multiplicity" assumption for surface adsorbates.

To be honest, I'm not very confident about radical adsorbates.

[alongd]

Some initial thoughts:

• If the only context for X is part of a bulk (rather than a few metal particles in a fluid), then I agree we should not care about the electronic state of X (set its u and p in the adjlist to either None or np.inf?).
• For cases like X=N-[O], seems like O should be more electronegative, and the structure suggested by RMG's resonance, X-N=O should perhaps be the representative / reactive structure.
• If a radical site is conjugated to X, we may allow its disappearance, otherwise, it should be preserved (?).

[rwest]

I permitted the different resonance structures of an adsorbate to have different multiplicities in rwest@d33df77
and it seems that nothing crashed.

[github-actions]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.