Make surface molecules from smiles using Argon workaround

rmgpy/molecule/adjlist.py

@@ -92,7 +92,7 @@ def check_partial_charge(atom):
         the theoretical one:
 
         """
-        if atom.symbol in {'X','L','R','e','H+','Li'}:
+        if atom.symbol in {'X','L','R','e','H+','Li'} or 'X' in [z.label for z in get_atomtype(atom, atom.bonds).generic]:
             return  # because we can't check it.
 
         valence = PeriodicSystem.valence_electrons[atom.symbol]

rmgpy/molecule/atomtype.py

@@ -270,7 +270,7 @@ def get_features(self):
     'Cl','Cl1s',
     'Br','Br1s',
     'I','I1s',
-    'F','F1s','X','Xv','Xo'])
+    'F','F1s','X','Xv','Xo','Pt','Ptv','Pto'])
 
 ATOMTYPES['Rx!H'] = AtomType(label='Rx!H', generic=['Rx'], specific=[
     'R!H',
@@ -286,20 +286,29 @@ def get_features(self):
     'Cl','Cl1s',
     'Br','Br1s',
     'I','I1s',
-    'F','F1s','X','Xv','Xo'])
+    'F','F1s','X','Xv','Xo','Pt','Ptv','Pto'])
 
 # Surface sites:
-ATOMTYPES['X'] = AtomType(label='X', generic=['Rx', 'Rx!H'], specific=['Xv', 'Xo'])
+ATOMTYPES['X'] = AtomType(label='X', generic=['Rx', 'Rx!H'], specific=['Xv', 'Xo', 'Pt', 'Ptv', 'Pto'],)
 
 # Vacant surface site:
-ATOMTYPES['Xv'] = AtomType('Xv', generic=['X','Rx', 'Rx!H'], specific=[],
-                             single=[0], all_double=[0], r_double=[], o_double=[], s_double=[], triple=[0], quadruple=[0],
+ATOMTYPES['Xv'] = AtomType('Xv', generic=['X','Rx', 'Rx!H'], specific=['Ptv'],
+                             single=[0], all_double=[0], r_double=[0], o_double=[0], s_double=[0], triple=[0], quadruple=[0],
                              benzene=[0], lone_pairs=[0])
 # Occupied surface site:
-ATOMTYPES['Xo'] = AtomType('Xo', generic=['X','Rx', 'Rx!H'], specific=[],
+ATOMTYPES['Xo'] = AtomType('Xo', generic=['X','Rx', 'Rx!H'], specific=['Pto'],
                              single=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], all_double=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], r_double=[], o_double=[], s_double=[], triple=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18],
                              quadruple=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], benzene=[0], lone_pairs=[0])
 
+ATOMTYPES['Pt'] = AtomType(label='Pt', generic=['Rx', 'Rx!H', 'X'], specific=['Ptv', 'Pto'])
+                            # single=[0], all_double=[0], r_double=[], o_double=[], s_double=[], triple=[0], quadruple=[0],
+                            # benzene=[0], lone_pairs=[0])
+ATOMTYPES['Ptv'] = AtomType(label='Ptv', generic=['Rx', 'Rx!H', 'X', 'Pt', 'Xv'], specific=[],
+                            single=[0], all_double=[0], r_double=[0], o_double=[0], s_double=[0], triple=[0], quadruple=[0],
+                            benzene=[0], lone_pairs=[0])
+ATOMTYPES['Pto'] = AtomType(label='Pto', generic=['Rx', 'Rx!H', 'X', 'Pt', 'Xo'], specific=[],
+                            single=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], all_double=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], r_double=[], o_double=[], s_double=[], triple=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18],
+                            quadruple=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18], benzene=[0], lone_pairs=[0])
 
 # Non-surface atomTypes, R being the most generic:
 ATOMTYPES['R'] = AtomType(label='R', generic=['Rx'], specific=[
@@ -847,8 +856,14 @@ def get_features(self):
 ATOMTYPES['F'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['F'], break_bond=['F'], increment_radical=['F'], decrement_radical=['F'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[])
 ATOMTYPES['F1s'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['F1s'], break_bond=['F1s'], increment_radical=['F1s'], decrement_radical=['F1s'], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[])
 
+ATOMTYPES['Pt'].set_actions(increment_bond=['Pt'], decrement_bond=['Pt'], form_bond=['Pt'], break_bond=['Pt'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[])
+ATOMTYPES['Ptv'].set_actions(increment_bond=[], decrement_bond=[], form_bond=['Pto'], break_bond=[], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[])
+ATOMTYPES['Pto'].set_actions(increment_bond=['Pto'], decrement_bond=['Pto'], form_bond=[], break_bond=['Ptv'], increment_radical=[], decrement_radical=[], increment_lone_pair=[], decrement_lone_pair=[], increment_charge=[], decrement_charge=[])
+
+
+
 # these are ordered in priority of picking if a more general atomtype is encountered
-allElements = ['H', 'C', 'O', 'N', 'S', 'P', 'Si', 'F', 'Cl', 'Br', 'I', 'Li', 'Ne', 'Ar', 'He', 'X', 'e', ]
+allElements = ['H', 'C', 'O', 'N', 'S', 'P', 'Si', 'F', 'Cl', 'Br', 'I', 'Li', 'Ne', 'Ar', 'He', 'X', 'e', 'Pt']
 # list of elements that do not have more specific atomTypes
 nonSpecifics = ['He', 'Ne', 'Ar', 'e']
 

rmgpy/molecule/element.py

@@ -123,14 +123,14 @@ class PeriodicSystem(object):
     isotopes of the same element may have slight different electronegativities, which is not reflected below
     """
     valences = {'H+':0, 'e': 0, 'H': 1, 'He': 0, 'C': 4, 'N': 3, 'O': 2, 'F': 1, 'Ne': 0,
-                'Si': 4, 'P': 3, 'S': 2, 'Cl': 1, 'Br': 1, 'Ar': 0, 'I': 1, 'X': 4, 'Li': 1}
+                'Si': 4, 'P': 3, 'S': 2, 'Cl': 1, 'Br': 1, 'Ar': 0, 'I': 1, 'X': 4, 'Li': 1, 'Pt': 4}
     valence_electrons = {'H+':0, 'e': 1, 'H': 1, 'He': 2, 'C': 4, 'N': 5, 'O': 6, 'F': 7, 'Ne': 8,
-                         'Si': 4, 'P': 5, 'S': 6, 'Cl': 7, 'Br': 7, 'Ar': 8, 'I': 7, 'X': 4, 'Li': 1}
+                         'Si': 4, 'P': 5, 'S': 6, 'Cl': 7, 'Br': 7, 'Ar': 8, 'I': 7, 'X': 4, 'Li': 1, 'Pt': 4}
     lone_pairs = {'H+':0, 'e': 0, 'H': 0, 'He': 1, 'C': 0, 'N': 1, 'O': 2, 'F': 3, 'Ne': 4,
-                  'Si': 0, 'P': 1, 'S': 2, 'Cl': 3, 'Br': 3, 'Ar': 4, 'I': 3, 'X': 0, 'Li': 0}
+                  'Si': 0, 'P': 1, 'S': 2, 'Cl': 3, 'Br': 3, 'Ar': 4, 'I': 3, 'X': 0, 'Li': 0, 'Pt': 0}
     electronegativity = {'H': 2.20, 'D': 2.20, 'T': 2.20, 'C': 2.55, 'C13': 2.55, 'N': 3.04, 'O': 3.44, 'O18': 3.44,
                          'F': 3.98, 'Si': 1.90, 'P': 2.19, 'S': 2.58, 'Cl': 3.16, 'Br': 2.96, 'I': 2.66, 'X': 0.0,
-                         'Li' : 0.98}
+                         'Li' : 0.98, 'Pt': 0}
 
 
 ################################################################################

rmgpy/molecule/molecule.py

@@ -38,6 +38,7 @@
 import itertools
 import logging
 import os
+import re
 from collections import OrderedDict, defaultdict
 from copy import deepcopy
 from urllib.parse import quote
@@ -432,7 +433,7 @@ def is_surface_site(self):
         """
         Return ``True`` if the atom represents a surface site or ``False`` if not.
         """
-        return self.symbol == 'X'
+        return self.symbol == 'X' or self.symbol in [z.label for z in ATOMTYPES['X'].specific]
 
     def is_bonded_to_surface(self):
         """
@@ -1205,9 +1206,11 @@ def contains_surface_site(self):
         Returns ``True`` iff the molecule contains an 'X' surface site.
         """
         cython.declare(atom=Atom)
+        cython.declare(z=AtomType)
         for atom in self.atoms:
-            if atom.symbol == 'X':
+            if atom.symbol == 'X' or atom.symbol in [z.label for z in ATOMTYPES['X'].specific]:
                 return True
+            # atom_type = get_atomtype(atom, atom.bonds)
         return False
 
     def number_of_surface_sites(self):
@@ -1862,6 +1865,30 @@ def from_smiles(self, smilesstr, backend='openbabel-first', raise_atomtype_excep
         single backend or try different backends in sequence. The available options for the ``backend``
         argument: 'openbabel-first'(default), 'rdkit-first', 'rdkit', or 'openbabel'.
         """
+        for surface_site_symbol in ['X', 'Pt']:
+            if surface_site_symbol in smilesstr:
+                assert 'Ar' not in smilesstr
+                self.from_smiles(smilesstr.replace(surface_site_symbol, 'Ar'))
+                lines = self.to_adjacency_list().split('\n')
+                for i, line in enumerate(lines):
+                    if 'Ar' in line:  # The adjacency list needs to use the identified 'X' for a site
+                        lines[i] = lines[i].replace('Ar', surface_site_symbol)
+                        # remove any extra electron pairs
+                        m = re.search(r'p[0-9]+', lines[i])  # searches for p0, p1, p2, p3, etc
+                        lines[i] = lines[i].replace(m[0], 'p0')
+                        m = re.search(r'u[0-9]+', lines[i])  # searches for u0, u1, u2, u3, etc
+                        lines[i] = lines[i].replace(m[0], 'u0')
+
+                        # remove any extra charge
+                        m = re.search(r'c[0-9+-]+', lines[i])  # searches for c0, c+2, c-1 etc
+                        lines[i] = lines[i].replace(m[0], 'c0')
+                adj_list = '\n'.join(lines)
+                self = self.from_adjacency_list(adj_list)
+                self._smiles = smilesstr
+                # but now we have to change the symbol back to 'Pt or 'X' for the smiles
+                # self.smiles = self.smiles.replace('X', surface_site_symbol)
+                return self
+
         translator.from_smiles(self, smilesstr, backend, raise_atomtype_exception=raise_atomtype_exception)
         return self
 

rmgpy/molecule/translator.py

@@ -139,7 +139,8 @@
     'I2': '[I][I]',
     'HI': 'I',
     'H': 'H+',
-    'e': 'e'
+    'e': 'e',
+
 }
 
 RADICAL_LOOKUPS = {
@@ -383,6 +384,13 @@ def _rdkit_translator(input_object, identifier_type, mol=None):
         output = from_rdkit_mol(mol, rdkitmol)
     elif isinstance(input_object, mm.Molecule):
         # We are converting from a molecule to a string identifier
+
+        # The to_rdkit_mol function in converter overwrites 'X' with 'Pt', so we need to keep track
+        # of whether this is a generic 'X' or specific 'Pt'
+        symbols = set([atom.element.symbol for atom in input_object.vertices])
+        assert 'X' not in symbols or 'Pt' not in symbols, 'X and Pt cannot be present at the same time'
+        generic_X = True if 'X' in symbols else False
+
         if identifier_type == 'smi':
             rdkitmol = to_rdkit_mol(input_object, sanitize=False)
         else:
@@ -401,6 +409,10 @@ def _rdkit_translator(input_object, identifier_type, mol=None):
                 output = Chem.MolToSmiles(rdkitmol, kekuleSmiles=True)
         else:
             raise ValueError('Identifier type {0} is not supported for writing using RDKit.'.format(identifier_type))
+
+        if generic_X:
+            output = output.replace('Pt', 'X')
+
     else:
         raise ValueError('Unexpected input format. Should be a Molecule or a string.')
 
@@ -439,6 +451,13 @@ def _openbabel_translator(input_object, identifier_type, mol=None):
         output = from_ob_mol(mol, obmol)
     elif isinstance(input_object, mm.Molecule):
         # We are converting from a Molecule to a string identifier
+
+        # The to_ob_mol function in converter overwrites 'X' with 'Pt', so we need to keep track
+        # of whether this is a generic 'X' or specific 'Pt'
+        symbols = set([atom.element.symbol for atom in input_object.vertices])
+        assert 'X' not in symbols or 'Pt' not in symbols, 'X and Pt cannot be present at the same time'
+        generic_X = True if 'X' in symbols else False
+
         if identifier_type == 'inchi':
             ob_conversion.SetOutFormat('inchi')
             ob_conversion.AddOption('w')
@@ -454,6 +473,9 @@ def _openbabel_translator(input_object, identifier_type, mol=None):
             raise ValueError('Unexpected identifier type {0}.'.format(identifier_type))
         obmol = to_ob_mol(input_object)
         output = ob_conversion.WriteString(obmol).strip()
+
+        if generic_X:
+            output = output.replace('Pt', 'X')
     else:
         raise ValueError('Unexpected input format. Should be a Molecule or a string.')
 

rmgpy/molecule/util.py

@@ -54,11 +54,6 @@ def get_element_count(obj):
                 else:
                     element_count[element] = int(count)
 
-        # For surface species, replace Pt with X again
-        if 'Pt' in element_count:
-            element_count['X'] = element_count['Pt']
-            del element_count['Pt']
-
         return element_count
 
     elif isinstance(obj, Molecule) or isinstance(obj, Fragment):

test/rmgpy/molecule/atomtypeTest.py

@@ -828,6 +828,71 @@ def setup_class(self):
         self.electron = Molecule().from_adjacency_list('''1 e u1 p0 c-1''')
         self.proton = Molecule().from_adjacency_list('''1 H u0 p0 c+1''')
 
+        self.X = Molecule().from_adjacency_list("""1 X u0 p0 c0""")
+        self.mol100 = Molecule().from_adjacency_list(
+            """1 H u0 p0 c0 {2,S}
+               2 X u0 p0 c0 {1,S}"""
+        )
+        self.mol101 = Molecule().from_adjacency_list(
+            """1 O u0 p2 c0 {2,D}
+               2 X u0 p0 c0 {1,D}"""
+        )
+        self.mol102 = Molecule().from_adjacency_list(
+            """1 N u0 p1 c0 {2,T}
+               2 X u0 p0 c0 {1,T}"""
+        )
+        self.mol103 = Molecule().from_adjacency_list(
+            """1 C u0 p0 c0 {2,Q}
+               2 X u0 p0 c0 {1,Q}"""
+        )
+        self.mol104 = Molecule().from_adjacency_list(  # bidentate
+            """1 N u0 p1 c0 {2,S} {5,D}
+               2 C u0 p0 c0 {1,S} {3,S} {4,D}
+               3 H u0 p0 c0 {2,S}
+               4 X u0 p0 c0 {2,D}
+               5 X u0 p0 c0 {1,D}"""
+        )
+        self.mol105 = Molecule().from_adjacency_list(  # vdW
+            """1 O u0 p2 c0 {2,S} {3,S}
+               2 H u0 p0 c0 {1,S}
+               3 H u0 p0 c0 {1,S}
+               4 X u0 p0 c0"""
+        )
+
+
+        self.Pt = Molecule().from_adjacency_list("""1 Pt u0 p0 c0""")
+        self.mol110 = Molecule().from_adjacency_list(
+            """1 H  u0 p0 c0 {2,S}
+               2 Pt u0 p0 c0 {1,S}"""
+        )
+        self.mol111 = Molecule().from_adjacency_list(
+            """1 O  u0 p2 c0 {2,D}
+               2 Pt u0 p0 c0 {1,D}"""
+        )
+        self.mol112 = Molecule().from_adjacency_list(
+            """1 N  u0 p1 c0 {2,T}
+               2 Pt u0 p0 c0 {1,T}"""
+        )
+        self.mol113 = Molecule().from_adjacency_list(
+            """1 C  u0 p0 c0 {2,Q}
+               2 Pt u0 p0 c0 {1,Q}"""
+        )
+        self.mol114 = Molecule().from_adjacency_list(  # bidentate
+            """1 N  u0 p1 c0 {2,S} {5,D}
+               2 C  u0 p0 c0 {1,S} {3,S} {4,D}
+               3 H  u0 p0 c0 {2,S}
+               4 Pt u0 p0 c0 {2,D}
+               5 Pt u0 p0 c0 {1,D}"""
+        )
+        self.mol115 = Molecule().from_adjacency_list(  # vdW
+            """1 O  u0 p2 c0 {2,S} {3,S}
+               2 H  u0 p0 c0 {1,S}
+               3 H  u0 p0 c0 {1,S}
+               4 Pt u0 p0 c0"""
+        )
+
+
+
     def atom_type(self, mol, atom_id):
         atom = mol.atoms[atom_id]
         atom_type = get_atomtype(atom, mol.get_bonds(atom))
@@ -836,6 +901,25 @@ def atom_type(self, mol, atom_id):
         else:
             return atom_type.label
 
+    def test_X_types(self):
+        assert self.atom_type(self.X, 0) == 'Xv'
+        assert self.atom_type(self.mol100, 1) == 'Xo'
+        assert self.atom_type(self.mol101, 1) == 'Xo'
+        assert self.atom_type(self.mol102, 1) == 'Xo'
+        assert self.atom_type(self.mol103, 1) == 'Xo'
+        assert self.atom_type(self.mol104, 3) == 'Xo'
+        assert self.atom_type(self.mol104, 4) == 'Xo'
+
+
+    def test_Pt_types(self):
+        assert self.atom_type(self.Pt, 0) == 'Ptv'
+        assert self.atom_type(self.mol110, 1) == 'Pto'
+        assert self.atom_type(self.mol111, 1) == 'Pto'
+        assert self.atom_type(self.mol112, 1) == 'Pto'
+        assert self.atom_type(self.mol113, 1) == 'Pto'
+        assert self.atom_type(self.mol114, 3) == 'Pto'
+        assert self.atom_type(self.mol114, 4) == 'Pto'
+
     def test_hydrogen_type(self):
         """
         Test that get_atomtype() returns the hydrogen atom type.

test/rmgpy/molecule/elementTest.py

@@ -44,6 +44,7 @@ class TestElement:
     def setup_class(self):
         self.element = rmgpy.molecule.element.C
         self.element_x = rmgpy.molecule.element.X
+        self.element_Pt = rmgpy.molecule.element.Pt
 
     def test_pickle(self):
         """
@@ -80,6 +81,7 @@ def test_get_element(self):
         assert rmgpy.molecule.element.get_element("C") is self.element
         assert rmgpy.molecule.element.get_element(0) is self.element_x
         assert rmgpy.molecule.element.get_element("X") is self.element_x
+        assert rmgpy.molecule.element.get_element("Pt") is self.element_Pt
 
     def test_get_element_isotope(self):
         """