diff --git a/.travis.yml b/.travis.yml
index d2ff12da4f..ff3c46baf9 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -9,6 +9,7 @@ before_install:
 #  - sudo apt-get install python-rdkit librdkit-dev librdkit1 rdkit-data
   - sudo apt-get install -qq python-numpy python-scipy python-matplotlib
   - cd ..
+
   - git clone https://github.com/GreenGroup/RMG-database.git
   - git clone https://github.com/GreenGroup/PyDAS.git
   - git clone https://github.com/GreenGroup/PyDQED.git
diff --git a/documentation/source/reference/molecule/adjlist.rst b/documentation/source/reference/molecule/adjlist.rst
index cbfa481d01..d1dd780bdd 100644
--- a/documentation/source/reference/molecule/adjlist.rst
+++ b/documentation/source/reference/molecule/adjlist.rst
@@ -7,6 +7,15 @@ Adjacency Lists
 
 .. module:: rmgpy.molecule.adjlist
 
+
+.. note::
+    The adjacency list syntax changed in July 2014.
+    The minimal requirement for most translations is to prefix the number
+    of unpaired electrons with the letter `u`.
+    The new syntax, however, allows much
+    greater flexibility, including definition of lone pairs, partial charges, 
+    wildcards, and molecule multiplicities.
+
 .. note::
     To quickly visualize any adjacency list, or to generate an adjacency list from
     other types of molecular representations such as SMILES, InChI, or even common
@@ -21,49 +30,100 @@ RMG -- but extended to allow for specification of extra semantic information.
 The first line of most adjacency lists is a unique identifier for the molecule
 or pattern the adjacency list represents. This is not strictly required, but
 is recommended in most cases. Generally the identifier should only use
-alphanumeric characters and the underscore, as if an identifer in many popular
+alphanumeric characters and the underscore, as if an identifier in many popular
 programming languages. However, strictly speaking any non-space ASCII character
 is allowed.
 
-After the identifier line, each subsequent line describes a single atom and its
+The subsequent lines may contain keyword-value pairs. Currently there is only
+one keyword, ``multiplicity``.
+
+For species or molecule declarations, the value after ``multiplicity`` defines 
+the spin multiplicity of the molecule. E.g. ``multiplicity 1`` for most ground state 
+closed shell species, ``multiplicity 2`` for most radical species, 
+and ``multiplicity 3`` for a triplet biradical.
+If the ``multiplicity`` line is not present then a value of 
+(1 + number of unpaired electrons) is assumed. 
+Thus, it can usually be omitted, but if present can be used to distinguish,
+for example, singlet CH2 from triplet CH2.
+
+If defining a Functional :class:`~rmgpy.molecule.Group`, then the value must be a list,
+which defines the multiplicities that will be matched by the group, eg. 
+``multiplicity [1,2,3,4,5]`` or, for a single value, ``multiplicity [1]``. 
+If the multiplicity line is omitted, then ``multiplicity [1,2,3,4,5]`` is assumed.
+
+After the identifier line and keyword-value lines,
+each subsequent line describes a single atom and its
 local bond structure. The format of these lines is a whitespace-delimited list
 with tokens ::
 
-    <number> [<label>] <element> <radicals> <bondlist>
+    <number> [<label>] <element> u<unpaired> [p<pairs>] [c<charge>] <bondlist>
 
 The first item is the number used to identify that atom. Any number may be used,
 though it is recommended to number the atoms sequentially starting from one.
 Next is an optional label used to tag that atom; this should be an
-asterisk followed by a unique number for the label, e.g. ``*1``. After that is
-the atom's element, indicated by its atomic symbol, followed by the number of
-radical electrons on the atom. The last set of tokens is the list of bonds.
+asterisk followed by a unique number for the label, e.g. ``*1``.
+In some cases (e.g. thermodynamics groups) there is only one labeled atom, and the label 
+is just an asterisk with no number: ``*``.
+
+After that is
+the atom's element or atom type, indicated by its atomic symbol, followed by 
+a sequence of tokens describing the electronic state of the atom:
+
+* ``u0`` number of **unpaired** electrons (eg. radicals)
+* ``p0`` number of lone **pairs** of electrons, common on oxygen and nitrogen.
+* ``c0`` formal **charge** on the atom, e.g. ``c-1`` (negatively charged),
+  ``c0``, ``c+1`` (positively charged)
+
+For :class:`~rmgpy.molecule.Molecule` definitions:
+The value must be a single integer (and for charge must have a + or - sign if not equal to 0)
+The number of unpaired electrons (i.e. radical electrons) is required, even if zero.
+The number of lone pairs and the formal charge are assumed to be zero if omitted.
+
+For :class:`~rmgpy.molecule.Group` definitions:
+The value can be an integer or a list of integers (with signs, for charges),
+ eg. ``u[0,1,2]`` or ``c[0,+1,+2,+3,+4]``, or may be a wildcard ``x`` 
+ which matches any valid value, 
+eg. ``px`` is the same as ``p[0,1,2,3,4]`` and ``cx`` is the same as
+``c[-4,-3,-2,-1,0,+1,+2,+3,+4]``. Lists must be enclosed is square brackets,
+and separated by commas, without spaces.
+If lone pairs or formal charges are omitted from a group definition,
+the wildcard is assumed.
+
+
+The last set of tokens is the list of bonds.
 To indicate a bond, place the number of the atom at the other end of the bond
 and the bond type within curly braces and separated by a comma, e.g. ``{2,S}``.
-Multiple bonds to the same atom should be separated by whitespace.
+Multiple bonds from the same atom should be separated by whitespace.
 
 .. note::
     You must take care to make sure each bond is listed on the lines of *both*
     atoms in the bond, and that these entries have the same bond type. RMG will
     raise an exception if it encounters such an invalid adjacency list.
 
+
 When writing a molecular substructure pattern, you may specify multiple 
-elements, radical counts, and bond types as a comma-separated list inside curly
-braces. For example, to specify any carbon or oxygen atom, use the syntax 
-``{C,O}``. Atom types may also be used as a shorthand. (Atom types can also be
+elements, radical counts, and bond types as a comma-separated list inside square
+brackets. For example, to specify any carbon or oxygen atom, use the syntax 
+``[C,O]``. For a single or double bond to atom 2, write ``{2,[S,D]}``.
+
+Atom types such as ``R!H`` or ``Cdd`` may also be used as a shorthand. (Atom types 
+like ``Cdd`` can also be
 used in full molecules, but this use is discouraged, as RMG can compute them
 automatically for full molecules.)
 
 Below is an example adjacency list, for 1,3-hexadiene, with the weakest bond in
 the molecule labeled with ``*1`` and ``*2``. Note that hydrogen atoms
-can be omitted if desired, as their presence is inferred::
+can be omitted if desired, as their presence is inferred, provided that unpaired
+electrons, lone pairs, and charges are all correctly defined::
 
     HXD13
-    1    C 0       {2,D}
-    2    C 0 {1,D} {3,S}
-    3    C 0 {2,S} {4,D}
-    4    C 0 {3,D} {5,S}
-    5 *1 C 0 {4,S} {6,S}
-    6 *2 C 0 {5,S}
+    multiplicity 1
+    1    C u0       {2,D}
+    2    C u0 {1,D} {3,S}
+    3    C u0 {2,S} {4,D}
+    4    C u0 {3,D} {5,S}
+    5 *1 C u0 {4,S} {6,S}
+    6 *2 C u0 {5,S}
     
 The allowed element types, radicals, and bonds are listed in the following table:
 
@@ -77,10 +137,10 @@ The allowed element types, radicals, and bonds are listed in the following table
  |                      | H        | Hydrogen atom       |
  |                      +----------+---------------------+
  |                      | S        | Sulfur atom         |
- +----------------------+----------+---------------------+
- | Nonreactive Elements | N        | Nitrogen atom       |
  |                      +----------+---------------------+
- |                      | Si       | Silicon atom        |
+ |                      | N        | Nitrogen atom       |
+ +----------------------+----------+---------------------+
+ | Nonreactive Elements | Si       | Silicon atom        |
  |                      +----------+---------------------+
  |                      | Cl       | Chlorine atom       |
  |                      +----------+---------------------+
@@ -89,18 +149,6 @@ The allowed element types, radicals, and bonds are listed in the following table
  |                      | Ar       | Argon atom          |
  |                      +----------+---------------------+
  +----------------------+----------+---------------------+
- | Free Electrons       | 0        | Non-radical         |
- |                      +----------+---------------------+
- |                      | 1        | Mono-radical        |
- |                      +----------+---------------------+
- |                      | 2        | Bi-radical          |
- |                      +----------+---------------------+
- |                      | 2T       | Triplet             |
- |                      +----------+---------------------+
- |                      | 2S       | Singlet             |
- |                      +----------+---------------------+
- |                      | 3        | Tri-radical         |
- +----------------------+----------+---------------------+
  | Chemical Bond        | S        | Single Bond         |
  |                      +----------+---------------------+
  |                      | D        | Double Bond         |
diff --git a/documentation/source/users/rmg/faq.rst b/documentation/source/users/rmg/faq.rst
index 25cf0e136b..1f880c4854 100644
--- a/documentation/source/users/rmg/faq.rst
+++ b/documentation/source/users/rmg/faq.rst
@@ -3,3 +3,37 @@
 **************************
 Frequently Asked Questions
 **************************
+
+
+Why can't my adjacency lists be read any more?
+==============================================
+
+The adjacency list syntax changed in July 2014.
+The minimal requirement for most translations is to prefix the number
+of unpaired electrons with the letter `u`.
+
+Example old syntax::
+
+    HXD13
+    1    C 0       {2,D}
+    2    C 0 {1,D} {3,S}
+    3    C 0 {2,S} {4,D}
+    4    C 0 {3,D} {5,S}
+    5 *1 C 0 {4,S} {6,S}
+    6 *2 C 0 {5,S}
+
+Example new syntax::
+
+    HXD13
+    1    C u0       {2,D}
+    2    C u0 {1,D} {3,S}
+    3    C u0 {2,S} {4,D}
+    4    C u0 {3,D} {5,S}
+    5 *1 C u0 {4,S} {6,S}
+    6 *2 C u0 {5,S}
+    
+The new syntax, however, allows much
+greater flexibility, including definition of lone pairs, partial charges, 
+wildcards, and molecule multiplicities, and was necessary to allow us to 
+add Nitrogen chemistry.
+See :ref:`rmgpy.molecule.adjlist` for details of the new syntax.
\ No newline at end of file
diff --git a/examples/rmg/1,3-hexadiene/input.py b/examples/rmg/1,3-hexadiene/input.py
index 28615b51a8..a77f17382f 100644
--- a/examples/rmg/1,3-hexadiene/input.py
+++ b/examples/rmg/1,3-hexadiene/input.py
@@ -8,6 +8,11 @@
     kineticsEstimator = 'rate rules',
 )
 
+# Constraints on generated species
+generatedSpeciesConstraints(
+    maximumRadicalElectrons = 2,
+)
+
 # List of species
 species(
     label='HXD13',
@@ -24,8 +29,8 @@
     reactive=True,
     structure=adjacencyList(
         """
-        1 H 0 {2,S}
-        2 H 0 {1,S}
+        1 H u0 p0 {2,S}
+        2 H u0 p0 {1,S}
         """),
 )
 species(
diff --git a/examples/rmg/c3h4/input.py b/examples/rmg/c3h4/input.py
index ffc50b4d3b..258b7c7f83 100644
--- a/examples/rmg/c3h4/input.py
+++ b/examples/rmg/c3h4/input.py
@@ -8,6 +8,10 @@
 	kineticsEstimator = 'rate rules',
 )
 
+generatedSpeciesConstraints(
+    maximumRadicalElectrons = 4,
+)
+
 # List of species
 species(
 	label='CH2',
diff --git a/examples/rmg/ch3no2/input.py b/examples/rmg/ch3no2/input.py
index 749cd7d7e5..f40eba031f 100644
--- a/examples/rmg/ch3no2/input.py
+++ b/examples/rmg/ch3no2/input.py
@@ -26,13 +26,13 @@
     reactive=True,
         structure=adjacencyList(
         """
-        1 C 0 0 {2,S} {3,S} {4,S} {5,S}
-        2 H 0 0 {1,S}
-        3 H 0 0 {1,S}
-        4 H 0 0 {1,S}
-        5 N 0 0 {1,S} {6,D} {7,S}
-        6 O 0 2 {5,D}
-        7 O 0 3 {5,S}
+        1 C u0 p0 {2,S} {3,S} {4,S} {5,S}
+        2 H u0 p0 {1,S}
+        3 H u0 p0 {1,S}
+        4 H u0 p0 {1,S}
+        5 N u0 p0 {1,S} {6,D} {7,S}
+        6 O u0 p2 {5,D}
+        7 O u0 p3 {5,S}
         """),
 )
 
@@ -41,8 +41,8 @@
     reactive=True,
         structure=adjacencyList(
         """
-        1 O 1 2 {2,S}
-        2 O 1 2 {1,S}
+        1 O u1 p2 {2,S}
+        2 O u1 p2 {1,S}
         """),
 )
 
@@ -51,8 +51,8 @@
     reactive=True,
         structure=adjacencyList(
         """
-        1 N 1 1 {2,T}
-        2 N 1 1 {1,T}
+        1 N u0 p1 {2,T}
+        2 N u0 p1 {1,T}
         """),
 )
 
diff --git a/examples/rmg/diesel/input.py b/examples/rmg/diesel/input.py
index 1b3a010db4..cfbd158f68 100644
--- a/examples/rmg/diesel/input.py
+++ b/examples/rmg/diesel/input.py
@@ -47,7 +47,7 @@
 species(
     label='O2',
     reactive=True,
-    structure=SMILES("O=O"),
+    structure=SMILES("[O][O]"),
 )
 
 # Reaction systems
diff --git a/examples/rmg/e85/input.py b/examples/rmg/e85/input.py
index 6ffb04c275..d1b478e443 100644
--- a/examples/rmg/e85/input.py
+++ b/examples/rmg/e85/input.py
@@ -8,11 +8,16 @@
     kineticsEstimator = 'rate rules',
 )
 
+# Constraints on generated species
+generatedSpeciesConstraints(
+    maximumRadicalElectrons = 2,
+)
+
 # List of species
 species(
     label='O2',     # oxygen
     reactive=True,
-    structure=SMILES("O=O"),
+    structure=SMILES("[O][O]"),
 )
 species(
     label='C8H18i', # isooctane
diff --git a/examples/rmg/liquid_phase/input.py b/examples/rmg/liquid_phase/input.py
index 1979dcc9d7..9ec6c98756 100644
--- a/examples/rmg/liquid_phase/input.py
+++ b/examples/rmg/liquid_phase/input.py
@@ -8,6 +8,11 @@
     kineticsEstimator = 'rate rules',
 )
 
+# Constraints on generated species
+generatedSpeciesConstraints(
+    maximumRadicalElectrons = 3,
+)
+
 # List of species
 species(
     label='octane',
diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py
index ccd3989236..07e68c26cc 100644
--- a/examples/rmg/methylformate/input.py
+++ b/examples/rmg/methylformate/input.py
@@ -29,8 +29,8 @@
     reactive=True,
     structure=adjacencyList(
         """
-        1     C     3 {2,S}
-        2     H     0 {1,S}
+        1 C u3 p0 {2,S}
+        2 H u0 p0 {1,S}
         """),
 )
 species(
@@ -46,7 +46,7 @@
 species(
     label='CO',
     reactive=True,
-    structure=SMILES("[C]=O"),
+    structure=SMILES("[C+]#[O-]"),
 )
 species(
     label='CO2',
@@ -157,4 +157,5 @@
     drawMolecules=False,
     generatePlots=False,
     saveConcentrationProfiles=False,
+    saveEdgeSpecies=True,
 )
diff --git a/examples/rmg/minimal_sensitivity/input.py b/examples/rmg/minimal_sensitivity/input.py
index c1c57f0936..85840e1eeb 100644
--- a/examples/rmg/minimal_sensitivity/input.py
+++ b/examples/rmg/minimal_sensitivity/input.py
@@ -8,6 +8,11 @@
     kineticsEstimator = 'rate rules',
 )
 
+# Constraints on generated species
+generatedSpeciesConstraints(
+    maximumRadicalElectrons = 2,
+)
+
 # List of species
 species(
     label='ethane',
diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py
index 1bf6fc0b1e..ae091512cc 100644
--- a/rmgpy/data/base.py
+++ b/rmgpy/data/base.py
@@ -1333,7 +1333,6 @@ def saveEntry(self, f, entry, name='entry'):
         
         f.write('{0}(\n'.format(name))
         f.write('    label = "{0}",\n'.format(entry.label))
-
         if isinstance(entry.item, Molecule):
             f.write('    molecule = \n')
             f.write('"""\n')
diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 511a047ac8..2e9f06aa21 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -841,7 +841,7 @@ def addKineticsRulesFromTrainingSet(self, thermoDatabase=None):
         reverse_entries = []
         for entry in entries:
             try:        
-                template = self.getReactionTemplate(deepcopy(entry.item))
+                template = self.getReactionTemplate(entry.item.copy())
             except UndeterminableKineticsError:
                 # Some entries might be stored in the reverse direction for
                 # this family; save them so we can try this
@@ -1141,172 +1141,79 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp
 
         
         # Generate other possible electronic states
-        electronicStructuresList1 = []
-        electronicStructuresList2 = []
-        
-        struct1 = productStructures[0]
-        struct1a = struct1.copy(True)
-        struct1a.updateAtomTypes()
-        electronicStructuresList1.append(struct1a)
-        atoms1 = struct1.getRadicalAtoms()
+        productStructuresList = []
+        totalSpin = [] # total spin times 2
         
-        for atom1 in atoms1:
+        # implement Angular Momentum Addition Theorem
+        if len(reactantStructures) == 1:
             
-            radical1 = atom1.radicalElectrons
-            spin1 = atom1.spinMultiplicity
+            totalSpin = [(reactantStructures[0].multiplicity-1.0)/2.0]
             
-            if atom1.label != '' and radical1 > 1 and radical1 < 4:
-                
-                if radical1 == 2 and spin1 == 3:
-                    atom1.setSpinMultiplicity(1)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
-                elif radical1 == 2 and spin1 == 1:
-                    atom1.setSpinMultiplicity(3)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
-                elif radical1 == 3 and spin1 == 4:
-                    atom1.setSpinMultiplicity(2)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
-                elif radical1 == 3 and spin1 == 2:
-                    atom1.setSpinMultiplicity(4)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
-                
-                for electronicStructures in electronicStructuresList1:
-                    if electronicStructures.isIsomorphic(struct1a):
-                        break
-                else:
-                    electronicStructuresList1.append(struct1a)
+        elif len(reactantStructures) == 2:
             
-            elif radical1 == 4:
+            spin1 = (reactantStructures[0].multiplicity-1.0)/2.0
+            spin2 = (reactantStructures[1].multiplicity-1.0)/2.0
+            
+            count = 0.0
+            
+            while (spin1+spin2-count) >= abs(spin1-spin2):
                 
-                if spin1 == 5:
-                    atom1.setSpinMultiplicity(3)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
+                totalSpin.append(spin1+spin2-count)
+                                
+                count += 1
+            
+        if len(productStructures) == 1:
+            
+            maxSpin1 = productStructures[0].getRadicalCount()/2.0
+            
+            count = 0.0
+            
+            while (maxSpin1-count) >= 0.0:
                 
-                    atom1.setSpinMultiplicity(1)
-                    struct1b = struct1.copy(True)
-                    struct1b.updateAtomTypes()
-                elif spin1 == 3:
-                    atom1.setSpinMultiplicity(5)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
+                if (maxSpin1-count) in totalSpin:
                 
-                    atom1.setSpinMultiplicity(1)
-                    struct1b = struct1.copy(True)
-                    struct1b.updateAtomTypes()
-                elif spin1 == 1:
-                    atom1.setSpinMultiplicity(5)
-                    struct1a = struct1.copy(True)
-                    struct1a.updateAtomTypes()
+                    struct = productStructures[0].copy(deep=True)
                 
-                    atom1.setSpinMultiplicity(3)
-                    struct1b = struct1.copy(True)
-                    struct1b.updateAtomTypes()
+                    struct.multiplicity = int((maxSpin1-count)*2.0+1.0)
                     
-                for electronicStructures in electronicStructuresList1:
-                    if electronicStructures.isIsomorphic(struct1a):
-                        break
-                else:
-                    electronicStructuresList1.append(struct1a)
+                    if not self.isMoleculeForbidden(struct):
+                        productStructuresList.append([struct])
                     
-                for electronicStructures in electronicStructuresList1:
-                    if electronicStructures.isIsomorphic(struct1b):
-                        break
-                else:
-                    electronicStructuresList1.append(struct1b)
-                            
-        if len(productStructures) == 2:
-        
-            struct2 = productStructures[1]
-            struct2a = struct2.copy(True)
-            struct2a.updateAtomTypes()
-            electronicStructuresList2.append(struct2a)
-            atoms2 = struct2.getRadicalAtoms()
-        
-            for atom2 in atoms2:
+                count += 1.0
+                    
+        elif len(productStructures) == 2:
             
-                radical2 = atom2.radicalElectrons
-                spin2 = atom2.spinMultiplicity
+            maxSpin1 = productStructures[0].getRadicalCount()/2.0
+            maxSpin2 = productStructures[1].getRadicalCount()/2.0
             
-                if atom2.label != '' and radical2 > 1 and radical2 < 4:
-                    
-                    if radical2 == 2 and spin2 == 3:
-                        atom2.setSpinMultiplicity(1)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
-                    elif radical2 == 2 and spin2 == 1:
-                        atom2.setSpinMultiplicity(3)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
-                    elif radical2 == 3 and spin2 == 4:
-                        atom2.setSpinMultiplicity(2)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
-                    elif radical2 == 3 and spin2 == 2:
-                        atom2.setSpinMultiplicity(4)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
-                
-                    for electronicStructures in electronicStructuresList2:
-                        if electronicStructures.isIsomorphic(struct2a):
-                            break
-                    else:
-                        electronicStructuresList2.append(struct2a)
+            count1 = 0.0
             
-                elif radical2 == 4:
+            while (maxSpin1-count1) >= 0.0:
                 
-                    if spin2 == 5:
-                        atom2.setSpinMultiplicity(3)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
+                count2 = 0.0
                 
-                        atom2.setSpinMultiplicity(1)
-                        struct2b = struct2.copy(True)
-                        struct2b.updateAtomTypes()
-                    elif spin2 == 3:
-                        atom2.setSpinMultiplicity(5)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
+                while (maxSpin2-count2) >= 0.0:
+                    
+                    count = 0.0
                 
-                        atom2.setSpinMultiplicity(1)
-                        struct2b = struct2.copy(True)
-                        struct2b.updateAtomTypes()
-                    elif spin2 == 1:
-                        atom2.setSpinMultiplicity(5)
-                        struct2a = struct2.copy(True)
-                        struct2a.updateAtomTypes()
+                    while (maxSpin1-count1+maxSpin2-count2-count) >= abs((maxSpin1-count1)-(maxSpin2-count2)):
+                        
+                        if (maxSpin1-count1+maxSpin2-count2-count) in totalSpin:
                 
-                        atom2.setSpinMultiplicity(3)
-                        struct2b = struct2.copy(True)
-                        struct2b.updateAtomTypes()
-                    
-                    for electronicStructures in electronicStructuresList2:
-                        if electronicStructures.isIsomorphic(struct2a):
-                            break
-                    else:
-                        electronicStructuresList2.append(struct2a)
+                            struct1 = productStructures[0].copy(deep=True)
+                            struct2 = productStructures[1].copy(deep=True)
+                
+                            struct1.multiplicity = int((maxSpin1-count1)*2.0+1.0)
+                            struct2.multiplicity = int((maxSpin2-count2)*2.0+1.0)
+                                
+                            if not self.isMoleculeForbidden(struct1) and not self.isMoleculeForbidden(struct2):
+                                productStructuresList.append([struct1,struct2])
+                        
+                        count += 1.0
+                        
+                    count2 += 1.0
                     
-                    for electronicStructures in electronicStructuresList2:
-                        if electronicStructures.isIsomorphic(struct2b):
-                            break
-                    else:
-                        electronicStructuresList2.append(struct2b)
-        
-        if len(productStructures) == 2:
-            
-            for structa in electronicStructuresList1:
-                for structb in electronicStructuresList2:
-                    if not (self.isMoleculeForbidden(structa) or self.isMoleculeForbidden(structb)):
-                        productStructuresList.append([structa,structb])
-        elif len(productStructures) == 1:
-            
-            for structa in electronicStructuresList1:
-                if not (self.isMoleculeForbidden(structa)):
-                    productStructuresList.append([structa])
+                count1 += 1.0
                     
         return productStructuresList
 
diff --git a/rmgpy/data/kinetics/groups.py b/rmgpy/data/kinetics/groups.py
index 3662c1f9b9..0460d820c8 100644
--- a/rmgpy/data/kinetics/groups.py
+++ b/rmgpy/data/kinetics/groups.py
@@ -159,10 +159,14 @@ def getReactionTemplate(self, reaction):
         # template is a list of the actual matched nodes
         # forwardTemplate is a list of the top level nodes that should be matched
         if len(template) != len(forwardTemplate):
+#            print 'len(template):', len(template)
+#            print 'len(forwardTemplate):', len(forwardTemplate)
             msg = 'Unable to find matching template for reaction {0} in reaction family {1}.'.format(str(reaction), str(self)) 
             msg += 'Trying to match {0} but matched {1}'.format(str(forwardTemplate),str(template))
+#            print 'reactants'
 #            for reactant in reaction.reactants:
 #                print reactant.toAdjacencyList() + '\n'
+#            print 'products'
 #            for product in reaction.products:
 #                print product.toAdjacencyList() + '\n'
             raise UndeterminableKineticsError(reaction, message=msg)
diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py
index d1c0a5c0ca..942dcce437 100644
--- a/rmgpy/data/kinetics/library.py
+++ b/rmgpy/data/kinetics/library.py
@@ -274,6 +274,7 @@ def loadEntry(self,
                   longDesc='',
                   ):
         
+        
         reactants = [Species(label=reactant1.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(reactant1)])]
         if reactant2 is not None: reactants.append(Species(label=reactant2.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(reactant2)]))
         if reactant3 is not None: reactants.append(Species(label=reactant3.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(reactant3)]))
diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py
index b11d6aa549..7103fa1812 100644
--- a/rmgpy/data/solvation.py
+++ b/rmgpy/data/solvation.py
@@ -726,6 +726,8 @@ def estimateSoluteViaGroupAdditivity(self, molecule):
             saturatedStruct.updateConnectivityValues()
             saturatedStruct.sortVertices()
             saturatedStruct.updateAtomTypes()
+            saturatedStruct.updateLonePairs()
+            saturatedStruct.multiplicity = 1
 
             # Get solute descriptor estimates for saturated form of structure
             soluteData = self.estimateSoluteViaGroupAdditivity(saturatedStruct)
diff --git a/rmgpy/data/statmech.py b/rmgpy/data/statmech.py
index cf4142f7d1..b8f37ef53b 100644
--- a/rmgpy/data/statmech.py
+++ b/rmgpy/data/statmech.py
@@ -131,10 +131,13 @@ def loadEntry(self,
                   shortDesc='',
                   longDesc='',
                   ):
+        
+        item = Molecule().fromAdjacencyList(molecule)
+        
         self.entries[label] = Entry(
             index = index,
             label = label,
-            item = Molecule().fromAdjacencyList(molecule),
+            item = item,
             data = statmech,
             reference = reference,
             referenceType = referenceType,
diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py
index c749b6575e..6d744aed56 100644
--- a/rmgpy/data/thermo.py
+++ b/rmgpy/data/thermo.py
@@ -54,8 +54,8 @@ def saveEntry(f, entry):
     """
     
     f.write('entry(\n')
-    f.write('    index = {0:d},\n'.format(entry.index))
-    f.write('    label = "{0}",\n'.format(entry.label))
+    f.write('    index        = {0:d},\n'.format(entry.index))
+    f.write('    label        = "{0}",\n'.format(entry.label))
 
     if isinstance(entry.item, Molecule):
         f.write('    molecule = \n')
@@ -233,7 +233,8 @@ def loadEntry(self,
         
         for entry in self.entries.values():
             if molecule.isIsomorphic(entry.item):
-                raise DatabaseError('Adjacency list of {0} matches that of existing molecule {1} in thermo library.  Please correct your library.'.format(label, entry.label))
+                if molecule.multiplicity == entry.item.multiplicity:
+                    raise DatabaseError('Adjacency list and multiplicity of {0} matches that of existing molecule {1} in thermo library.  Please correct your library.'.format(label, entry.label))
         
         self.entries[label] = Entry(
             index = index,
@@ -810,6 +811,8 @@ def estimateRadicalThermoViaHBI(self, molecule, stableThermoEstimator ):
         saturatedStruct.updateConnectivityValues()
         saturatedStruct.sortVertices()
         saturatedStruct.updateAtomTypes()
+        saturatedStruct.updateLonePairs()
+        saturatedStruct.multiplicity = 1
         
         # Get thermo estimate for saturated form of structure
         thermoData = stableThermoEstimator(saturatedStruct)
diff --git a/rmgpy/data/thermoTest.py b/rmgpy/data/thermoTest.py
index 374ff2de7c..5f80b3b27f 100644
--- a/rmgpy/data/thermoTest.py
+++ b/rmgpy/data/thermoTest.py
@@ -19,8 +19,8 @@ class TestThermoDatabase(unittest.TestCase):
     # Only load these once to save time
     database = ThermoDatabase()
     database.load(os.path.join(settings['database.directory'], 'thermo'))
-    oldDatabase = ThermoDatabase()
-    oldDatabase.loadOld(os.path.join(settings['database.directory'], '../output/RMG_database'))
+#    oldDatabase = ThermoDatabase()
+#    oldDatabase.loadOld(os.path.join(settings['database.directory'], '../output/RMG_database'))
     
     
     def setUp(self):
@@ -29,7 +29,7 @@ def setUp(self):
         """
         
         self.database = self.__class__.database
-        self.oldDatabase = self.__class__.oldDatabase
+#        self.oldDatabase = self.__class__.oldDatabase
 
         self.Tlist = [300, 400, 500, 600, 800, 1000, 1500]
         
@@ -65,10 +65,11 @@ def testNewThermoGeneration(self):
         
         for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases:
             Cplist = [Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500]
-            species = Species(molecule=[Molecule(SMILES=smiles)])
+            molecule=Molecule(SMILES=smiles)
+            species = Species(molecule=molecule)
             species.generateResonanceIsomers()
             species.molecule[0]
-            thermoData = self.database.getThermoDataFromGroups(Species(molecule=[species.molecule[0]]))
+            thermoData = self.database.getThermoDataFromGroups(species)
             molecule = species.molecule[0]
             for mol in species.molecule[1:]:
                 thermoData0 = self.database.getAllThermoData(Species(molecule=[mol]))[0][0]
@@ -92,7 +93,8 @@ def testSymmetryNumberGeneration(self):
         to select the stablest resonance isomer.
         """
         for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases:
-            species = Species(molecule=[Molecule(SMILES=smiles)])
+            molecule=Molecule(SMILES=smiles)
+            species = Species(molecule=molecule)
             species.generateResonanceIsomers()
             thermoData = self.database.getThermoDataFromGroups(Species(molecule=[species.molecule[0]]))
             # pick the molecule with lowest H298
@@ -107,30 +109,30 @@ def testSymmetryNumberGeneration(self):
                     molecule = mol
             self.assertEqual(molecule.calculateSymmetryNumber(), symm, msg="Symmetry number error for {0}".format(smiles))
 
-    @work_in_progress
-    def testOldThermoGeneration(self):
-        """
-        Test that the old ThermoDatabase generates relatively accurate thermo data.
-        """
-        for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases:
-            Cplist = [Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500]
-            species = Species(molecule=[Molecule(SMILES=smiles)])
-            species.generateResonanceIsomers()
-            thermoData = self.oldDatabase.getThermoData(Species(molecule=[species.molecule[0]]))
-            molecule = species.molecule[0]
-            for mol in species.molecule[1:]:
-                thermoData0 = self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[0][0]
-                for data in self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[1:]:
-                    if data.getEnthalpy(298) < thermoData0.getEnthalpy(298):
-                        thermoData0 = data
-                if thermoData0.getEnthalpy(298) < thermoData.getEnthalpy(298):
-                    thermoData = thermoData0
-                    molecule = mol
-            
-            self.assertAlmostEqual(H298, thermoData.getEnthalpy(298) / 4184, places=1, msg="H298 error for {0}".format(smiles))
-            self.assertAlmostEqual(S298, thermoData.getEntropy(298) / 4.184, places=1, msg="S298 error for {0}".format(smiles))
-            for T, Cp in zip(self.Tlist, Cplist):
-                self.assertAlmostEqual(Cp, thermoData.getHeatCapacity(T) / 4.184, places=1, msg="Cp{1} error for {0}".format(smiles, T))
+#    @work_in_progress
+#    def testOldThermoGeneration(self):
+#        """
+#        Test that the old ThermoDatabase generates relatively accurate thermo data.
+#        """
+#        for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases:
+#            Cplist = [Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500]
+#            species = Species(molecule=[Molecule(SMILES=smiles)])
+#            species.generateResonanceIsomers()
+#            thermoData = self.oldDatabase.getThermoData(Species(molecule=[species.molecule[0]]))
+#            molecule = species.molecule[0]
+#            for mol in species.molecule[1:]:
+#                thermoData0 = self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[0][0]
+#                for data in self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[1:]:
+#                    if data.getEnthalpy(298) < thermoData0.getEnthalpy(298):
+#                        thermoData0 = data
+#                if thermoData0.getEnthalpy(298) < thermoData.getEnthalpy(298):
+#                    thermoData = thermoData0
+#                    molecule = mol
+#            
+#            self.assertAlmostEqual(H298, thermoData.getEnthalpy(298) / 4184, places=1, msg="H298 error for {0}".format(smiles))
+#            self.assertAlmostEqual(S298, thermoData.getEntropy(298) / 4.184, places=1, msg="S298 error for {0}".format(smiles))
+#            for T, Cp in zip(self.Tlist, Cplist):
+#                self.assertAlmostEqual(Cp, thermoData.getHeatCapacity(T) / 4.184, places=1, msg="Cp{1} error for {0}".format(smiles, T))
 
 class TestThermoDatabaseAromatics(TestThermoDatabase):
     """
diff --git a/rmgpy/data/transport.py b/rmgpy/data/transport.py
index 8c1b2d1c73..72b65d94a0 100644
--- a/rmgpy/data/transport.py
+++ b/rmgpy/data/transport.py
@@ -52,8 +52,8 @@ def saveEntry(f, entry):
     database to the file object `f`.
     """
     f.write('entry(\n')
-    f.write('    index = {0:d},\n'.format(entry.index))
-    f.write('    label = "{0}",\n'.format(entry.label))
+    f.write('    index        = {0:d},\n'.format(entry.index))
+    f.write('    label        = "{0}",\n'.format(entry.label))
 
     if isinstance(entry.item, Molecule):
         f.write('    molecule = \n')
@@ -142,10 +142,13 @@ def loadEntry(self,
                   shortDesc='',
                   longDesc='',
                   ):
+        
+        item = Molecule().fromAdjacencyList(molecule)
+        
         self.entries[label] = Entry(
             index = index,
             label = label,
-            item = Molecule().fromAdjacencyList(molecule),
+            item = item,
             data = transport,
             reference = reference,
             referenceType = referenceType,
diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py
index 01cc58da1e..2f899fa979 100644
--- a/rmgpy/molecule/adjlist.py
+++ b/rmgpy/molecule/adjlist.py
@@ -32,6 +32,7 @@
 adjacency list format used by Reaction Mechanism Generator (RMG).
 """
 
+import logging
 import re
 from .molecule import Atom, Bond
 from .group import GroupAtom, GroupBond
@@ -48,17 +49,18 @@ class InvalidAdjacencyListError(Exception):
 
 ################################################################################
 
-def fromAdjacencyList(adjlist, group=False, saturateH=False):
+def fromOldAdjacencyList(adjlist, group=False, saturateH=False):
     """
-    Convert a string adjacency list `adjlist` into a set of :class:`Atom` and
-    :class:`Bond` objects.
+    Convert a pre-June-2014 string adjacency list `adjlist` into a set of :class:`Atom` and
+    :class:`Bond` objects. 
+    It can read both "old style" that existed for years, an the "intermediate style" that
+    existed for a few months in 2014, with the extra column of integers for lone pairs.
     """
     atoms = []
     atomdict = {}
     bonds = {}
     
     try:
-        
         adjlist = adjlist.strip()
         lines = adjlist.splitlines()
         if adjlist == '' or len(lines) == 0:
@@ -71,6 +73,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
                 raise InvalidAdjacencyListError('No atoms specified in adjacency list.')
         
         mistake1 = re.compile('\{[^}]*\s+[^}]*\}')
+        atomicMultiplicities = {} # these are no longer stored on atoms, so we make a separate dictionary
         # Iterate over the remaining lines, generating Atom or GroupAtom objects
         for line in lines:
 
@@ -111,7 +114,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
             index += 1
             
             # Next is the electron state
-            radicalElectrons = []; spinMultiplicity = []
+            radicalElectrons = []; atomSpinMultiplicity = []
             elecState = data[index].upper()
             if elecState[0] == '{':
                 elecState = elecState[1:-1].split(',')
@@ -119,62 +122,419 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
                 elecState = [elecState]
             for e in elecState:
                 if e == '0':
-                    radicalElectrons.append(0); spinMultiplicity.append(1)
+                    radicalElectrons.append(0); atomSpinMultiplicity.append(1)
                 elif e == '1':
-                    radicalElectrons.append(1); spinMultiplicity.append(2)
+                    radicalElectrons.append(1); atomSpinMultiplicity.append(2)
                 elif e == '2':
-                    radicalElectrons.append(2); spinMultiplicity.append(1)
-                    radicalElectrons.append(2); spinMultiplicity.append(3)
+                    radicalElectrons.append(2); atomSpinMultiplicity.append(1)
+                    radicalElectrons.append(2); atomSpinMultiplicity.append(3)
                 elif e == '2S':
-                    radicalElectrons.append(2); spinMultiplicity.append(1)
+                    radicalElectrons.append(2); atomSpinMultiplicity.append(1)
                 elif e == '2T':
-                    radicalElectrons.append(2); spinMultiplicity.append(3)
+                    radicalElectrons.append(2); atomSpinMultiplicity.append(3)
                 elif e == '3':
-                    radicalElectrons.append(3); spinMultiplicity.append(4)
+                    radicalElectrons.append(3); atomSpinMultiplicity.append(4)
                 elif e == '3D':
-                    radicalElectrons.append(3); spinMultiplicity.append(2)
+                    radicalElectrons.append(3); atomSpinMultiplicity.append(2)
                 elif e == '3Q':
-                    radicalElectrons.append(3); spinMultiplicity.append(4)
+                    radicalElectrons.append(3); atomSpinMultiplicity.append(4)
                 elif e == '4':
-                    radicalElectrons.append(4); spinMultiplicity.append(5)
+                    radicalElectrons.append(4); atomSpinMultiplicity.append(5)
                 elif e == '4S':
-                    radicalElectrons.append(4); spinMultiplicity.append(1)
+                    radicalElectrons.append(4); atomSpinMultiplicity.append(1)
                 elif e == '4T':
-                    radicalElectrons.append(4); spinMultiplicity.append(3)
+                    radicalElectrons.append(4); atomSpinMultiplicity.append(3)
                 elif e == '4V':
-                    radicalElectrons.append(4); spinMultiplicity.append(5)
+                    radicalElectrons.append(4); atomSpinMultiplicity.append(5)
                 elif e == 'X':
                     radicalElectrons.extend([0,1,2,2])
-                    spinMultiplicity.extend([1,2,1,3])
+                    atomSpinMultiplicity.extend([1,2,1,3])
             index += 1
             
             # Next number defines the number of lone electron pairs (if provided)
-            lonePairElectrons = -1
-         #   if not group and len(data) > index:
+            lonePairsOfElectrons = -1
             if len(data) > index:
                 lpState = data[index]
-                if lpState[0] != '{':
+                if lpState[0] == '{':
+                    # this is the start of the chemical bonds - no lone pair info was provided
+                    lonePairsOfElectrons = -1
+                else:
                     if lpState == '0':
-                        lonePairElectrons = 0
+                        lonePairsOfElectrons = 0
                     if lpState == '1':
-                        lonePairElectrons = 1
+                        lonePairsOfElectrons = 1
                     if lpState == '2':
-                        lonePairElectrons = 2
+                        lonePairsOfElectrons = 2
                     if lpState == '3':
-                        lonePairElectrons = 3
+                        lonePairsOfElectrons = 3
                     if lpState == '4':
-                        lonePairElectrons = 4
+                        lonePairsOfElectrons = 4
                     index += 1
+            else: # no bonds or lone pair info provided.
+                lonePairsOfElectrons = -1
+            
+            # Create a new atom based on the above information
+            if group:
+                # charge currently not allowed
+                atom = GroupAtom(atomType=atomType,
+                                 radicalElectrons=sorted(set(radicalElectrons)),
+                                 charge=[0],
+                                 label=label,
+                                 lonePairs=(None if lonePairsOfElectrons==-1 else [lonePairsOfElectrons])
+                                 )
+            else:
+                atom = Atom(element=atomType[0],
+                            radicalElectrons=radicalElectrons[0],
+                            charge=0,
+                            label=label,
+                            lonePairs=lonePairsOfElectrons
+                            )
+
+            atomicMultiplicities[atom] = atomSpinMultiplicity
+            # Add the atom to the list
+            atoms.append(atom)
+            atomdict[aid] = atom
+            
+            # Process list of bonds
+            bonds[aid] = {}
+            for datum in data[index:]:
+
+                # Sometimes commas are used to delimit bonds in the bond list,
+                # so strip them just in case
+                datum = datum.strip(',')
+                
+                aid2, comma, order = datum[1:-1].partition(',')
+                aid2 = int(aid2)
+                if aid == aid2:
+                    raise InvalidAdjacencyListError('Attempted to create a bond between atom {0:d} and itself.'.format(aid))
+                
+                if order[0] == '{':
+                    order = order[1:-1].split(',')
                 else:
-                    lonePairElectrons = -1
+                    order = [order]
+
+                bonds[aid][aid2] = order
+
+        if group:
+            multiplicity = None
+        else:
+            multiplicity = 1
+            for atom in atoms:
+                multiplicity += max(atomicMultiplicities[atom]) - 1
+
+        # Check consistency using bonddict
+        for atom1 in bonds:
+            for atom2 in bonds[atom1]:
+                if atom2 not in bonds:
+                    raise InvalidAdjacencyListError('Atom {0:d} not in bond dictionary.'.format(atom2))
+                elif atom1 not in bonds[atom2]:
+                    raise InvalidAdjacencyListError('Found bond between {0:d} and {1:d}, but not the reverse.'.format(atom1, atom2))
+                elif bonds[atom1][atom2] != bonds[atom2][atom1]:
+                    raise InvalidAdjacencyListError('Found bonds between {0:d} and {1:d}, but of different orders "{2}" and "{3}".'.format(atom1, atom2, bonds[atom1][atom2], bonds[atom2][atom1]))
+
+        # Convert bonddict to use Atom[group] and Bond[group] objects
+        atomkeys = atomdict.keys()
+        atomkeys.sort()
+        for aid1 in atomkeys:
+            atomkeys2 = bonds[aid1].keys()
+            atomkeys2.sort()
+            for aid2 in atomkeys2:
+                if aid1 < aid2:
+                    atom1 = atomdict[aid1]
+                    atom2 = atomdict[aid2]
+                    order = bonds[aid1][aid2]
+                    if group:
+                        bond = GroupBond(atom1, atom2, order)
+                    elif len(order) == 1:
+                        bond = Bond(atom1, atom2, order[0])
+                    else:
+                        raise InvalidAdjacencyListError('Multiple bond orders specified for an atom in a Molecule.')
+                    atom1.edges[atom2] = bond
+                    atom2.edges[atom1] = bond
+        
+        if saturateH and not group:
+            # Add explicit hydrogen atoms to complete structure if desired
+            valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'Cl': 1, 'He': 0, 'Ne': 0, 'Ar': 0}
+            orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
+            newAtoms = []
+            for atom in atoms:
+                try:
+                    valence = valences[atom.symbol]
+                except KeyError:
+                    raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol))
+                radical = atom.radicalElectrons
+                order = 0
+                for atom2, bond in atom.bonds.items():
+                    order += orders[bond.order]
+                count = valence - radical - int(order)
+                for i in range(count):
+                    a = Atom(element='H', radicalElectrons=0, charge=0, label='')
+                    b = Bond(atom, a, 'S')
+                    newAtoms.append(a)
+                    atom.bonds[a] = b
+                    a.bonds[atom] = b
+            atoms.extend(newAtoms)
+        
+        # Calculate the number of lone pair electrons requiring molecule with all hydrogen atoms present
+        if not group and lonePairsOfElectrons == -1:
+            valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0, 'Cl': 1}
+            orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
+            for atom in atoms:
+                try:
+                    valence = valences[atom.symbol]
+                except KeyError:
+                    raise InvalidAdjacencyListError('Cannot calculate lone pairs: Unknown valence for atom "{0}".'.format(atom.symbol))
+                radical = atom.radicalElectrons
+                order = 0
+                for atom2, bond in atom.bonds.items():
+                    order += orders[bond.order]
+                lonePairs = (1 if atom.isHydrogen() else 4) - order - radical
+                charge = (2 if atom.isHydrogen() else 8) - valence - order - radical - 2 * lonePairs
+                atom.setLonePairs(lonePairs)
+                atom.updateCharge()
+                "It is pointless that we calculate this twice, but as we do, we may as well check they agree"
+                assert charge == atom.charge, "Charge calculations didn't agree!"
+
+        elif not group:
+            for atom in atoms:
+                atom.updateCharge()
+                    
+    except InvalidAdjacencyListError:
+        logging.error("Troublesome adjacency list:\n" + adjlist)
+        raise
+    
+    return atoms, multiplicity
+
+###############################
+
+re_IntermediateAdjList = re.compile('^\s*(\d)\s+' +  # atom number digit
+                          '(?P<label>\*\d*\s+)?' +  # optional label eg * or *2
+                          '(?P<atomtype>\{?[A-Z]\S*)\s+' +  # atomtype eg R!H or {Cb,Cd}
+                          '(?P<radicals>X|\d[STDQV]?|\{?\d[^}]*\})\s+' +  #radicals eg. X or 2T or {1,2,2T}
+                          '(?P<lonepairs>\d)' +  # lone pairs eg. 0
+                          '(?P<bonds>(\s+\{\d+\,(?:[SDTB]|\{.+?\})\},?)*)' +  # bonds, eg {2,S} {4,{S,D}}
+                          '\s*$')  # the end!
+
+re_OldAdjList = re.compile('^\s*(\d)\s+' +  # atom number digit
+                          '(?P<label>\*\d*\s+)?' +  # optional label eg * or *2
+                          '(?P<atomtype>\{?[A-Z]\S*)\s+' +  # atomtype eg R!H or {Cb,Cd}
+                          '(?P<radicals>X|\d[STDQV]?|\{?\d[^}]*\})' +  #radicals eg. X or 2T or {1,2,2T}
+                          '(?P<bonds>(\s+\{\d+\,(?:[SDTB]|\{.+?\})\},?)*)' +  # bonds, eg {2,S} {4,{S,D}}
+                          '\s*$')  # the end!
+
+def fromAdjacencyList(adjlist, group=False, saturateH=False):
+    """
+    Convert a string adjacency list `adjlist` into a set of :class:`Atom` and
+    :class:`Bond` objects.
+    """
+    atoms = []
+    atomdict = {}
+    bonds = {}
+    multiplicity = None
+    
+    try:
+        adjlist = adjlist.strip()
+        lines = adjlist.splitlines()
+        if adjlist == '' or len(lines) == 0:
+            raise InvalidAdjacencyListError('Empty adjacency list.')
+
+        # Detect old-style adjacency lists by looking at the last line's syntax
+        lastLine = lines[-1].strip()
+        while not lastLine:  # Remove any empty lines from the end
+            lines.pop()
+            lastLine = lines[-1].strip()
+        if re_IntermediateAdjList.match(lastLine):
+            logging.debug("Adjacency list line '{0}' looks like an intermediate style adjacency list".format(lastLine))
+            return fromOldAdjacencyList(adjlist, group=group, saturateH=saturateH)
+        if re_OldAdjList.match(lastLine):
+            logging.debug("Adjacency list line '{0}' looks like an old style adjacency list".format(lastLine))
+            if not group:
+                logging.debug("Will assume implicit H atoms")
+            return fromOldAdjacencyList(adjlist, group=group, saturateH=(not group))
+
+        # Interpret the first line if it contains a label
+        if len(lines[0].split()) == 1:
+            label = lines.pop(0)
+            if len(lines) == 0:
+                raise InvalidAdjacencyListError('No atoms specified in adjacency list.')
+            
+        # Interpret the second line if it contains a multiplicity
+        if lines[0].split()[0] == 'multiplicity':
+            line = lines.pop(0)
+            if group:
+                match = re.match('\s*multiplicity\s+\[\d(?:,\s*\d)*\]\s*$', line)
+                assert match, "Invalid multiplicity line '{0}'. Should be a list like 'multiplicity [1,2,3]'".format(line)
+                multiplicities = match.group(1).split(',')
+                multiplicity = [int(i) for i in multiplicities]
+            else:
+                match = re.match('\s*multiplicity\s+\d+\s*$', line)
+                assert match, "Invalid multiplicity line '{0}'. Should be an integer like 'multiplicity 2'".format(line)
+                multiplicity = int(line.split()[1])
+            if len(lines) == 0:
+                raise InvalidAdjacencyListError('No atoms specified in adjacency list.')
+        
+        mistake1 = re.compile('\{[^}]*\s+[^}]*\}')
+        # Iterate over the remaining lines, generating Atom or GroupAtom objects
+        for line in lines:
+
+            # Sometimes people put spaces after commas, which messes up the
+            # parse-by-whitespace. Examples include '[Cd, Ct]'.
+            if mistake1.search(line):
+                raise InvalidAdjacencyListError(
+                    "Shouldn't have spaces inside braces: {0}".format(mistake1.search(line).group())
+                    )
+
+            # Sometimes commas are used to delimit bonds in the bond list,
+            # so replace them just in case
+            line = line.replace('},{', '} {')
+            
+            data = line.split()
+
+            # Skip if blank line
+            if len(data) == 0: continue
+
+            # First item is index for atom
+            # Sometimes these have a trailing period (as if in a numbered list),
+            # so remove it just in case
+            aid = int(data[0].strip('.'))
+
+            # If second item starts with '*', then atom is labeled
+            label = ''; index = 1
+            if data[1][0] == '*':
+                label = data[1]
+                index += 1
+
+            # Next is the element or functional group element
+            # A list can be specified with the {,} syntax
+            atomType = data[index]
+            if atomType[0] == '[':
+                atomType = atomType[1:-1].split(',')
+            else:
+                atomType = [atomType]
+            index += 1
+            
+            # Next the number of unpaired electrons
+            unpairedElectrons = []
+            uState = data[index]
+            if uState[0] == 'u':
+                if uState[1] == '[':
+                    uState = uState[2:-1].split(',')
+                else:
+                    uState = [uState[1]]
+                for u in uState:
+                    if u == '0':
+                        unpairedElectrons.append(0)
+                    elif u == '1':
+                        unpairedElectrons.append(1)
+                    elif u == '2':
+                        unpairedElectrons.append(2)
+                    elif u == '3':
+                        unpairedElectrons.append(3)
+                    elif u == '4':
+                        unpairedElectrons.append(4)
+                    elif u == 'x':
+                        unpairedElectrons.append(0)
+                        unpairedElectrons.append(1)
+                        unpairedElectrons.append(2)
+                        unpairedElectrons.append(3)
+                        unpairedElectrons.append(4)
+                    else:
+                        raise InvalidAdjacencyListError('Number of unpaired electrons not recognized.')
+                index += 1
+            else:
+                raise InvalidAdjacencyListError('Number of unpaired electrons not defined.')
+            
+            # Next the number of lone electron pairs (if provided)
+            lonePairs = []
+            if len(data) > index:
+                lpState = data[index]
+                if lpState[0] == 'p':
+                    if lpState[1] == '[':
+                        lpState = lpState[2:-1].split(',')
+                    else:
+                        lpState = [lpState[1]]
+                    for l in lpState:
+                        if l == '0':
+                            lonePairs.append(0)
+                        elif l == '1':
+                            lonePairs.append(1)
+                        elif l == '2':
+                            lonePairs.append(2)
+                        elif l == '3':
+                            lonePairs.append(3)
+                        elif l == '4':
+                            lonePairs.append(4)
+                        elif l == 'x':
+                            lonePairs.append(0)
+                            lonePairs.append(1)
+                            lonePairs.append(2)
+                            lonePairs.append(3)
+                            lonePairs.append(4)
+                        else:
+                            raise InvalidAdjacencyListError('Number of lone electron pairs not recognized.')
+                    index += 1
+                else:
+                    if group:
+                        lonePairs.append(None)
+                    else:
+                        lonePairs.append(0)
+            else:
+                if group:
+                    lonePairs.append(None)
+                else:
+                    lonePairs.append(0)
+                
+            # Next the number of partial charges (if provided)
+            partialCharges = []
+            if len(data) > index:
+                eState = data[index]
+                if eState[0] == 'c':
+                    if eState[1] == '[':
+                        eState = eState[2:-1].split(',')
+                    else:
+                        eState = [eState[1:]]
+                    for e in eState:
+                        if e == '0':
+                            partialCharges.append(0)
+                        elif e == '+1':
+                            partialCharges.append(1)
+                        elif e == '+2':
+                            partialCharges.append(2)
+                        elif e == '+3':
+                            partialCharges.append(3)
+                        elif e == '+4':
+                            partialCharges.append(4)
+                        elif e == '-1':
+                            partialCharges.append(-1)
+                        elif e == '-2':
+                            partialCharges.append(-2)
+                        elif e == '-3':
+                            partialCharges.append(-3)
+                        elif e == '-4':
+                            partialCharges.append(-4)
+                        elif e == 'x':
+                            partialCharges.append(0)
+                            partialCharges.append(1)
+                            partialCharges.append(2)
+                            partialCharges.append(3)
+                            partialCharges.append(4)
+                            partialCharges.append(-1)
+                            partialCharges.append(-2)
+                            partialCharges.append(-3)
+                            partialCharges.append(-4)
+                        else:
+                            raise InvalidAdjacencyListError('Number of partial charges not recognized.')
+                    index += 1
+                else:
+                    partialCharges.append(0)
             else:
-                lonePairElectrons = -1
+                partialCharges.append(0)
             
             # Create a new atom based on the above information
             if group:
-                atom = GroupAtom(atomType, radicalElectrons, spinMultiplicity, [0 for e in radicalElectrons], label, [lonePairElectrons])
+                atom = GroupAtom(atomType, unpairedElectrons, partialCharges, label, lonePairs)
             else:
-                atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label, lonePairElectrons)
+                atom = Atom(atomType[0], unpairedElectrons[0], partialCharges[0], label, lonePairs[0])
 
             # Add the atom to the list
             atoms.append(atom)
@@ -193,7 +553,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
                 if aid == aid2:
                     raise InvalidAdjacencyListError('Attempted to create a bond between atom {0:d} and itself.'.format(aid))
                 
-                if order[0] == '{':
+                if order[0] == '[':
                     order = order[1:-1].split(',')
                 else:
                     order = [order]
@@ -233,21 +593,25 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
         if saturateH:
             # Add explicit hydrogen atoms to complete structure if desired
             if not group:
-                valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'Cl': 1, 'He': 0, 'Ne': 0, 'Ar': 0}
+                orbitals = {'H': 1, 'C': 4, 'O': 4, 'N': 4, 'S': 4, 'Si': 4, 'Cl': 1, 'He': 0, 'Ne': 0, 'Ar': 0}
                 orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
                 newAtoms = []
                 for atom in atoms:
                     try:
-                        valence = valences[atom.symbol]
+                        orbital = orbitals[atom.symbol]
                     except KeyError:
-                        raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol))
-                    radical = atom.radicalElectrons
+                        raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown orbital for atom "{0}".'.format(atom.symbol))
+                    
                     order = 0
                     for atom2, bond in atom.bonds.items():
                         order += orders[bond.order]
-                    count = valence - radical - int(order)
+                    count = orbital - atom.radicalElectrons - atom.lonePairs - int(order)
+                    
+                    if count < 0:
+                        raise InvalidAdjacencyListError('Incorrect electron configuration on atom.')
+                        
                     for i in range(count):
-                        a = Atom('H', 0, 1, 0, '')
+                        a = Atom(element='H', radicalElectrons=0, charge=0, label='', lonePairs=0)
                         b = Bond(atom, a, 'S')
                         newAtoms.append(a)
                         atom.bonds[a] = b
@@ -255,7 +619,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
                 atoms.extend(newAtoms)
         
         # Calculate the number of lone pair electrons requiring molecule with all hydrogen atoms present
-        if not group and lonePairElectrons == -1:
+        if not group and lonePairs is not None:
             valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0, 'Cl': 1}
             orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
             for atom in atoms:
@@ -286,44 +650,32 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False):
             for atom in atoms:
                 atom.updateCharge()
                     
-    except InvalidAdjacencyListError:
-        print adjlist
-        raise
+    except InvalidAdjacencyListError, e:
+        logging.debug("Assuming this adjacency list is old-style because {1}:\n{0}".format(adjlist, e.message))
+        logging.debug("Though it should have been detected by the regex test")
+        return fromOldAdjacencyList(adjlist, group=group, saturateH=saturateH)
     
-    return atoms
-
-################################################################################
-
-def getElectronState(radicalElectrons, spinMultiplicity):
-    """
-    Return the electron state corresponding to the given number of radical
-    electrons `radicalElectrons` and spin multiplicity `spinMultiplicity`.
-    Raise a :class:`ValueError` if the electron state cannot be determined.
-    """
-    electronState = ''
-    if radicalElectrons == 0: 
-        electronState = '0'
-    elif radicalElectrons == 1:
-        electronState = '1'
-    elif radicalElectrons == 2 and spinMultiplicity == 1:
-        electronState = '2S'
-    elif radicalElectrons == 2 and spinMultiplicity == 3: 
-        electronState = '2T'
-    elif radicalElectrons == 3 and spinMultiplicity == 2: 
-        electronState = '3D'
-    elif radicalElectrons == 3 and spinMultiplicity == 4: 
-        electronState = '3Q'
-    elif radicalElectrons == 4 and spinMultiplicity == 1: 
-        electronState = '4S'
-    elif radicalElectrons == 4 and spinMultiplicity == 3: 
-        electronState = '4T'
-    elif radicalElectrons == 4 and spinMultiplicity == 5: 
-        electronState = '4V'
+    if not group:
+        nRad = 0
+        for atom in atoms:
+            nRad += atom.radicalElectrons
+        if multiplicity == None:
+            multiplicity = nRad + 1
+        n = 0
+        while (nRad + 1 - n*2) > 0:
+            if (nRad + 1 - n*2) == multiplicity:
+                break
+            n=n+1
+        else:
+            print adjlist
+            raise InvalidAdjacencyListError('Multiplicity not in agreement with total number of radicals.')
+            
+        return atoms, multiplicity
     else:
-        raise ValueError('Unable to determine electron state for {0:d} radical electrons with spin multiplicity of {1:d}.'.format(radicalElectrons, spinMultiplicity))
-    return electronState
+        return atoms, multiplicity
+
 
-def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePairs=False):
+def toAdjacencyList(atoms, multiplicity, label=None, group=False, removeH=False, removeLonePairs=False):
     """
     Convert a chemical graph defined by a list of `atoms` into a string
     adjacency list.
@@ -337,6 +689,16 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai
         pass
 
     if label: adjlist += label + '\n'
+    
+    if group:
+        if multiplicity is not None:
+            assert isinstance(multiplicity, list), "Functional group should have a list of possible multiplicities"
+            if multiplicity != [1,2,3,4,5]:
+                adjlist += 'multiplicity [{0!s}]\n'.format(','.join(str(i) for i in multiplicity))
+    else:
+        assert isinstance(multiplicity, int), "Molecule should have an integer multiplicity"
+        if multiplicity != 1 or any( atom.radicalElectrons for atom in atoms ):
+            adjlist += 'multiplicity {0!r}\n'.format(multiplicity)
 
     # Determine the numbers to use for each atom
     atomNumbers = {}
@@ -349,37 +711,52 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai
     atomLabels = dict([(atom, '{0}'.format(atom.label)) for atom in atomNumbers])
     
     atomTypes = {}
-    atomElectronStates = {}
+    atomUnpairedElectrons = {}
     atomLonePairs = {}
+    atomCharge = {}
     if group:
         for atom in atomNumbers:
             # Atom type(s)
-            if len(atom.atomType) == 1:
+            if len(atom.atomType) == 1: 
                 atomTypes[atom] = atom.atomType[0].label
             else:
-                atomTypes[atom] = '{{{0}}}'.format(','.join([a.label for a in atom.atomType]))
-            # Electron state(s)
+                atomTypes[atom] = '[{0}]'.format(','.join([a.label for a in atom.atomType]))
+            # Unpaired Electron(s)
             if len(atom.radicalElectrons) == 1: 
-                atomElectronStates[atom] = getElectronState(atom.radicalElectrons[0], atom.spinMultiplicity[0])
+                atomUnpairedElectrons[atom] = str(atom.radicalElectrons[0])
+            elif set(atom.radicalElectrons) == set(range(5)):
+                atomUnpairedElectrons[atom] = 'x'  # the wildcard represents [0,1,2,3,4]
             else:
-                atomElectronStates[atom] = '{{{0}}}'.format(','.join([getElectronState(radical, spin) for radical, spin in zip(atom.radicalElectrons, atom.spinMultiplicity)]))  
+                atomUnpairedElectrons[atom] = '[{0}]'.format(','.join([str(radical) for radical in atom.radicalElectrons]))
+
+            # Lone Electron Pair(s)
+            if atom.lonePairs is None or not atom.lonePairs: # if None or []
+                atomLonePairs[atom] = None
+            elif len(atom.lonePairs) == 1: 
+                atomLonePairs[atom] = str(atom.lonePairs[0])
+            elif set(atom.lonePairs) == set(range(5)):
+                atomLonePairs[atom] = 'x'  # the wildcard represents [0,1,2,3,4]
+            else:
+                atomLonePairs[atom] = '[{0}]'.format(','.join([str(pair) for pair in atom.lonePairs]))
     else:
         for atom in atomNumbers:
             # Atom type
             atomTypes[atom] = '{0}'.format(atom.element.symbol)
-            # Electron state(s)
-            atomElectronStates[atom] = '{0}'.format(getElectronState(atom.radicalElectrons, atom.spinMultiplicity))    
-            if not removeLonePairs:
-                # Lone Pair(s)
-                atomLonePairs[atom] = atom.lonePairs
+            # Unpaired Electron(s)
+            atomUnpairedElectrons[atom] = '{0}'.format(atom.radicalElectrons)
+            # Lone Electron Pair(s)
+            atomLonePairs[atom] = atom.lonePairs
+            # Partial Charge(s)
+            atomCharge[atom] = atom.charge
     
     # Determine field widths
     atomNumberWidth = max([len(s) for s in atomNumbers.values()]) + 1
     atomLabelWidth = max([len(s) for s in atomLabels.values()])
     if atomLabelWidth > 0: atomLabelWidth += 1
     atomTypeWidth = max([len(s) for s in atomTypes.values()]) + 1
-    atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()])
-    atomLonePairWidth = 2
+    atomUnpairedElectronsWidth = max([len(s) for s in atomUnpairedElectrons.values()])
+    atomLonePairWidth = 1
+    atomChargeWidth = 1
     
     # Assemble the adjacency list
     for atom in atoms:
@@ -391,11 +768,19 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai
         adjlist += '{0:<{1:d}}'.format(atomLabels[atom], atomLabelWidth)
         # Atom type(s)
         adjlist += '{0:<{1:d}}'.format(atomTypes[atom], atomTypeWidth)
-        # Electron state(s)
-        adjlist += '{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth)
-        if group == False and not removeLonePairs:
-            # Lone Pair(s)
-            adjlist += '{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth)
+        # Unpaired Electron(s)
+        adjlist += 'u{0:<{1:d}}'.format(atomUnpairedElectrons[atom], atomUnpairedElectronsWidth)
+        # Lone Electron Pair(s)
+        if str(atomLonePairs[atom]) != 'None':
+            adjlist += ' p{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth)
+        if not group:
+            # Partial Charge(s)
+            if atomCharge[atom] > 0:
+                adjlist += ' c+{0:>{1:d}}'.format(atomCharge[atom], atomChargeWidth)
+            elif atomCharge[atom] < 0:
+                adjlist += ' c{0:>{1:d}}'.format(atomCharge[atom], atomChargeWidth)
+            else:
+                adjlist += ' c{0:>{1:d}} '.format(atomCharge[atom], atomChargeWidth)
         
         # Bonds list
         atoms2 = atom.bonds.keys()
@@ -413,7 +798,7 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai
                 if len(bond.order) == 1:
                     adjlist += bond.order[0]
                 else:
-                    adjlist += '{{{0}}}'.format(','.join(bond.order))
+                    adjlist += '[{0}]'.format(','.join(bond.order))
             else:
                 adjlist += bond.order
             adjlist += '}'
diff --git a/rmgpy/molecule/adjlistTest.py b/rmgpy/molecule/adjlistTest.py
new file mode 100644
index 0000000000..2265ccc635
--- /dev/null
+++ b/rmgpy/molecule/adjlistTest.py
@@ -0,0 +1,402 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+import unittest
+from external.wip import work_in_progress
+from rmgpy.molecule.molecule import *
+from rmgpy.molecule.group import Group
+from rmgpy.molecule.element import getElement, elementList
+
+import logging
+logging.basicConfig(level=logging.DEBUG)
+
+################################################################################
+
+class TestGroupAdjLists(unittest.TestCase):
+    """
+    Contains adjacency list unit tests of the Graph class.
+    """
+    def setUp(self):
+        pass
+
+    def testFromOldAdjacencyList1(self):
+        """
+        adjlist: Test the Group.fromAdjacencyList() method on an old style adjacency list.
+        """
+        adjlist = """
+1 *2 {Cs,Cd} 0 {2,{S,D}} {3,S}
+2 *1 {Os,Od}  0   {1,{S,D}}
+3    R!H     {0,1} {1,S}
+            """
+        group = Group().fromAdjacencyList(adjlist)
+
+        atom1, atom2, atom3 = group.atoms
+        self.assertTrue(group.hasBond(atom1, atom2))
+        self.assertTrue(group.hasBond(atom1, atom3))
+        self.assertFalse(group.hasBond(atom2, atom3))
+        bond12 = atom1.bonds[atom2]
+        bond13 = atom1.bonds[atom3]
+
+        self.assertTrue(atom1.label == '*2')
+        self.assertTrue(atom1.atomType[0].label in ['Cs', 'Cd'])
+        self.assertTrue(atom1.atomType[1].label in ['Cs', 'Cd'])
+        self.assertTrue(atom1.radicalElectrons == [0])
+
+        self.assertTrue(atom2.label == '*1')
+        self.assertTrue(atom2.atomType[0].label in ['Os', 'Od'])
+        self.assertTrue(atom2.atomType[1].label in ['Os', 'Od'])
+        self.assertTrue(atom2.radicalElectrons == [0])
+
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.atomType[0].label == 'R!H')
+        self.assertTrue(atom3.radicalElectrons == [0, 1])
+
+        self.assertTrue(bond12.order == ['S', 'D'])
+        self.assertTrue(bond13.order == ['S'])
+
+
+    def testFromAdjacencyList(self):
+        """
+        adjlist: Test the Group.fromAdjacencyList() method.
+        """
+        adjlist = """
+1 *2 [Cs,Cd] u0 {2,[S,D]} {3,S}
+2 *1 [Os,Od] u0 {1,[S,D]}
+3    R!H     u0 {1,S}
+            """
+        group = Group().fromAdjacencyList(adjlist)
+        
+        atom1, atom2, atom3 = group.atoms
+        self.assertTrue(group.hasBond(atom1, atom2))
+        self.assertTrue(group.hasBond(atom1, atom3))
+        self.assertFalse(group.hasBond(atom2, atom3))
+        bond12 = atom1.bonds[atom2]
+        bond13 = atom1.bonds[atom3]
+           
+        self.assertTrue(atom1.label == '*2')
+        self.assertTrue(atom1.atomType[0].label in ['Cs', 'Cd'])
+        self.assertTrue(atom1.atomType[1].label in ['Cs', 'Cd'])
+        self.assertTrue(atom1.radicalElectrons == [0])
+
+        self.assertTrue(atom2.label == '*1')
+        self.assertTrue(atom2.atomType[0].label in ['Os', 'Od'])
+        self.assertTrue(atom2.atomType[1].label in ['Os', 'Od'])
+        self.assertTrue(atom2.radicalElectrons == [0])
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.atomType[0].label == 'R!H')
+        self.assertTrue(atom3.radicalElectrons == [0])
+
+        self.assertTrue(bond12.order == ['S', 'D'])
+        self.assertTrue(bond13.order == ['S'])
+
+    def testToAdjacencyList(self):
+        """
+        adjlist: Test the Group.toAdjacencyList() method.
+        """
+        adjlist = """
+1 *2 [Cs,Cd] u0 {2,[S,D]} {3,S}
+2 *1 [Os,Od] u0 {1,[S,D]}
+3    R!H     u0 {1,S}
+            """
+        group = Group().fromAdjacencyList(adjlist)
+        adjlist2 = group.toAdjacencyList()
+
+        self.assertEqual(adjlist.strip(), adjlist2.strip())
+
+
+class TestMoleculeAdjLists(unittest.TestCase):
+    """
+    adjlist: Contains adjacency list unit tests of the Molecule class.
+    """
+    
+    def setUp(self):
+        pass
+
+    def testFromAdjacencyList1(self):
+        """
+        adjlist: Test the Molecule.fromAdjacencyList() method 1.
+        """
+        # molecule 1
+        adjlist = """
+1 *1 C u1 p0 c0  {2,S} {3,S} {4,S}
+2    H u0 p0 c0  {1,S}
+3    H u0 p0 c0  {1,S}
+4 *2 N u0 p0 c+1 {1,S} {5,S} {6,D}
+5    O u0 p3 c-1 {4,S}
+6    O u0 p2 c0  {4,D}
+            """
+        molecule = Molecule().fromAdjacencyList(adjlist)
+        
+        self.assertTrue(molecule.multiplicity == 2)
+        
+        atom1 = molecule.atoms[0]
+        atom2 = molecule.atoms[3]
+        atom3 = molecule.atoms[4]
+        atom4 = molecule.atoms[5]
+        self.assertTrue(molecule.hasBond(atom2, atom1))
+        self.assertTrue(molecule.hasBond(atom2, atom3))
+        self.assertTrue(molecule.hasBond(atom2, atom4))
+        self.assertFalse(molecule.hasBond(atom1, atom3))
+        self.assertFalse(molecule.hasBond(atom1, atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+        
+    def testFromAdjacencyList2(self):
+        """
+        adjlist: Test the Molecule.fromAdjacencyList() method 2.
+        """
+        # molecule 2
+        adjlist = """
+1 *1 C u1 p0 {2,S} {3,S} {4,S}
+2    H u0 p0 {1,S}
+3    H u0 p0 {1,S}
+4 *2 N u0 p0 {1,S} {5,S} {6,D}
+5    O u0 p3 {4,S}
+6    O u0 p2 {4,D}
+            """
+        molecule = Molecule().fromAdjacencyList(adjlist)
+        
+        self.assertTrue(molecule.multiplicity == 2)
+        
+        atom1 = molecule.atoms[0]
+        atom2 = molecule.atoms[3]
+        atom3 = molecule.atoms[4]
+        atom4 = molecule.atoms[5]
+        self.assertTrue(molecule.hasBond(atom2, atom1))
+        self.assertTrue(molecule.hasBond(atom2, atom3))
+        self.assertTrue(molecule.hasBond(atom2, atom4))
+        self.assertFalse(molecule.hasBond(atom1, atom3))
+        self.assertFalse(molecule.hasBond(atom1, atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+    def testFromAdjacencyList3(self):
+        """
+        adjlist: Test the Molecule.fromAdjacencyList() method 3.
+        """
+        # molecule 3
+        adjlist = """
+1 *1 C u1 {2,S} {3,S} {4,S}
+2    H u0 {1,S}
+3    H u0 {1,S}
+4 *2 N u0 {1,S} {5,S} {6,D}
+5    O u0 {4,S}
+6    O u0 {4,D}
+            """
+        molecule = Molecule().fromAdjacencyList(adjlist)
+        
+        self.assertTrue(molecule.multiplicity == 2)
+        
+        atom1 = molecule.atoms[0]
+        atom2 = molecule.atoms[3]
+        atom3 = molecule.atoms[4]
+        atom4 = molecule.atoms[5]
+        self.assertTrue(molecule.hasBond(atom2, atom1))
+        self.assertTrue(molecule.hasBond(atom2, atom3))
+        self.assertTrue(molecule.hasBond(atom2, atom4))
+        self.assertFalse(molecule.hasBond(atom1, atom3))
+        self.assertFalse(molecule.hasBond(atom1, atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+        
+    def testFromAdjacencyList4(self):
+        """
+        adjlist: Test the Molecule.fromAdjacencyList() method 4.
+        """
+        # molecule 4
+        adjlist = """
+1 *1 C u1 p0 {2,S}
+2 *2 N u0 p0 {1,S} {3,S} {4,D}
+3    O u0 p3 {2,S}
+4    O u0 p2 {2,D}
+            """
+        molecule = Molecule().fromAdjacencyList(adjlist, saturateH=True)
+        
+        self.assertTrue(molecule.multiplicity == 2)
+        
+        atom1 = molecule.atoms[0]
+        atom2 = molecule.atoms[1]
+        atom3 = molecule.atoms[2]
+        atom4 = molecule.atoms[3]
+        self.assertTrue(molecule.hasBond(atom2, atom1))
+        self.assertTrue(molecule.hasBond(atom2, atom3))
+        self.assertTrue(molecule.hasBond(atom2, atom4))
+        self.assertFalse(molecule.hasBond(atom1, atom3))
+        self.assertFalse(molecule.hasBond(atom1, atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+
+
+    def testToAdjacencyList(self):
+        """
+        adjlist: Test the Molecule.toAdjacencyList() method.
+        """
+        adjlist = """
+1 *1 C u1 p0 c0  {2,S} {3,S} {4,S}
+2    H u0 p0 c0  {1,S}
+3    H u0 p0 c0  {1,S}
+4 *2 N u0 p0 c+1 {1,S} {5,S} {6,D}
+5    O u0 p3 c-1 {4,S}
+6    O u0 p2 c0  {4,D}
+            """
+        molecule = Molecule().fromAdjacencyList(adjlist)
+        adjlist_1 = molecule.toAdjacencyList(removeH=False)
+        newMolecule = Molecule().fromAdjacencyList(adjlist_1)
+        self.assertTrue(molecule.isIsomorphic(newMolecule))
+        
+        #self.assertEqual(adjlist_1.strip(), adjlist.strip())
+        
+    def testFromIntermediateAdjacencyList1(self):
+        """
+        Test we can read an intermediate style adjacency list with implicit hydrogens 1
+        """
+        adjList = """
+        1 O 0 0
+        """  # should be Water
+        molecule = Molecule().fromAdjacencyList(adjList, saturateH=True)  # only works with saturateH=True
+        self.assertEqual(molecule.getFormula(), 'H2O')
+
+    def testFromOldAdjacencyList1(self):
+        """
+        Test we can read an old style adjacency list with implicit hydrogens 1
+        """
+        adjList = """
+        1 O 0 
+        """  # should be Water
+        molecule = Molecule().fromAdjacencyList(adjList)  # only works with saturateH=True
+        self.assertEqual(molecule.getFormula(), 'H2O')
+
+
+    def testAdjacencyList(self):
+        """
+        adjlist: Check the adjacency list read/write functions for a full molecule.
+        """
+        molecule1 = Molecule().fromAdjacencyList("""
+        1  C u0 {2,D} {7,S} {8,S}
+        2  C u0 {1,D} {3,S} {9,S}
+        3  C u0 {2,S} {4,D} {10,S}
+        4  C u0 {3,D} {5,S} {11,S}
+        5  C u1 {4,S} {6,S} {12,S}
+        6  C u0 {5,S} {13,S} {14,S} {15,S}
+        7  H u0 {1,S}
+        8  H u0 {1,S}
+        9  H u0 {2,S}
+        10 H u0 {3,S}
+        11 H u0 {4,S}
+        12 H u0 {5,S}
+        13 H u0 {6,S}
+        14 H u0 {6,S}
+        15 H u0 {6,S}
+        """)
+        molecule2 = Molecule().fromSMILES('C=CC=C[CH]C')
+        self.assertTrue(molecule1.isIsomorphic(molecule2))
+        self.assertTrue(molecule2.isIsomorphic(molecule1))
+
+################################################################################
+
+if __name__ == '__main__':
+
+    unittest.main(testRunner=unittest.TextTestRunner(verbosity=3))
diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py
index 899fd51ae4..ddfc141abe 100644
--- a/rmgpy/molecule/atomtype.py
+++ b/rmgpy/molecule/atomtype.py
@@ -157,6 +157,12 @@ def isSpecificCaseOf(self, other):
 
 ################################################################################
 
+
+"""
+Note: function to read adjacency lists assumes that all atom types begin
+with a capital letter [A-Z]
+"""
+
 atomTypes = {}
 atomTypes['R']    = AtomType(label='R', generic=[], specific=[
     'R!H',
diff --git a/rmgpy/molecule/atomtypeTest.py b/rmgpy/molecule/atomtypeTest.py
index 97de266933..3750ed469f 100644
--- a/rmgpy/molecule/atomtypeTest.py
+++ b/rmgpy/molecule/atomtypeTest.py
@@ -112,75 +112,75 @@ def setUp(self):
         self.mol1 = Molecule().fromSMILES('COC(=O)CC=C=CC#C')
         # self.mol2 = Molecule().fromSMILES('c1ccccc1')
         ## the fromSMILES method currently Kekulizes, so to test Benzene we use fromAdjacencyList
-        self.mol2 = Molecule().fromAdjacencyList('''1  C 0 0 {2,B} {6,B} {7,S}
-                                                    2  C 0 0 {1,B} {3,B} {8,S}
-                                                    3  C 0 0 {2,B} {4,B} {9,S}
-                                                    4  C 0 0 {3,B} {5,B} {10,S}
-                                                    5  C 0 0 {4,B} {6,B} {11,S}
-                                                    6  C 0 0 {1,B} {5,B} {12,S}
-                                                    7  H 0 0 {1,S}
-                                                    8  H 0 0 {2,S}
-                                                    9  H 0 0 {3,S}
-                                                    10 H 0 0 {4,S}
-                                                    11 H 0 0 {5,S}
-                                                    12 H 0 0 {6,S}''')
+        self.mol2 = Molecule().fromAdjacencyList('''1  C u0 p0 {2,B} {6,B} {7,S}
+                                                    2  C u0 p0 {1,B} {3,B} {8,S}
+                                                    3  C u0 p0 {2,B} {4,B} {9,S}
+                                                    4  C u0 p0 {3,B} {5,B} {10,S}
+                                                    5  C u0 p0 {4,B} {6,B} {11,S}
+                                                    6  C u0 p0 {1,B} {5,B} {12,S}
+                                                    7  H u0 p0 {1,S}
+                                                    8  H u0 p0 {2,S}
+                                                    9  H u0 p0 {3,S}
+                                                    10 H u0 p0 {4,S}
+                                                    11 H u0 p0 {5,S}
+                                                    12 H u0 p0 {6,S}''')
         self.mol3 = Molecule().fromSMILES('[H]')
         self.mol4 = Molecule().fromSMILES(
                                 'O=[Si][Si][Si]=[Si]=[Si][Si]#[Si]SS=S')
-        self.mol5 = Molecule().fromAdjacencyList('''1 H 0 0 {3,S}
-                                                    2 H 0 0 {3,S}
-                                                    3 N 0 0 {1,S} {2,S} {4,D}
-                                                    4 N 0 2 {3,D}''')
+        self.mol5 = Molecule().fromAdjacencyList('''1 H u0 p0 {3,S}
+                                                    2 H u0 p0 {3,S}
+                                                    3 N u0 p0 {1,S} {2,S} {4,D}
+                                                    4 N u0 p2 {3,D}''')
         self.mol6 = Molecule().fromSMILES('[Ar]')
         self.mol7 = Molecule().fromSMILES('[He]')
         self.mol8 = Molecule().fromSMILES('[Ne]')
-        self.mol9 = Molecule().fromAdjacencyList('''1 N 0 1 {2,S} {3,S} {4,S}
-                                                    2 H 0 0 {1,S}
-                                                    3 H 0 0 {1,S}
-                                                    4 H 0 0 {1,S}''')
+        self.mol9 = Molecule().fromAdjacencyList('''1 N u0 p1 {2,S} {3,S} {4,S}
+                                                    2 H u0 p0 {1,S}
+                                                    3 H u0 p0 {1,S}
+                                                    4 H u0 p0 {1,S}''')
         
-        self.mol10 = Molecule().fromAdjacencyList('''1 N 1 1 {2,S} {3,S}
-                                                     2 H 0 0 {1,S}
-                                                     3 H 0 0 {1,S}''')
+        self.mol10 = Molecule().fromAdjacencyList('''1 N u1 p1 {2,S} {3,S}
+                                                     2 H u0 p0 {1,S}
+                                                     3 H u0 p0 {1,S}''')
         
-        self.mol11 = Molecule().fromAdjacencyList('''1 N 2 1 {2,S}
-                                                     2 H 0 0 {1,S}''')
+        self.mol11 = Molecule().fromAdjacencyList('''1 N u2 p1 {2,S}
+                                                     2 H u0 p0 {1,S}''')
         
-        self.mol12 = Molecule().fromAdjacencyList('''1 N 0 1 {2,T}
-                                                     2 C 1 0 {1,T}''')
+        self.mol12 = Molecule().fromAdjacencyList('''1 N u0 p1 {2,T}
+                                                     2 C u1 p0 {1,T}''')
         
-        self.mol13 = Molecule().fromAdjacencyList('''1 N 0 0 {2,S} {3,S} {4,S} {5,S}
-                                                     2 H 0 0 {1,S}
-                                                     3 H 0 0 {1,S}
-                                                     4 H 0 0 {1,S}
-                                                     5 O 0 3 {1,S}''')
+        self.mol13 = Molecule().fromAdjacencyList('''1 N u0 p0 {2,S} {3,S} {4,S} {5,S}
+                                                     2 H u0 p0 {1,S}
+                                                     3 H u0 p0 {1,S}
+                                                     4 H u0 p0 {1,S}
+                                                     5 O u0 p3 {1,S}''')
         
-        self.mol14 = Molecule().fromAdjacencyList('''1 N 0 2 {2,D}
-                                                     2 N 0 0 {1,D} {3,D}
-                                                     3 O 0 2 {2,D}''')
+        self.mol14 = Molecule().fromAdjacencyList('''1 N u0 p2 {2,D}
+                                                     2 N u0 p0 {1,D} {3,D}
+                                                     3 O u0 p2 {2,D}''')
         
-        self.mol15 = Molecule().fromAdjacencyList('''1 N 0 1 {2,T}
-                                                     2 N 0 0 {1,T} {3,S}
-                                                     3 O 0 3 {2,S}''')
+        self.mol15 = Molecule().fromAdjacencyList('''1 N u0 p1 {2,T}
+                                                     2 N u0 p0 {1,T} {3,S}
+                                                     3 O u0 p3 {2,S}''')
         
-        self.mol16 = Molecule().fromAdjacencyList('''1 N 0 1 {2,D} {3,S}
-                                                     2 O 0 2 {1,D}
-                                                     3 O 1 2 {1,S}''')
+        self.mol16 = Molecule().fromAdjacencyList('''1 N u0 p1 {2,D} {3,S}
+                                                     2 O u0 p2 {1,D}
+                                                     3 O u1 p2 {1,S}''')
         
-        self.mol17 = Molecule().fromAdjacencyList('''1 N 1 1 {2,D}
-                                                     2 O 0 2 {1,D}''')
+        self.mol17 = Molecule().fromAdjacencyList('''1 N u1 p1 {2,D}
+                                                     2 O u0 p2 {1,D}''')
         
-        self.mol18 = Molecule().fromAdjacencyList('''1  N 0 0 {2,B} {6,B} {7,S}
-                                                     2  C 0 0 {1,B} {3,B} {8,S}
-                                                     3  C 0 0 {2,B} {4,B} {9,S}
-                                                     4  C 0 0 {3,B} {5,B} {10,S}
-                                                     5  C 0 0 {4,B} {6,B} {11,S}
-                                                     6  N 0 1 {1,B} {5,B}
-                                                     7  O 0 3 {1,S}
-                                                     8  H 0 0 {2,S}
-                                                     9  H 0 0 {3,S}
-                                                     10 H 0 0 {4,S}
-                                                     11 H 0 0 {5,S}''')
+        self.mol18 = Molecule().fromAdjacencyList('''1  N u0 p0 {2,B} {6,B} {7,S}
+                                                     2  C u0 p0 {1,B} {3,B} {8,S}
+                                                     3  C u0 p0 {2,B} {4,B} {9,S}
+                                                     4  C u0 p0 {3,B} {5,B} {10,S}
+                                                     5  C u0 p0 {4,B} {6,B} {11,S}
+                                                     6  N u0 p1 {1,B} {5,B}
+                                                     7  O u0 p3 {1,S}
+                                                     8  H u0 p0 {2,S}
+                                                     9  H u0 p0 {3,S}
+                                                     10 H u0 p0 {4,S}
+                                                     11 H u0 p0 {5,S}''')
         
     
     def atomType(self, mol, atomID):
diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py
index 9b2bd57c70..354fb46e54 100644
--- a/rmgpy/molecule/draw.py
+++ b/rmgpy/molecule/draw.py
@@ -166,7 +166,7 @@ def draw(self, molecule, format, path=None):
         self.__findRingGroups()
         
         # Handle carbon monoxide special case
-        if self.molecule.toSMILES() == 'C#O':
+        if self.molecule.getFormula == 'CO':
             # RDKit does not accept atom type Ot
             self.molecule.removeAtom(self.molecule.atoms[-1])
             self.symbols = ['CO']
diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd
index efdc8bdde5..87378d5eb3 100644
--- a/rmgpy/molecule/group.pxd
+++ b/rmgpy/molecule/group.pxd
@@ -33,7 +33,6 @@ cdef class GroupAtom(Vertex):
 
     cdef public list atomType
     cdef public list radicalElectrons
-    cdef public list spinMultiplicity
     cdef public list charge
     cdef public str label
     cdef public list lonePairs
@@ -80,6 +79,8 @@ cdef class GroupBond(Edge):
 
 cdef class Group(Graph):
 
+    cdef public list multiplicity
+
     # These read-only attribues act as a "fingerprint" for accelerating
     # subgraph isomorphism checks
     cdef public short carbonCount
diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py
index e9eea060ae..03b5df26d7 100644
--- a/rmgpy/molecule/group.py
+++ b/rmgpy/molecule/group.py
@@ -61,7 +61,6 @@ class GroupAtom(Vertex):
     =================== =================== ====================================
     `atomType`          ``list``            The allowed atom types (as :class:`AtomType` objects)
     `radicalElectrons`  ``list``            The allowed numbers of radical electrons (as short integers)
-    `spinMultiplicity`  ``list``            The allowed spin multiplicities (as short integers)
     `charge`            ``list``            The allowed formal charges (as short integers)
     `label`             ``str``             A string label that can be used to tag individual atoms
     `lonePairs`         ``list``            The number of lone electron pairs
@@ -69,19 +68,18 @@ class GroupAtom(Vertex):
 
     Each list represents a logical OR construct, i.e. an atom will match the
     group if it matches *any* item in the list. However, the
-    `radicalElectrons`, `spinMultiplicity`, and `charge` attributes are linked
+    `radicalElectrons`, and `charge` attributes are linked
     such that an atom must match values from the same index in each of these in
     order to match.
     """
 
-    def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label='', lonePairs=None):
+    def __init__(self, atomType=None, radicalElectrons=None, charge=None, label='', lonePairs=None):
         Vertex.__init__(self)
         self.atomType = atomType or []
         for index in range(len(self.atomType)):
             if isinstance(self.atomType[index], str):
                 self.atomType[index] = atomTypes[self.atomType[index]]
         self.radicalElectrons = radicalElectrons or []
-        self.spinMultiplicity = spinMultiplicity or []
         self.charge = charge or []
         self.label = label
         self.lonePairs = lonePairs or []
@@ -100,7 +98,7 @@ def __reduce__(self):
         atomType = self.atomType
         if atomType is not None:
             atomType = [a.label for a in atomType]
-        return (GroupAtom, (atomType, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label, self.lonePairs), d)
+        return (GroupAtom, (atomType, self.radicalElectrons, self.charge, self.label, self.lonePairs), d)
 
     def __setstate__(self, d):
         """
@@ -132,7 +130,7 @@ def copy(self):
         Return a deep copy of the :class:`GroupAtom` object. Modifying the
         attributes of the copy will not affect the original.
         """
-        return GroupAtom(self.atomType[:], self.radicalElectrons[:], self.spinMultiplicity[:], self.charge[:], self.label)
+        return GroupAtom(self.atomType[:], self.radicalElectrons[:], self.charge[:], self.label)
 
     def __changeBond(self, order):
         """
@@ -191,15 +189,12 @@ def __gainRadical(self, radical):
         where `radical` specifies the number of radical electrons to add.
         """
         radicalElectrons = []
-        spinMultiplicity = []
         if any([len(atomType.incrementRadical) == 0 for atomType in self.atomType]):
             raise ActionError('Unable to update GroupAtom due to GAIN_RADICAL action: Unknown atom type produced from set "{0}".'.format(self.atomType))
-        for electron, spin in zip(self.radicalElectrons, self.spinMultiplicity):
+        for electron in self.radicalElectrons:
             radicalElectrons.append(electron + radical)
-            spinMultiplicity.append(spin + radical)
-        # Set the new radical electron counts and spin multiplicities
+        # Set the new radical electron counts
         self.radicalElectrons = radicalElectrons
-        self.spinMultiplicity = spinMultiplicity
 
     def __loseRadical(self, radical):
         """
@@ -207,28 +202,17 @@ def __loseRadical(self, radical):
         where `radical` specifies the number of radical electrons to remove.
         """
         radicalElectrons = []
-        spinMultiplicity = []
         pairs = set()
         if any([len(atomType.decrementRadical) == 0 for atomType in self.atomType]):
             raise ActionError('Unable to update GroupAtom due to LOSE_RADICAL action: Unknown atom type produced from set "{0}".'.format(self.atomType))
-        for electron, spin in zip(self.radicalElectrons, self.spinMultiplicity):
+        for electron in self.radicalElectrons:
             electron = electron - radical
             if electron < 0:
-                raise ActionError('Unable to update GroupAtom due to LOSE_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons))
-            spin = spin - radical
-            if spin <= 0:
-                spin += 2
-            
-            pair = (electron,spin)
-            if pair in pairs:
-                continue # with next electron,spin pair, so we don't get redundant answers. Otherwise....
-            pairs.add(pair)
+                raise ActionError('Unable to update GroupAtom due to LOSE_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons))    
             radicalElectrons.append(electron)
-            spinMultiplicity.append(spin)
             
-        # Set the new radical electron counts and spin multiplicities
+        # Set the new radical electron counts
         self.radicalElectrons = radicalElectrons
-        self.spinMultiplicity = spinMultiplicity
         
     def __gainPair(self, pair):
         """
@@ -310,16 +294,17 @@ def equivalent(self, other):
             else:
                 return False
         # Each free radical electron state in self must have an equivalent in other (and vice versa)
-        for radical1, spin1 in zip(self.radicalElectrons, self.spinMultiplicity):
-            for radical2, spin2 in zip(other.radicalElectrons, other.spinMultiplicity):
-                if radical1 == radical2 and spin1 == spin2: break
+        for radical1 in self.radicalElectrons:
+            for radical2  in other.radicalElectrons:
+                if radical1 == radical2: break
             else:
                 return False
-        for radical1, spin1 in zip(other.radicalElectrons, other.spinMultiplicity):
-            for radical2, spin2 in zip(self.radicalElectrons, self.spinMultiplicity):
-                if radical1 == radical2 and spin1 == spin2: break
+        for radical1 in other.radicalElectrons:
+            for radical2 in self.radicalElectrons:
+                if radical1 == radical2: break
             else:
                 return False
+        #TODO: Lone pairs??
         #Each charge in self must have an equivalent in other (and vice versa)
         for charge1 in self.charge:
             for charge2 in other.charge:
@@ -355,11 +340,12 @@ def isSpecificCaseOf(self, other):
             else:
                 return False
         # Each free radical electron state in self must have an equivalent in other (and vice versa)
-        for radical1, spin1 in zip(self.radicalElectrons, self.spinMultiplicity): # all these must match
-            for radical2, spin2 in zip(other.radicalElectrons, other.spinMultiplicity): # can match any of these
-                if radical1 == radical2 and spin1 == spin2: break
+        for radical1 in self.radicalElectrons:
+            for radical2 in other.radicalElectrons:
+                if radical1 == radical2: break
             else:
                 return False
+        # TODO: Lone pairs??
         #Each charge in self must have an equivalent in other
         for charge1 in self.charge:
             for charge2 in other.charge:
@@ -514,6 +500,7 @@ class Group(Graph):
 
     def __init__(self, atoms=None):
         Graph.__init__(self, atoms)
+        self.multiplicity = [1,2,3,4,5]
         self.updateConnectivityValues()
         self.updateFingerprint()
     
@@ -671,7 +658,9 @@ def fromAdjacencyList(self, adjlist):
         ``False``.
         """
         from .adjlist import fromAdjacencyList
-        self.vertices = fromAdjacencyList(adjlist, group=True)
+        self.vertices, multiplicity = fromAdjacencyList(adjlist, group=True)
+        if multiplicity is not None:
+            self.multiplicity = multiplicity
         self.updateConnectivityValues()
         self.updateFingerprint()
         return self
@@ -681,7 +670,7 @@ def toAdjacencyList(self, label=''):
         Convert the molecular structure to a string adjacency list.
         """
         from .adjlist import toAdjacencyList
-        return toAdjacencyList(self.vertices, label='', group=True)
+        return toAdjacencyList(self.vertices, multiplicity=self.multiplicity, label='', group=True)
 
     def updateFingerprint(self):
         """
diff --git a/rmgpy/molecule/groupTest.py b/rmgpy/molecule/groupTest.py
index 5c12cbea2d..a62d623f93 100644
--- a/rmgpy/molecule/groupTest.py
+++ b/rmgpy/molecule/groupTest.py
@@ -18,7 +18,7 @@ def setUp(self):
         """
         A method called before each unit test in this class.
         """
-        self.atom = GroupAtom(atomType=[atomTypes['Cd']], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+        self.atom = GroupAtom(atomType=[atomTypes['Cd']], radicalElectrons=[1], charge=[0], label='*1')
     
     def testApplyActionBreakBond(self):
         """
@@ -26,7 +26,7 @@ def testApplyActionBreakBond(self):
         """
         action = ['BREAK_BOND', '*1', 'S', '*2']
         for label, atomType in atomTypes.iteritems():
-            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1')
             atom = atom0.copy()
             try:
                 atom.applyAction(action)
@@ -34,7 +34,6 @@ def testApplyActionBreakBond(self):
                 for a in atomType.breakBond:
                     self.assertTrue(a in atom.atomType)
                 self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-                self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
                 self.assertEqual(atom0.charge, atom.charge)
                 self.assertEqual(atom0.label, atom.label)
             except ActionError:
@@ -46,7 +45,7 @@ def testApplyActionFormBond(self):
         """
         action = ['FORM_BOND', '*1', 'S', '*2']
         for label, atomType in atomTypes.iteritems():
-            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1')
             atom = atom0.copy()
             try:
                 atom.applyAction(action)
@@ -54,7 +53,6 @@ def testApplyActionFormBond(self):
                 for a in atomType.formBond:
                     self.assertTrue(a in atom.atomType)
                 self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-                self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
                 self.assertEqual(atom0.charge, atom.charge)
                 self.assertEqual(atom0.label, atom.label)
             except ActionError:
@@ -66,7 +64,7 @@ def testApplyActionIncrementBond(self):
         """
         action = ['CHANGE_BOND', '*1', 1, '*2']
         for label, atomType in atomTypes.iteritems():
-            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1')
             atom = atom0.copy()
             try:
                 atom.applyAction(action)
@@ -74,7 +72,6 @@ def testApplyActionIncrementBond(self):
                 for a in atomType.incrementBond:
                     self.assertTrue(a in atom.atomType)
                 self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-                self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
                 self.assertEqual(atom0.charge, atom.charge)
                 self.assertEqual(atom0.label, atom.label)
             except ActionError:
@@ -86,7 +83,7 @@ def testApplyActionDecrementBond(self):
         """
         action = ['CHANGE_BOND', '*1', -1, '*2']
         for label, atomType in atomTypes.iteritems():
-            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1')
             atom = atom0.copy()
             try:
                 atom.applyAction(action)
@@ -94,7 +91,6 @@ def testApplyActionDecrementBond(self):
                 for a in atomType.decrementBond:
                     self.assertTrue(a in atom.atomType)
                 self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-                self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
                 self.assertEqual(atom0.charge, atom.charge)
                 self.assertEqual(atom0.label, atom.label)
             except ActionError:
@@ -106,7 +102,7 @@ def testApplyActionGainRadical(self):
         """
         action = ['GAIN_RADICAL', '*1', 1]
         for label, atomType in atomTypes.iteritems():
-            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1')
             atom = atom0.copy()
             try:
                 atom.applyAction(action)
@@ -114,7 +110,6 @@ def testApplyActionGainRadical(self):
                 for a in atomType.incrementRadical:
                     self.assertTrue(a in atom.atomType, "GAIN_RADICAL on {0} gave {1} not {2}".format(atomType, atom.atomType, atomType.incrementRadical))
                 self.assertEqual(atom0.radicalElectrons, [r - 1 for r in atom.radicalElectrons])
-                self.assertEqual(atom0.spinMultiplicity, [s - 1 for s in atom.spinMultiplicity])
                 self.assertEqual(atom0.charge, atom.charge)
                 self.assertEqual(atom0.label, atom.label)
             except ActionError:
@@ -126,7 +121,7 @@ def testApplyActionLoseRadical(self):
         """
         action = ['LOSE_RADICAL', '*1', 1]
         for label, atomType in atomTypes.iteritems():
-            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+            atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1')
             atom = atom0.copy()
             try:
                 atom.applyAction(action)
@@ -134,7 +129,6 @@ def testApplyActionLoseRadical(self):
                 for a in atomType.incrementRadical:
                     self.assertTrue(a in atom.atomType, "LOSE_RADICAL on {0} gave {1} not {2}".format(atomType, atom.atomType, atomType.decrementRadical))
                 self.assertEqual(atom0.radicalElectrons, [r + 1 for r in atom.radicalElectrons])
-                self.assertEqual(atom0.spinMultiplicity, [s + 1 for s in atom.spinMultiplicity])
                 self.assertEqual(atom0.charge, atom.charge)
                 self.assertEqual(atom0.label, atom.label)
             except ActionError:
@@ -146,8 +140,8 @@ def testEquivalent(self):
         """
         for label1, atomType1 in atomTypes.iteritems():
             for label2, atomType2 in atomTypes.iteritems():
-                atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
-                atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+                atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], charge=[0], label='*1')
+                atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], charge=[0], label='*1')
                 if label1 == label2 or atomType2 in atomType1.generic or atomType1 in atomType2.generic:
                     self.assertTrue(atom1.equivalent(atom2), '{0!s} is not equivalent to {1!s}'.format(atom1, atom2))
                     self.assertTrue(atom2.equivalent(atom1), '{0!s} is not equivalent to {1!s}'.format(atom2, atom1))
@@ -157,7 +151,7 @@ def testEquivalent(self):
             # Now see if charge and radical count are checked properly
             for charge in range(3):
                 for radicals in range(2):
-                    atom3 = GroupAtom(atomType=[atomType1], radicalElectrons=[radicals], spinMultiplicity=[radicals + 1], charge=[charge], label='*1')
+                    atom3 = GroupAtom(atomType=[atomType1], radicalElectrons=[radicals], charge=[charge], label='*1')
                     if radicals == 1 and charge == 0:
                         self.assertTrue(atom1.equivalent(atom3), '{0!s} is not equivalent to {1!s}'.format(atom1, atom3))
                         self.assertTrue(atom1.equivalent(atom3), '{0!s} is not equivalent to {1!s}'.format(atom3, atom1))
@@ -171,11 +165,11 @@ def testIsSpecificCaseOf(self):
         """
         for label1, atomType1 in atomTypes.iteritems():
             for label2, atomType2 in atomTypes.iteritems():
-                atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
-                atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1')
+                atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], charge=[0], label='*1')
+                atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], charge=[0], label='*1')
                 # And make more generic types of these two atoms
-                atom1gen = GroupAtom(atomType=[atomType1], radicalElectrons=[0, 1], spinMultiplicity=[1, 2], charge=[0, 1], label='*1')
-                atom2gen = GroupAtom(atomType=[atomType2], radicalElectrons=[0, 1], spinMultiplicity=[1, 2], charge=[0, 1], label='*1')
+                atom1gen = GroupAtom(atomType=[atomType1], radicalElectrons=[0, 1], charge=[0, 1], label='*1')
+                atom2gen = GroupAtom(atomType=[atomType2], radicalElectrons=[0, 1], charge=[0, 1], label='*1')
                 if label1 == label2 or atomType2 in atomType1.generic:
                     self.assertTrue(atom1.isSpecificCaseOf(atom2), '{0!s} is not a specific case of {1!s}'.format(atom1, atom2))
                     self.assertTrue(atom1.isSpecificCaseOf(atom2gen), '{0!s} is not a specific case of {1!s}'.format(atom1, atom2gen))
@@ -193,7 +187,6 @@ def testCopy(self):
         self.assertEqual(len(self.atom.atomType), len(atom.atomType))
         self.assertEqual(self.atom.atomType[0].label, atom.atomType[0].label)
         self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons)
-        self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity)
         self.assertEqual(self.atom.charge, atom.charge)
         self.assertEqual(self.atom.label, atom.label)
     
@@ -207,7 +200,6 @@ def testPickle(self):
         self.assertEqual(len(self.atom.atomType), len(atom.atomType))
         self.assertEqual(self.atom.atomType[0].label, atom.atomType[0].label)
         self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons)
-        self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity)
         self.assertEqual(self.atom.charge, atom.charge)
         self.assertEqual(self.atom.label, atom.label)
 
@@ -362,9 +354,9 @@ class TestGroup(unittest.TestCase):
 
     def setUp(self):
         self.adjlist = """
-1 *2 {Cs,Cd} 0 {2,{S,D}} {3,S}
-2 *1 {Os,Od} 0 {1,{S,D}}
-3    R!H     0 {1,S}
+1 *2 [Cs,Cd] u0 {2,[S,D]} {3,S}
+2 *1 [Os,Od] u0 {1,[S,D]}
+3    R!H     u0 {1,S}
             """
         self.group = Group().fromAdjacencyList(self.adjlist)
         
@@ -454,9 +446,9 @@ def testIsIsomorphic(self):
         Test the Group.isIsomorphic() method.
         """
         adjlist = """
-1  *1 {Os,Od} 0 {3,{S,D}}
-2     R!H 0 {3,S}
-3  *2 {Cs,Cd} 0 {1,{S,D}} {2,S}
+1  *1 [Os,Od] u0 {3,[S,D]}
+2     R!H     u0 {3,S}
+3  *2 [Cs,Cd] u0 {1,[S,D]} {2,S}
             """
         group = Group().fromAdjacencyList(adjlist)
         self.assertTrue(self.group.isIsomorphic(group))
@@ -467,9 +459,9 @@ def testFindIsomorphism(self):
         Test the Group.findIsomorphism() method.
         """
         adjlist = """
-1  *1 {Os,Od} 0 {3,{S,D}}
-2     R!H 0 {3,S}
-3  *2 {Cs,Cd} 0 {1,{S,D}} {2,S}
+1  *1 [Os,Od] u0 {3,[S,D]}
+2     R!H     u0 {3,S}
+3  *2 [Cs,Cd] u0 {1,[S,D]} {2,S}
             """
         group = Group().fromAdjacencyList(adjlist)
         result = self.group.findIsomorphism(group)
@@ -491,7 +483,7 @@ def testIsSubgraphIsomorphic(self):
         Test the Group.isSubgraphIsomorphic() method.
         """
         adjlist = """
-1  *1 {Cs,Cd} 0
+1  *1 [Cs,Cd] u0
             """
         group = Group().fromAdjacencyList(adjlist)
         self.assertTrue(self.group.isSubgraphIsomorphic(group))
@@ -502,7 +494,7 @@ def testFindSubgraphIsomorphisms(self):
         Test the Group.findSubgraphIsomorphisms() method.
         """
         adjlist = """
-1  *1 {Cs,Cd} 0
+1  *1 [Cs,Cd] u0
             """
         group = Group().fromAdjacencyList(adjlist)
         result = self.group.findSubgraphIsomorphisms(group)
diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd
index 3e474908b8..ee47e73907 100644
--- a/rmgpy/molecule/molecule.pxd
+++ b/rmgpy/molecule/molecule.pxd
@@ -38,7 +38,6 @@ cdef class Atom(Vertex):
 
     cdef public Element element
     cdef public short radicalElectrons
-    cdef public short spinMultiplicity
     cdef public short charge
     cdef public str label
     cdef public AtomType atomType
@@ -103,6 +102,7 @@ cdef class Molecule(Graph):
 
     cdef public bint implicitHydrogens
     cdef public int symmetryNumber
+    cdef public int multiplicity
     cdef public object rdMol
     cdef public int rdMolConfId
     cdef str _fingerprint
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 7e81b424df..fcf4247ccd 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -65,7 +65,7 @@
                  'He': '[He]',
                  'CH4O': 'CO',
                  'CO2': 'O=C=O',
-                 'CO': 'C#O',
+                 'CO': '[C+]#[O-]',
                  'C2H4': 'C=C',
                  'O2': 'O=O'
              }
@@ -97,7 +97,6 @@ class Atom(Vertex):
     `atomType`          :class:`AtomType`   The :ref:`atom type <atom-types>`
     `element`           :class:`Element`    The chemical element the atom represents
     `radicalElectrons`  ``short``           The number of radical electrons
-    `spinMultiplicity`  ``short``           The spin multiplicity of the atom
     `charge`            ``short``           The formal charge of the atom
     `label`             ``str``             A string label that can be used to tag individual atoms
     `coords`            ``numpy array``     The (x,y,z) coordinates in Angstrom
@@ -109,14 +108,13 @@ class Atom(Vertex):
     e.g. ``atom.symbol`` instead of ``atom.element.symbol``.
     """
 
-    def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label='', lonePairs=0, coords=numpy.array([])):
+    def __init__(self, element=None, radicalElectrons=0, charge=0, label='', lonePairs=-100, coords=numpy.array([])):
         Vertex.__init__(self)
         if isinstance(element, str):
             self.element = elements.__dict__[element]
         else:
             self.element = element
         self.radicalElectrons = radicalElectrons
-        self.spinMultiplicity = spinMultiplicity
         self.charge = charge
         self.label = label
         self.atomType = None
@@ -150,8 +148,9 @@ def __reduce__(self):
             'connectivity3': self.connectivity3,
             'sortingLabel': self.sortingLabel,
             'atomType': self.atomType.label if self.atomType else None,
+            'lonePairs': self.lonePairs,
         }
-        return (Atom, (self.element.symbol, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label), d)
+        return (Atom, (self.element.symbol, self.radicalElectrons, self.charge, self.label), d)
 
     def __setstate__(self, d):
         """
@@ -163,6 +162,7 @@ def __setstate__(self, d):
         self.connectivity3 = d['connectivity3']
         self.sortingLabel = d['sortingLabel']
         self.atomType = atomTypes[d['atomType']] if d['atomType'] else None
+        self.lonePairs = d['lonePairs']
     
     @property
     def mass(self): return self.element.mass
@@ -189,7 +189,6 @@ def equivalent(self, other):
             atom = other
             return (self.element is atom.element and
                 self.radicalElectrons == atom.radicalElectrons and
-                self.spinMultiplicity == atom.spinMultiplicity and
                 self.charge == atom.charge)
         elif isinstance(other, GroupAtom):
             cython.declare(a=AtomType, radical=cython.short, spin=cython.short, charge=cython.short)
@@ -198,8 +197,8 @@ def equivalent(self, other):
                 if self.atomType.equivalent(a): break
             else:
                 return False
-            for radical, spin in zip(ap.radicalElectrons, ap.spinMultiplicity):
-                if self.radicalElectrons == radical and self.spinMultiplicity == spin: break
+            for radical in ap.radicalElectrons:
+                if self.radicalElectrons == radical: break
             else:
                 return False
             for charge in ap.charge:
@@ -219,7 +218,7 @@ def isSpecificCaseOf(self, other):
         if isinstance(other, Atom):
             return self.equivalent(other)
         elif isinstance(other, GroupAtom):
-            cython.declare(atom=GroupAtom, a=AtomType, radical=cython.short, spin=cython.short, charge=cython.short, index=cython.int)
+            cython.declare(atom=GroupAtom, a=AtomType, radical=cython.short, charge=cython.short)
             atom = other
             if self.atomType is None:
                 return False
@@ -227,12 +226,11 @@ def isSpecificCaseOf(self, other):
                 if self.atomType.isSpecificCaseOf(a): break
             else:
                 return False
-            for index in range(len(atom.radicalElectrons)):
-                radical = atom.radicalElectrons[index]
-                spin = atom.spinMultiplicity[index]
-                if self.radicalElectrons == radical and self.spinMultiplicity == spin: break
+            for radical in atom.radicalElectrons:
+                if self.radicalElectrons == radical: break
             else:
                 return False
+            #TODO: lone pairs and charge
 # until we have charges and lone pairs in the group values we neglect them here
 #            for charge in atom.charge:
 #                if self.charge == charge: break
@@ -252,7 +250,6 @@ def copy(self):
         a.resetConnectivityValues()
         a.element = self.element
         a.radicalElectrons = self.radicalElectrons
-        a.spinMultiplicity = self.spinMultiplicity
         a.charge = self.charge
         a.label = self.label
         a.atomType = self.atomType
@@ -300,16 +297,10 @@ def incrementRadical(self):
         Update the atom pattern as a result of applying a GAIN_RADICAL action,
         where `radical` specifies the number of radical electrons to add.
         """
-        # Set the new radical electron counts and spin multiplicities
+        # Set the new radical electron count
         self.radicalElectrons += 1
         if self.radicalElectrons <= 0:
             raise ActionError('Unable to update Atom due to GAIN_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons))
-        elif self.radicalElectrons == 1:
-            self.spinMultiplicity = 2
-        elif self.radicalElectrons == 2:
-            self.spinMultiplicity = 3   # Assume this always results in the triplet, as they tend to be more stable than the singlet (though there are exceptions!)
-        else:
-            self.spinMultiplicity = self.radicalElectrons + 1
 
     def decrementRadical(self):
         """
@@ -317,18 +308,10 @@ def decrementRadical(self):
         where `radical` specifies the number of radical electrons to remove.
         """
         cython.declare(radicalElectrons=cython.short)
-        # Set the new radical electron counts and spin multiplicities
+        # Set the new radical electron count
         radicalElectrons = self.radicalElectrons = self.radicalElectrons - 1
         if radicalElectrons  < 0:
             raise ActionError('Unable to update Atom due to LOSE_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons))
-        elif radicalElectrons == 0:
-            self.spinMultiplicity = 1
-        elif radicalElectrons == 1:
-            self.spinMultiplicity = 2
-        elif radicalElectrons == 2:
-            self.spinMultiplicity = 3   # Assume this always results in the triplet, as they tend to be more stable than the singlet (though there are exceptions!)
-        else:
-            self.spinMultiplicity = self.radicalElectrons + 1
 
     def setLonePairs(self,lonePairs):
         """
@@ -361,6 +344,10 @@ def decrementLonePairs(self):
         self.updateCharge()
         
     def updateCharge(self):
+        """
+        Update self.charge, according to the valence, and the 
+        number and types of bonds, radicals, and lone pairs.
+        """
         valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0, 'Cl': 1}
         orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
         valence = valences[self.symbol]
@@ -401,6 +388,7 @@ def setSpinMultiplicity(self,spinMultiplicity):
         """
         Set the spin multiplicity.
         """
+        raise NotImplementedError("I thought multiplicity was now a molecule attribute not atom?")
         # Set the spin multiplicity
         self.spinMultiplicity = spinMultiplicity
         if self.spinMultiplicity < 0:
@@ -592,19 +580,23 @@ class Molecule(Graph):
     Attribute               Type        Description
     ======================= =========== ========================================
     `symmetryNumber`        ``int``     The (estimated) external + internal symmetry number of the molecule
+    `multiplicity`          ``int``     The multiplicity of this species, multiplicity = 2*total_spin+1
     ======================= =========== ========================================
 
     A new molecule object can be easily instantiated by passing the `SMILES` or
     `InChI` string representing the molecular structure.
     """
 
-    def __init__(self, atoms=None, symmetry=1, SMILES='', InChI='', SMARTS = ''):
+    def __init__(self, atoms=None, symmetry=1, multiplicity=-187, SMILES='', InChI='', SMARTS=''):
         Graph.__init__(self, atoms)
         self.symmetryNumber = symmetry
+        self.multiplicity = multiplicity
         self._fingerprint = None
         if SMILES != '': self.fromSMILES(SMILES)
         elif InChI != '': self.fromInChI(InChI)
         elif SMARTS != '': self.fromSMARTS(SMARTS)
+        if multiplicity != -187:  # it was set explicitly, so re-set it (fromSMILES etc may have changed it)
+            self.multiplicity = multiplicity
     
     def __str__(self):
         """
@@ -616,13 +608,17 @@ def __repr__(self):
         """
         Return a representation that can be used to reconstruct the object.
         """
+        cython.declare(multiplicity=cython.int)
+        multiplicity = self.multiplicity
+        if multiplicity != self.getRadicalCount() + 1:
+            return 'Molecule(SMILES="{0}", multiplicity={1:d})'.format(self.toSMILES(), multiplicity)
         return 'Molecule(SMILES="{0}")'.format(self.toSMILES())
 
     def __reduce__(self):
         """
         A helper function used when pickling an object.
         """
-        return (Molecule, (self.vertices, self.symmetryNumber))
+        return (Molecule, (self.vertices, self.symmetryNumber, self.multiplicity))
 
     def __getAtoms(self): return self.vertices
     def __setAtoms(self, atoms): self.vertices = atoms
@@ -788,6 +784,7 @@ def copy(self, deep=False):
         other = cython.declare(Molecule)
         g = Graph.copy(self, deep)
         other = Molecule(g.vertices)
+        other.multiplicity = self.multiplicity
         return other
 
     def merge(self, other):
@@ -946,6 +943,7 @@ def isIsomorphic(self, other, initialMap=None):
         mapping from `self` to `other` (i.e. the atoms of `self` are the keys,
         while the atoms of `other` are the values). The `other` parameter must
         be a :class:`Molecule` object, or a :class:`TypeError` is raised.
+        Also ensures multiplicities are also equal.
         """
         # It only makes sense to compare a Molecule to a Molecule for full
         # isomorphism, so raise an exception if this is not what was requested
@@ -956,6 +954,9 @@ def isIsomorphic(self, other, initialMap=None):
         # sufficient!) condition for the associated molecules to be isomorphic
         if self.getFingerprint() != other.getFingerprint():
             return False
+        # check multiplicity
+        if self.multiplicity != other.multiplicity:
+            return False
         # Do the full isomorphism comparison
         result = Graph.isIsomorphic(self, other, initialMap)
         return result
@@ -1014,7 +1015,8 @@ def isSubgraphIsomorphic(self, other, initialMap=None):
             carbonCount < group.carbonCount or
             nitrogenCount < group.nitrogenCount or
             oxygenCount < group.oxygenCount or
-            sulfurCount < group.sulfurCount):
+            sulfurCount < group.sulfurCount or
+            self.multiplicity not in group.multiplicity):
             return False
 
         # Do the isomorphism comparison
@@ -1102,7 +1104,22 @@ def fromSMILES(self, smilesstr):
         # Special handling of helium
         if smilesstr == '[He]':
             # RDKit improperly handles helium and returns it in a triplet state
-            self.fromAdjacencyList('1 He 0 1')
+            self.fromAdjacencyList(
+            """
+            He
+            multiplicity 1
+            1 He u0 p1
+            """)
+            return self
+        elif smilesstr == '[C+]#[O-]':
+            #  carbon monoxide
+            self.fromAdjacencyList(
+            """
+            CO
+            multiplicity 1
+            1 C u0 p1 {2,T}
+            2 O u0 p1 {1,T}
+            """)
             return self
         
         else:
@@ -1130,7 +1147,7 @@ def fromRDKitMol(self, rdkitmol):
         """
         # Below are the declared variables for cythonizing the module
         cython.declare(i=cython.int)
-        cython.declare(radicalElectrons=cython.int, spinMultiplicity=cython.int, charge=cython.int, lonePairs=cython.int)
+        cython.declare(radicalElectrons=cython.int, charge=cython.int, lonePairs=cython.int)
         cython.declare(atom=Atom, atom1=Atom, atom2=Atom, bond=Bond)
         
         self.vertices = []
@@ -1139,26 +1156,19 @@ def fromRDKitMol(self, rdkitmol):
         rdkitmol = Chem.AddHs(rdkitmol)
         Chem.rdmolops.Kekulize(rdkitmol, clearAromaticFlags=True)
         
-        # iterate though atoms in rdkitmol
+        # iterate through atoms in rdkitmol
         for i in range(rdkitmol.GetNumAtoms()):
             rdkitatom = rdkitmol.GetAtomWithIdx(i)
             
             # Use atomic number as key for element
             number = rdkitatom.GetAtomicNum()
             element = elements.getElement(number)
-            
-            # Process spin multiplicity
-            radicalElectrons = rdkitatom.GetNumRadicalElectrons()
-            
-            # Assume this is always true
-            # There are cases where 2 radicalElectrons is a singlet, but
-            # the triplet is often more stable, 
-            spinMultiplicity = radicalElectrons + 1
                 
             # Process charge
             charge = rdkitatom.GetFormalCharge()
+            radicalElectrons = rdkitatom.GetNumRadicalElectrons()
             
-            atom = Atom(element, radicalElectrons, spinMultiplicity, charge, '', 0)
+            atom = Atom(element, radicalElectrons, charge, '', 0)
             self.vertices.append(atom)
             
             # Add bonds by iterating again through atoms
@@ -1183,9 +1193,13 @@ def fromRDKitMol(self, rdkitmol):
         self.updateLonePairs()
         self.updateAtomTypes()
         
+        # Assume this is always true
+        # There are cases where 2 radicalElectrons is a singlet, but
+        # the triplet is often more stable, 
+        self.multiplicity = self.getRadicalCount() + 1
+        
         return self
 
-    newStyleAdjMatcher = re.compile('\s*\d+\s+(\*\d*\s+)?[A-Za-z]+\s+\S+\s+\d').match
     def fromAdjacencyList(self, adjlist, saturateH=False):
         """
         Convert a string adjacency list `adjlist` to a molecular structure.
@@ -1193,14 +1207,17 @@ def fromAdjacencyList(self, adjlist, saturateH=False):
         ``False``.
         """
         from .adjlist import fromAdjacencyList
-        if not self.newStyleAdjMatcher(adjlist[adjlist.find('1 '):]):
-            logging.warning("There could be an old-style adjacency list. Please check library and input before proceeding!")
-            print adjlist
-            
-        self.vertices = fromAdjacencyList(adjlist, False, saturateH=saturateH)
+        
+        self.vertices, self.multiplicity = fromAdjacencyList(adjlist, group=False, saturateH=saturateH)
         self.updateConnectivityValues()
         self.updateAtomTypes()
         
+        # Check if multiplicity is possible
+        n_rad = self.getRadicalCount() 
+        multiplicity = self.multiplicity
+        if not (n_rad + 1 == multiplicity or n_rad - 1 == multiplicity or n_rad - 3 == multiplicity or n_rad - 5 == multiplicity):
+            raise ValueError('Impossible multiplicity for molecule\n{0}\n multiplicity = {1} and number of unpaired electrons = {2}'.format(self.toAdjacencyList(),multiplicity,n_rad))
+
         return self
         
     def fromXYZ(self, atomicNums, coordinates):
@@ -1450,7 +1467,7 @@ def toAdjacencyList(self, label='', removeH=False, removeLonePairs=False):
         Convert the molecular structure to a string adjacency list.
         """
         from .adjlist import toAdjacencyList
-        result = toAdjacencyList(self.vertices, label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs)
+        result = toAdjacencyList(self.vertices, self.multiplicity,  label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs)
         return result
 
     def isLinear(self):
@@ -1883,38 +1900,38 @@ def getURL(self):
     
     def isBiradicalSinglet(self):
         """
-        Return ``True`` if the molecule is a 1-centered biradical in singlet state,
+        Return ``True`` if the molecule is a 1-centered biradical, and the molecule is in singlet state,
         or ``False`` otherwise.
         """
         cython.declare(atom=Atom)
         for atom in self.vertices:
             if atom.radicalElectrons == 2:
-                if atom.spinMultiplicity == 1:
+                if self.multiplicity == 1:
                     return True
         return False
     
     def isBiradicalTriplet(self):
         """
-        Return ``True`` if the molecule is a 1-centered biradical in triplet state,
+        Return ``True`` if the molecule is a 1-centered biradical, and the molecule is in triplet state,
         or ``False`` otherwise.
         """
         cython.declare(atom=Atom)
         for atom in self.vertices:
             if atom.radicalElectrons == 2:
-                if atom.spinMultiplicity == 3:
+                if self.multiplicity == 3:
                     return True
         return False
     
     def changeTripletSinglet(self):
         """
-        If the molecule is a 1-centered biradical in triplet state,
+        If the molecule is a 1-centered biradical, and the molecule is in a triplet state,
         change it to singlet state.
         """
         cython.declare(atom=Atom)
         for atom in self.vertices:
             if atom.radicalElectrons == 2:
-                if atom.spinMultiplicity == 3:
-                    atom.spinMultiplicity = 1
+                if self.multiplicity == 3:
+                    self.multiplicity = 1
                     
     def getRadicalAtoms(self):
         """
diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py
index fcb085f5e0..bfd6130ab6 100644
--- a/rmgpy/molecule/moleculeTest.py
+++ b/rmgpy/molecule/moleculeTest.py
@@ -18,7 +18,7 @@ def setUp(self):
         """
         A method called before each unit test in this class.
         """
-        self.atom = Atom(element=getElement('C'), radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+        self.atom = Atom(element=getElement('C'), radicalElectrons=1, charge=0, label='*1', lonePairs=0)
     
     def testMass(self):
         """
@@ -43,7 +43,7 @@ def testIsHydrogen(self):
         Test the Atom.isHydrogen() method.
         """
         for element in elementList:
-            atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             if element.symbol == 'H':
                 self.assertTrue(atom.isHydrogen())
             else:
@@ -54,7 +54,7 @@ def testIsNonHydrogen(self):
         Test the Atom.isNonHydrogen() method.
         """
         for element in elementList:
-            atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             if element.symbol == 'H':
                 self.assertFalse(atom.isNonHydrogen())
             else:
@@ -65,7 +65,7 @@ def testIsCarbon(self):
         Test the Atom.isCarbon() method.
         """
         for element in elementList:
-            atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             if element.symbol == 'C':
                 self.assertTrue(atom.isCarbon())
             else:
@@ -76,7 +76,7 @@ def testIsOxygen(self):
         Test the Atom.isOxygen() method.
         """
         for element in elementList:
-            atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=2)
             if element.symbol == 'O':
                 self.assertTrue(atom.isOxygen())
             else:
@@ -87,20 +87,16 @@ def testIncrementRadical(self):
         Test the Atom.incrementRadical() method.
         """
         radicalElectrons = self.atom.radicalElectrons
-        spinMultiplicity = self.atom.spinMultiplicity
         self.atom.incrementRadical()
         self.assertEqual(self.atom.radicalElectrons, radicalElectrons + 1)
-        self.assertEqual(self.atom.spinMultiplicity, spinMultiplicity + 1)
     
     def testDecrementRadical(self):
         """
         Test the Atom.decrementRadical() method.
         """
         radicalElectrons = self.atom.radicalElectrons
-        spinMultiplicity = self.atom.spinMultiplicity
         self.atom.decrementRadical()
         self.assertEqual(self.atom.radicalElectrons, radicalElectrons - 1)
-        self.assertEqual(self.atom.spinMultiplicity, spinMultiplicity - 1)
            
     def testApplyActionBreakBond(self):
         """
@@ -108,12 +104,11 @@ def testApplyActionBreakBond(self):
         """
         action = ['BREAK_BOND', '*1', 'S', '*2']
         for element in elementList:
-            atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             atom = atom0.copy()
             atom.applyAction(action)
             self.assertEqual(atom0.element, atom.element)
             self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-            self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
             self.assertEqual(atom0.charge, atom.charge)
             self.assertEqual(atom0.label, atom.label)
     
@@ -123,12 +118,11 @@ def testApplyActionFormBond(self):
         """
         action = ['FORM_BOND', '*1', 'S', '*2']
         for element in elementList:
-            atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             atom = atom0.copy()
             atom.applyAction(action)
             self.assertEqual(atom0.element, atom.element)
             self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-            self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
             self.assertEqual(atom0.charge, atom.charge)
             self.assertEqual(atom0.label, atom.label)
     
@@ -138,12 +132,11 @@ def testApplyActionIncrementBond(self):
         """
         action = ['CHANGE_BOND', '*1', 1, '*2']
         for element in elementList:
-            atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             atom = atom0.copy()
             atom.applyAction(action)
             self.assertEqual(atom0.element, atom.element)
             self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-            self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
             self.assertEqual(atom0.charge, atom.charge)
             self.assertEqual(atom0.label, atom.label)
     
@@ -153,12 +146,11 @@ def testApplyActionDecrementBond(self):
         """
         action = ['CHANGE_BOND', '*1', -1, '*2']
         for element in elementList:
-            atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             atom = atom0.copy()
             atom.applyAction(action)
             self.assertEqual(atom0.element, atom.element)
             self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons)
-            self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity)
             self.assertEqual(atom0.charge, atom.charge)
             self.assertEqual(atom0.label, atom.label)
     
@@ -168,12 +160,11 @@ def testApplyActionGainRadical(self):
         """
         action = ['GAIN_RADICAL', '*1', 1]
         for element in elementList:
-            atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             atom = atom0.copy()
             atom.applyAction(action)
             self.assertEqual(atom0.element, atom.element)
             self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons - 1)
-            self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity - 1)
             self.assertEqual(atom0.charge, atom.charge)
             self.assertEqual(atom0.label, atom.label)
     
@@ -183,12 +174,11 @@ def testApplyActionLoseRadical(self):
         """
         action = ['LOSE_RADICAL', '*1', 1]
         for element in elementList:
-            atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+            atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
             atom = atom0.copy()
             atom.applyAction(action)
             self.assertEqual(atom0.element, atom.element)
             self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons + 1)
-            self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity + 1)
             self.assertEqual(atom0.charge, atom.charge)
             self.assertEqual(atom0.label, atom.label)
     
@@ -198,8 +188,8 @@ def testEquivalent(self):
         """
         for index1, element1 in enumerate(elementList[0:10]):
             for index2, element2 in enumerate(elementList[0:10]):
-                atom1 = Atom(element=element1, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
-                atom2 = Atom(element=element2, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+                atom1 = Atom(element=element1, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
+                atom2 = Atom(element=element2, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
                 if index1 == index2:
                     self.assertTrue(atom1.equivalent(atom2))
                     self.assertTrue(atom2.equivalent(atom1))
@@ -213,8 +203,8 @@ def testIsSpecificCaseOf(self):
         """
         for index1, element1 in enumerate(elementList[0:10]):
             for index2, element2 in enumerate(elementList[0:10]):
-                atom1 = Atom(element=element1, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
-                atom2 = Atom(element=element2, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1')
+                atom1 = Atom(element=element1, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
+                atom2 = Atom(element=element2, radicalElectrons=1, charge=0, label='*1', lonePairs=0)
                 if index1 == index2:
                     self.assertTrue(atom1.isSpecificCaseOf(atom2))
                 else:
@@ -228,7 +218,6 @@ def testCopy(self):
         self.assertEqual(self.atom.element.symbol, atom.element.symbol)
         self.assertEqual(self.atom.atomType, atom.atomType)
         self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons)
-        self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity)
         self.assertEqual(self.atom.charge, atom.charge)
         self.assertEqual(self.atom.label, atom.label)
     
@@ -242,7 +231,6 @@ def testPickle(self):
         self.assertEqual(self.atom.element.symbol, atom.element.symbol)
         self.assertEqual(self.atom.atomType, atom.atomType)
         self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons)
-        self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity)
         self.assertEqual(self.atom.charge, atom.charge)
         self.assertEqual(self.atom.label, atom.label)
         
@@ -468,45 +456,73 @@ class TestMolecule(unittest.TestCase):
     """
     
     def setUp(self):
-        self.adjlist = """
-1 *2 C 1 0 {2,D} {3,S}
-2 *1 O 0 2 {1,D}
-3    C 0 0 {1,S} {4,S} {5,S} {6,S}
-4    H 0 0 {3,S}
-5    H 0 0 {3,S}
-6    H 0 0 {3,S}
+        self.adjlist_1 = """
+1 *1 C u1 p0 c0  {2,S} {3,S} {4,S}
+2    H u0 p0 c0  {1,S}
+3    H u0 p0 c0  {1,S}
+4 *2 N u0 p0 c+1 {1,S} {5,S} {6,D}
+5    O u0 p3 c-1 {4,S}
+6    O u0 p2 c0  {4,D}
             """
-        self.molecule = Molecule().fromAdjacencyList(self.adjlist)
+        self.molecule = [Molecule().fromAdjacencyList(self.adjlist_1)]
+        
+        self.adjlist_2 = """
+1 *1 C u1 p0 {2,S} {3,S} {4,S}
+2    H u0 p0 {1,S}
+3    H u0 p0 {1,S}
+4 *2 N u0 p0 {1,S} {5,S} {6,D}
+5    O u0 p3 {4,S}
+6    O u0 p2 {4,D}
+            """
+        self.molecule.append(Molecule().fromAdjacencyList(self.adjlist_2))
+        
+        self.adjlist_3 = """
+1 *1 C u1 {2,S} {3,S} {4,S}
+2    H u0 {1,S}
+3    H u0 {1,S}
+4 *2 N u0 {1,S} {5,S} {6,D}
+5    O u0 {4,S}
+6    O u0 {4,D}
+            """
+        self.molecule.append(Molecule().fromAdjacencyList(self.adjlist_3))
+        
+        self.adjlist_4 = """
+1 *1 C u1 p0 {2,S}
+2 *2 N u0 p0 {1,S} {3,S} {4,D}
+3    O u0 p3 {2,S}
+4    O u0 p2 {2,D}
+            """
+        self.molecule.append(Molecule().fromAdjacencyList(self.adjlist_4,saturateH=True))
         
     def testClearLabeledAtoms(self):
         """
         Test the Molecule.clearLabeledAtoms() method.
         """
-        self.molecule.clearLabeledAtoms()
-        for atom in self.molecule.atoms:
+        self.molecule[0].clearLabeledAtoms()
+        for atom in self.molecule[0].atoms:
             self.assertEqual(atom.label, '')
 
     def testContainsLabeledAtom(self):
         """
         Test the Molecule.containsLabeledAtom() method.
         """
-        for atom in self.molecule.atoms:
+        for atom in self.molecule[0].atoms:
             if atom.label != '':
-                self.assertTrue(self.molecule.containsLabeledAtom(atom.label))
-        self.assertFalse(self.molecule.containsLabeledAtom('*3'))
-        self.assertFalse(self.molecule.containsLabeledAtom('*4'))
-        self.assertFalse(self.molecule.containsLabeledAtom('*5'))
-        self.assertFalse(self.molecule.containsLabeledAtom('*6'))
+                self.assertTrue(self.molecule[0].containsLabeledAtom(atom.label))
+        self.assertFalse(self.molecule[0].containsLabeledAtom('*3'))
+        self.assertFalse(self.molecule[0].containsLabeledAtom('*4'))
+        self.assertFalse(self.molecule[0].containsLabeledAtom('*5'))
+        self.assertFalse(self.molecule[0].containsLabeledAtom('*6'))
         
     def testGetLabeledAtom(self):
         """
         Test the Molecule.getLabeledAtom() method.
         """
-        for atom in self.molecule.atoms:
+        for atom in self.molecule[0].atoms:
             if atom.label != '':
-                self.assertEqual(atom, self.molecule.getLabeledAtom(atom.label))
+                self.assertEqual(atom, self.molecule[0].getLabeledAtom(atom.label))
         try:
-            self.molecule.getLabeledAtom('*3')
+            self.molecule[0].getLabeledAtom('*3')
             self.fail('Unexpected successful return from Molecule.getLabeledAtom() with invalid atom label.')
         except ValueError:
             pass
@@ -515,8 +531,8 @@ def testGetLabeledAtoms(self):
         """
         Test the Molecule.getLabeledAtoms() method.
         """
-        labeled = self.molecule.getLabeledAtoms()
-        for atom in self.molecule.atoms:
+        labeled = self.molecule[0].getLabeledAtoms()
+        for atom in self.molecule[0].atoms:
             if atom.label != '':
                 self.assertTrue(atom.label in labeled)
                 self.assertTrue(atom in labeled.values())
@@ -528,67 +544,223 @@ def testGetFormula(self):
         """
         Test the Molecule.getLabeledAtoms() method.
         """
-        self.assertEqual(self.molecule.getFormula(), 'C2H3O')
+        self.assertEqual(self.molecule[0].getFormula(), 'CH2NO2')
+        self.assertEqual(self.molecule[1].getFormula(), 'CH2NO2')
+        self.assertEqual(self.molecule[2].getFormula(), 'CH2NO2')
+        self.assertEqual(self.molecule[3].getFormula(), 'CH2NO2')
 
 
     def testRadicalCount(self):
         """
         Test the Molecule.getRadicalCount() method.
         """
-        self.assertEqual( self.molecule.getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule.atoms]) )
+        self.assertEqual( self.molecule[0].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[0].atoms]) )
+        self.assertEqual( self.molecule[1].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[1].atoms]) )
+        self.assertEqual( self.molecule[2].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[2].atoms]) )
+        self.assertEqual( self.molecule[3].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[3].atoms]) )
         
     def testGetMolecularWeight(self):
         """
         Test the Molecule.getMolecularWeight() method.
         """
-        self.assertAlmostEqual(self.molecule.getMolecularWeight() * 1000, 43.04, 2)
+        self.assertAlmostEqual(self.molecule[0].getMolecularWeight() * 1000, 60.03, 2)
+        self.assertAlmostEqual(self.molecule[1].getMolecularWeight() * 1000, 60.03, 2)
+        self.assertAlmostEqual(self.molecule[2].getMolecularWeight() * 1000, 60.03, 2)
+        self.assertAlmostEqual(self.molecule[3].getMolecularWeight() * 1000, 60.03, 2)
 
     def testFromAdjacencyList(self):
         """
         Test the Molecule.fromAdjacencyList() method.
         """
-        atom1, atom2, atom3 = self.molecule.atoms[0:3]
-        self.assertTrue(self.molecule.hasBond(atom1,atom2))
-        self.assertTrue(self.molecule.hasBond(atom1,atom3))
-        self.assertFalse(self.molecule.hasBond(atom2,atom3))
-        bond12 = atom1.bonds[atom2]
-        bond13 = atom1.bonds[atom3]
+        
+        # molecule 1
+        
+        self.assertTrue(self.molecule[0].multiplicity == 2)
+        
+        atom1 = self.molecule[0].atoms[0]
+        atom2 = self.molecule[0].atoms[3]
+        atom3 = self.molecule[0].atoms[4]
+        atom4 = self.molecule[0].atoms[5]
+        self.assertTrue(self.molecule[0].hasBond(atom2,atom1))
+        self.assertTrue(self.molecule[0].hasBond(atom2,atom3))
+        self.assertTrue(self.molecule[0].hasBond(atom2,atom4))
+        self.assertFalse(self.molecule[0].hasBond(atom1,atom3))
+        self.assertFalse(self.molecule[0].hasBond(atom1,atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+        
+        # molecule 2
+        
+        self.assertTrue(self.molecule[1].multiplicity == 2)
+        
+        atom1 = self.molecule[1].atoms[0]
+        atom2 = self.molecule[1].atoms[3]
+        atom3 = self.molecule[1].atoms[4]
+        atom4 = self.molecule[1].atoms[5]
+        self.assertTrue(self.molecule[1].hasBond(atom2,atom1))
+        self.assertTrue(self.molecule[1].hasBond(atom2,atom3))
+        self.assertTrue(self.molecule[1].hasBond(atom2,atom4))
+        self.assertFalse(self.molecule[1].hasBond(atom1,atom3))
+        self.assertFalse(self.molecule[1].hasBond(atom1,atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+        
+        # molecule 3
+        
+        self.assertTrue(self.molecule[2].multiplicity == 2)
+        
+        atom1 = self.molecule[2].atoms[0]
+        atom2 = self.molecule[2].atoms[3]
+        atom3 = self.molecule[2].atoms[4]
+        atom4 = self.molecule[2].atoms[5]
+        self.assertTrue(self.molecule[2].hasBond(atom2,atom1))
+        self.assertTrue(self.molecule[2].hasBond(atom2,atom3))
+        self.assertTrue(self.molecule[2].hasBond(atom2,atom4))
+        self.assertFalse(self.molecule[2].hasBond(atom1,atom3))
+        self.assertFalse(self.molecule[2].hasBond(atom1,atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
+           
+        self.assertTrue(atom1.label == '*1')
+        self.assertTrue(atom1.element.symbol == 'C')
+        self.assertTrue(atom1.radicalElectrons == 1)
+        self.assertTrue(atom1.charge == 0)
+        
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
+        self.assertTrue(atom2.radicalElectrons == 0)
+        self.assertTrue(atom2.charge == 1)
+        
+        self.assertTrue(atom3.label == '')
+        self.assertTrue(atom3.element.symbol == 'O')
+        self.assertTrue(atom3.radicalElectrons == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
+
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+        
+        # molecule 4
+        
+        self.assertTrue(self.molecule[3].multiplicity == 2)
+        
+        atom1 = self.molecule[3].atoms[0]
+        atom2 = self.molecule[3].atoms[1]
+        atom3 = self.molecule[3].atoms[2]
+        atom4 = self.molecule[3].atoms[3]
+        self.assertTrue(self.molecule[3].hasBond(atom2,atom1))
+        self.assertTrue(self.molecule[3].hasBond(atom2,atom3))
+        self.assertTrue(self.molecule[3].hasBond(atom2,atom4))
+        self.assertFalse(self.molecule[3].hasBond(atom1,atom3))
+        self.assertFalse(self.molecule[3].hasBond(atom1,atom4))
+        bond21 = atom2.bonds[atom1]
+        bond23 = atom2.bonds[atom3]
+        bond24 = atom2.bonds[atom4]
            
-        self.assertTrue(atom1.label == '*2')
+        self.assertTrue(atom1.label == '*1')
         self.assertTrue(atom1.element.symbol == 'C')
         self.assertTrue(atom1.radicalElectrons == 1)
-        self.assertTrue(atom1.spinMultiplicity == 2)
         self.assertTrue(atom1.charge == 0)
         
-        self.assertTrue(atom2.label == '*1')
-        self.assertTrue(atom2.element.symbol == 'O')
+        self.assertTrue(atom2.label == '*2')
+        self.assertTrue(atom2.element.symbol == 'N')
         self.assertTrue(atom2.radicalElectrons == 0)
-        self.assertTrue(atom2.spinMultiplicity == 1)
-        self.assertTrue(atom2.charge == 0)
+        self.assertTrue(atom2.charge == 1)
         
         self.assertTrue(atom3.label == '')
-        self.assertTrue(atom3.element.symbol == 'C')
+        self.assertTrue(atom3.element.symbol == 'O')
         self.assertTrue(atom3.radicalElectrons == 0)
-        self.assertTrue(atom3.spinMultiplicity == 1)
-        self.assertTrue(atom3.charge == 0)
+        self.assertTrue(atom3.charge == -1)
+        
+        self.assertTrue(atom4.label == '')
+        self.assertTrue(atom4.element.symbol == 'O')
+        self.assertTrue(atom4.radicalElectrons == 0)
+        self.assertTrue(atom4.charge == 0)
 
-        self.assertTrue(bond12.isDouble())
-        self.assertTrue(bond13.isSingle())
+        self.assertTrue(bond21.isSingle())
+        self.assertTrue(bond23.isSingle())
+        self.assertTrue(bond24.isDouble())
+        
+        
 
     def testToAdjacencyList(self):
         """
         Test the Molecule.toAdjacencyList() method.
         """
-        adjlist = self.molecule.toAdjacencyList(removeH=False)
-        self.assertEqual(adjlist.strip(), self.adjlist.strip())
+        adjlist_1 = self.molecule[0].toAdjacencyList(removeH=False)
+        newMolecule = Molecule().fromAdjacencyList(adjlist_1)
+        self.assertTrue(self.molecule[0].isIsomorphic(newMolecule))
         
-    def testFromOldAdjacencyList(self):
-        """
-        Test we can read things with implicit hydrogens.
-        """
-        adjList = "1 O 0" # should be Water
-        molecule = Molecule().fromAdjacencyList(adjList, saturateH=True) # only works with saturateH=True
-        self.assertEqual(molecule.getFormula(),'H2O')
+        #self.assertEqual(adjlist_1.strip(), self.adjlist_1.strip())
+        
+#    def testFromOldAdjacencyList(self):
+#        """
+#        Test we can read things with implicit hydrogens.
+#        """
+#        adjList = """
+#        1 O 0 
+#        """ # should be Water
+#        molecule = Molecule().fromAdjacencyList(adjList, saturateH=True) # only works with saturateH=True
+#        self.assertEqual(molecule.getFormula(),'H2O')
 
     def testIsomorphism(self):
         """
@@ -605,8 +777,8 @@ def testSubgraphIsomorphism(self):
         """
         molecule = Molecule().fromSMILES('C=CC=C[CH]C')
         group = Group().fromAdjacencyList("""
-        1 Cd 0 {2,D}
-        2 Cd 0 {1,D}
+        1 Cd u0 {2,D}
+        2 Cd u0 {1,D}
         """)
 
         self.assertTrue(molecule.isSubgraphIsomorphic(group))
@@ -621,30 +793,30 @@ def testSubgraphIsomorphism(self):
     def testSubgraphIsomorphismAgain(self):
         molecule = Molecule()
         molecule.fromAdjacencyList("""
-        1 * C 0 {2,D} {7,S} {8,S}
-        2 C 0 {1,D} {3,S} {9,S}
-        3 C 0 {2,S} {4,D} {10,S}
-        4 C 0 {3,D} {5,S} {11,S}
-        5 C 0 {4,S} {6,S} {12,S} {13,S}
-        6 C 0 {5,S} {14,S} {15,S} {16,S}
-        7 H 0 {1,S}
-        8 H 0 {1,S}
-        9 H 0 {2,S}
-        10 H 0 {3,S}
-        11 H 0 {4,S}
-        12 H 0 {5,S}
-        13 H 0 {5,S}
-        14 H 0 {6,S}
-        15 H 0 {6,S}
-        16 H 0 {6,S}
+        1 * C u0 {2,D} {7,S} {8,S}
+        2   C u0 {1,D} {3,S} {9,S}
+        3   C u0 {2,S} {4,D} {10,S}
+        4   C u0 {3,D} {5,S} {11,S}
+        5   C u0 {4,S} {6,S} {12,S} {13,S}
+        6   C u0 {5,S} {14,S} {15,S} {16,S}
+        7   H u0 {1,S}
+        8   H u0 {1,S}
+        9   H u0 {2,S}
+        10  H u0 {3,S}
+        11  H u0 {4,S}
+        12  H u0 {5,S}
+        13  H u0 {5,S}
+        14  H u0 {6,S}
+        15  H u0 {6,S}
+        16  H u0 {6,S}
         """)
 
         group = Group()
         group.fromAdjacencyList("""
-        1 * C 0 {2,D} {3,S} {4,S}
-        2   C 0 {1,D}
-        3   H 0 {1,S}
-        4   H 0 {1,S}
+        1 * C u0 {2,D} {3,S} {4,S}
+        2   C u0 {1,D}
+        3   H u0 {1,S}
+        4   H u0 {1,S}
         """)
 
         labeled1 = molecule.getLabeledAtoms().values()[0]
@@ -665,22 +837,21 @@ def testSubgraphIsomorphismAgain(self):
     def testSubgraphIsomorphismManyLabels(self):
         molecule = Molecule() # specific case (species)
         molecule.fromAdjacencyList("""
-1 *1 C  1 {2,S} {3,S} {4,S}
-2    C  0 {1,S} {3,S} {5,S} {6,S}
-3    C  0 {1,S} {2,S} {7,S} {8,S}
-4    H  0 {1,S}
-5    H  0 {2,S}
-6    H  0 {2,S}
-7    H  0 {3,S}
-8    H  0 {3,S}
+1 *1 C  u1 {2,S} {3,S} {4,S}
+2    C  u0 {1,S} {3,S} {5,S} {6,S}
+3    C  u0 {1,S} {2,S} {7,S} {8,S}
+4    H  u0 {1,S}
+5    H  u0 {2,S}
+6    H  u0 {2,S}
+7    H  u0 {3,S}
+8    H  u0 {3,S}
         """)
-        print molecule.toAdjacencyList()
 
         group = Group() # general case (functional group)
         group.fromAdjacencyList("""
-1 *1 C   1 {2,S}, {3,S}
-2    R!H 0 {1,S}
-3    R!H 0 {1,S}
+1 *1 C   u1 {2,S}, {3,S}
+2    R!H u0 {1,S}
+3    R!H u0 {1,S}
         """)
 
         labeled1 = molecule.getLabeledAtoms()
@@ -703,21 +874,21 @@ def testAdjacencyList(self):
         Check the adjacency list read/write functions for a full molecule.
         """
         molecule1 = Molecule().fromAdjacencyList("""
-        1 C 0       {2,D} {7,S} {8,S}
-        2 C 0 {1,D} {3,S} {9,S}
-        3 C 0 {2,S} {4,D} {10,S}
-        4 C 0 {3,D} {5,S} {11,S}
-        5 C 1 {4,S} {6,S} {12,S}
-        6 C 0 {5,S} {13,S} {14,S} {15,S}
-        7 H 0 {1,S}
-        8 H 0 {1,S}
-        9 H 0 {2,S}
-        10 H 0 {3,S}
-        11 H 0 {4,S}
-        12 H 0 {5,S}
-        13 H 0 {6,S}
-        14 H 0 {6,S}
-        15 H 0 {6,S}
+        1  C u0 {2,D} {7,S} {8,S}
+        2  C u0 {1,D} {3,S} {9,S}
+        3  C u0 {2,S} {4,D} {10,S}
+        4  C u0 {3,D} {5,S} {11,S}
+        5  C u1 {4,S} {6,S} {12,S}
+        6  C u0 {5,S} {13,S} {14,S} {15,S}
+        7  H u0 {1,S}
+        8  H u0 {1,S}
+        9  H u0 {2,S}
+        10 H u0 {3,S}
+        11 H u0 {4,S}
+        12 H u0 {5,S}
+        13 H u0 {6,S}
+        14 H u0 {6,S}
+        15 H u0 {6,S}
         """)
         molecule2 = Molecule().fromSMILES('C=CC=C[CH]C')
         self.assertTrue(molecule1.isIsomorphic(molecule2))
@@ -806,7 +977,7 @@ def testRadicalCH(self):
         """
         molecule = Molecule().fromSMILES('[CH]')
         self.assertEqual(molecule.atoms[0].radicalElectrons, 3)
-        self.assertEqual(molecule.atoms[0].spinMultiplicity, 4)
+        self.assertEqual(molecule.multiplicity, 4)
         self.assertEqual(molecule.getRadicalCount(), 3)
 
     def testRadicalCH2(self):
@@ -815,7 +986,7 @@ def testRadicalCH2(self):
         """
         molecule = Molecule().fromSMILES('[CH2]')
         self.assertEqual(molecule.atoms[0].radicalElectrons, 2)
-        self.assertEqual(molecule.atoms[0].spinMultiplicity, 3)
+        self.assertEqual(molecule.multiplicity, 3)
         self.assertEqual(molecule.getRadicalCount(), 2)
         
     def testRadicalCH2CH2CH2(self):
@@ -830,7 +1001,7 @@ def testSMILES(self):
         Test that we can generate a few SMILES strings as expected
         """
         import rmgpy.molecule
-        test_strings =['CO', '[C]', '[CH]', 'OO', '[H][H]', '[H]', '[He]', '[O]', 'O', '[CH3]', 'C', '[OH]', 'CCC', 'CC', 'N#N', '[O]O', 'C[CH2]', '[Ar]', 'CCCC','O=C=O','CN']
+        test_strings =['[C+]#[O-]', '[C]', '[CH]', 'OO', '[H][H]', '[H]', '[He]', '[O]', 'O', '[CH3]', 'C', '[OH]', 'CCC', 'CC', 'N#N', '[O]O', 'C[CH2]', '[Ar]', 'CCCC','O=C=O','[C]#N']
         for s in test_strings:
             molecule = Molecule(SMILES=s)
             self.assertEqual(s,molecule.toSMILES())
@@ -848,8 +1019,8 @@ def testAugmentedInChI(self):
         Test that Augmented InChI generation is printing the /mult layer
         """
         mol = Molecule().fromAdjacencyList("""
-            1     C     1 {2,S}
-            2     C     1 {1,S}
+            1     C     u1 {2,S}
+            2     C     u1 {1,S}
         """)
         
         self.assertEqual(mol.toAugmentedInChI(), 'InChI=1S/C2H4/c1-2/h1-2H2/mult3')
@@ -859,8 +1030,8 @@ def testAugmentedInChIKey(self):
         Test that Augmented InChI Key generation is printing the mult layer
         """
         mol = Molecule().fromAdjacencyList("""
-            1     C     1 {2,S}
-            2     C     1 {1,S}
+            1     C     u1 {2,S}
+            2     C     u1 {1,S}
         """)
         
         self.assertEqual(mol.toAugmentedInChIKey(), 'VGGSQFUCUMXWEO-UHFFFAOYSAmult3')
diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py
index 4af083e9dc..72487d789a 100644
--- a/rmgpy/qm/molecule.py
+++ b/rmgpy/qm/molecule.py
@@ -32,18 +32,15 @@ class Geometry:
     `settings`          :class:`QMSettings`     Settings for QM calculations
     `uniqueID`          ``str``                 A short ID such as an augmented InChI Key
     `molecule`          :class:`Molecule`       RMG Molecule object
-    `multiplicity`      ``int``                 The multiplicity of the molecule
     `uniqueIDlong`      ``str``                 A long, truly unique ID such as an augmented InChI
     =================== ======================= ====================================
     
     """
-    def __init__(self, settings, uniqueID, molecule, multiplicity, uniqueIDlong=None):
+    def __init__(self, settings, uniqueID, molecule, uniqueIDlong=None):
         self.settings = settings
         #: A short unique ID such as an augmented InChI Key.
         self.uniqueID = uniqueID
         self.molecule = molecule
-        #: The multiplicity, eg. the number of free radicals plus one.
-        self.multiplicity = multiplicity
         if uniqueIDlong is None:
             self.uniqueIDlong = uniqueID
         else:
@@ -236,8 +233,7 @@ def createGeometry(self):
         """
         Creates self.geometry with RDKit geometries
         """
-        multiplicity = sum([i.radicalElectrons for i in self.molecule.atoms]) + 1
-        self.geometry = Geometry(self.settings, self.uniqueID, self.molecule, multiplicity, uniqueIDlong=self.uniqueIDlong)
+        self.geometry = Geometry(self.settings, self.uniqueID, self.molecule, uniqueIDlong=self.uniqueIDlong)
         self.geometry.generateRDKitGeometries()
         return self.geometry
     
diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py
index e26f1e74bd..ce359859b4 100644
--- a/rmgpy/qm/mopac.py
+++ b/rmgpy/qm/mopac.py
@@ -256,7 +256,7 @@ def inputFileKeywords(self, attempt):
         if attempt > self.scriptAttempts:
             attempt -= self.scriptAttempts
         
-        multiplicity_keys = self.multiplicityKeywords[self.geometry.multiplicity]
+        multiplicity_keys = self.multiplicityKeywords[self.geometry.molecule.multiplicity]
 
         top_keys = "{method} {mult} {top}".format(
                 method = self.pm_method,
@@ -270,7 +270,7 @@ def inputFileKeywords(self, attempt):
                 )
         polar_keys = "oldgeo {polar} nosym precise {method} {mult}".format(
                 method = self.pm_method,
-                polar = ('polar' if self.geometry.multiplicity == 1 else 'static'),
+                polar = ('polar' if self.geometry.molecule.multiplicity == 1 else 'static'),
                 mult = multiplicity_keys,
                 )
 
diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd
index 8eeb90797f..230852e40a 100644
--- a/rmgpy/reaction.pxd
+++ b/rmgpy/reaction.pxd
@@ -100,6 +100,8 @@ cdef class Reaction:
     
     cpdef generatePairs(self)
     
+    cpdef copy(self)
+    
 ################################################################################
 
 cdef class ReactionModel:
diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py
index 83f81fceb2..bad5b2ee1c 100755
--- a/rmgpy/reaction.py
+++ b/rmgpy/reaction.py
@@ -45,6 +45,7 @@
 import logging
 import re
 import os.path
+from copy import copy, deepcopy
 
 import rmgpy.constants as constants
 from rmgpy.molecule.molecule import Molecule, Atom
@@ -980,6 +981,31 @@ def generate3dTS(self, reactants, products):
                         zCoord[k] = dirVec[k].z*lenVec[k]
             reactionAxis = [sum(xCoord), sum(yCoord), sum(zCoord)]
             products[i].reactionAxis = reactionAxis
+            
+    def copy(self):
+        """
+        Create a deep copy of the current reaction.
+        """
+        
+        cython.declare(other=Reaction)
+
+        other = Reaction.__new__(Reaction)
+        other.index = self.index
+        other.label = self.label
+        other.reactants = []
+        for reactant in self.reactants:
+            other.reactants.append(reactant.copy(deep=True))
+        other.products = []
+        for product in self.products:
+            other.products.append(product.copy(deep=True))
+        other.kinetics = deepcopy(self.kinetics)
+        other.reversible = self.reversible
+        other.transitionState = deepcopy(self.transitionState)
+        other.duplicate = self.duplicate
+        other.degeneracy = self.degeneracy
+        other.pairs = deepcopy(self.pairs)
+        
+        return other
 
                 
 ################################################################################
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 1e252a604b..5d9ceca3cf 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -61,6 +61,9 @@
 from pdep import PDepReaction, PDepNetwork, PressureDependenceError
 # generateThermoDataFromQM under the Species class imports the qm package
 
+#: This dictionary is used to add multiplicity to species label
+_multiplicity_labels = {1:'S',2:'D',3:'T',4:'Q',5:'V',}
+
 
 ################################################################################
 
@@ -126,7 +129,6 @@ def generateThermoData(self, database, thermoClass=NASA, quantumMechanics=None):
                     thermo = []
                     for molecule in self.molecule:
                         molecule.clearLabeledAtoms()
-                        molecule.updateAtomTypes()
                         tdata = database.thermo.estimateRadicalThermoViaHBI(molecule, quantumMechanics.getThermoData)
                         if tdata is not None:
                             thermo.append(tdata)
@@ -151,7 +153,7 @@ def generateThermoData(self, database, thermoClass=NASA, quantumMechanics=None):
                 if thermo0 is not None:
                     # Write the QM molecule thermo to a library so that can be used in future RMG jobs.
                     quantumMechanics.database.loadEntry(index = len(quantumMechanics.database.entries) + 1,
-                                                        label = molecule.toSMILES(),
+                                                        label = molecule.toSMILES() + '_({0})'.format(_multiplicity_labels[molecule.multiplicity]),
                                                         molecule = molecule.toAdjacencyList(),
                                                         thermo = thermo0,
                                                         shortDesc = thermo0.comment
diff --git a/rmgpy/species.py b/rmgpy/species.py
index 4fdf54da93..82c29eaa1a 100755
--- a/rmgpy/species.py
+++ b/rmgpy/species.py
@@ -49,6 +49,9 @@
 import rmgpy.quantity as quantity
 from rmgpy.molecule import Molecule
 
+#: This dictionary is used to add multiplicity to species label
+_multiplicity_labels = {1:'S',2:'D',3:'T',4:'Q',5:'V',}
+                           
 ################################################################################
 
 class SpeciesError(Exception):
@@ -104,6 +107,14 @@ def __init__(self, index=-1, label='', thermo=None, conformer=None,
         self.polarizability = polarizability
         self.Zrot = Zrot
         self.energyTransferModel = energyTransferModel
+        
+        # Check multiplicity of each molecule is the same
+        if molecule is not None and len(molecule)>1:
+            mult = molecule[0].multiplicity
+            for m in molecule[1:]:
+                if mult != m.multiplicity:
+                    raise SpeciesError('Multiplicities of molecules in species {species} do not match.'.format(species=label))
+
 
     def __repr__(self):
         """
@@ -189,10 +200,10 @@ def isIsomorphic(self, other):
                 if molecule.isIsomorphic(other):
                     return True
         elif isinstance(other, Species):
-            for molecule1 in self.molecule:
-                for molecule2 in other.molecule:
-                    if molecule1.isIsomorphic(molecule2):
-                        return True
+                for molecule1 in self.molecule:
+                    for molecule2 in other.molecule:
+                        if molecule1.isIsomorphic(molecule2):
+                            return True
         else:
             raise ValueError('Unexpected value "{0!r}" for other parameter; should be a Molecule or Species object.'.format(other))
         return False
diff --git a/rmgpy/speciesTest.py b/rmgpy/speciesTest.py
index 5ea45667f5..3961a71a1e 100644
--- a/rmgpy/speciesTest.py
+++ b/rmgpy/speciesTest.py
@@ -61,6 +61,7 @@ def testPickle(self):
         species = cPickle.loads(cPickle.dumps(self.species,-1))
         self.assertEqual(self.species.index, species.index)
         self.assertEqual(self.species.label, species.label)
+        self.assertEqual(self.species.molecule[0].multiplicity, species.molecule[0].multiplicity)
         self.assertEqual(self.species.thermo.H298.value_si, species.thermo.H298.value_si)
         self.assertEqual(self.species.thermo.H298.units, species.thermo.H298.units)
         self.assertEqual(len(self.species.conformer.modes), len(species.conformer.modes))
@@ -86,6 +87,7 @@ def testOutput(self):
         exec('species = {0!r}'.format(self.species))
         self.assertEqual(self.species.index, species.index)
         self.assertEqual(self.species.label, species.label)
+        self.assertEqual(self.species.molecule[0].multiplicity, species.molecule[0].multiplicity)
         self.assertEqual(self.species.thermo.H298.value_si, species.thermo.H298.value_si)
         self.assertEqual(self.species.thermo.H298.units, species.thermo.H298.units)
         self.assertEqual(len(self.species.conformer.modes), len(species.conformer.modes))
diff --git a/rmgpy/transportDataTest.py b/rmgpy/transportDataTest.py
index 12981e4c2e..e979a9659e 100644
--- a/rmgpy/transportDataTest.py
+++ b/rmgpy/transportDataTest.py
@@ -252,7 +252,8 @@ def testJoback(self):
             ['benzene', 'c1ccccc1', None, None, None],
             ]
         for name, smiles, sigma, epsilon, comment in self.testCases:
-            species = Species(molecule=[Molecule(SMILES=smiles)])
+            molecule=Molecule(SMILES=smiles)
+            species = Species(molecule=[molecule])
             transportData, blank, blank2 = self.transportdb.getTransportPropertiesViaGroupEstimates(species)
             # check Joback worked.
             # If we don't know what to expect, don't check (just make sure we didn't crash)
@@ -266,18 +267,18 @@ def testJoback(self):
     def testJobackOnBenzeneBonds(self):
         "Test Joback doesn't crash on Cb desription of beneze"
         adjlist = """
-                    1  C 0 0 {2,D} {6,S} {7,S}
-                    2  C 0 0 {1,D} {3,S} {8,S}
-                    3  C 0 0 {2,S} {4,D} {9,S}
-                    4  C 0 0 {3,D} {5,S} {10,S}
-                    5  C 0 0 {4,S} {6,D} {11,S}
-                    6  C 0 0 {1,S} {5,D} {12,S}
-                    7  H 0 0 {1,S}
-                    8  H 0 0 {2,S}
-                    9  H 0 0 {3,S}
-                    10 H 0 0 {4,S}
-                    11 H 0 0 {5,S}
-                    12 H 0 0 {6,S}
+                    1  C u0 p0 {2,D} {6,S} {7,S}
+                    2  C u0 p0 {1,D} {3,S} {8,S}
+                    3  C u0 p0 {2,S} {4,D} {9,S}
+                    4  C u0 p0 {3,D} {5,S} {10,S}
+                    5  C u0 p0 {4,S} {6,D} {11,S}
+                    6  C u0 p0 {1,S} {5,D} {12,S}
+                    7  H u0 p0 {1,S}
+                    8  H u0 p0 {2,S}
+                    9  H u0 p0 {3,S}
+                    10 H u0 p0 {4,S}
+                    11 H u0 p0 {5,S}
+                    12 H u0 p0 {6,S}
                     """
         m = Molecule().fromAdjacencyList(adjlist)
         species = Species(molecule=[m])