diff --git a/.benchmarks/Darwin-CPython-3.12-64bit/0001_latest.json b/.benchmarks/Darwin-CPython-3.12-64bit/0001_latest.json
index 79f05f6..c94dedd 100644
--- a/.benchmarks/Darwin-CPython-3.12-64bit/0001_latest.json
+++ b/.benchmarks/Darwin-CPython-3.12-64bit/0001_latest.json
@@ -11,7 +11,7 @@
             "main",
             "Apr 10 2025 22:19:24"
         ],
-        "release": "25.1.0",
+        "release": "25.4.0",
         "system": "Darwin",
         "cpu": {
             "python_version": "3.12.10.final.0 (64 bit)",
@@ -29,11 +29,11 @@
         }
     },
     "commit_info": {
-        "id": "659c1303a77acd9517a592564d2b16e15534ff26",
-        "time": "2025-11-30T15:37:16-05:00",
-        "author_time": "2025-11-30T15:37:16-05:00",
+        "id": "878eec2b40e5bb093bd8d1e091a728b8a8b25aea",
+        "time": "2026-05-20T14:05:03-04:00",
+        "author_time": "2026-05-20T14:05:03-04:00",
         "dirty": true,
-        "project": "ChemPy",
+        "project": "chempy",
         "branch": "master"
     },
     "benchmarks": [
@@ -53,131 +53,26 @@
                 "warmup": false
             },
             "stats": {
-                "min": 0.0002426248975098133,
-                "max": 0.00026537501253187656,
-                "mean": 0.000249016797170043,
-                "stddev": 9.422936996810157e-06,
+                "min": 0.0003945000935345888,
+                "max": 0.00044408394023776054,
+                "mean": 0.00041614188812673093,
+                "stddev": 2.196570902992079e-05,
                 "rounds": 5,
-                "median": 0.0002448339946568012,
-                "iqr": 9.437382686883211e-06,
-                "q1": 0.00024337507784366608,
-                "q3": 0.0002528124605305493,
+                "median": 0.00040483311749994755,
+                "iqr": 3.705290146172047e-05,
+                "q1": 0.00040028116200119257,
+                "q3": 0.00043733406346291304,
                 "iqr_outliers": 0,
                 "stddev_outliers": 1,
                 "outliers": "1;0",
-                "ld15iqr": 0.0002426248975098133,
-                "hd15iqr": 0.00026537501253187656,
-                "ops": 4015.79335757476,
-                "total": 0.001245083985850215,
+                "ld15iqr": 0.0003945000935345888,
+                "hd15iqr": 0.00044408394023776054,
+                "ops": 2403.0265362170467,
+                "total": 0.0020807094406336546,
                 "iterations": 1
             }
-        },
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_ethane_rotors",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_ethane_rotors",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.7457870095968246e-05,
-                "max": 0.0021958330180495977,
-                "mean": 5.593497004731008e-05,
-                "stddev": 4.265468175877866e-05,
-                "rounds": 11817,
-                "median": 4.970794543623924e-05,
-                "iqr": 2.3748725652694702e-06,
-                "q1": 4.9042049795389175e-05,
-                "q3": 5.1416922360658646e-05,
-                "iqr_outliers": 1759,
-                "stddev_outliers": 502,
-                "outliers": "502;1759",
-                "ld15iqr": 4.7457870095968246e-05,
-                "hd15iqr": 5.4999953135848045e-05,
-                "ops": 17877.903557545396,
-                "total": 0.6609835410490632,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_density_of_states_ilt",
-            "fullname": "unittest/benchmarksTest.py::test_bench_density_of_states_ilt",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.03818258293904364,
-                "max": 0.03883983311243355,
-                "mean": 0.03844411661848426,
-                "stddev": 0.0002482778623857746,
-                "rounds": 5,
-                "median": 0.038354666903615,
-                "iqr": 0.00028325035236775875,
-                "q1": 0.03830248938174918,
-                "q3": 0.03858573973411694,
-                "iqr_outliers": 0,
-                "stddev_outliers": 2,
-                "outliers": "2;0",
-                "ld15iqr": 0.03818258293904364,
-                "hd15iqr": 0.03883983311243355,
-                "ops": 26.01178250300051,
-                "total": 0.1922205830924213,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_states_construction",
-            "fullname": "unittest/benchmarksTest.py::test_bench_states_construction",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.978966899216175e-07,
-                "max": 4.279147833585739e-06,
-                "mean": 5.243171563162537e-07,
-                "stddev": 4.060613257512666e-08,
-                "rounds": 80542,
-                "median": 5.146022886037826e-07,
-                "iqr": 1.2503005564212757e-08,
-                "q1": 5.103996954858303e-07,
-                "q3": 5.229027010500431e-07,
-                "iqr_outliers": 8745,
-                "stddev_outliers": 4937,
-                "outliers": "4937;8745",
-                "ld15iqr": 4.978966899216175e-07,
-                "hd15iqr": 5.416572093963623e-07,
-                "ops": 1907242.5686502378,
-                "total": 0.04222955240402371,
-                "iterations": 20
-            }
         }
     ],
-    "datetime": "2025-11-30T20:52:28.842699+00:00",
+    "datetime": "2026-05-20T18:07:09.210269+00:00",
     "version": "5.2.3"
-}
+}
\ No newline at end of file
diff --git a/.benchmarks/Darwin-CPython-3.12-64bit/0002_latest.json b/.benchmarks/Darwin-CPython-3.12-64bit/0002_latest.json
deleted file mode 100644
index fe612ff..0000000
--- a/.benchmarks/Darwin-CPython-3.12-64bit/0002_latest.json
+++ /dev/null
@@ -1,183 +0,0 @@
-{
-    "machine_info": {
-        "node": "Georges-Mini",
-        "processor": "arm",
-        "machine": "arm64",
-        "python_compiler": "Clang 18.1.8 ",
-        "python_implementation": "CPython",
-        "python_implementation_version": "3.12.10",
-        "python_version": "3.12.10",
-        "python_build": [
-            "main",
-            "Apr 10 2025 22:19:24"
-        ],
-        "release": "25.1.0",
-        "system": "Darwin",
-        "cpu": {
-            "python_version": "3.12.10.final.0 (64 bit)",
-            "cpuinfo_version": [
-                9,
-                0,
-                0
-            ],
-            "cpuinfo_version_string": "9.0.0",
-            "arch": "ARM_8",
-            "bits": 64,
-            "count": 10,
-            "arch_string_raw": "arm64",
-            "brand_raw": "Apple M4"
-        }
-    },
-    "commit_info": {
-        "id": "659c1303a77acd9517a592564d2b16e15534ff26",
-        "time": "2025-11-30T15:37:16-05:00",
-        "author_time": "2025-11-30T15:37:16-05:00",
-        "dirty": true,
-        "project": "ChemPy",
-        "branch": "master"
-    },
-    "benchmarks": [
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_benzene",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_benzene",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.0002381661906838417,
-                "max": 0.00027883402071893215,
-                "mean": 0.000252208299934864,
-                "stddev": 1.6047700560055338e-05,
-                "rounds": 5,
-                "median": 0.0002476251684129238,
-                "iqr": 1.9103987142443657e-05,
-                "q1": 0.00024122907780110836,
-                "q3": 0.000260333064943552,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.0002381661906838417,
-                "hd15iqr": 0.00027883402071893215,
-                "ops": 3964.9765699949708,
-                "total": 0.0012610414996743202,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_ethane_rotors",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_ethane_rotors",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.729209467768669e-05,
-                "max": 0.0001233331859111786,
-                "mean": 4.990232637113832e-05,
-                "stddev": 3.2006128975859647e-06,
-                "rounds": 13491,
-                "median": 4.924996756017208e-05,
-                "iqr": 1.0421499609947205e-06,
-                "q1": 4.870793782174587e-05,
-                "q3": 4.975008778274059e-05,
-                "iqr_outliers": 1434,
-                "stddev_outliers": 930,
-                "outliers": "930;1434",
-                "ld15iqr": 4.729209467768669e-05,
-                "hd15iqr": 5.13328704982996e-05,
-                "ops": 20039.14592202987,
-                "total": 0.673232285073027,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_density_of_states_ilt",
-            "fullname": "unittest/benchmarksTest.py::test_bench_density_of_states_ilt",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.03864358412101865,
-                "max": 0.0394084588624537,
-                "mean": 0.039007241977378725,
-                "stddev": 0.0003365920352392336,
-                "rounds": 5,
-                "median": 0.03888758295215666,
-                "iqr": 0.0005912806373089552,
-                "q1": 0.0387470840360038,
-                "q3": 0.03933836467331275,
-                "iqr_outliers": 0,
-                "stddev_outliers": 2,
-                "outliers": "2;0",
-                "ld15iqr": 0.03864358412101865,
-                "hd15iqr": 0.0394084588624537,
-                "ops": 25.636265198650165,
-                "total": 0.19503620988689363,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_states_construction",
-            "fullname": "unittest/benchmarksTest.py::test_bench_states_construction",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.979432560503483e-07,
-                "max": 2.8499995823949576e-05,
-                "mean": 5.381731284998154e-07,
-                "stddev": 1.2095207186923548e-07,
-                "rounds": 82474,
-                "median": 5.270587280392646e-07,
-                "iqr": 1.040752977132793e-08,
-                "q1": 5.228910595178605e-07,
-                "q3": 5.332985892891884e-07,
-                "iqr_outliers": 10686,
-                "stddev_outliers": 1611,
-                "outliers": "1611;10686",
-                "ld15iqr": 5.082925781607628e-07,
-                "hd15iqr": 5.499925464391708e-07,
-                "ops": 1858138.1102909208,
-                "total": 0.04438529059989378,
-                "iterations": 20
-            }
-        }
-    ],
-    "datetime": "2025-11-30T20:53:42.147668+00:00",
-    "version": "5.2.3"
-}
diff --git a/.benchmarks/Darwin-CPython-3.12-64bit/0003_latest.json b/.benchmarks/Darwin-CPython-3.12-64bit/0003_latest.json
deleted file mode 100644
index 2441ad6..0000000
--- a/.benchmarks/Darwin-CPython-3.12-64bit/0003_latest.json
+++ /dev/null
@@ -1,183 +0,0 @@
-{
-    "machine_info": {
-        "node": "Georges-Mini",
-        "processor": "arm",
-        "machine": "arm64",
-        "python_compiler": "Clang 18.1.8 ",
-        "python_implementation": "CPython",
-        "python_implementation_version": "3.12.10",
-        "python_version": "3.12.10",
-        "python_build": [
-            "main",
-            "Apr 10 2025 22:19:24"
-        ],
-        "release": "25.1.0",
-        "system": "Darwin",
-        "cpu": {
-            "python_version": "3.12.10.final.0 (64 bit)",
-            "cpuinfo_version": [
-                9,
-                0,
-                0
-            ],
-            "cpuinfo_version_string": "9.0.0",
-            "arch": "ARM_8",
-            "bits": 64,
-            "count": 10,
-            "arch_string_raw": "arm64",
-            "brand_raw": "Apple M4"
-        }
-    },
-    "commit_info": {
-        "id": "659c1303a77acd9517a592564d2b16e15534ff26",
-        "time": "2025-11-30T15:37:16-05:00",
-        "author_time": "2025-11-30T15:37:16-05:00",
-        "dirty": true,
-        "project": "ChemPy",
-        "branch": "master"
-    },
-    "benchmarks": [
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_benzene",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_benzene",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.00022962503135204315,
-                "max": 0.0002653328701853752,
-                "mean": 0.00024175834842026233,
-                "stddev": 1.4386789822774536e-05,
-                "rounds": 5,
-                "median": 0.00024091685190796852,
-                "iqr": 1.767714275047183e-05,
-                "q1": 0.00023037503706291318,
-                "q3": 0.000248052179813385,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.00022962503135204315,
-                "hd15iqr": 0.0002653328701853752,
-                "ops": 4136.36181143016,
-                "total": 0.0012087917421013117,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_ethane_rotors",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_ethane_rotors",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.854216240346432e-05,
-                "max": 0.000122792087495327,
-                "mean": 5.1181247923011944e-05,
-                "stddev": 3.397949587214006e-06,
-                "rounds": 11841,
-                "median": 5.0416914746165276e-05,
-                "iqr": 1.0421499609947205e-06,
-                "q1": 4.9916794523596764e-05,
-                "q3": 5.0958944484591484e-05,
-                "iqr_outliers": 1296,
-                "stddev_outliers": 865,
-                "outliers": "865;1296",
-                "ld15iqr": 4.854216240346432e-05,
-                "hd15iqr": 5.2540795877575874e-05,
-                "ops": 19538.40597056609,
-                "total": 0.6060371566563845,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_density_of_states_ilt",
-            "fullname": "unittest/benchmarksTest.py::test_bench_density_of_states_ilt",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.03916516713798046,
-                "max": 0.04022625018842518,
-                "mean": 0.03967965845949948,
-                "stddev": 0.0004099150461965831,
-                "rounds": 5,
-                "median": 0.03960829204879701,
-                "iqr": 0.00060533348005265,
-                "q1": 0.039395729021634907,
-                "q3": 0.040001062501687557,
-                "iqr_outliers": 0,
-                "stddev_outliers": 2,
-                "outliers": "2;0",
-                "ld15iqr": 0.03916516713798046,
-                "hd15iqr": 0.04022625018842518,
-                "ops": 25.201829824737207,
-                "total": 0.1983982922974974,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_states_construction",
-            "fullname": "unittest/benchmarksTest.py::test_bench_states_construction",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.958012141287326e-07,
-                "max": 3.7853955291211604e-06,
-                "mean": 5.251837434976344e-07,
-                "stddev": 4.427463320219111e-08,
-                "rounds": 82191,
-                "median": 5.166977643966674e-07,
-                "iqr": 1.4598481357097583e-08,
-                "q1": 5.103996954858303e-07,
-                "q3": 5.249981768429279e-07,
-                "iqr_outliers": 7921,
-                "stddev_outliers": 3390,
-                "outliers": "3390;7921",
-                "ld15iqr": 4.958012141287326e-07,
-                "hd15iqr": 5.47897070646286e-07,
-                "ops": 1904095.4949217774,
-                "total": 0.04316537706181407,
-                "iterations": 20
-            }
-        }
-    ],
-    "datetime": "2025-11-30T20:59:14.332285+00:00",
-    "version": "5.2.3"
-}
diff --git a/.benchmarks/Darwin-CPython-3.12-64bit/0004_latest.json b/.benchmarks/Darwin-CPython-3.12-64bit/0004_latest.json
deleted file mode 100644
index ee45745..0000000
--- a/.benchmarks/Darwin-CPython-3.12-64bit/0004_latest.json
+++ /dev/null
@@ -1,183 +0,0 @@
-{
-    "machine_info": {
-        "node": "Georges-Mini",
-        "processor": "arm",
-        "machine": "arm64",
-        "python_compiler": "Clang 18.1.8 ",
-        "python_implementation": "CPython",
-        "python_implementation_version": "3.12.10",
-        "python_version": "3.12.10",
-        "python_build": [
-            "main",
-            "Apr 10 2025 22:19:24"
-        ],
-        "release": "25.1.0",
-        "system": "Darwin",
-        "cpu": {
-            "python_version": "3.12.10.final.0 (64 bit)",
-            "cpuinfo_version": [
-                9,
-                0,
-                0
-            ],
-            "cpuinfo_version_string": "9.0.0",
-            "arch": "ARM_8",
-            "bits": 64,
-            "count": 10,
-            "arch_string_raw": "arm64",
-            "brand_raw": "Apple M4"
-        }
-    },
-    "commit_info": {
-        "id": "659c1303a77acd9517a592564d2b16e15534ff26",
-        "time": "2025-11-30T15:37:16-05:00",
-        "author_time": "2025-11-30T15:37:16-05:00",
-        "dirty": true,
-        "project": "ChemPy",
-        "branch": "master"
-    },
-    "benchmarks": [
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_benzene",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_benzene",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.00023558316752314568,
-                "max": 0.0002960828132927418,
-                "mean": 0.00025211661122739317,
-                "stddev": 2.5213789042751517e-05,
-                "rounds": 5,
-                "median": 0.00024266703985631466,
-                "iqr": 2.4936278350651264e-05,
-                "q1": 0.00023633387172594666,
-                "q3": 0.00026127015007659793,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.00023558316752314568,
-                "hd15iqr": 0.0002960828132927418,
-                "ops": 3966.4185359768444,
-                "total": 0.0012605830561369658,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_ethane_rotors",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_ethane_rotors",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.8625050112605095e-05,
-                "max": 0.0026821668725460768,
-                "mean": 5.8376059826447295e-05,
-                "stddev": 5.2342820895647595e-05,
-                "rounds": 12904,
-                "median": 5.1083043217658997e-05,
-                "iqr": 2.874992787837982e-06,
-                "q1": 5.02921175211668e-05,
-                "q3": 5.3167110309004784e-05,
-                "iqr_outliers": 2071,
-                "stddev_outliers": 481,
-                "outliers": "481;2071",
-                "ld15iqr": 4.8625050112605095e-05,
-                "hd15iqr": 5.7499855756759644e-05,
-                "ops": 17130.309975921835,
-                "total": 0.7532846760004759,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_density_of_states_ilt",
-            "fullname": "unittest/benchmarksTest.py::test_bench_density_of_states_ilt",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.03824570798315108,
-                "max": 0.038739207899197936,
-                "mean": 0.038405808387324214,
-                "stddev": 0.0001997663437607365,
-                "rounds": 5,
-                "median": 0.03832020913250744,
-                "iqr": 0.00023793673608452082,
-                "q1": 0.038275614962913096,
-                "q3": 0.03851355169899762,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.03824570798315108,
-                "hd15iqr": 0.038739207899197936,
-                "ops": 26.037728197645453,
-                "total": 0.19202904193662107,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_states_construction",
-            "fullname": "unittest/benchmarksTest.py::test_bench_states_construction",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 5.00003807246685e-07,
-                "max": 4.391651600599289e-06,
-                "mean": 5.282627402545919e-07,
-                "stddev": 3.937120054251051e-08,
-                "rounds": 82475,
-                "median": 5.229027010500431e-07,
-                "iqr": 1.040752977132793e-08,
-                "q1": 5.187466740608215e-07,
-                "q3": 5.291542038321495e-07,
-                "iqr_outliers": 6764,
-                "stddev_outliers": 4044,
-                "outliers": "4044;6764",
-                "ld15iqr": 5.04148192703724e-07,
-                "hd15iqr": 5.457899533212185e-07,
-                "ops": 1892997.4117009619,
-                "total": 0.043568469502497466,
-                "iterations": 20
-            }
-        }
-    ],
-    "datetime": "2025-11-30T21:01:24.508702+00:00",
-    "version": "5.2.3"
-}
diff --git a/.benchmarks/Darwin-CPython-3.12-64bit/0005_latest.json b/.benchmarks/Darwin-CPython-3.12-64bit/0005_latest.json
deleted file mode 100644
index aa0e986..0000000
--- a/.benchmarks/Darwin-CPython-3.12-64bit/0005_latest.json
+++ /dev/null
@@ -1,183 +0,0 @@
-{
-    "machine_info": {
-        "node": "Georges-Mini",
-        "processor": "arm",
-        "machine": "arm64",
-        "python_compiler": "Clang 18.1.8 ",
-        "python_implementation": "CPython",
-        "python_implementation_version": "3.12.10",
-        "python_version": "3.12.10",
-        "python_build": [
-            "main",
-            "Apr 10 2025 22:19:24"
-        ],
-        "release": "25.1.0",
-        "system": "Darwin",
-        "cpu": {
-            "python_version": "3.12.10.final.0 (64 bit)",
-            "cpuinfo_version": [
-                9,
-                0,
-                0
-            ],
-            "cpuinfo_version_string": "9.0.0",
-            "arch": "ARM_8",
-            "bits": 64,
-            "count": 10,
-            "arch_string_raw": "arm64",
-            "brand_raw": "Apple M4"
-        }
-    },
-    "commit_info": {
-        "id": "659c1303a77acd9517a592564d2b16e15534ff26",
-        "time": "2025-11-30T15:37:16-05:00",
-        "author_time": "2025-11-30T15:37:16-05:00",
-        "dirty": true,
-        "project": "ChemPy",
-        "branch": "master"
-    },
-    "benchmarks": [
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_benzene",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_benzene",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.0002497918903827667,
-                "max": 0.0002986670006066561,
-                "mean": 0.00026836679317057135,
-                "stddev": 2.0308224101156007e-05,
-                "rounds": 5,
-                "median": 0.0002666250802576542,
-                "iqr": 3.159424522891641e-05,
-                "q1": 0.000250291486736387,
-                "q3": 0.0002818857319653034,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.0002497918903827667,
-                "hd15iqr": 0.0002986670006066561,
-                "ops": 3726.2434304396584,
-                "total": 0.0013418339658528566,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_ethane_rotors",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_ethane_rotors",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.7582900151610374e-05,
-                "max": 0.00011350004933774471,
-                "mean": 5.0070155716427545e-05,
-                "stddev": 3.071057377016809e-06,
-                "rounds": 8993,
-                "median": 4.929094575345516e-05,
-                "iqr": 8.33999365568161e-07,
-                "q1": 4.891608841717243e-05,
-                "q3": 4.975008778274059e-05,
-                "iqr_outliers": 1214,
-                "stddev_outliers": 733,
-                "outliers": "733;1214",
-                "ld15iqr": 4.7666020691394806e-05,
-                "hd15iqr": 5.1040900871157646e-05,
-                "ops": 19971.97703285571,
-                "total": 0.4502809103578329,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_density_of_states_ilt",
-            "fullname": "unittest/benchmarksTest.py::test_bench_density_of_states_ilt",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.03873725002631545,
-                "max": 0.039327667094767094,
-                "mean": 0.038956049783155325,
-                "stddev": 0.00023144431210765445,
-                "rounds": 5,
-                "median": 0.03895545797422528,
-                "iqr": 0.0002846981515176594,
-                "q1": 0.03877571824705228,
-                "q3": 0.03906041639856994,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.03873725002631545,
-                "hd15iqr": 0.039327667094767094,
-                "ops": 25.66995384712754,
-                "total": 0.1947802489157766,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_states_construction",
-            "fullname": "unittest/benchmarksTest.py::test_bench_states_construction",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 5.00003807246685e-07,
-                "max": 4.854146391153336e-06,
-                "mean": 5.320958957982121e-07,
-                "stddev": 5.2786756999192506e-08,
-                "rounds": 80809,
-                "median": 5.228910595178605e-07,
-                "iqr": 1.2549571692943594e-08,
-                "q1": 5.166511982679367e-07,
-                "q3": 5.292007699608803e-07,
-                "iqr_outliers": 8575,
-                "stddev_outliers": 3113,
-                "outliers": "3113;8575",
-                "ld15iqr": 5.00003807246685e-07,
-                "hd15iqr": 5.499925464391708e-07,
-                "ops": 1879360.4835080935,
-                "total": 0.04299813724355772,
-                "iterations": 20
-            }
-        }
-    ],
-    "datetime": "2025-11-30T21:08:01.357121+00:00",
-    "version": "5.2.3"
-}
diff --git a/.benchmarks/Darwin-CPython-3.12-64bit/0006_latest.json b/.benchmarks/Darwin-CPython-3.12-64bit/0006_latest.json
deleted file mode 100644
index 5b02d38..0000000
--- a/.benchmarks/Darwin-CPython-3.12-64bit/0006_latest.json
+++ /dev/null
@@ -1,183 +0,0 @@
-{
-    "machine_info": {
-        "node": "Georges-Mini",
-        "processor": "arm",
-        "machine": "arm64",
-        "python_compiler": "Clang 18.1.8 ",
-        "python_implementation": "CPython",
-        "python_implementation_version": "3.12.10",
-        "python_version": "3.12.10",
-        "python_build": [
-            "main",
-            "Apr 10 2025 22:19:24"
-        ],
-        "release": "25.1.0",
-        "system": "Darwin",
-        "cpu": {
-            "python_version": "3.12.10.final.0 (64 bit)",
-            "cpuinfo_version": [
-                9,
-                0,
-                0
-            ],
-            "cpuinfo_version_string": "9.0.0",
-            "arch": "ARM_8",
-            "bits": 64,
-            "count": 10,
-            "arch_string_raw": "arm64",
-            "brand_raw": "Apple M4"
-        }
-    },
-    "commit_info": {
-        "id": "659c1303a77acd9517a592564d2b16e15534ff26",
-        "time": "2025-11-30T15:37:16-05:00",
-        "author_time": "2025-11-30T15:37:16-05:00",
-        "dirty": true,
-        "project": "ChemPy",
-        "branch": "master"
-    },
-    "benchmarks": [
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_benzene",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_benzene",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.00023345905356109142,
-                "max": 0.0002743750810623169,
-                "mean": 0.00024532503448426723,
-                "stddev": 1.6727842485121804e-05,
-                "rounds": 5,
-                "median": 0.0002387501299381256,
-                "iqr": 1.6510195564478636e-05,
-                "q1": 0.00023523950949311256,
-                "q3": 0.0002517497050575912,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.00023345905356109142,
-                "hd15iqr": 0.0002743750810623169,
-                "ops": 4076.224842288284,
-                "total": 0.0012266251724213362,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "molecule",
-            "name": "test_bench_molecule_from_smiles_ethane_rotors",
-            "fullname": "unittest/benchmarksTest.py::test_bench_molecule_from_smiles_ethane_rotors",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.812492989003658e-05,
-                "max": 0.00010950001887977123,
-                "mean": 5.0903510526966985e-05,
-                "stddev": 3.158549425730003e-06,
-                "rounds": 13282,
-                "median": 5.0083035603165627e-05,
-                "iqr": 1.2081582099199295e-06,
-                "q1": 4.9582915380597115e-05,
-                "q3": 5.0791073590517044e-05,
-                "iqr_outliers": 1549,
-                "stddev_outliers": 1128,
-                "outliers": "1128;1549",
-                "ld15iqr": 4.812492989003658e-05,
-                "hd15iqr": 5.262484773993492e-05,
-                "ops": 19645.010523787612,
-                "total": 0.6761004268191755,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_density_of_states_ilt",
-            "fullname": "unittest/benchmarksTest.py::test_bench_density_of_states_ilt",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 0.038504458963871,
-                "max": 0.03934745816513896,
-                "mean": 0.03882397529669106,
-                "stddev": 0.00032009935441353494,
-                "rounds": 5,
-                "median": 0.038748042192310095,
-                "iqr": 0.00035987450974062085,
-                "q1": 0.03862152132205665,
-                "q3": 0.03898139583179727,
-                "iqr_outliers": 0,
-                "stddev_outliers": 1,
-                "outliers": "1;0",
-                "ld15iqr": 0.038504458963871,
-                "hd15iqr": 0.03934745816513896,
-                "ops": 25.757279937410978,
-                "total": 0.1941198764834553,
-                "iterations": 1
-            }
-        },
-        {
-            "group": "states",
-            "name": "test_bench_states_construction",
-            "fullname": "unittest/benchmarksTest.py::test_bench_states_construction",
-            "params": null,
-            "param": null,
-            "extra_info": {},
-            "options": {
-                "disable_gc": false,
-                "timer": "perf_counter",
-                "min_rounds": 5,
-                "max_time": 1.0,
-                "min_time": 5e-06,
-                "warmup": false
-            },
-            "stats": {
-                "min": 4.937406629323959e-07,
-                "max": 3.8000056520104407e-06,
-                "mean": 5.247793831464644e-07,
-                "stddev": 4.439141982743334e-08,
-                "rounds": 81914,
-                "median": 5.166977643966674e-07,
-                "iqr": 1.040752977132793e-08,
-                "q1": 5.124951712787152e-07,
-                "q3": 5.229027010500431e-07,
-                "iqr_outliers": 7398,
-                "stddev_outliers": 3940,
-                "outliers": "3940;7398",
-                "ld15iqr": 4.978966899216175e-07,
-                "hd15iqr": 5.395500920712948e-07,
-                "ops": 1905562.6652179337,
-                "total": 0.04298677839105949,
-                "iterations": 20
-            }
-        }
-    ],
-    "datetime": "2025-11-30T21:11:14.075478+00:00",
-    "version": "5.2.3"
-}
diff --git a/.github/workflows/benchmarks.yml b/.github/workflows/benchmarks.yml
deleted file mode 100644
index ebd8fa4..0000000
--- a/.github/workflows/benchmarks.yml
+++ /dev/null
@@ -1,54 +0,0 @@
-name: Benchmarks
-
-on:
-  push:
-    branches: [ master ]
-    paths:
-      - 'chempy/**'
-      - 'benchmarks/**'
-      - 'unittest/benchmarksTest.py'
-      - 'pyproject.toml'
-      - '.github/workflows/benchmarks.yml'
-  workflow_dispatch:  # Manual trigger
-
-permissions:
-  contents: read
-
-jobs:
-  benchmark:
-    runs-on: ubuntu-latest
-    # Note: Cython 3.2.2 compilation currently fails due to compatibility issues
-    # with the codebase (designed for Cython 2.x). Running pure Python benchmarks only.
-    # To compare with Cython performance, compile locally with Cython<3.
-
-    steps:
-    - uses: actions/checkout@v4
-
-    - name: Set up Python 3.12
-      uses: actions/setup-python@v5
-      with:
-        python-version: '3.12'
-        cache: 'pip'
-
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip setuptools wheel
-        pip install numpy scipy openbabel-wheel pytest pytest-benchmark
-        pip install -e ".[dev]"
-
-    - name: Run legacy benchmarks
-      run: |
-        pytest -q unittest/benchmarksTest.py --benchmark-only --benchmark-min-rounds=1 --benchmark-autosave
-
-    - name: Run new benchmarks (pure Python)
-      run: |
-        pytest benchmarks/ --benchmark-only --benchmark-json=benchmark-python.json
-
-    - name: Upload benchmark results
-      uses: actions/upload-artifact@v4
-      with:
-        name: benchmark-results-python-py3.12
-        path: |
-          benchmark-python.json
-          .benchmarks/**
-        if-no-files-found: ignore
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 72b1b98..4f2eb9a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -2,167 +2,31 @@ name: CI
 
 on:
   push:
-    branches: [ master ]
+    branches: [ master, rust-conversion ]
   pull_request:
-    branches: [ master ]
+    branches: [ master, rust-conversion ]
 
 permissions:
   contents: read
-  actions: read
 
 jobs:
-  test-and-type:
-    runs-on: macos-latest
+  test:
+    name: Test and Lint
+    runs-on: ubuntu-latest
     steps:
       - name: Checkout
         uses: actions/checkout@v4
 
-      - name: Set up Python
-        uses: actions/setup-python@v5
+      - name: Install Rust
+        uses: dtolnay/rust-toolchain@stable
         with:
-          python-version: '3.12'
+          components: rustfmt, clippy
 
-      - name: Install dependencies
-        run: |
-          python -m pip install --upgrade pip
-          pip install -e .
-          pip install pytest pytest-benchmark mypy
-          # Open Babel wheel for macOS/Linux
-          pip install openbabel-wheel || true
+      - name: Check formatting
+        run: cargo fmt -- --check
 
-      - name: Type check (strict for key modules)
-        run: |
-          mypy chempy/graph.py chempy/molecule.py --show-error-codes --check-untyped-defs
+      - name: Run clippy
+        run: cargo clippy -- -D warnings
 
-      - name: Run tests and save benchmarks
-        run: |
-          # Create a dummy benchmark that has no dependencies to verify pytest-benchmark is working
-          echo "def test_simple_bench(benchmark): benchmark(lambda: sum(range(100)))" > benchmarks/test_simple_bench.py
-          
-          # Print environment info for debugging
-          python --version
-          pip list | grep -E "pytest|benchmark|openbabel"
-          python -c "import openbabel; print('OpenBabel OK') if hasattr(openbabel, 'OBConversion') else print('OpenBabel partial')" || echo "OpenBabel import failed"
-          
-          # Run pytest with explicit config and multiple output formats
-          # We use --benchmark-json to force a file creation we can definitely find
-          python -m pytest -v tests/ unittest/ benchmarks/ --benchmark-autosave --benchmark-json=benchmarks-result.json
-
-      - name: List benchmark files (diagnostic)
-        run: |
-          echo "Checking for .benchmarks directory:"
-          ls -R .benchmarks || echo ".benchmarks directory not found"
-          echo "Checking for explicit JSON output:"
-          ls -l benchmarks-result.json || echo "benchmarks-result.json not found"
-          echo "Searching for any JSON files generated:"
-          find . -name "*.json" -not -path "./.git/*"
-
-      - name: Upload benchmark artifacts
-        if: always()
-        uses: actions/upload-artifact@v4
-        with:
-          name: benchmarks-json
-          path: |
-            .benchmarks/
-            benchmarks-result.json
-          if-no-files-found: warn
-
-  compare-artifacts:
-    runs-on: macos-latest
-    needs: test-and-type
-    env:
-      REGRESS_THRESHOLD_OPS: "-10.0"
-      REGRESS_THRESHOLD_MEAN: "-10.0"
-      REGRESS_THRESHOLD_MEDIAN: "-10.0"
-      REGRESS_FILTER_REGEX: ""
-    steps:
-      - name: Checkout
-        uses: actions/checkout@v4
-
-      - name: Set up Python
-        uses: actions/setup-python@v5
-        with:
-          python-version: '3.12'
-
-      - name: Install dependencies
-        run: |
-          python -m pip install --upgrade pip
-          pip install -e .
-          pip install pytest pytest-benchmark
-
-      - name: Download benchmark artifacts
-        uses: actions/download-artifact@v4
-        with:
-          name: benchmarks-json
-          path: .benchmarks
-      - name: Generate CSV and JSON comparisons
-        shell: bash
-        run: |
-          # Latest summary JSON for 'states' benchmarks
-          python3 scripts/compare_benchmarks.py --latest 1 --regex 'states_' --output json --save benchmark_states.json || true
-          # Two-run comparison CSV for molecule benchmarks
-          python3 scripts/compare_benchmarks.py --latest 2 --group molecule --output csv --save benchmark_molecule.csv || true
-
-      - name: Fail on significant performance regressions
-        shell: bash
-        run: |
-          python - <<'PY'
-          import csv, sys, os, re
-          try:
-            with open('benchmark_molecule.csv', newline='') as f:
-              reader = csv.DictReader(f)
-              worst_ops = 0.0
-              worst_mean = 0.0
-              worst_median = 0.0
-              rows = list(reader)
-              # Optional name regex filter
-              pattern = os.environ.get('REGRESS_FILTER_REGEX', '')
-              if pattern:
-                try:
-                  rx = re.compile(pattern)
-                  rows = [r for r in rows if rx.search(r.get('name',''))]
-                except re.error:
-                  print(f"Invalid REGRESS_FILTER_REGEX: {pattern}")
-              for r in rows:
-                for key, target in (
-                  ('ops_delta', 'ops'),
-                  ('mean_delta', 'mean'),
-                  ('median_delta', 'median'),
-                ):
-                  s = r.get(key, '')
-                  if s.endswith('%'):
-                    val = float(s.strip('%'))
-                    if target == 'ops' and val < worst_ops:
-                      worst_ops = val
-                    if target == 'mean' and val < worst_mean:
-                      worst_mean = val
-                    if target == 'median' and val < worst_median:
-                      worst_median = val
-              def get_thr(name, default):
-                try:
-                  return float(os.environ.get(name, str(default)))
-                except Exception:
-                  return default
-              threshold_ops = get_thr('REGRESS_THRESHOLD_OPS', -10.0)
-              threshold_mean = get_thr('REGRESS_THRESHOLD_MEAN', -10.0)
-              threshold_median = get_thr('REGRESS_THRESHOLD_MEDIAN', -10.0)
-              msg = (
-                f"Worst deltas — ops: {worst_ops:.2f}%, mean: {worst_mean:.2f}%, median: {worst_median:.2f}%"
-              )
-              print(msg)
-              if worst_ops < threshold_ops or worst_mean < threshold_mean or worst_median < threshold_median:
-                print("Regression detected beyond thresholds.")
-                sys.exit(1)
-              print("No regressions beyond thresholds.")
-          except FileNotFoundError:
-            print('No benchmark_molecule.csv found; skipping regression check')
-          PY
-
-      - name: Upload comparison artifacts
-        uses: actions/upload-artifact@v4
-        with:
-          name: benchmark-comparisons
-          path: |
-            benchmark_states.json
-            benchmark_molecule.csv
-          if-no-files-found: ignore
+      - name: Run tests
+        run: cargo test
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
deleted file mode 100644
index 5d8d41f..0000000
--- a/.github/workflows/docs.yml
+++ /dev/null
@@ -1,44 +0,0 @@
-name: Deploy Documentation
-
-on:
-  push:
-    branches:
-      - master
-    paths:
-      - 'chempy/**'
-      - 'documentation/**'
-      - '.github/workflows/docs.yml'
-  workflow_dispatch:
-
-permissions:
-  contents: write
-
-jobs:
-  build-and-deploy:
-    runs-on: ubuntu-latest
-    steps:
-      - name: Checkout
-        uses: actions/checkout@v4
-        with:
-          fetch-depth: 0
-
-      - name: Set up Python
-        uses: actions/setup-python@v5
-        with:
-          python-version: '3.12'
-
-      - name: Install dependencies
-        run: |
-          python -m pip install --upgrade pip
-          pip install .[docs]
-
-      - name: Build documentation
-        run: |
-          make docs
-
-      - name: Deploy to GitHub Pages
-        uses: JamesIves/github-pages-deploy-action@v4
-        with:
-          folder: documentation/build/html
-          branch: gh-pages
-          clean: true
diff --git a/.github/workflows/lint-and-test.yml b/.github/workflows/lint-and-test.yml
deleted file mode 100644
index 268186e..0000000
--- a/.github/workflows/lint-and-test.yml
+++ /dev/null
@@ -1,62 +0,0 @@
-name: Lint and Test
-
-on:
-  workflow_dispatch:
-  push:
-    branches: [ "master" ]
-  pull_request:
-    branches: [ "master" ]
-
-permissions:
-  contents: read
-
-jobs:
-  build:
-    runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python-version: ["3.12", "3.13"]
-
-    steps:
-      - name: Checkout
-        uses: actions/checkout@v4
-
-      - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v5
-        with:
-          python-version: ${{ matrix.python-version }}
-
-      - name: Install dependencies
-        run: |
-          python -m pip install --upgrade pip
-          python -m pip install -e .[dev,test]
-          python -m pip install openbabel-wheel
-
-      - name: Check formatting with black (line length 120)
-        run: |
-          python -m pip install black
-          python -m black --version
-          python -m black --check --line-length=120 chempy unittest tests
-
-      - name: Lint with flake8
-        run: |
-          python -m pip install flake8
-          python -m flake8
-
-      - name: Run tests
-        run: |
-          pytest -v --cov=chempy --cov-report=xml
-
-      - name: Upload coverage to Codecov
-        uses: codecov/codecov-action@v4
-        with:
-          token: ${{ secrets.CODECOV_TOKEN }}
-          files: ./coverage.xml
-          flags: unit
-          fail_ci_if_error: true
-
-      - name: Mypy type check
-        run: |
-          python -m pip install mypy
-          mypy --version
-          mypy chempy
diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml
deleted file mode 100644
index 60754d8..0000000
--- a/.github/workflows/pre-commit.yml
+++ /dev/null
@@ -1,26 +0,0 @@
-name: Pre-commit
-
-on:
-  pull_request:
-    branches: [ master, main, develop ]
-
-permissions:
-  contents: read
-
-jobs:
-  pre-commit:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v4
-      - name: Set up Python 3.12
-        uses: actions/setup-python@v5
-        with:
-          python-version: '3.12'
-          cache: 'pip'
-      - name: Install pre-commit and deps
-        run: |
-          python -m pip install --upgrade pip
-          pip install -e ".[dev]"
-          pre-commit install
-      - name: Run pre-commit on all files
-        run: pre-commit run --all-files --show-diff-on-failure
diff --git a/.github/workflows/smoke.yml b/.github/workflows/smoke.yml
deleted file mode 100644
index 2f3c092..0000000
--- a/.github/workflows/smoke.yml
+++ /dev/null
@@ -1,27 +0,0 @@
-name: Smoke Tests
-
-on:
-  pull_request:
-    branches: [ master, main, develop ]
-
-permissions:
-  contents: read
-
-jobs:
-  smoke:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v4
-      - name: Set up Python 3.12
-        uses: actions/setup-python@v5
-        with:
-          python-version: '3.12'
-          cache: 'pip'
-      - name: Install dependencies
-        run: |
-          python -m pip install --upgrade pip setuptools wheel
-          pip install numpy scipy
-          pip install -e ".[dev]"
-      - name: Run fast subset
-        run: |
-          pytest -q unittest/geometryTest.py unittest/graphTest.py unittest/moleculeTest.py unittest/reactionTest.py unittest/statesTest.py unittest/thermoTest.py
diff --git a/.github/workflows/stubs.yml b/.github/workflows/stubs.yml
deleted file mode 100644
index 6d948e4..0000000
--- a/.github/workflows/stubs.yml
+++ /dev/null
@@ -1,26 +0,0 @@
-name: Stub Type Check
-
-on:
-  pull_request:
-    branches: [ master, main, develop ]
-
-permissions:
-  contents: read
-
-jobs:
-  mypy-stubs:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v4
-      - name: Set up Python 3.12
-        uses: actions/setup-python@v5
-        with:
-          python-version: '3.12'
-          cache: 'pip'
-      - name: Install mypy and package
-        run: |
-          python -m pip install --upgrade pip
-          pip install -e ".[dev]"
-      - name: Run mypy (strict) on stubs
-        run: |
-          mypy --strict chempy/ext/*.pyi chempy/io/*.pyi
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
deleted file mode 100644
index fd22f07..0000000
--- a/.github/workflows/tests.yml
+++ /dev/null
@@ -1,114 +0,0 @@
-name: Tests
-
-on:
-  push:
-    branches: [ master, main, develop ]
-  pull_request:
-    branches: [ master, main, develop ]
-  schedule:
-    # Run tests daily at 2 AM UTC to catch any dependency issues
-    - cron: '0 2 * * *'
-
-permissions:
-  contents: read
-
-jobs:
-  test:
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        os: [ubuntu-latest, macos-latest]
-        python-version: ['3.8', '3.9', '3.10', '3.11', '3.12', '3.13']
-        exclude:
-          - os: windows-latest
-            python-version: '3.8'
-          - os: windows-latest
-            python-version: '3.9'
-          - os: windows-latest
-            python-version: '3.10'
-          - os: windows-latest
-            python-version: '3.11'
-          - os: windows-latest
-            python-version: '3.13'
-
-    steps:
-    - uses: actions/checkout@v4
-
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v5
-      with:
-        python-version: ${{ matrix.python-version }}
-        cache: 'pip'
-
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip setuptools wheel
-        pip install numpy scipy
-        pip install -e ".[dev]"
-    - name: Cache pip wheels
-      uses: actions/cache@v4
-      with:
-        path: ~/.cache/pip
-        key: pip-${{ runner.os }}-${{ matrix.python-version }}-${{ hashFiles('pyproject.toml') }}
-        restore-keys: |
-          pip-${{ runner.os }}-${{ matrix.python-version }}-
-
-    - name: Install OpenBabel
-      run: pip install openbabel-wheel
-      continue-on-error: true
-
-    - name: Run tests with pytest
-      run: |
-        if [[ "${{ matrix.os }}" == "ubuntu-latest" && "${{ matrix.python-version }}" == "3.12" ]]; then
-          pytest -v --cov=chempy --cov-report=xml
-        else
-          pytest -q
-        fi
-
-    - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v4
-      if: matrix.os == 'ubuntu-latest' && matrix.python-version == '3.12'
-      with:
-        files: ./coverage.xml
-        flags: unittests
-        name: codecov-umbrella
-        fail_ci_if_error: false
-        token: ${{ secrets.CODECOV_TOKEN }}
-
-  quality:
-    runs-on: ubuntu-latest
-    strategy:
-      matrix:
-        python-version: ['3.12', '3.13']
-
-    steps:
-    - uses: actions/checkout@v4
-
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v5
-      with:
-        python-version: ${{ matrix.python-version }}
-        cache: 'pip'
-
-    - name: Install dependencies
-      run: |
-        python -m pip install --upgrade pip
-        pip install -e ".[dev]"
-
-    - name: Check formatting with black
-      run: black --check chempy unittest
-      continue-on-error: true
-
-    - name: Check import sorting with isort
-      run: isort --check-only chempy unittest
-      continue-on-error: true
-
-    - name: Lint with flake8
-      run: |
-        flake8 chempy unittest --max-line-length=100 --extend-ignore=E203,W503,E501
-      continue-on-error: true
-
-    - name: Type check with mypy
-      run: mypy chempy --ignore-missing-imports
-      continue-on-error: true
diff --git a/.gitignore b/.gitignore
index 054c1d8..615a4f7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -83,3 +83,12 @@ make.inc
 .cache/
 *.egg-info
 .pytest_cache/
+
+
+# Added by cargo
+
+/target
+
+# Rust
+target/
+Cargo.lock
diff --git a/Cargo.toml b/Cargo.toml
new file mode 100644
index 0000000..33da035
--- /dev/null
+++ b/Cargo.toml
@@ -0,0 +1,13 @@
+[package]
+name = "chempy"
+version = "0.1.0"
+edition = "2024"
+
+[dependencies]
+
+[dev-dependencies]
+criterion = "0.5"
+
+[[bench]]
+name = "chempy_benchmarks"
+harness = false
diff --git a/README.md b/README.md
index 636d965..20d0058 100644
--- a/README.md
+++ b/README.md
@@ -1,256 +1,51 @@
-# ChemPy Toolkit
+# ChemPy (Rust)
 
-[![Python 3.8+](https://img.shields.io/badge/python-3.8%2B-blue.svg)](https://www.python.org/downloads/)
-[![Python 3.13](https://img.shields.io/badge/python-3.13-brightgreen.svg)](https://www.python.org/downloads/)
 [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
-[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
-[![Checked with mypy](http://www.mypy-lang.org/static/mypy_badge.svg)](http://mypy-lang.org/)
-[![Lint & Test](https://github.com/elkins/ChemPy/actions/workflows/lint-and-test.yml/badge.svg?branch=master)](https://github.com/elkins/ChemPy/actions/workflows/lint-and-test.yml)
-[![Tests](https://github.com/elkins/ChemPy/actions/workflows/tests.yml/badge.svg?branch=master)](https://github.com/elkins/ChemPy/actions/workflows/tests.yml)
-[![Codecov](https://codecov.io/gh/elkins/ChemPy/branch/master/graph/badge.svg)](https://codecov.io/gh/elkins/ChemPy)
-[![PEP 561 Compliant](https://img.shields.io/badge/pep-561-blue.svg)](https://www.python.org/dev/peps/pep-0561/)
-[![Benchmarks](https://github.com/elkins/ChemPy/actions/workflows/benchmarks.yml/badge.svg?branch=master)](https://github.com/elkins/ChemPy/actions/workflows/benchmarks.yml?query=branch%3Amaster)
 
 > [!NOTE]
-> **Project Status & RMG-Py Integration**
-> This repository represents the foundational toolkit originally developed by Joshua W. Allen. While this standalone version remains a stable and lightweight library for molecular graphs, thermodynamics, and kinetics, it has effectively been integrated into and superseded by the **[RMG-Py (Reaction Mechanism Generator)](https://github.com/ReactionMechanismGenerator/RMG-Py)** ecosystem.
->
-> Most active development, including new thermodynamics models and kinetics solvers, now takes place within the RMG-Py project. If you are looking for the most feature-rich and actively maintained version of these tools, we recommend using [RMG-Py](https://github.com/ReactionMechanismGenerator/RMG-Py).
+> **Project Evolution**
+> This project has been converted from its original Python implementation to a **Pure Rust crate**. 
+> While the Python version has been integrated into the [RMG-Py](https://github.com/ReactionMechanismGenerator/RMG-Py) ecosystem, this repository now serves as a high-performance Rust port of the foundational ChemPy toolkit.
 
-**ChemPy Toolkit** is a free, open-source Python toolkit for chemistry, chemical engineering, and materials science applications, with a focus on molecular structures, thermodynamics, and chemical kinetics.
-
-> [!IMPORTANT]
-> **Naming & Installation Notice**
-> This project is the **ChemPy Toolkit** (distribution name: `chempy-toolkit`), originally developed by Joshua W. Allen as part of the [RMG](https://rmgpy.github.io/) ecosystem.
->
-> It is **distinct** from the general-purpose `chempy` package on PyPI by Björn Dahlgren.
-> - To install this toolkit, use: `pip install chempy-toolkit`
-> - Once installed, it is imported as: `import chempy`
-
-## Quick Links
-
-- 📖 **[Documentation](https://elkins.github.io/ChemPy)** - Full documentation and API reference
-- 🐛 **[Issue Tracker](https://github.com/elkins/ChemPy/issues)** - Report bugs and request features
-- 📝 **[Contributing](CONTRIBUTING.md)** - How to contribute to ChemPy
-- 📋 **[Changelog](CHANGELOG.md)** - Version history and release notes
-- 🔐 **[Security](SECURITY.md)** - Security policy and vulnerability reporting
-- 🔧 **[TODO](TODO.md)** - Future improvements and known issues
-- 👨‍💻 **[Developer Docs](docs/)** - Development guides and technical documentation
+**ChemPy (Rust)** is a high-performance toolkit for chemistry, chemical engineering, and materials science applications, with a focus on molecular structures, thermodynamics, and chemical kinetics.
 
 ## Features
+- **Fast Graph Isomorphism:** Efficient VF2-based molecular graph comparison.
+- **Thermodynamic Models:** Support for NASA polynomials and Wilhoit models.
+- **Kinetics:** Arrhenius models and rate coefficient calculations.
+- **States:** Partition function and density of states calculations.
 
-- Python 3.13 support: Updated and tested on latest Python.
-- Open Babel 3.x integration: Modern molecular format handling.
-- Type hints (PEP 561): Full type annotation coverage with `py.typed`.
-- Test suite: All tests passing; legacy and modern suites maintained.
-- Code quality: `black`, `isort`, `flake8`, and mypy checks.
-- GitHub Actions CI/CD: Automated linting and testing across Python 3.8–3.13.
-- NumPy compatibility: Addressed array-to-scalar deprecations.
-- Modern packaging: PEP 517/518 with `pyproject.toml`.
-
-## Platform Support
-
-**Windows:** Experimental. Unit tests are not run on Windows in CI due to persistent failures and lack of a Windows development environment. Use at your own risk.
-
-If you are able to help improve Windows compatibility, contributions and fixes are very welcome!
-
-**macOS and Linux:** Fully supported and tested in CI.
-## Installation
-
-### Requirements
-
-
-Note: Features such as SMILES parsing and certain rotor-counting utilities depend on Open Babel. On macOS/Linux, install `openbabel-wheel` to enable these features. Windows support for Open Babel is currently experimental.
-
-### Optional Dependencies
-
-
-### Quick Start
-
-Install via pip:
-
-```bash
-pip install chempy-toolkit
+## Getting Started
+Add this to your `Cargo.toml`:
+```toml
+[dependencies]
+chempy = { git = "https://github.com/elkins/ChemPy.git", branch = "rust-conversion" }
 ```
 
-Or install from source with development dependencies:
-
+## Development
+Run tests with:
 ```bash
-git clone https://github.com/elkins/ChemPy.git
-cd ChemPy
-pip install -e ".[dev]"
-make build
+cargo test
 ```
 
+## Python Comparison (Legacy)
+The original Python implementation is preserved in the `python/` directory for behavioral and performance comparison.
 
-### Setup Development Environment
-
+To run the original Python tests:
 ```bash
-# Install with all development tools
-pip install -e ".[dev,docs]"
-
-# Install pre-commit hooks for automatic quality checks
-pre-commit install
-
-# Build Cython extensions for performance
-make build
+cd python
+pip install -e .
+pytest unittest/
 ```
 
-### Running Tests
-
-```bash
-# Run all tests
-make test
-
-# Run with coverage report
-make test-cov
-
-# Run tests in parallel
-make test-fast
-
-# Run specific test file
-pytest unittest/moleculeTest.py -v
-```
-
-### Benchmarking
-
-ChemPy includes a small benchmark suite using `pytest-benchmark` to track performance of key hot-paths (SMILES parsing, rotor counting, density-of-states ILT, etc.).
-
-Run locally:
-
+To run original Python benchmarks:
 ```bash
+cd python
 pytest unittest/benchmarksTest.py --benchmark-only
 ```
 
-Compare two runs (e.g., branch vs. main):
-
-```bash
-# On main
-pytest unittest/benchmarksTest.py --benchmark-only --benchmark-save=main
-
-# On your branch
-pytest unittest/benchmarksTest.py --benchmark-only --benchmark-save=feature
-
-# Compare
-pytest unittest/benchmarksTest.py --benchmark-only --benchmark-compare
-```
-
-CI runs a quick benchmark job on Ubuntu/Python 3.12 and uploads JSON results as an artifact for trend tracking.
-
-Latest CI benchmark artifacts (master):
-
-- https://github.com/elkins/ChemPy/actions/workflows/benchmarks.yml?query=branch%3Amaster
-
-Open the most recent run, then download the artifact named `benchmark-results-<OS>-py312`.
-
-### Code Quality
-
-```bash
-# Format code automatically
-make format
-
-# Check code style
-make lint
-
-# Type checking
-make type-check
-
-# All quality checks
-make check
-```
-
-### Building Documentation
-
-```bash
-make docs
-cd documentation
-open build/html/index.html
-```
-
-## Project Structure
-
-```
-chempy/
-├── constants.py          # Physical constants (SI units)
-├── element.py            # Element properties and data
-├── molecule.py           # Molecular structure representation
-├── reaction.py           # Chemical reactions
-├── kinetics.py           # Reaction kinetics modeling
-├── thermo.py             # Thermodynamic calculations
-├── species.py            # Chemical species definitions
-├── geometry.py           # Geometric calculations
-├── graph.py              # Graph-based molecular analysis
-├── pattern.py            # Pattern matching for molecules
-├── states.py             # Physical/chemical state variables
-├── py.typed              # PEP 561 type hint marker
-└── ext/                  # Extensions
-    ├── molecule_draw.py  # Molecular visualization
-    └── thermo_converter.py  # Thermodynamics converters
-
-tests/                     # Modern test directory
-unittest/                  # Legacy test suite
-docs/                      # Documentation
-documentation/             # Sphinx documentation source
-```
-
-## Documentation
-
- The Sphinx docs homepage includes a Codecov badge; see `documentation/build/html/index.html` after building.
- The contents page also shows the badge for quick visibility.
-
-## Manual CI
-
-- Purpose: Run lint (`flake8`) and tests (`pytest`) without pushing.
-- Trigger: Go to `Actions` → select `Lint & Test` → `Run workflow`.
-- Branch: Choose a branch and optionally a specific commit SHA.
-- Outputs: See job logs; test results appear inline in the workflow run.
-
-## Citation
-
-If you use ChemPy in your research, please cite:
-
-```bibtex
-@software{chempy2024,
-  title={ChemPy: A Chemistry Toolkit for Python},
-  author={Allen, Joshua W.},
-  year={2024},
-  url={https://github.com/elkins/ChemPy}
-}
-```
-
 ## License
-
-
-## Troubleshooting
-
-See the [Developer Documentation](docs/DEVELOPMENT.md) for detailed troubleshooting, including:
-- Coverage uploads and Codecov configuration
-- Type checking with mypy
-- Lint and formatting tools
-- CI debugging
-
-## License
-
 ChemPy is licensed under the MIT License - see [LICENSE](LICENSE) for details.
 
-## Related Projects
-
-- [RMG](https://rmgpy.github.io/) - Reaction Mechanism Generator
-- [Cantera](https://cantera.org/) - Chemical kinetics and thermodynamics
-- [OpenBabel](http://openbabel.org/) - Molecular structures and formats
-- [GAMESS](https://www.msg.chem.iastate.edu/gamess/) - Quantum chemistry
-
-## Support
-
-For questions and discussions:
-- Open an [issue](https://github.com/elkins/ChemPy/issues)
-- Read the [documentation](https://chempy.readthedocs.io)
-- Browse existing issues or propose enhancements via the Issue Tracker
-
 ## Acknowledgments
-
-ChemPy was originally developed by Joshua W. Allen and is maintained by the open-source community.
-
----
-
-Made with ❤️ for the chemistry community
+ChemPy was originally developed by Joshua W. Allen in Python and has been ported to Rust to ensure its long-term performance and maintainability.
diff --git a/benches/chempy_benchmarks.rs b/benches/chempy_benchmarks.rs
new file mode 100644
index 0000000..64074c9
--- /dev/null
+++ b/benches/chempy_benchmarks.rs
@@ -0,0 +1,42 @@
+use chempy::element;
+use chempy::kinetics::{ArrheniusModel, KineticsModel};
+use chempy::molecule::{Atom, Bond, BondOrder, Molecule};
+use criterion::{Criterion, black_box, criterion_group, criterion_main};
+
+fn setup_benzene() -> Molecule {
+    let mut benzene = Molecule::new();
+    let carbons: Vec<usize> = (0..6)
+        .map(|_| benzene.add_atom(Atom::new(&element::C)))
+        .collect();
+    for i in 0..6 {
+        let order = if i % 2 == 0 {
+            BondOrder::Double
+        } else {
+            BondOrder::Single
+        };
+        benzene.add_bond(carbons[i], carbons[(i + 1) % 6], Bond::new(order));
+    }
+    benzene
+}
+
+fn bench_isomorphism(c: &mut Criterion) {
+    let benzene1 = setup_benzene();
+    let benzene2 = setup_benzene();
+
+    c.bench_function("isomorphism_benzene", |b| {
+        b.iter(|| benzene1.is_isomorphic(black_box(&benzene2)))
+    });
+}
+
+fn bench_kinetics(c: &mut Criterion) {
+    let model = ArrheniusModel::new(1.0e10, 0.5, 50000.0, 1.0);
+    let t = 1000.0;
+    let p = 1.0e5;
+
+    c.bench_function("kinetics_arrhenius", |b| {
+        b.iter(|| model.get_rate_coefficient(black_box(t), black_box(p)))
+    });
+}
+
+criterion_group!(benches, bench_isomorphism, bench_kinetics);
+criterion_main!(benches);
diff --git a/.pre-commit-config.yaml b/python/.pre-commit-config.yaml
similarity index 100%
rename from .pre-commit-config.yaml
rename to python/.pre-commit-config.yaml
diff --git a/.python-version b/python/.python-version
similarity index 100%
rename from .python-version
rename to python/.python-version
diff --git a/MANIFEST.in b/python/MANIFEST.in
similarity index 100%
rename from MANIFEST.in
rename to python/MANIFEST.in
diff --git a/Makefile b/python/Makefile
similarity index 100%
rename from Makefile
rename to python/Makefile
diff --git a/benchmarks/README.md b/python/benchmarks/README.md
similarity index 100%
rename from benchmarks/README.md
rename to python/benchmarks/README.md
diff --git a/benchmarks/__init__.py b/python/benchmarks/__init__.py
similarity index 100%
rename from benchmarks/__init__.py
rename to python/benchmarks/__init__.py
diff --git a/benchmarks/benchmark_graph.py b/python/benchmarks/benchmark_graph.py
similarity index 100%
rename from benchmarks/benchmark_graph.py
rename to python/benchmarks/benchmark_graph.py
diff --git a/benchmarks/benchmark_kinetics.py b/python/benchmarks/benchmark_kinetics.py
similarity index 100%
rename from benchmarks/benchmark_kinetics.py
rename to python/benchmarks/benchmark_kinetics.py
diff --git a/benchmarks/compare_benchmarks.py b/python/benchmarks/compare_benchmarks.py
similarity index 100%
rename from benchmarks/compare_benchmarks.py
rename to python/benchmarks/compare_benchmarks.py
diff --git a/benchmarks/conftest.py b/python/benchmarks/conftest.py
similarity index 100%
rename from benchmarks/conftest.py
rename to python/benchmarks/conftest.py
diff --git a/chempy/__init__.py b/python/chempy/__init__.py
similarity index 100%
rename from chempy/__init__.py
rename to python/chempy/__init__.py
diff --git a/chempy/_cython_compat.py b/python/chempy/_cython_compat.py
similarity index 100%
rename from chempy/_cython_compat.py
rename to python/chempy/_cython_compat.py
diff --git a/chempy/constants.py b/python/chempy/constants.py
similarity index 100%
rename from chempy/constants.py
rename to python/chempy/constants.py
diff --git a/chempy/element.pxd b/python/chempy/element.pxd
similarity index 100%
rename from chempy/element.pxd
rename to python/chempy/element.pxd
diff --git a/chempy/element.py b/python/chempy/element.py
similarity index 100%
rename from chempy/element.py
rename to python/chempy/element.py
diff --git a/chempy/exception.py b/python/chempy/exception.py
similarity index 100%
rename from chempy/exception.py
rename to python/chempy/exception.py
diff --git a/chempy/ext/__init__.py b/python/chempy/ext/__init__.py
similarity index 100%
rename from chempy/ext/__init__.py
rename to python/chempy/ext/__init__.py
diff --git a/chempy/ext/molecule_draw.py b/python/chempy/ext/molecule_draw.py
similarity index 100%
rename from chempy/ext/molecule_draw.py
rename to python/chempy/ext/molecule_draw.py
diff --git a/chempy/ext/molecule_draw.pyi b/python/chempy/ext/molecule_draw.pyi
similarity index 100%
rename from chempy/ext/molecule_draw.pyi
rename to python/chempy/ext/molecule_draw.pyi
diff --git a/chempy/ext/thermo_converter.pxd b/python/chempy/ext/thermo_converter.pxd
similarity index 100%
rename from chempy/ext/thermo_converter.pxd
rename to python/chempy/ext/thermo_converter.pxd
diff --git a/chempy/ext/thermo_converter.py b/python/chempy/ext/thermo_converter.py
similarity index 100%
rename from chempy/ext/thermo_converter.py
rename to python/chempy/ext/thermo_converter.py
diff --git a/chempy/ext/thermo_converter.pyi b/python/chempy/ext/thermo_converter.pyi
similarity index 100%
rename from chempy/ext/thermo_converter.pyi
rename to python/chempy/ext/thermo_converter.pyi
diff --git a/chempy/geometry.pxd b/python/chempy/geometry.pxd
similarity index 100%
rename from chempy/geometry.pxd
rename to python/chempy/geometry.pxd
diff --git a/chempy/geometry.py b/python/chempy/geometry.py
similarity index 100%
rename from chempy/geometry.py
rename to python/chempy/geometry.py
diff --git a/chempy/graph.pxd b/python/chempy/graph.pxd
similarity index 100%
rename from chempy/graph.pxd
rename to python/chempy/graph.pxd
diff --git a/chempy/graph.py b/python/chempy/graph.py
similarity index 100%
rename from chempy/graph.py
rename to python/chempy/graph.py
diff --git a/chempy/io/__init__.py b/python/chempy/io/__init__.py
similarity index 100%
rename from chempy/io/__init__.py
rename to python/chempy/io/__init__.py
diff --git a/chempy/io/gaussian.py b/python/chempy/io/gaussian.py
similarity index 100%
rename from chempy/io/gaussian.py
rename to python/chempy/io/gaussian.py
diff --git a/chempy/io/gaussian.pyi b/python/chempy/io/gaussian.pyi
similarity index 100%
rename from chempy/io/gaussian.pyi
rename to python/chempy/io/gaussian.pyi
diff --git a/chempy/kinetics.pxd b/python/chempy/kinetics.pxd
similarity index 100%
rename from chempy/kinetics.pxd
rename to python/chempy/kinetics.pxd
diff --git a/chempy/kinetics.py b/python/chempy/kinetics.py
similarity index 100%
rename from chempy/kinetics.py
rename to python/chempy/kinetics.py
diff --git a/chempy/molecule.pxd b/python/chempy/molecule.pxd
similarity index 100%
rename from chempy/molecule.pxd
rename to python/chempy/molecule.pxd
diff --git a/chempy/molecule.py b/python/chempy/molecule.py
similarity index 100%
rename from chempy/molecule.py
rename to python/chempy/molecule.py
diff --git a/chempy/pattern.pxd b/python/chempy/pattern.pxd
similarity index 100%
rename from chempy/pattern.pxd
rename to python/chempy/pattern.pxd
diff --git a/chempy/pattern.py b/python/chempy/pattern.py
similarity index 100%
rename from chempy/pattern.py
rename to python/chempy/pattern.py
diff --git a/chempy/py.typed b/python/chempy/py.typed
similarity index 100%
rename from chempy/py.typed
rename to python/chempy/py.typed
diff --git a/chempy/reaction.pxd b/python/chempy/reaction.pxd
similarity index 100%
rename from chempy/reaction.pxd
rename to python/chempy/reaction.pxd
diff --git a/chempy/reaction.py b/python/chempy/reaction.py
similarity index 100%
rename from chempy/reaction.py
rename to python/chempy/reaction.py
diff --git a/chempy/species.pxd b/python/chempy/species.pxd
similarity index 100%
rename from chempy/species.pxd
rename to python/chempy/species.pxd
diff --git a/chempy/species.py b/python/chempy/species.py
similarity index 100%
rename from chempy/species.py
rename to python/chempy/species.py
diff --git a/chempy/states.pxd b/python/chempy/states.pxd
similarity index 100%
rename from chempy/states.pxd
rename to python/chempy/states.pxd
diff --git a/chempy/states.py b/python/chempy/states.py
similarity index 100%
rename from chempy/states.py
rename to python/chempy/states.py
diff --git a/chempy/thermo.pxd b/python/chempy/thermo.pxd
similarity index 100%
rename from chempy/thermo.pxd
rename to python/chempy/thermo.pxd
diff --git a/chempy/thermo.py b/python/chempy/thermo.py
similarity index 100%
rename from chempy/thermo.py
rename to python/chempy/thermo.py
diff --git a/docs/.gitkeep b/python/docs/.gitkeep
similarity index 100%
rename from docs/.gitkeep
rename to python/docs/.gitkeep
diff --git a/docs/DEVELOPMENT.md b/python/docs/DEVELOPMENT.md
similarity index 100%
rename from docs/DEVELOPMENT.md
rename to python/docs/DEVELOPMENT.md
diff --git a/docs/README.md b/python/docs/README.md
similarity index 100%
rename from docs/README.md
rename to python/docs/README.md
diff --git a/docs/STRUCTURE.md b/python/docs/STRUCTURE.md
similarity index 100%
rename from docs/STRUCTURE.md
rename to python/docs/STRUCTURE.md
diff --git a/docs/TYPE_HINTS.md b/python/docs/TYPE_HINTS.md
similarity index 100%
rename from docs/TYPE_HINTS.md
rename to python/docs/TYPE_HINTS.md
diff --git a/docs/__init__.py b/python/docs/__init__.py
similarity index 100%
rename from docs/__init__.py
rename to python/docs/__init__.py
diff --git a/docs/conf.py b/python/docs/conf.py
similarity index 100%
rename from docs/conf.py
rename to python/docs/conf.py
diff --git a/documentation/Makefile b/python/documentation/Makefile
similarity index 100%
rename from documentation/Makefile
rename to python/documentation/Makefile
diff --git a/documentation/make.bat b/python/documentation/make.bat
similarity index 100%
rename from documentation/make.bat
rename to python/documentation/make.bat
diff --git a/documentation/source/_static/chempy_logo.png b/python/documentation/source/_static/chempy_logo.png
similarity index 100%
rename from documentation/source/_static/chempy_logo.png
rename to python/documentation/source/_static/chempy_logo.png
diff --git a/documentation/source/_static/chempy_logo.svg b/python/documentation/source/_static/chempy_logo.svg
similarity index 100%
rename from documentation/source/_static/chempy_logo.svg
rename to python/documentation/source/_static/chempy_logo.svg
diff --git a/documentation/source/_static/default.css b/python/documentation/source/_static/default.css
similarity index 100%
rename from documentation/source/_static/default.css
rename to python/documentation/source/_static/default.css
diff --git a/documentation/source/_templates/index.html b/python/documentation/source/_templates/index.html
similarity index 100%
rename from documentation/source/_templates/index.html
rename to python/documentation/source/_templates/index.html
diff --git a/documentation/source/_templates/indexsidebar.html b/python/documentation/source/_templates/indexsidebar.html
similarity index 100%
rename from documentation/source/_templates/indexsidebar.html
rename to python/documentation/source/_templates/indexsidebar.html
diff --git a/documentation/source/_templates/layout.html b/python/documentation/source/_templates/layout.html
similarity index 100%
rename from documentation/source/_templates/layout.html
rename to python/documentation/source/_templates/layout.html
diff --git a/documentation/source/conf.py b/python/documentation/source/conf.py
similarity index 100%
rename from documentation/source/conf.py
rename to python/documentation/source/conf.py
diff --git a/documentation/source/constants.rst b/python/documentation/source/constants.rst
similarity index 100%
rename from documentation/source/constants.rst
rename to python/documentation/source/constants.rst
diff --git a/documentation/source/contents.rst b/python/documentation/source/contents.rst
similarity index 100%
rename from documentation/source/contents.rst
rename to python/documentation/source/contents.rst
diff --git a/documentation/source/element.rst b/python/documentation/source/element.rst
similarity index 100%
rename from documentation/source/element.rst
rename to python/documentation/source/element.rst
diff --git a/documentation/source/exception.rst b/python/documentation/source/exception.rst
similarity index 100%
rename from documentation/source/exception.rst
rename to python/documentation/source/exception.rst
diff --git a/documentation/source/geometry.rst b/python/documentation/source/geometry.rst
similarity index 100%
rename from documentation/source/geometry.rst
rename to python/documentation/source/geometry.rst
diff --git a/documentation/source/graph.rst b/python/documentation/source/graph.rst
similarity index 100%
rename from documentation/source/graph.rst
rename to python/documentation/source/graph.rst
diff --git a/documentation/source/introduction.rst b/python/documentation/source/introduction.rst
similarity index 100%
rename from documentation/source/introduction.rst
rename to python/documentation/source/introduction.rst
diff --git a/documentation/source/kinetics.rst b/python/documentation/source/kinetics.rst
similarity index 100%
rename from documentation/source/kinetics.rst
rename to python/documentation/source/kinetics.rst
diff --git a/documentation/source/molecule.rst b/python/documentation/source/molecule.rst
similarity index 100%
rename from documentation/source/molecule.rst
rename to python/documentation/source/molecule.rst
diff --git a/documentation/source/pattern.rst b/python/documentation/source/pattern.rst
similarity index 100%
rename from documentation/source/pattern.rst
rename to python/documentation/source/pattern.rst
diff --git a/documentation/source/reaction.rst b/python/documentation/source/reaction.rst
similarity index 100%
rename from documentation/source/reaction.rst
rename to python/documentation/source/reaction.rst
diff --git a/documentation/source/species.rst b/python/documentation/source/species.rst
similarity index 100%
rename from documentation/source/species.rst
rename to python/documentation/source/species.rst
diff --git a/documentation/source/states.rst b/python/documentation/source/states.rst
similarity index 100%
rename from documentation/source/states.rst
rename to python/documentation/source/states.rst
diff --git a/documentation/source/thermo.rst b/python/documentation/source/thermo.rst
similarity index 100%
rename from documentation/source/thermo.rst
rename to python/documentation/source/thermo.rst
diff --git a/pyproject.toml b/python/pyproject.toml
similarity index 100%
rename from pyproject.toml
rename to python/pyproject.toml
diff --git a/scripts/compare_benchmarks.py b/python/scripts/compare_benchmarks.py
similarity index 100%
rename from scripts/compare_benchmarks.py
rename to python/scripts/compare_benchmarks.py
diff --git a/setup.cfg b/python/setup.cfg
similarity index 100%
rename from setup.cfg
rename to python/setup.cfg
diff --git a/setup.py b/python/setup.py
similarity index 100%
rename from setup.py
rename to python/setup.py
diff --git a/tests/__init__.py b/python/tests/__init__.py
similarity index 100%
rename from tests/__init__.py
rename to python/tests/__init__.py
diff --git a/tests/conftest.py b/python/tests/conftest.py
similarity index 100%
rename from tests/conftest.py
rename to python/tests/conftest.py
diff --git a/tests/test_constants.py b/python/tests/test_constants.py
similarity index 100%
rename from tests/test_constants.py
rename to python/tests/test_constants.py
diff --git a/tests/test_element.py b/python/tests/test_element.py
similarity index 100%
rename from tests/test_element.py
rename to python/tests/test_element.py
diff --git a/tests/test_graph_iso.py b/python/tests/test_graph_iso.py
similarity index 100%
rename from tests/test_graph_iso.py
rename to python/tests/test_graph_iso.py
diff --git a/tests/test_kinetics_models.py b/python/tests/test_kinetics_models.py
similarity index 100%
rename from tests/test_kinetics_models.py
rename to python/tests/test_kinetics_models.py
diff --git a/tests/test_kinetics_smoke.py b/python/tests/test_kinetics_smoke.py
similarity index 100%
rename from tests/test_kinetics_smoke.py
rename to python/tests/test_kinetics_smoke.py
diff --git a/tests/test_molecule_min.py b/python/tests/test_molecule_min.py
similarity index 100%
rename from tests/test_molecule_min.py
rename to python/tests/test_molecule_min.py
diff --git a/tests/test_reaction_smoke.py b/python/tests/test_reaction_smoke.py
similarity index 100%
rename from tests/test_reaction_smoke.py
rename to python/tests/test_reaction_smoke.py
diff --git a/tests/test_species_smoke.py b/python/tests/test_species_smoke.py
similarity index 100%
rename from tests/test_species_smoke.py
rename to python/tests/test_species_smoke.py
diff --git a/tests/test_states_smoke.py b/python/tests/test_states_smoke.py
similarity index 100%
rename from tests/test_states_smoke.py
rename to python/tests/test_states_smoke.py
diff --git a/tests/test_thermo_models.py b/python/tests/test_thermo_models.py
similarity index 100%
rename from tests/test_thermo_models.py
rename to python/tests/test_thermo_models.py
diff --git a/tests/test_thermo_smoke.py b/python/tests/test_thermo_smoke.py
similarity index 100%
rename from tests/test_thermo_smoke.py
rename to python/tests/test_thermo_smoke.py
diff --git a/tests/test_tst_smoke.py b/python/tests/test_tst_smoke.py
similarity index 100%
rename from tests/test_tst_smoke.py
rename to python/tests/test_tst_smoke.py
diff --git a/tox.ini b/python/tox.ini
similarity index 100%
rename from tox.ini
rename to python/tox.ini
diff --git a/unittest/benchmarksTest.py b/python/unittest/benchmarksTest.py
similarity index 100%
rename from unittest/benchmarksTest.py
rename to python/unittest/benchmarksTest.py
diff --git a/unittest/conftest.py b/python/unittest/conftest.py
similarity index 100%
rename from unittest/conftest.py
rename to python/unittest/conftest.py
diff --git a/unittest/ethylene.log b/python/unittest/ethylene.log
similarity index 100%
rename from unittest/ethylene.log
rename to python/unittest/ethylene.log
diff --git a/unittest/gaussianTest.py b/python/unittest/gaussianTest.py
similarity index 100%
rename from unittest/gaussianTest.py
rename to python/unittest/gaussianTest.py
diff --git a/unittest/geometryTest.py b/python/unittest/geometryTest.py
similarity index 100%
rename from unittest/geometryTest.py
rename to python/unittest/geometryTest.py
diff --git a/unittest/graphTest.py b/python/unittest/graphTest.py
similarity index 100%
rename from unittest/graphTest.py
rename to python/unittest/graphTest.py
diff --git a/unittest/moleculeTest.py b/python/unittest/moleculeTest.py
similarity index 100%
rename from unittest/moleculeTest.py
rename to python/unittest/moleculeTest.py
diff --git a/unittest/oxygen.log b/python/unittest/oxygen.log
similarity index 100%
rename from unittest/oxygen.log
rename to python/unittest/oxygen.log
diff --git a/unittest/reactionTest.py b/python/unittest/reactionTest.py
similarity index 100%
rename from unittest/reactionTest.py
rename to python/unittest/reactionTest.py
diff --git a/unittest/statesTest.py b/python/unittest/statesTest.py
similarity index 100%
rename from unittest/statesTest.py
rename to python/unittest/statesTest.py
diff --git a/unittest/test.py b/python/unittest/test.py
similarity index 100%
rename from unittest/test.py
rename to python/unittest/test.py
diff --git a/unittest/thermoTest.py b/python/unittest/thermoTest.py
similarity index 100%
rename from unittest/thermoTest.py
rename to python/unittest/thermoTest.py
diff --git a/src/constants.rs b/src/constants.rs
new file mode 100644
index 0000000..17b73eb
--- /dev/null
+++ b/src/constants.rs
@@ -0,0 +1,23 @@
+//! Physical constants used throughout ChemPy.
+//! All constants are in SI units (m, s, kg, mol, etc.).
+//! Values are taken from NIST.
+
+use std::f64::consts::PI;
+
+/// The Avogadro constant (particles/mol)
+pub const NA: f64 = 6.02214179e23;
+
+/// The Boltzmann constant (J/K)
+pub const KB: f64 = 1.3806504e-23;
+
+/// The gas law constant (J/(mol·K))
+pub const R: f64 = 8.314472;
+
+/// The Planck constant (J·s)
+pub const H: f64 = 6.62606896e-34;
+
+/// The speed of light in a vacuum (m/s)
+pub const C: f64 = 299792458.0;
+
+/// pi (dimensionless)
+pub const PI_CONST: f64 = PI;
diff --git a/src/element.rs b/src/element.rs
new file mode 100644
index 0000000..a5f63ac
--- /dev/null
+++ b/src/element.rs
@@ -0,0 +1,745 @@
+use std::fmt;
+
+/// A chemical element.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct Element {
+    /// The atomic number of the element
+    pub number: u16,
+    /// The symbol used for the element
+    pub symbol: &'static str,
+    /// The IUPAC name of the element
+    pub name: &'static str,
+    /// The mass of the element in kg/mol
+    pub mass: f64,
+}
+
+impl fmt::Display for Element {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "{}", self.symbol)
+    }
+}
+
+/// Return the `Element` object with attributes defined by the given parameters.
+/// Only the parameters explicitly given will be used.
+pub fn get_element(number: u16, symbol: &str) -> Option<&'static Element> {
+    ELEMENT_LIST
+        .iter()
+        .find(|e| (number == 0 || e.number == number) && (symbol.is_empty() || e.symbol == symbol))
+}
+
+// Period 1
+pub const H: Element = Element {
+    number: 1,
+    symbol: "H",
+    name: "hydrogen",
+    mass: 0.00100794,
+};
+pub const HE: Element = Element {
+    number: 2,
+    symbol: "He",
+    name: "helium",
+    mass: 0.004002602,
+};
+
+// Period 2
+pub const LI: Element = Element {
+    number: 3,
+    symbol: "Li",
+    name: "lithium",
+    mass: 0.006941,
+};
+pub const BE: Element = Element {
+    number: 4,
+    symbol: "Be",
+    name: "beryllium",
+    mass: 0.009012182,
+};
+pub const B: Element = Element {
+    number: 5,
+    symbol: "B",
+    name: "boron",
+    mass: 0.010811,
+};
+pub const C: Element = Element {
+    number: 6,
+    symbol: "C",
+    name: "carbon",
+    mass: 0.0120107,
+};
+pub const N: Element = Element {
+    number: 7,
+    symbol: "N",
+    name: "nitrogen",
+    mass: 0.01400674,
+};
+pub const O: Element = Element {
+    number: 8,
+    symbol: "O",
+    name: "oxygen",
+    mass: 0.0159994,
+};
+pub const F: Element = Element {
+    number: 9,
+    symbol: "F",
+    name: "fluorine",
+    mass: 0.018998403,
+};
+pub const NE: Element = Element {
+    number: 10,
+    symbol: "Ne",
+    name: "neon",
+    mass: 0.0201797,
+};
+
+// Period 3
+pub const NA: Element = Element {
+    number: 11,
+    symbol: "Na",
+    name: "sodium",
+    mass: 0.022989770,
+};
+pub const MG: Element = Element {
+    number: 12,
+    symbol: "Mg",
+    name: "magnesium",
+    mass: 0.0243050,
+};
+pub const AL: Element = Element {
+    number: 13,
+    symbol: "Al",
+    name: "aluminium",
+    mass: 0.026981538,
+};
+pub const SI: Element = Element {
+    number: 14,
+    symbol: "Si",
+    name: "silicon",
+    mass: 0.0280855,
+};
+pub const P: Element = Element {
+    number: 15,
+    symbol: "P",
+    name: "phosphorus",
+    mass: 0.030973761,
+};
+pub const S: Element = Element {
+    number: 16,
+    symbol: "S",
+    name: "sulfur",
+    mass: 0.032065,
+};
+pub const CL: Element = Element {
+    number: 17,
+    symbol: "Cl",
+    name: "chlorine",
+    mass: 0.035453,
+};
+pub const AR: Element = Element {
+    number: 18,
+    symbol: "Ar",
+    name: "argon",
+    mass: 0.039348,
+};
+
+// Period 4
+pub const K: Element = Element {
+    number: 19,
+    symbol: "K",
+    name: "potassium",
+    mass: 0.0390983,
+};
+pub const CA: Element = Element {
+    number: 20,
+    symbol: "Ca",
+    name: "calcium",
+    mass: 0.040078,
+};
+pub const SC: Element = Element {
+    number: 21,
+    symbol: "Sc",
+    name: "scandium",
+    mass: 0.044955910,
+};
+pub const TI: Element = Element {
+    number: 22,
+    symbol: "Ti",
+    name: "titanium",
+    mass: 0.047867,
+};
+pub const V: Element = Element {
+    number: 23,
+    symbol: "V",
+    name: "vanadium",
+    mass: 0.0509415,
+};
+pub const CR: Element = Element {
+    number: 24,
+    symbol: "Cr",
+    name: "chromium",
+    mass: 0.0519961,
+};
+pub const MN: Element = Element {
+    number: 25,
+    symbol: "Mn",
+    name: "manganese",
+    mass: 0.054938049,
+};
+pub const FE: Element = Element {
+    number: 26,
+    symbol: "Fe",
+    name: "iron",
+    mass: 0.055845,
+};
+pub const CO: Element = Element {
+    number: 27,
+    symbol: "Co",
+    name: "cobalt",
+    mass: 0.058933200,
+};
+pub const NI: Element = Element {
+    number: 28,
+    symbol: "Ni",
+    name: "nickel",
+    mass: 0.0586934,
+};
+pub const CU: Element = Element {
+    number: 29,
+    symbol: "Cu",
+    name: "copper",
+    mass: 0.063546,
+};
+pub const ZN: Element = Element {
+    number: 30,
+    symbol: "Zn",
+    name: "zinc",
+    mass: 0.065409,
+};
+pub const GA: Element = Element {
+    number: 31,
+    symbol: "Ga",
+    name: "gallium",
+    mass: 0.069723,
+};
+pub const GE: Element = Element {
+    number: 32,
+    symbol: "Ge",
+    name: "germanium",
+    mass: 0.07264,
+};
+pub const AS: Element = Element {
+    number: 33,
+    symbol: "As",
+    name: "arsenic",
+    mass: 0.07492160,
+};
+pub const SE: Element = Element {
+    number: 34,
+    symbol: "Se",
+    name: "selenium",
+    mass: 0.07896,
+};
+pub const BR: Element = Element {
+    number: 35,
+    symbol: "Br",
+    name: "bromine",
+    mass: 0.079904,
+};
+pub const KR: Element = Element {
+    number: 36,
+    symbol: "Kr",
+    name: "krypton",
+    mass: 0.083798,
+};
+
+// Period 5
+pub const RB: Element = Element {
+    number: 37,
+    symbol: "Rb",
+    name: "rubidium",
+    mass: 0.0854678,
+};
+pub const SR: Element = Element {
+    number: 38,
+    symbol: "Sr",
+    name: "strontium",
+    mass: 0.08762,
+};
+pub const Y: Element = Element {
+    number: 39,
+    symbol: "Y",
+    name: "yttrium",
+    mass: 0.08890585,
+};
+pub const ZR: Element = Element {
+    number: 40,
+    symbol: "Zr",
+    name: "zirconium",
+    mass: 0.091224,
+};
+pub const NB: Element = Element {
+    number: 41,
+    symbol: "Nb",
+    name: "niobium",
+    mass: 0.09290638,
+};
+pub const MO: Element = Element {
+    number: 42,
+    symbol: "Mo",
+    name: "molybdenum",
+    mass: 0.09594,
+};
+pub const TC: Element = Element {
+    number: 43,
+    symbol: "Tc",
+    name: "technetium",
+    mass: 0.098,
+};
+pub const RU: Element = Element {
+    number: 44,
+    symbol: "Ru",
+    name: "ruthenium",
+    mass: 0.10107,
+};
+pub const RH: Element = Element {
+    number: 45,
+    symbol: "Rh",
+    name: "rhodium",
+    mass: 0.10290550,
+};
+pub const PD: Element = Element {
+    number: 46,
+    symbol: "Pd",
+    name: "palladium",
+    mass: 0.10642,
+};
+pub const AG: Element = Element {
+    number: 47,
+    symbol: "Ag",
+    name: "silver",
+    mass: 0.1078682,
+};
+pub const CD: Element = Element {
+    number: 48,
+    symbol: "Cd",
+    name: "cadmium",
+    mass: 0.112411,
+};
+pub const IN: Element = Element {
+    number: 49,
+    symbol: "In",
+    name: "indium",
+    mass: 0.114818,
+};
+pub const SN: Element = Element {
+    number: 50,
+    symbol: "Sn",
+    name: "tin",
+    mass: 0.118710,
+};
+pub const SB: Element = Element {
+    number: 51,
+    symbol: "Sb",
+    name: "antimony",
+    mass: 0.121760,
+};
+pub const TE: Element = Element {
+    number: 52,
+    symbol: "Te",
+    name: "tellurium",
+    mass: 0.12760,
+};
+pub const I: Element = Element {
+    number: 53,
+    symbol: "I",
+    name: "iodine",
+    mass: 0.12690447,
+};
+pub const XE: Element = Element {
+    number: 54,
+    symbol: "Xe",
+    name: "xenon",
+    mass: 0.131293,
+};
+
+// Period 6
+pub const CS: Element = Element {
+    number: 55,
+    symbol: "Cs",
+    name: "caesium",
+    mass: 0.13290545,
+};
+pub const BA: Element = Element {
+    number: 56,
+    symbol: "Ba",
+    name: "barium",
+    mass: 0.137327,
+};
+pub const LA: Element = Element {
+    number: 57,
+    symbol: "La",
+    name: "lanthanum",
+    mass: 0.1389055,
+};
+pub const CE: Element = Element {
+    number: 58,
+    symbol: "Ce",
+    name: "cerium",
+    mass: 0.140116,
+};
+pub const PR: Element = Element {
+    number: 59,
+    symbol: "Pr",
+    name: "praesodymium",
+    mass: 0.14090765,
+};
+pub const ND: Element = Element {
+    number: 60,
+    symbol: "Nd",
+    name: "neodymium",
+    mass: 0.14424,
+};
+pub const PM: Element = Element {
+    number: 61,
+    symbol: "Pm",
+    name: "promethium",
+    mass: 0.145,
+};
+pub const SM: Element = Element {
+    number: 62,
+    symbol: "Sm",
+    name: "samarium",
+    mass: 0.15036,
+};
+pub const EU: Element = Element {
+    number: 63,
+    symbol: "Eu",
+    name: "europium",
+    mass: 0.151964,
+};
+pub const GD: Element = Element {
+    number: 64,
+    symbol: "Gd",
+    name: "gadolinium",
+    mass: 0.15725,
+};
+pub const TB: Element = Element {
+    number: 65,
+    symbol: "Tb",
+    name: "terbium",
+    mass: 0.15892534,
+};
+pub const DY: Element = Element {
+    number: 66,
+    symbol: "Dy",
+    name: "dysprosium",
+    mass: 0.162500,
+};
+pub const HO: Element = Element {
+    number: 67,
+    symbol: "Ho",
+    name: "holmium",
+    mass: 0.16493032,
+};
+pub const ER: Element = Element {
+    number: 68,
+    symbol: "Er",
+    name: "erbium",
+    mass: 0.167259,
+};
+pub const TM: Element = Element {
+    number: 69,
+    symbol: "Tm",
+    name: "thulium",
+    mass: 0.16893421,
+};
+pub const YB: Element = Element {
+    number: 70,
+    symbol: "Yb",
+    name: "ytterbium",
+    mass: 0.17304,
+};
+pub const LU: Element = Element {
+    number: 71,
+    symbol: "Lu",
+    name: "lutetium",
+    mass: 0.174967,
+};
+pub const HF: Element = Element {
+    number: 72,
+    symbol: "Hf",
+    name: "hafnium",
+    mass: 0.17849,
+};
+pub const TA: Element = Element {
+    number: 73,
+    symbol: "Ta",
+    name: "tantalum",
+    mass: 0.1809479,
+};
+pub const W: Element = Element {
+    number: 74,
+    symbol: "W",
+    name: "tungsten",
+    mass: 0.18384,
+};
+pub const RE: Element = Element {
+    number: 75,
+    symbol: "Re",
+    name: "rhenium",
+    mass: 0.186207,
+};
+pub const OS: Element = Element {
+    number: 76,
+    symbol: "Os",
+    name: "osmium",
+    mass: 0.19023,
+};
+pub const IR: Element = Element {
+    number: 77,
+    symbol: "Ir",
+    name: "iridium",
+    mass: 0.192217,
+};
+pub const PT: Element = Element {
+    number: 78,
+    symbol: "Pt",
+    name: "platinum",
+    mass: 0.195078,
+};
+pub const AU: Element = Element {
+    number: 79,
+    symbol: "Au",
+    name: "gold",
+    mass: 0.19696655,
+};
+pub const HG: Element = Element {
+    number: 80,
+    symbol: "Hg",
+    name: "mercury",
+    mass: 0.20059,
+};
+pub const TL: Element = Element {
+    number: 81,
+    symbol: "Tl",
+    name: "thallium",
+    mass: 0.2043833,
+};
+pub const PB: Element = Element {
+    number: 82,
+    symbol: "Pb",
+    name: "lead",
+    mass: 0.2072,
+};
+pub const BI: Element = Element {
+    number: 83,
+    symbol: "Bi",
+    name: "bismuth",
+    mass: 0.20898038,
+};
+pub const PO: Element = Element {
+    number: 84,
+    symbol: "Po",
+    name: "polonium",
+    mass: 0.209,
+};
+pub const AT: Element = Element {
+    number: 85,
+    symbol: "At",
+    name: "astatine",
+    mass: 0.210,
+};
+pub const RN: Element = Element {
+    number: 86,
+    symbol: "Rn",
+    name: "radon",
+    mass: 0.222,
+};
+
+// Period 7
+pub const FR: Element = Element {
+    number: 87,
+    symbol: "Fr",
+    name: "francium",
+    mass: 0.223,
+};
+pub const RA: Element = Element {
+    number: 88,
+    symbol: "Ra",
+    name: "radium",
+    mass: 0.226,
+};
+pub const AC: Element = Element {
+    number: 89,
+    symbol: "Ac",
+    name: "actinum",
+    mass: 0.227,
+};
+pub const TH: Element = Element {
+    number: 90,
+    symbol: "Th",
+    name: "thorium",
+    mass: 0.2320381,
+};
+pub const PA: Element = Element {
+    number: 91,
+    symbol: "Pa",
+    name: "protactinum",
+    mass: 0.23103588,
+};
+pub const U: Element = Element {
+    number: 92,
+    symbol: "U",
+    name: "uranium",
+    mass: 0.23802891,
+};
+pub const NP: Element = Element {
+    number: 93,
+    symbol: "Np",
+    name: "neptunium",
+    mass: 0.237,
+};
+pub const PU: Element = Element {
+    number: 94,
+    symbol: "Pu",
+    name: "plutonium",
+    mass: 0.244,
+};
+pub const AM: Element = Element {
+    number: 95,
+    symbol: "Am",
+    name: "americium",
+    mass: 0.243,
+};
+pub const CM: Element = Element {
+    number: 96,
+    symbol: "Cm",
+    name: "curium",
+    mass: 0.247,
+};
+pub const BK: Element = Element {
+    number: 97,
+    symbol: "Bk",
+    name: "berkelium",
+    mass: 0.247,
+};
+pub const CF: Element = Element {
+    number: 98,
+    symbol: "Cf",
+    name: "californium",
+    mass: 0.251,
+};
+pub const ES: Element = Element {
+    number: 99,
+    symbol: "Es",
+    name: "einsteinium",
+    mass: 0.252,
+};
+pub const FM: Element = Element {
+    number: 100,
+    symbol: "Fm",
+    name: "fermium",
+    mass: 0.257,
+};
+pub const MD: Element = Element {
+    number: 101,
+    symbol: "Md",
+    name: "mendelevium",
+    mass: 0.258,
+};
+pub const NO: Element = Element {
+    number: 102,
+    symbol: "No",
+    name: "nobelium",
+    mass: 0.259,
+};
+pub const LR: Element = Element {
+    number: 103,
+    symbol: "Lr",
+    name: "lawrencium",
+    mass: 0.262,
+};
+pub const RF: Element = Element {
+    number: 104,
+    symbol: "Rf",
+    name: "rutherfordium",
+    mass: 0.261,
+};
+pub const DB: Element = Element {
+    number: 105,
+    symbol: "Db",
+    name: "dubnium",
+    mass: 0.262,
+};
+pub const SG: Element = Element {
+    number: 106,
+    symbol: "Sg",
+    name: "seaborgium",
+    mass: 0.266,
+};
+pub const BH: Element = Element {
+    number: 107,
+    symbol: "Bh",
+    name: "bohrium",
+    mass: 0.264,
+};
+pub const HS: Element = Element {
+    number: 108,
+    symbol: "Hs",
+    name: "hassium",
+    mass: 0.277,
+};
+pub const MT: Element = Element {
+    number: 109,
+    symbol: "Mt",
+    name: "meitnerium",
+    mass: 0.268,
+};
+pub const DS: Element = Element {
+    number: 110,
+    symbol: "Ds",
+    name: "darmstadtium",
+    mass: 0.281,
+};
+pub const RG: Element = Element {
+    number: 111,
+    symbol: "Rg",
+    name: "roentgenium",
+    mass: 0.272,
+};
+pub const CN: Element = Element {
+    number: 112,
+    symbol: "Cn",
+    name: "copernicum",
+    mass: 0.285,
+};
+
+pub const ELEMENT_LIST: [Element; 112] = [
+    H, HE, LI, BE, B, C, N, O, F, NE, NA, MG, AL, SI, P, S, CL, AR, K, CA, SC, TI, V, CR, MN, FE,
+    CO, NI, CU, ZN, GA, GE, AS, SE, BR, KR, RB, SR, Y, ZR, NB, MO, TC, RU, RH, PD, AG, CD, IN, SN,
+    SB, TE, I, XE, CS, BA, LA, CE, PR, ND, PM, SM, EU, GD, TB, DY, HO, ER, TM, YB, LU, HF, TA, W,
+    RE, OS, IR, PT, AU, HG, TL, PB, BI, PO, AT, RN, FR, RA, AC, TH, PA, U, NP, PU, AM, CM, BK, CF,
+    ES, FM, MD, NO, LR, RF, DB, SG, BH, HS, MT, DS, RG, CN,
+];
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_get_element() {
+        let carbon = get_element(6, "").unwrap();
+        assert_eq!(carbon.symbol, "C");
+        assert_eq!(carbon.name, "carbon");
+
+        let hydrogen = get_element(0, "H").unwrap();
+        assert_eq!(hydrogen.number, 1);
+
+        let none = get_element(999, "");
+        assert!(none.is_none());
+    }
+
+    #[test]
+    fn test_display() {
+        assert_eq!(format!("{}", C), "C");
+    }
+}
diff --git a/src/graph.rs b/src/graph.rs
new file mode 100644
index 0000000..9eabcb9
--- /dev/null
+++ b/src/graph.rs
@@ -0,0 +1,592 @@
+use std::collections::HashMap;
+
+/// A base trait for vertices in a graph.
+pub trait Vertex: Clone {
+    fn equivalent(&self, _other: &Self) -> bool {
+        true
+    }
+
+    fn is_specific_case_of(&self, _other: &Self) -> bool {
+        true
+    }
+}
+
+/// A base trait for edges in a graph.
+pub trait Edge: Clone {
+    fn equivalent(&self, _other: &Self) -> bool {
+        true
+    }
+
+    fn is_specific_case_of(&self, _other: &Self) -> bool {
+        true
+    }
+}
+
+pub trait HasConnectivity {
+    fn connectivity1(&self) -> i32;
+    fn set_connectivity1(&mut self, value: i32);
+    fn connectivity2(&self) -> i32;
+    fn set_connectivity2(&mut self, value: i32);
+    fn connectivity3(&self) -> i32;
+    fn set_connectivity3(&mut self, value: i32);
+}
+
+/// A simple implementation of a vertex for generic graphs.
+#[derive(Debug, Clone, Default)]
+pub struct BaseVertex {
+    pub connectivity1: i32,
+    pub connectivity2: i32,
+    pub connectivity3: i32,
+    pub sorting_label: i32,
+}
+
+impl Vertex for BaseVertex {}
+
+impl HasConnectivity for BaseVertex {
+    fn connectivity1(&self) -> i32 {
+        self.connectivity1
+    }
+    fn set_connectivity1(&mut self, value: i32) {
+        self.connectivity1 = value;
+    }
+    fn connectivity2(&self) -> i32 {
+        self.connectivity2
+    }
+    fn set_connectivity2(&mut self, value: i32) {
+        self.connectivity2 = value;
+    }
+    fn connectivity3(&self) -> i32 {
+        self.connectivity3
+    }
+    fn set_connectivity3(&mut self, value: i32) {
+        self.connectivity3 = value;
+    }
+}
+
+/// A simple implementation of an edge for generic graphs.
+#[derive(Debug, Clone, Default)]
+pub struct BaseEdge {}
+
+impl Edge for BaseEdge {}
+
+/// A graph data type.
+#[derive(Debug, Clone)]
+pub struct Graph<V: Vertex, E: Edge> {
+    pub vertices: Vec<V>,
+    pub edges: Vec<HashMap<usize, E>>,
+}
+
+impl<V: Vertex, E: Edge> Graph<V, E> {
+    pub fn new() -> Self {
+        Graph {
+            vertices: Vec::new(),
+            edges: Vec::new(),
+        }
+    }
+
+    pub fn add_vertex(&mut self, vertex: V) -> usize {
+        let index = self.vertices.len();
+        self.vertices.push(vertex);
+        self.edges.push(HashMap::new());
+        index
+    }
+
+    pub fn add_edge(&mut self, v1: usize, v2: usize, edge: E) {
+        self.edges[v1].insert(v2, edge.clone());
+        self.edges[v2].insert(v1, edge);
+    }
+
+    pub fn get_edge(&self, v1: usize, v2: usize) -> Option<&E> {
+        self.edges.get(v1)?.get(&v2)
+    }
+
+    pub fn has_edge(&self, v1: usize, v2: usize) -> bool {
+        self.edges.get(v1).is_some_and(|adj| adj.contains_key(&v2))
+    }
+
+    pub fn remove_vertex(&mut self, index: usize) {
+        if index >= self.vertices.len() {
+            return;
+        }
+
+        // Remove all edges connected to this vertex
+        self.edges.remove(index);
+        self.vertices.remove(index);
+
+        // Update remaining edges to reflect new indices
+        for adj in self.edges.iter_mut() {
+            let mut new_adj = HashMap::new();
+            for (&neighbor, edge) in adj.iter() {
+                if neighbor == index {
+                    continue;
+                }
+                let new_neighbor = if neighbor > index {
+                    neighbor - 1
+                } else {
+                    neighbor
+                };
+                new_adj.insert(new_neighbor, edge.clone());
+            }
+            *adj = new_adj;
+        }
+    }
+
+    pub fn remove_edge(&mut self, v1: usize, v2: usize) {
+        if let Some(adj1) = self.edges.get_mut(v1) {
+            adj1.remove(&v2);
+        }
+        if let Some(adj2) = self.edges.get_mut(v2) {
+            adj2.remove(&v1);
+        }
+    }
+
+    pub fn update_connectivity_values(&mut self)
+    where
+        V: HasConnectivity,
+    {
+        for i in 0..self.vertices.len() {
+            self.vertices[i].set_connectivity1(self.edges[i].len() as i32);
+        }
+
+        for i in 0..self.vertices.len() {
+            let mut cv2 = 0;
+            for &neighbor in self.edges[i].keys() {
+                cv2 += self.vertices[neighbor].connectivity1();
+            }
+            self.vertices[i].set_connectivity2(cv2);
+        }
+
+        for i in 0..self.vertices.len() {
+            let mut cv3 = 0;
+            for &neighbor in self.edges[i].keys() {
+                cv3 += self.vertices[neighbor].connectivity2();
+            }
+            self.vertices[i].set_connectivity3(cv3);
+        }
+    }
+
+    pub fn split(&self) -> Vec<Graph<V, E>> {
+        let mut components = Vec::new();
+        let mut visited = vec![false; self.vertices.len()];
+
+        for i in 0..self.vertices.len() {
+            if !visited[i] {
+                let mut component_indices = Vec::new();
+                let mut stack = vec![i];
+                visited[i] = true;
+
+                while let Some(u) = stack.pop() {
+                    component_indices.push(u);
+                    for &v in self.edges[u].keys() {
+                        if !visited[v] {
+                            visited[v] = true;
+                            stack.push(v);
+                        }
+                    }
+                }
+
+                // Sort indices to maintain order and help with mapping
+                component_indices.sort();
+                let mut new_graph = Graph::new();
+                let mut old_to_new = HashMap::new();
+
+                for &old_idx in &component_indices {
+                    let new_idx = new_graph.add_vertex(self.vertices[old_idx].clone());
+                    old_to_new.insert(old_idx, new_idx);
+                }
+
+                for &old_u in &component_indices {
+                    for (&old_v, edge) in &self.edges[old_u] {
+                        if old_u < old_v {
+                            new_graph.add_edge(
+                                *old_to_new.get(&old_u).unwrap(),
+                                *old_to_new.get(&old_v).unwrap(),
+                                edge.clone(),
+                            );
+                        }
+                    }
+                }
+                components.push(new_graph);
+            }
+        }
+        components
+    }
+
+    pub fn merge(&self, other: &Graph<V, E>) -> Graph<V, E> {
+        let mut new_graph = self.clone();
+        let mut old_to_new = HashMap::new();
+
+        for vertex in &other.vertices {
+            let new_idx = new_graph.add_vertex(vertex.clone());
+            old_to_new.insert(old_to_new.len(), new_idx);
+        }
+
+        for (u_idx, adj) in other.edges.iter().enumerate() {
+            for (&v_idx, edge) in adj {
+                if u_idx < v_idx {
+                    new_graph.add_edge(
+                        *old_to_new.get(&u_idx).unwrap(),
+                        *old_to_new.get(&v_idx).unwrap(),
+                        edge.clone(),
+                    );
+                }
+            }
+        }
+        new_graph
+    }
+
+    pub fn is_cyclic(&self) -> bool {
+        let mut visited = vec![false; self.vertices.len()];
+        for i in 0..self.vertices.len() {
+            if !visited[i] && self.has_cycle_from(i, None, &mut visited) {
+                return true;
+            }
+        }
+        false
+    }
+
+    fn has_cycle_from(&self, u: usize, parent: Option<usize>, visited: &mut Vec<bool>) -> bool {
+        visited[u] = true;
+        for &v in self.edges[u].keys() {
+            if Some(v) == parent {
+                continue;
+            }
+            if visited[v] || self.has_cycle_from(v, Some(u), visited) {
+                return true;
+            }
+        }
+        false
+    }
+
+    pub fn is_vertex_in_cycle(&self, u: usize) -> bool {
+        // A vertex is in a cycle if it can reach itself without using the same edge twice
+        for &v in self.edges[u].keys() {
+            if self.can_reach(v, u, Some(u)) {
+                return true;
+            }
+        }
+        false
+    }
+
+    fn can_reach(&self, start: usize, target: usize, forbidden_parent: Option<usize>) -> bool {
+        let mut visited = vec![false; self.vertices.len()];
+        if let Some(p) = forbidden_parent {
+            visited[p] = true;
+        }
+        let mut stack = vec![start];
+        visited[start] = true;
+
+        while let Some(u) = stack.pop() {
+            if u == target {
+                return true;
+            }
+            for &v in self.edges[u].keys() {
+                if !visited[v] {
+                    visited[v] = true;
+                    stack.push(v);
+                }
+            }
+        }
+        false
+    }
+
+    pub fn is_isomorphic(&self, other: &Graph<V, E>) -> bool {
+        if self.vertices.len() != other.vertices.len() {
+            return false;
+        }
+        if self.vertices.is_empty() {
+            return true;
+        }
+
+        let mut mapping = HashMap::new();
+        let mut reverse_mapping = HashMap::new();
+        self.vf2_match(other, &mut mapping, &mut reverse_mapping, 0, false)
+    }
+
+    pub fn is_subgraph_isomorphic(&self, other: &Graph<V, E>) -> bool {
+        if self.vertices.len() < other.vertices.len() {
+            return false;
+        }
+        if other.vertices.is_empty() {
+            return true;
+        }
+
+        let mut mapping = HashMap::new();
+        let mut reverse_mapping = HashMap::new();
+        // VF2 for subgraph: swap self and other?
+        // Actually, Python's self.isSubgraphIsomorphic(other) checks if 'other' is in 'self'
+        other.vf2_match(self, &mut mapping, &mut reverse_mapping, 0, true)
+    }
+
+    pub fn find_subgraph_isomorphisms(&self, other: &Graph<V, E>) -> Vec<HashMap<usize, usize>> {
+        let mut mappings = Vec::new();
+        if self.vertices.len() < other.vertices.len() {
+            return mappings;
+        }
+        let mut mapping = HashMap::new();
+        let mut reverse_mapping = HashMap::new();
+        other.vf2_all_matches(
+            self,
+            &mut mapping,
+            &mut reverse_mapping,
+            0,
+            true,
+            &mut mappings,
+        );
+        mappings
+    }
+
+    fn vf2_match(
+        &self,
+        other: &Graph<V, E>,
+        mapping: &mut HashMap<usize, usize>,
+        reverse_mapping: &mut HashMap<usize, usize>,
+        depth: usize,
+        subgraph: bool,
+    ) -> bool {
+        if depth == self.vertices.len() {
+            return true;
+        }
+
+        let v1 = depth;
+        for v2 in 0..other.vertices.len() {
+            if !reverse_mapping.contains_key(&v2)
+                && self.is_feasible(v1, v2, other, mapping, subgraph)
+            {
+                mapping.insert(v1, v2);
+                reverse_mapping.insert(v2, v1);
+
+                if self.vf2_match(other, mapping, reverse_mapping, depth + 1, subgraph) {
+                    return true;
+                }
+
+                mapping.remove(&v1);
+                reverse_mapping.remove(&v2);
+            }
+        }
+        false
+    }
+
+    fn vf2_all_matches(
+        &self,
+        other: &Graph<V, E>,
+        mapping: &mut HashMap<usize, usize>,
+        reverse_mapping: &mut HashMap<usize, usize>,
+        depth: usize,
+        subgraph: bool,
+        mappings: &mut Vec<HashMap<usize, usize>>,
+    ) {
+        if depth == self.vertices.len() {
+            mappings.push(mapping.clone());
+            return;
+        }
+
+        let v1 = depth;
+        for v2 in 0..other.vertices.len() {
+            if !reverse_mapping.contains_key(&v2)
+                && self.is_feasible(v1, v2, other, mapping, subgraph)
+            {
+                mapping.insert(v1, v2);
+                reverse_mapping.insert(v2, v1);
+
+                self.vf2_all_matches(
+                    other,
+                    mapping,
+                    reverse_mapping,
+                    depth + 1,
+                    subgraph,
+                    mappings,
+                );
+
+                mapping.remove(&v1);
+                reverse_mapping.remove(&v2);
+            }
+        }
+    }
+
+    fn is_feasible(
+        &self,
+        v1: usize,
+        v2: usize,
+        other: &Graph<V, E>,
+        mapping: &HashMap<usize, usize>,
+        subgraph: bool,
+    ) -> bool {
+        // Semantic check
+        if !self.vertices[v1].equivalent(&other.vertices[v2]) {
+            return false;
+        }
+
+        // Structural check
+        for (&neighbor1, edge1) in &self.edges[v1] {
+            if let Some(&neighbor2_mapped) = mapping.get(&neighbor1) {
+                if let Some(edge2) = other.get_edge(v2, neighbor2_mapped) {
+                    if !edge1.equivalent(edge2) {
+                        return false;
+                    }
+                } else {
+                    return false;
+                }
+            }
+        }
+
+        // Degree check
+        if subgraph {
+            if self.edges[v1].len() > other.edges[v2].len() {
+                return false;
+            }
+        } else if self.edges[v1].len() != other.edges[v2].len() {
+            return false;
+        }
+
+        true
+    }
+}
+
+impl<V: Vertex, E: Edge> Default for Graph<V, E> {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_graph_basic() {
+        let mut g = Graph::<BaseVertex, BaseEdge>::new();
+        let v1 = g.add_vertex(BaseVertex::default());
+        let v2 = g.add_vertex(BaseVertex::default());
+        g.add_edge(v1, v2, BaseEdge::default());
+
+        assert_eq!(g.vertices.len(), 2);
+        assert!(g.has_edge(v1, v2));
+        assert!(g.has_edge(v2, v1));
+        assert!(g.get_edge(v1, v2).is_some());
+    }
+
+    #[test]
+    fn test_remove_vertex() {
+        let mut g = Graph::<BaseVertex, BaseEdge>::new();
+        let v1 = g.add_vertex(BaseVertex::default());
+        let v2 = g.add_vertex(BaseVertex::default());
+        let v3 = g.add_vertex(BaseVertex::default());
+        g.add_edge(v1, v2, BaseEdge::default());
+        g.add_edge(v2, v3, BaseEdge::default());
+
+        g.remove_vertex(v1); // v1 is gone, v2 becomes 0, v3 becomes 1
+        assert_eq!(g.vertices.len(), 2);
+        assert!(g.has_edge(0, 1));
+    }
+
+    #[test]
+    fn test_isomorphism() {
+        let mut g1 = Graph::<BaseVertex, BaseEdge>::new();
+        let v1 = g1.add_vertex(BaseVertex::default());
+        let v2 = g1.add_vertex(BaseVertex::default());
+        g1.add_edge(v1, v2, BaseEdge::default());
+
+        let mut g2 = Graph::<BaseVertex, BaseEdge>::new();
+        let v3 = g2.add_vertex(BaseVertex::default());
+        let v4 = g2.add_vertex(BaseVertex::default());
+        g2.add_edge(v3, v4, BaseEdge::default());
+
+        assert!(g1.is_isomorphic(&g2));
+
+        let mut g3 = Graph::<BaseVertex, BaseEdge>::new();
+        g3.add_vertex(BaseVertex::default());
+        g3.add_vertex(BaseVertex::default());
+        // No edge
+        assert!(!g1.is_isomorphic(&g3));
+    }
+
+    #[test]
+    fn test_connectivity_values() {
+        // 0-1-2-3-4
+        //   |
+        //   5
+        let mut g = Graph::<BaseVertex, BaseEdge>::new();
+        let vertices: Vec<usize> = (0..6)
+            .map(|_| g.add_vertex(BaseVertex::default()))
+            .collect();
+        g.add_edge(vertices[0], vertices[1], BaseEdge::default());
+        g.add_edge(vertices[1], vertices[2], BaseEdge::default());
+        g.add_edge(vertices[2], vertices[3], BaseEdge::default());
+        g.add_edge(vertices[3], vertices[4], BaseEdge::default());
+        g.add_edge(vertices[1], vertices[5], BaseEdge::default());
+
+        g.update_connectivity_values();
+
+        let expected_cv1 = [1, 3, 2, 2, 1, 1];
+        let expected_cv2 = [3, 4, 5, 3, 2, 3];
+        let expected_cv3 = [4, 11, 7, 7, 3, 4];
+
+        for i in 0..6 {
+            assert_eq!(g.vertices[i].connectivity1, expected_cv1[i]);
+            assert_eq!(g.vertices[i].connectivity2, expected_cv2[i]);
+            assert_eq!(g.vertices[i].connectivity3, expected_cv3[i]);
+        }
+    }
+
+    #[test]
+    fn test_split() {
+        let mut g = Graph::<BaseVertex, BaseEdge>::new();
+        let v: Vec<usize> = (0..6)
+            .map(|_| g.add_vertex(BaseVertex::default()))
+            .collect();
+        g.add_edge(v[0], v[1], BaseEdge::default());
+        g.add_edge(v[1], v[2], BaseEdge::default());
+        g.add_edge(v[2], v[3], BaseEdge::default());
+        g.add_edge(v[4], v[5], BaseEdge::default());
+
+        let components = g.split();
+        assert_eq!(components.len(), 2);
+        let lens: Vec<usize> = components.iter().map(|c| c.vertices.len()).collect();
+        assert!(lens.contains(&4));
+        assert!(lens.contains(&2));
+    }
+
+    #[test]
+    fn test_merge() {
+        let mut g1 = Graph::<BaseVertex, BaseEdge>::new();
+        let v1: Vec<usize> = (0..4)
+            .map(|_| g1.add_vertex(BaseVertex::default()))
+            .collect();
+        g1.add_edge(v1[0], v1[1], BaseEdge::default());
+
+        let mut g2 = Graph::<BaseVertex, BaseEdge>::new();
+        let v2: Vec<usize> = (0..3)
+            .map(|_| g2.add_vertex(BaseVertex::default()))
+            .collect();
+        g2.add_edge(v2[0], v2[1], BaseEdge::default());
+
+        let g = g1.merge(&g2);
+        assert_eq!(g.vertices.len(), 7);
+    }
+
+    #[test]
+    fn test_subgraph_isomorphism() {
+        let mut g1 = Graph::<BaseVertex, BaseEdge>::new();
+        let v1: Vec<usize> = (0..6)
+            .map(|_| g1.add_vertex(BaseVertex::default()))
+            .collect();
+        // Path graph 0-1-2-3-4-5
+        for i in 0..5 {
+            g1.add_edge(v1[i], v1[i + 1], BaseEdge::default());
+        }
+
+        let mut g2 = Graph::<BaseVertex, BaseEdge>::new();
+        let v2: Vec<usize> = (0..2)
+            .map(|_| g2.add_vertex(BaseVertex::default()))
+            .collect();
+        g2.add_edge(v2[0], v2[1], BaseEdge::default());
+
+        assert!(g1.is_subgraph_isomorphic(&g2));
+        let mappings = g1.find_subgraph_isomorphisms(&g2);
+        // A single edge (g2) can be mapped to any of the 5 edges in the path (g1)
+        // Each edge can be mapped in 2 directions.
+        // 5 edges * 2 directions = 10 mappings.
+        assert_eq!(mappings.len(), 10);
+    }
+}
diff --git a/src/kinetics.rs b/src/kinetics.rs
new file mode 100644
index 0000000..160f72f
--- /dev/null
+++ b/src/kinetics.rs
@@ -0,0 +1,56 @@
+use crate::constants;
+
+pub trait KineticsModel {
+    fn get_rate_coefficient(&self, t: f64, p: f64) -> f64;
+
+    fn is_temperature_valid(&self, t: f64, t_min: f64, t_max: f64) -> bool {
+        t >= t_min && t <= t_max
+    }
+
+    fn is_pressure_valid(&self, p: f64, p_min: f64, p_max: f64) -> bool {
+        p >= p_min && p <= p_max
+    }
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct ArrheniusModel {
+    pub a: f64,
+    pub n: f64,
+    pub ea: f64,
+    pub t0: f64,
+    pub t_min: f64,
+    pub t_max: f64,
+}
+
+impl ArrheniusModel {
+    pub fn new(a: f64, n: f64, ea: f64, t0: f64) -> Self {
+        ArrheniusModel {
+            a,
+            n,
+            ea,
+            t0,
+            t_min: 0.0,
+            t_max: 1.0e10,
+        }
+    }
+}
+
+impl KineticsModel for ArrheniusModel {
+    fn get_rate_coefficient(&self, t: f64, _p: f64) -> f64 {
+        self.a * (t / self.t0).powf(self.n) * (-self.ea / constants::R / t).exp()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_arrhenius_rate_coefficient() {
+        let model = ArrheniusModel::new(1.0e10, 0.0, 50000.0, 1.0);
+        let t = 1000.0;
+        let k_expected = 1.0e10 * (-50000.0 / (constants::R * t)).exp();
+        let k_actual = model.get_rate_coefficient(t, 1.0e5);
+        assert!((k_actual - k_expected).abs() < 1e-10 * k_expected);
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
new file mode 100644
index 0000000..e5ce860
--- /dev/null
+++ b/src/lib.rs
@@ -0,0 +1,9 @@
+pub mod constants;
+pub mod element;
+pub mod graph;
+pub mod kinetics;
+pub mod molecule;
+pub mod reaction;
+pub mod species;
+pub mod states;
+pub mod thermo;
diff --git a/src/molecule.rs b/src/molecule.rs
new file mode 100644
index 0000000..5ef90df
--- /dev/null
+++ b/src/molecule.rs
@@ -0,0 +1,383 @@
+use crate::element::Element;
+use crate::graph::{Edge, Graph, Vertex};
+use std::fmt;
+
+/// An atom.
+#[derive(Debug, Clone, PartialEq)]
+pub struct Atom {
+    pub element: &'static Element,
+    pub radical_electrons: i16,
+    pub spin_multiplicity: i16,
+    pub implicit_hydrogens: i16,
+    pub charge: i16,
+    pub label: String,
+    pub connectivity1: i32,
+    pub connectivity2: i32,
+    pub connectivity3: i32,
+}
+
+impl Atom {
+    pub fn new(element: &'static Element) -> Self {
+        Atom {
+            element,
+            radical_electrons: 0,
+            spin_multiplicity: 1,
+            implicit_hydrogens: 0,
+            charge: 0,
+            label: String::new(),
+            connectivity1: 0,
+            connectivity2: 0,
+            connectivity3: 0,
+        }
+    }
+}
+
+impl crate::graph::HasConnectivity for Atom {
+    fn connectivity1(&self) -> i32 {
+        self.connectivity1
+    }
+    fn set_connectivity1(&mut self, value: i32) {
+        self.connectivity1 = value;
+    }
+    fn connectivity2(&self) -> i32 {
+        self.connectivity2
+    }
+    fn set_connectivity2(&mut self, value: i32) {
+        self.connectivity2 = value;
+    }
+    fn connectivity3(&self) -> i32 {
+        self.connectivity3
+    }
+    fn set_connectivity3(&mut self, value: i32) {
+        self.connectivity3 = value;
+    }
+}
+
+impl Vertex for Atom {
+    fn equivalent(&self, other: &Self) -> bool {
+        self.element == other.element
+            && self.radical_electrons == other.radical_electrons
+            && self.spin_multiplicity == other.spin_multiplicity
+            && self.implicit_hydrogens == other.implicit_hydrogens
+            && self.charge == other.charge
+    }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum BondOrder {
+    Single,
+    Double,
+    Triple,
+    Benzene,
+}
+
+impl fmt::Display for BondOrder {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            BondOrder::Single => write!(f, "S"),
+            BondOrder::Double => write!(f, "D"),
+            BondOrder::Triple => write!(f, "T"),
+            BondOrder::Benzene => write!(f, "B"),
+        }
+    }
+}
+
+/// A chemical bond.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub struct Bond {
+    pub order: BondOrder,
+}
+
+impl Bond {
+    pub fn new(order: BondOrder) -> Self {
+        Bond { order }
+    }
+}
+
+impl Edge for Bond {
+    fn equivalent(&self, other: &Self) -> bool {
+        self.order == other.order
+    }
+}
+
+/// A representation of a molecular structure.
+#[derive(Debug, Clone, Default)]
+pub struct Molecule {
+    pub graph: Graph<Atom, Bond>,
+}
+
+impl Molecule {
+    pub fn new() -> Self {
+        Molecule {
+            graph: Graph::new(),
+        }
+    }
+
+    pub fn add_atom(&mut self, atom: Atom) -> usize {
+        self.graph.add_vertex(atom)
+    }
+
+    pub fn add_bond(&mut self, v1: usize, v2: usize, bond: Bond) {
+        self.graph.add_edge(v1, v2, bond);
+    }
+
+    pub fn get_atom(&self, index: usize) -> Option<&Atom> {
+        self.graph.vertices.get(index)
+    }
+
+    pub fn get_bond(&self, v1: usize, v2: usize) -> Option<&Bond> {
+        self.graph.get_edge(v1, v2)
+    }
+
+    pub fn to_adjacency_list(&self) -> String {
+        let mut result = String::new();
+        for (i, atom) in self.graph.vertices.iter().enumerate() {
+            let mut line = format!(
+                "{} {} {}",
+                i + 1,
+                atom.element.symbol,
+                atom.radical_electrons
+            );
+            let mut neighbors: Vec<_> = self.graph.edges[i].keys().collect();
+            neighbors.sort();
+            for &neighbor in neighbors {
+                let bond = self.get_bond(i, neighbor).unwrap();
+                line.push_str(&format!(" {{{},{}}}", neighbor + 1, bond.order));
+            }
+            result.push_str(&line);
+            result.push('\n');
+        }
+        result
+    }
+
+    pub fn from_adjacency_list(&mut self, adj_list: &str) {
+        use crate::element::get_element;
+        self.graph = Graph::new();
+        let lines: Vec<&str> = adj_list.trim().lines().collect();
+        let mut bond_info = Vec::new();
+
+        for line in &lines {
+            let parts: Vec<&str> = line.split_whitespace().collect();
+            if parts.len() < 3 {
+                continue;
+            }
+            let symbol = parts[1];
+            let radical = parts[2].parse::<i16>().unwrap_or(0);
+            let element = get_element(0, symbol).expect("Unknown element");
+            let mut atom = Atom::new(element);
+            atom.radical_electrons = radical;
+            self.add_atom(atom);
+
+            // Collect bonds to add after all atoms are created
+            for &part in &parts[3..] {
+                if part.starts_with('{') && part.ends_with('}') {
+                    let inner = &part[1..part.len() - 1];
+                    let bond_parts: Vec<&str> = inner.split(',').collect();
+                    if bond_parts.len() == 2 {
+                        let target_idx = bond_parts[0].parse::<usize>().unwrap() - 1;
+                        let order_str = bond_parts[1];
+                        let order = match order_str {
+                            "S" => BondOrder::Single,
+                            "D" => BondOrder::Double,
+                            "T" => BondOrder::Triple,
+                            "B" => BondOrder::Benzene,
+                            _ => BondOrder::Single,
+                        };
+                        bond_info.push((self.graph.vertices.len() - 1, target_idx, order));
+                    }
+                }
+            }
+        }
+
+        for (v1, v2, order) in bond_info {
+            if v1 < v2 {
+                self.add_bond(v1, v2, Bond::new(order));
+            }
+        }
+    }
+
+    pub fn is_isomorphic(&self, other: &Molecule) -> bool {
+        self.graph.is_isomorphic(&other.graph)
+    }
+
+    pub fn is_cyclic(&self) -> bool {
+        self.graph.is_cyclic()
+    }
+
+    pub fn is_linear(&self) -> bool {
+        let atom_count = self.graph.vertices.len();
+
+        if atom_count <= 1 {
+            return false;
+        }
+        if atom_count == 2 {
+            return true;
+        }
+        if self.is_cyclic() {
+            return false;
+        }
+
+        // A molecule is linear if all atoms have degree <= 2
+        for adj in &self.graph.edges {
+            if adj.len() > 2 {
+                return false;
+            }
+        }
+
+        // Check for specific linear bond patterns:
+        // 1. All double bonds (e.g., O=C=O)
+        let mut all_double = true;
+        for adj in &self.graph.edges {
+            for bond in adj.values() {
+                if bond.order != BondOrder::Double {
+                    all_double = false;
+                    break;
+                }
+            }
+        }
+        if all_double {
+            return true;
+        }
+
+        // 2. Alternating single and triple bonds (e.g., H-C#C-H)
+        let mut single_triple = true;
+        for adj in &self.graph.edges {
+            for bond in adj.values() {
+                if bond.order != BondOrder::Single && bond.order != BondOrder::Triple {
+                    single_triple = false;
+                    break;
+                }
+            }
+        }
+        if single_triple {
+            // Need at least one triple bond for this to be definitely linear in this simplified model
+            let mut has_triple = false;
+            for adj in &self.graph.edges {
+                for bond in adj.values() {
+                    if bond.order == BondOrder::Triple {
+                        has_triple = true;
+                        break;
+                    }
+                }
+            }
+            if has_triple {
+                return true;
+            }
+        }
+
+        false
+    }
+
+    pub fn is_subgraph_isomorphic(&self, other: &Molecule) -> bool {
+        self.graph.is_subgraph_isomorphic(&other.graph)
+    }
+
+    pub fn find_subgraph_isomorphisms(
+        &self,
+        other: &Molecule,
+    ) -> Vec<std::collections::HashMap<usize, usize>> {
+        self.graph.find_subgraph_isomorphisms(&other.graph)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::element;
+
+    #[test]
+    fn test_molecule_basic() {
+        let mut mol = Molecule::new();
+        let c1 = mol.add_atom(Atom::new(&element::C));
+        let c2 = mol.add_atom(Atom::new(&element::C));
+        mol.add_bond(c1, c2, Bond::new(BondOrder::Single));
+
+        assert_eq!(mol.graph.vertices.len(), 2);
+        assert_eq!(mol.get_atom(c1).unwrap().element.symbol, "C");
+        assert_eq!(mol.get_bond(c1, c2).unwrap().order, BondOrder::Single);
+    }
+
+    #[test]
+    fn test_atom_equivalence() {
+        let a1 = Atom::new(&element::C);
+        let a2 = Atom::new(&element::C);
+        let a3 = Atom::new(&element::O);
+
+        assert!(a1.equivalent(&a2));
+        assert!(!a1.equivalent(&a3));
+    }
+
+    #[test]
+    fn test_from_adjacency_list() {
+        let adj_list = "
+        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 C 1 {4,S} {6,S}
+        6 C 0 {5,S}
+        ";
+        let mut mol = Molecule::new();
+        mol.from_adjacency_list(adj_list);
+
+        assert_eq!(mol.graph.vertices.len(), 6);
+        assert_eq!(mol.get_atom(4).unwrap().radical_electrons, 1);
+        assert_eq!(mol.get_bond(0, 1).unwrap().order, BondOrder::Double);
+    }
+
+    #[test]
+    fn test_subgraph_isomorphism() {
+        let mut mol = Molecule::new();
+        mol.from_adjacency_list(
+            "
+        1 C 0 {2,D}
+        2 C 0 {1,D} {3,S}
+        3 C 0 {2,S}
+        ",
+        );
+
+        let mut pattern = Molecule::new();
+        pattern.from_adjacency_list(
+            "
+        1 C 0 {2,D}
+        2 C 0 {1,D}
+        ",
+        );
+
+        assert!(mol.is_subgraph_isomorphic(&pattern));
+        let mappings = mol.find_subgraph_isomorphisms(&pattern);
+        assert_eq!(mappings.len(), 2); // C=C can be mapped in 2 ways
+    }
+
+    #[test]
+    fn test_is_linear() {
+        let mut mol = Molecule::new();
+        mol.from_adjacency_list(
+            "
+        1 O 0 {2,D}
+        2 O 0 {1,D}
+        ",
+        );
+        assert!(mol.is_linear());
+
+        let mut mol2 = Molecule::new();
+        mol2.from_adjacency_list(
+            "
+        1 O 0 {2,D}
+        2 C 0 {1,D} {3,D}
+        3 O 0 {2,D}
+        ",
+        );
+        assert!(mol2.is_linear());
+
+        let mut mol3 = Molecule::new();
+        mol3.from_adjacency_list(
+            "
+        1 C 0 {2,S} {3,S}
+        2 H 0 {1,S}
+        3 H 0 {1,S}
+        ",
+        );
+        assert!(!mol3.is_linear());
+    }
+}
diff --git a/src/reaction.rs b/src/reaction.rs
new file mode 100644
index 0000000..76a73fd
--- /dev/null
+++ b/src/reaction.rs
@@ -0,0 +1,97 @@
+use crate::constants;
+use crate::kinetics::KineticsModel;
+use crate::species::{Species, TransitionState};
+use std::sync::Arc;
+
+pub struct Reaction {
+    pub index: i32,
+    pub reactants: Vec<Arc<Species>>,
+    pub products: Vec<Arc<Species>>,
+    pub kinetics: Option<Box<dyn KineticsModel>>,
+    pub reversible: bool,
+    pub transition_state: Option<TransitionState>,
+}
+
+impl Reaction {
+    pub fn new(reactants: Vec<Arc<Species>>, products: Vec<Arc<Species>>) -> Self {
+        Reaction {
+            index: -1,
+            reactants,
+            products,
+            kinetics: None,
+            reversible: true,
+            transition_state: None,
+        }
+    }
+
+    pub fn get_enthalpy_of_reaction(&self, t: f64) -> f64 {
+        let mut dh_rxn = 0.0;
+        for reactant in &self.reactants {
+            if let Some(thermo) = &reactant.thermo {
+                dh_rxn -= thermo.get_enthalpy(t);
+            }
+        }
+        for product in &self.products {
+            if let Some(thermo) = &product.thermo {
+                dh_rxn += thermo.get_enthalpy(t);
+            }
+        }
+        dh_rxn
+    }
+
+    pub fn get_entropy_of_reaction(&self, t: f64) -> f64 {
+        let mut ds_rxn = 0.0;
+        for reactant in &self.reactants {
+            if let Some(thermo) = &reactant.thermo {
+                ds_rxn -= thermo.get_entropy(t);
+            }
+        }
+        for product in &self.products {
+            if let Some(thermo) = &product.thermo {
+                ds_rxn += thermo.get_entropy(t);
+            }
+        }
+        ds_rxn
+    }
+
+    pub fn get_free_energy_of_reaction(&self, t: f64) -> f64 {
+        let mut dg_rxn = 0.0;
+        for reactant in &self.reactants {
+            if let Some(thermo) = &reactant.thermo {
+                dg_rxn -= thermo.get_free_energy(t);
+            }
+        }
+        for product in &self.products {
+            if let Some(thermo) = &product.thermo {
+                dg_rxn += thermo.get_free_energy(t);
+            }
+        }
+        dg_rxn
+    }
+
+    pub fn get_equilibrium_constant(&self, t: f64, k_type: &str) -> f64 {
+        let dg_rxn = self.get_free_energy_of_reaction(t);
+        let mut k = (-dg_rxn / constants::R / t).exp();
+
+        let p0 = 1.0e5;
+        match k_type {
+            "Kc" => {
+                let c0 = p0 / constants::R / t;
+                k *= c0.powi(self.products.len() as i32 - self.reactants.len() as i32);
+            }
+            "Kp" => {
+                k *= p0.powi(self.products.len() as i32 - self.reactants.len() as i32);
+            }
+            _ => {}
+        }
+        k
+    }
+
+    pub fn get_rate_coefficient(&self, t: f64, p: f64) -> f64 {
+        if let Some(kinetics) = &self.kinetics {
+            kinetics.get_rate_coefficient(t, p)
+        } else {
+            0.0
+        }
+    }
+}
diff --git a/src/species.rs b/src/species.rs
new file mode 100644
index 0000000..2ada807
--- /dev/null
+++ b/src/species.rs
@@ -0,0 +1,50 @@
+use crate::molecule::Molecule;
+use crate::states::StatesModel;
+use crate::thermo::ThermoModel;
+
+/// A chemical species.
+pub struct Species {
+    pub index: i32,
+    pub label: String,
+    pub thermo: Option<Box<dyn ThermoModel>>,
+    pub states: Option<StatesModel>,
+    pub molecule: Vec<Molecule>,
+    pub e0: f64,
+    pub molecular_weight: f64,
+    pub reactive: bool,
+}
+
+impl Species {
+    pub fn new(label: &str) -> Self {
+        Species {
+            index: -1,
+            label: label.to_string(),
+            thermo: None,
+            states: None,
+            molecule: Vec::new(),
+            e0: 0.0,
+            molecular_weight: 0.0,
+            reactive: true,
+        }
+    }
+}
+
+pub struct TransitionState {
+    pub label: String,
+    pub states: Option<StatesModel>,
+    pub e0: f64,
+    pub frequency: f64,
+    pub degeneracy: i32,
+}
+
+impl TransitionState {
+    pub fn new(label: &str) -> Self {
+        TransitionState {
+            label: label.to_string(),
+            states: None,
+            e0: 0.0,
+            frequency: 0.0,
+            degeneracy: 1,
+        }
+    }
+}
diff --git a/src/states.rs b/src/states.rs
new file mode 100644
index 0000000..4abd246
--- /dev/null
+++ b/src/states.rs
@@ -0,0 +1,243 @@
+use crate::constants;
+
+pub trait Mode {
+    fn get_partition_function(&self, t: f64) -> f64;
+    fn get_heat_capacity(&self, t: f64) -> f64;
+    fn get_enthalpy(&self, t: f64) -> f64;
+    fn get_entropy(&self, t: f64) -> f64;
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct Translation {
+    pub mass: f64, // kg/mol
+}
+
+impl Translation {
+    pub fn new(mass: f64) -> Self {
+        Translation { mass }
+    }
+}
+
+impl Mode for Translation {
+    fn get_partition_function(&self, t: f64) -> f64 {
+        let qt = ((2.0 * constants::PI_CONST * self.mass / constants::NA)
+            / (constants::H * constants::H))
+            .powf(1.5)
+            / 1.0e5;
+        qt * (constants::KB * t).powf(2.5)
+    }
+
+    fn get_heat_capacity(&self, _t: f64) -> f64 {
+        1.5 * constants::R
+    }
+
+    fn get_enthalpy(&self, t: f64) -> f64 {
+        1.5 * constants::R * t
+    }
+
+    fn get_entropy(&self, t: f64) -> f64 {
+        (self.get_partition_function(t).ln() + 1.5 + 1.0) * constants::R
+    }
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct RigidRotor {
+    pub linear: bool,
+    pub inertia: Vec<f64>, // kg*m^2
+    pub symmetry: i32,
+}
+
+impl RigidRotor {
+    pub fn new(linear: bool, inertia: Vec<f64>, symmetry: i32) -> Self {
+        RigidRotor {
+            linear,
+            inertia,
+            symmetry,
+        }
+    }
+}
+
+impl Mode for RigidRotor {
+    fn get_partition_function(&self, t: f64) -> f64 {
+        if self.linear {
+            let inertia = if !self.inertia.is_empty() {
+                self.inertia[0]
+            } else {
+                0.0
+            };
+            if inertia == 0.0 {
+                return 0.0;
+            }
+            constants::KB * t
+                / (self.symmetry as f64 * constants::H * constants::H
+                    / (8.0 * constants::PI_CONST * constants::PI_CONST * inertia))
+        } else {
+            if self.inertia.len() < 3 || self.inertia.contains(&0.0) {
+                return 0.0;
+            }
+            let mut theta = (constants::KB * t).powf(1.5)
+                * (8.0 * constants::PI_CONST.powi(2) / constants::H.powi(2)).powf(1.5);
+            theta *= (self.inertia[0] * self.inertia[1] * self.inertia[2]).sqrt();
+            theta *= constants::PI_CONST.sqrt() / self.symmetry as f64;
+            theta
+        }
+    }
+
+    fn get_heat_capacity(&self, _t: f64) -> f64 {
+        if self.linear {
+            constants::R
+        } else {
+            1.5 * constants::R
+        }
+    }
+
+    fn get_enthalpy(&self, t: f64) -> f64 {
+        if self.linear {
+            constants::R * t
+        } else {
+            1.5 * constants::R * t
+        }
+    }
+
+    fn get_entropy(&self, t: f64) -> f64 {
+        if self.linear {
+            (self.get_partition_function(t).ln() + 1.0) * constants::R
+        } else {
+            (self.get_partition_function(t).ln() + 1.5) * constants::R
+        }
+    }
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct HarmonicOscillator {
+    pub frequencies: Vec<f64>, // cm^-1
+}
+
+impl HarmonicOscillator {
+    pub fn new(frequencies: Vec<f64>) -> Self {
+        HarmonicOscillator { frequencies }
+    }
+}
+
+impl Mode for HarmonicOscillator {
+    fn get_partition_function(&self, t: f64) -> f64 {
+        let mut q = 1.0;
+        for &freq in &self.frequencies {
+            q /= 1.0 - (-freq / (0.695039 * t)).exp();
+        }
+        q
+    }
+
+    fn get_heat_capacity(&self, t: f64) -> f64 {
+        let mut cv = 0.0;
+        for &freq in &self.frequencies {
+            let x = freq / (0.695039 * t);
+            let exp_x = x.exp();
+            let one_minus_exp_x = 1.0 - exp_x;
+            cv += x * x * exp_x / (one_minus_exp_x * one_minus_exp_x);
+        }
+        cv * constants::R
+    }
+
+    fn get_enthalpy(&self, t: f64) -> f64 {
+        let mut h = 0.0;
+        for &freq in &self.frequencies {
+            let x = freq / (0.695039 * t);
+            h += x / (x.exp() - 1.0);
+        }
+        h * constants::R * t
+    }
+
+    fn get_entropy(&self, t: f64) -> f64 {
+        let mut s = self.get_partition_function(t).ln();
+        for &freq in &self.frequencies {
+            let x = freq / (0.695039 * t);
+            s += x / (x.exp() - 1.0);
+        }
+        s * constants::R
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_ethylene_modes() {
+        let t = 298.15;
+        let trans = Translation::new(0.02803);
+        let rot = RigidRotor::new(false, vec![5.6952e-47, 2.7758e-46, 3.3454e-46], 1);
+        let vib = HarmonicOscillator::new(vec![
+            834.50, 973.31, 975.37, 1067.1, 1238.5, 1379.5, 1472.3, 1691.3, 3121.6, 3136.7, 3192.5,
+            3221.0,
+        ]);
+
+        // Partition functions
+        assert!((trans.get_partition_function(t) / 1.01325 / 5.83338e6 - 1.0).abs() < 1e-3);
+        assert!((rot.get_partition_function(t) / 2.59622e3 - 1.0).abs() < 1e-3);
+        assert!((vib.get_partition_function(t) / 1.0481e0 - 1.0).abs() < 1e-3);
+
+        // Heat capacities (converted from cal/mol*K in original test to J/mol*K)
+        // Original used / 4.184 / 2.981, etc.
+        assert!((trans.get_heat_capacity(t) / (4.184 * 2.981) - 1.0).abs() < 1e-3);
+        assert!((rot.get_heat_capacity(t) / (4.184 * 2.981) - 1.0).abs() < 1e-3);
+    }
+
+    #[test]
+    fn test_oxygen_modes() {
+        let t = 298.15;
+        let trans = Translation::new(0.03199);
+        let rot = RigidRotor::new(true, vec![1.9271e-46], 2);
+        let vib = HarmonicOscillator::new(vec![1637.9]);
+
+        assert!((trans.get_partition_function(t) / 1.01325 / 7.11169e6 - 1.0).abs() < 1e-3);
+        assert!((rot.get_partition_function(t) / 7.13316e1 - 1.0).abs() < 1e-3);
+        assert!((vib.get_partition_function(t) / 1.000037e0 - 1.0).abs() < 1e-3);
+    }
+}
+
+pub struct StatesModel {
+    pub modes: Vec<Box<dyn Mode>>,
+    pub spin_multiplicity: i32,
+}
+
+impl StatesModel {
+    pub fn new(modes: Vec<Box<dyn Mode>>, spin_multiplicity: i32) -> Self {
+        StatesModel {
+            modes,
+            spin_multiplicity,
+        }
+    }
+
+    pub fn get_partition_function(&self, t: f64) -> f64 {
+        let mut q = 1.0;
+        for mode in &self.modes {
+            q *= mode.get_partition_function(t);
+        }
+        q * self.spin_multiplicity as f64
+    }
+
+    pub fn get_heat_capacity(&self, t: f64) -> f64 {
+        let mut cp = constants::R;
+        for mode in &self.modes {
+            cp += mode.get_heat_capacity(t);
+        }
+        cp
+    }
+
+    pub fn get_enthalpy(&self, t: f64) -> f64 {
+        let mut h = constants::R * t;
+        for mode in &self.modes {
+            h += mode.get_enthalpy(t);
+        }
+        h
+    }
+
+    pub fn get_entropy(&self, t: f64) -> f64 {
+        let mut s = 0.0;
+        for mode in &self.modes {
+            s += mode.get_entropy(t);
+        }
+        s
+    }
+}
diff --git a/src/thermo.rs b/src/thermo.rs
new file mode 100644
index 0000000..1fe8619
--- /dev/null
+++ b/src/thermo.rs
@@ -0,0 +1,180 @@
+use crate::constants;
+
+pub trait ThermoModel {
+    fn get_heat_capacity(&self, t: f64) -> f64;
+    fn get_enthalpy(&self, t: f64) -> f64;
+    fn get_entropy(&self, t: f64) -> f64;
+
+    fn get_free_energy(&self, t: f64) -> f64 {
+        self.get_enthalpy(t) - t * self.get_entropy(t)
+    }
+
+    fn is_temperature_valid(&self, t: f64, t_min: f64, t_max: f64) -> bool {
+        t >= t_min && t <= t_max
+    }
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct NASAPolynomial {
+    pub t_min: f64,
+    pub t_max: f64,
+    pub coeffs: [f64; 7],
+}
+
+impl NASAPolynomial {
+    pub fn new(t_min: f64, t_max: f64, coeffs: [f64; 7]) -> Self {
+        NASAPolynomial {
+            t_min,
+            t_max,
+            coeffs,
+        }
+    }
+}
+
+impl ThermoModel for NASAPolynomial {
+    fn get_heat_capacity(&self, t: f64) -> f64 {
+        let [c0, c1, c2, c3, c4, _, _] = self.coeffs;
+        (c0 + t * (c1 + t * (c2 + t * (c3 + c4 * t)))) * constants::R
+    }
+
+    fn get_enthalpy(&self, t: f64) -> f64 {
+        let [c0, c1, c2, c3, c4, c5, _] = self.coeffs;
+        let t2 = t * t;
+        let t4 = t2 * t2;
+        (c0 + c1 * t / 2.0 + c2 * t2 / 3.0 + c3 * t2 * t / 4.0 + c4 * t4 / 5.0 + c5 / t)
+            * constants::R
+            * t
+    }
+
+    fn get_entropy(&self, t: f64) -> f64 {
+        let [c0, c1, c2, c3, c4, _, c6] = self.coeffs;
+        let t2 = t * t;
+        let t4 = t2 * t2;
+        (c0 * t.ln() + c1 * t + c2 * t2 / 2.0 + c3 * t2 * t / 3.0 + c4 * t4 / 4.0 + c6)
+            * constants::R
+    }
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct WilhoitModel {
+    pub cp0: f64,
+    pub cp_inf: f64,
+    pub b: f64,
+    pub a0: f64,
+    pub a1: f64,
+    pub a2: f64,
+    pub a3: f64,
+    pub h0: f64,
+    pub s0: f64,
+}
+
+impl WilhoitModel {
+    #[allow(clippy::too_many_arguments)]
+    pub fn new(
+        cp0: f64,
+        cp_inf: f64,
+        a0: f64,
+        a1: f64,
+        a2: f64,
+        a3: f64,
+        h0: f64,
+        s0: f64,
+        b: f64,
+    ) -> Self {
+        WilhoitModel {
+            cp0,
+            cp_inf,
+            b,
+            a0,
+            a1,
+            a2,
+            a3,
+            h0,
+            s0,
+        }
+    }
+}
+
+impl ThermoModel for WilhoitModel {
+    fn get_heat_capacity(&self, t: f64) -> f64 {
+        let y = t / (t + self.b);
+        self.cp0
+            + (self.cp_inf - self.cp0)
+                * y
+                * y
+                * (1.0 + (y - 1.0) * (self.a0 + y * (self.a1 + y * (self.a2 + y * self.a3))))
+    }
+
+    fn get_enthalpy(&self, t: f64) -> f64 {
+        let y = t / (t + self.b);
+        let y2 = y * y;
+        let log_b_plus_t = (self.b + t).ln();
+        self.h0 + self.cp0 * t
+            - (self.cp_inf - self.cp0)
+                * t
+                * (y2
+                    * ((3.0 * self.a0 + self.a1 + self.a2 + self.a3) / 6.0
+                        + (4.0 * self.a1 + self.a2 + self.a3) * y / 12.0
+                        + (5.0 * self.a2 + self.a3) * y2 / 20.0
+                        + self.a3 * y2 * y / 5.0)
+                    + (2.0 + self.a0 + self.a1 + self.a2 + self.a3)
+                        * (y / 2.0 - 1.0 + (1.0 / y - 1.0) * log_b_plus_t))
+    }
+
+    fn get_entropy(&self, t: f64) -> f64 {
+        let y = t / (t + self.b);
+        self.s0 + self.cp_inf * t.ln()
+            - (self.cp_inf - self.cp0)
+                * (y.ln()
+                    + y * (1.0
+                        + y * (self.a0 / 2.0
+                            + y * (self.a1 / 3.0 + y * (self.a2 / 4.0 + y * self.a3 / 5.0)))))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::constants;
+
+    #[test]
+    fn test_wilhoit_regression() {
+        let wilhoit = WilhoitModel::new(
+            4.0 * constants::R,
+            21.0 * constants::R,
+            -3.95,
+            9.26,
+            -15.6,
+            8.55,
+            -6.151e04,
+            -790.2,
+            500.0,
+        );
+
+        let t_list = [
+            200.0, 400.0, 600.0, 800.0, 1000.0, 1200.0, 1400.0, 1600.0, 1800.0, 2000.0,
+        ];
+        let cp_expected = [
+            64.398, 94.765, 116.464, 131.392, 141.658, 148.830, 153.948, 157.683, 160.469, 162.589,
+        ];
+        let h_expected = [
+            -166312.0, -150244.0, -128990.0, -104110.0, -76742.9, -47652.6, -17347.1, 13834.8,
+            45663.0, 77978.1,
+        ];
+        let s_expected = [
+            287.421, 341.892, 384.685, 420.369, 450.861, 477.360, 500.708, 521.521, 540.262,
+            557.284,
+        ];
+
+        for i in 0..t_list.len() {
+            let t = t_list[i];
+            let cp = wilhoit.get_heat_capacity(t);
+            let h = wilhoit.get_enthalpy(t);
+            let s = wilhoit.get_entropy(t);
+
+            assert!((cp - cp_expected[i]).abs() / cp_expected[i] < 1e-3);
+            assert!((h - h_expected[i]).abs() / h_expected[i].abs() < 1e-3);
+            assert!((s - s_expected[i]).abs() / s_expected[i] < 1e-3);
+        }
+    }
+}