diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index 129ab5a5a2d8cf191f44b78550d017e2dfb7853f..11f0dc706bc420707672c7e6f1995db41935a376 100644
--- a/pycgtool/bondset.py
+++ b/pycgtool/bondset.py
@@ -226,26 +226,29 @@ class BondSet:
                     e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, mol),)
                     raise
 
-    def write_itp(self, filename, mapping):
+    def write_itp(self, filename, mapping, exclude=set()):
         """
         Output a GROMACS .itp file containing atoms/beads and bonded terms.
 
         :param filename: Name of output file
         :param mapping: AA->CG Mapping from which to collect bead properties
+        :param exclude: Set of molecule names to be excluded from itp
         """
         self._populate_atom_numbers(mapping)
         backup_file(filename, verbose=True)
 
-        def write_bond_angle_dih(bonds, section_header, itp, fconst=True):
+        def write_bond_angle_dih(bonds, section_header, itp, fconst=True, multiplicity=None):
             if bonds:
                 print("\n[ {0:s} ]".format(section_header), file=itp)
             for bond in bonds:
                 # Factor is usually 1, unless doing correction
                 eqm = bond.eqm * self._empirical_correction_factor
                 line = " ".join(["{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers])
-                line += " {0:12.5f}".format(eqm)
+                line += " {0:4d} {1:12.5f}".format(1, eqm)
                 if fconst:
                     line += " {0:12.5f}".format(bond.fconst)
+                if multiplicity is not None:
+                    line += " {0:4d}".format(multiplicity)
                 print(line, file=itp)
 
         with open(filename, "w") as itp:
@@ -258,7 +261,7 @@ class BondSet:
             print(header, file=itp)
 
             # Print molecule
-            for mol in self._molecules:
+            for mol in filter(lambda mol: mol not in exclude, self._molecules):
                 print("\n[ moleculetype ]", file=itp)
                 print("{0:4s} {1:4d}".format(mol, 1), file=itp)
 
@@ -271,7 +274,7 @@ class BondSet:
 
                 write_bond_angle_dih(self.get_bond_lengths(mol), "bonds", itp)
                 write_bond_angle_dih(self.get_bond_angles(mol), "angles", itp)
-                write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", itp)
+                write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", itp, multiplicity=1)
                 write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", itp, fconst=False)
 
     def apply(self, frame):
diff --git a/test/data/sugar_out.itp b/test/data/sugar_out.itp
new file mode 100644
index 0000000000000000000000000000000000000000..1fc051bbd25ecb080da05cba82d823b88d49b704
--- /dev/null
+++ b/test/data/sugar_out.itp
@@ -0,0 +1,43 @@
+; 
+; Topology prepared automatically using PyCGTOOL 
+; James Graham <J.A.Graham@soton.ac.uk> 2016 
+; University of Southampton 
+; https://github.com/jag1g13/pycgtool 
+;
+
+[ moleculetype ]
+ALLA    1
+
+[ atoms ]
+   1 P3      1 ALLA C1      1    0.000
+   2 P3      1 ALLA C2      2    0.000
+   3 P3      1 ALLA C3      3    0.000
+   4 P3      1 ALLA C4      4    0.000
+   5 P2      1 ALLA C5      5    0.000
+   6 P4      1 ALLA O5      6    0.000
+
+[ bonds ]
+   2    3    1      0.21998  35162.10022
+   3    4    1      0.21648  65545.10313
+   4    5    1      0.26454  19022.17121
+   5    6    1      0.20452  53457.20659
+
+[ angles ]
+   1    2    3    1     77.69557    383.16060
+   2    3    4    1    115.59599    262.66648
+   3    4    5    1    108.51038    698.56915
+   4    5    6    1     81.23827   1098.41706
+   5    6    1    1    146.49898     48.08599
+   6    1    2    1     99.84045   1571.73513
+
+[ dihedrals ]
+   1    2    3    4    1    -85.45004    234.10405    1
+   2    3    4    5    1     65.21378    117.46132    1
+   3    4    5    6    1    -22.73796     94.20175    1
+   4    5    6    1    1     55.18087     37.92107    1
+   5    6    1    2    1    -96.52349     40.73340    1
+   6    1    2    3    1     73.82477    259.20126    1
+
+[ constraints ]
+   1    2    1      0.22209
+   6    1    1      0.17931
diff --git a/test/test_bondset.py b/test/test_bondset.py
index faf529ef398fbdf779a25ec372f46c75574f6c39..5028dc810c13f64c326ca95ec9733c0ec052cc97 100644
--- a/test/test_bondset.py
+++ b/test/test_bondset.py
@@ -1,4 +1,5 @@
 import unittest
+import filecmp
 
 from pycgtool.bondset import BondSet
 from pycgtool.frame import Frame
@@ -103,3 +104,18 @@ class BondSetTest(unittest.TestCase):
         for bond in bondset.get_bond_lengths("ETH", True):
             self.assertAlmostEqual(1., bond.eqm)
             self.assertEqual(float("inf"), bond.fconst)
+
+    def test_full_itp_sugar(self):
+        measure = BondSet("test/data/sugar.bnd", DummyOptions)
+        frame = Frame("test/data/sugar.gro", xtc="test/data/sugar.xtc")
+        mapping = Mapping("test/data/sugar.map", DummyOptions)
+        cgframe = mapping.apply(frame)
+
+        while frame.next_frame():
+            cgframe = mapping.apply(frame, cgframe=cgframe, exclude={"SOL"})
+            measure.apply(cgframe)
+
+        measure.boltzmann_invert()
+        measure.write_itp("sugar_out.itp", mapping, exclude={"SOL"})
+
+        self.assertTrue(filecmp.cmp("sugar_out.itp", "test/data/sugar_out.itp"))