diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..79a6816
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,31 @@
+name: CI
+
+on:
+ push:
+ branches: [master]
+ pull_request:
+
+jobs:
+ test:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ - name: Install package
+ run: pip install -e ".[dev]"
+ - name: Run fast tests (units, VV equivalence, ideal-gas limit, energy conservation)
+ run: pytest tests/ --ignore=tests/test_vs_c_oracle.py -v
+
+ test-vs-oracle:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ - name: Install package
+ run: pip install -e ".[dev]"
+ - name: Build C oracle and run oracle-agreement tests
+ run: pytest tests/test_vs_c_oracle.py -v
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
new file mode 100644
index 0000000..5eee236
--- /dev/null
+++ b/.github/workflows/release.yml
@@ -0,0 +1,72 @@
+name: Publish to PyPI
+
+# Builds the pure-Python wheel/sdist and publishes via PyPI Trusted Publishing
+# (OIDC) -- no long-lived API token is stored in this repo. Publishing itself
+# is authorized by whoever configures the "noblegasmd" trusted publisher on
+# PyPI/TestPyPI to point at this repo + workflow; this workflow cannot publish
+# until that's done.
+#
+# Trigger: push a tag matching v*.*.* (e.g. v0.1.0) to publish to PyPI, or run
+# manually via workflow_dispatch with target=testpypi for a TestPyPI dry run.
+
+on:
+ push:
+ tags:
+ - "v*.*.*"
+ workflow_dispatch:
+ inputs:
+ target:
+ description: "Where to publish"
+ required: true
+ default: "testpypi"
+ type: choice
+ options:
+ - testpypi
+ - pypi
+
+jobs:
+ build:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.12"
+ - name: Build sdist and wheel
+ run: |
+ python -m pip install --upgrade build
+ python -m build
+ - uses: actions/upload-artifact@v4
+ with:
+ name: dist
+ path: dist/
+
+ publish-testpypi:
+ needs: build
+ if: github.event_name == 'workflow_dispatch' && inputs.target == 'testpypi'
+ runs-on: ubuntu-latest
+ environment: testpypi
+ permissions:
+ id-token: write # required for OIDC trusted publishing
+ steps:
+ - uses: actions/download-artifact@v4
+ with:
+ name: dist
+ path: dist/
+ - uses: pypa/gh-action-pypi-publish@release/v1
+ with:
+ repository-url: https://test.pypi.org/legacy/
+
+ publish-pypi:
+ needs: build
+ if: github.event_name == 'push' || (github.event_name == 'workflow_dispatch' && inputs.target == 'pypi')
+ runs-on: ubuntu-latest
+ environment: pypi
+ permissions:
+ id-token: write # required for OIDC trusted publishing
+ steps:
+ - uses: actions/download-artifact@v4
+ with:
+ name: dist
+ path: dist/
+ - uses: pypa/gh-action-pypi-publish@release/v1
diff --git a/.gitignore b/.gitignore
index ce2f69c..44c2c3e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,13 @@
MD.exe
+*.egg-info/
+dist/
+build/
+__pycache__/
+*.pyc
+.pytest_cache/
+.ipynb_checkpoints/
+tests/oracle/MD_oracle.exe
+tests/oracle/*_traj.xyz
+tests/oracle/*_output.txt
+tests/oracle/*_average.txt
+.venv/
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..9cecc1d
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ {one line to give the program's name and a brief idea of what it does.}
+ Copyright (C) {year} {name of author}
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see .
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ {project} Copyright (C) {year} {fullname}
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+.
diff --git a/MD_python_port_spec.md b/MD_python_port_spec.md
new file mode 100644
index 0000000..460f096
--- /dev/null
+++ b/MD_python_port_spec.md
@@ -0,0 +1,240 @@
+# Handoff Spec: Pure-Python (numba) Port of the NVE Lennard-Jones MD Engine
+
+**Purpose.** Reimplement the real-gas molecular-dynamics engine currently written in C
+(`MD.cpp`, Foley/Sweet/Akinfenwa, GPLv3) as a small, `pip install`-able Python package.
+The package must reproduce the published quasi-isotherms (Figure 4 of the *J. Chem. Educ.*
+article) and give students a clean programmatic API for parameter sweeps and numerical
+experiments — runnable in Google Colab with zero build step.
+
+**Source of truth.** `MD.cpp` is authoritative for all physics and all numeric constants.
+Where this spec and the C disagree, **the C wins** — transcribe constants directly from the
+source rather than from this document, and cross-check against the reference values quoted here.
+
+**The golden rule of validation.** MD is chaotic and the port uses a different RNG, so
+trajectories will *not* match step-for-step and are not expected to. Validate on
+**statistical/thermodynamic averages** (Z, P, T) and on invariants (energy conservation),
+never on trajectory identity.
+
+---
+
+## 1. Physics specification (must match `MD.cpp` exactly)
+
+### Ensemble & integrator
+- Pure **NVE** (microcanonical). No thermostat or barostat anywhere in the main loop.
+- **Velocity Verlet**, timestep `dt`:
+ 1. `r += v·dt + 0.5·a·dt²`
+ 2. `v += 0.5·a·dt` (first half-kick)
+ 3. recompute accelerations at new `r`
+ 4. `v += 0.5·a·dt` (second half-kick)
+ 5. apply walls (below)
+- **Optimization note:** the reference recomputes forces *twice* per step (once redundantly at
+ the top of `VelocityVerlet`, once mid-step). A standard velocity Verlet that caches the
+ previous step's second force computation is numerically identical and ~2× faster. Implement
+ the cached single-force-per-step version, but assert equivalence against the reference during
+ testing.
+
+### Forces (no cutoff)
+- Full **O(N²)** pairwise Lennard-Jones over all distinct pairs `i 0 departs
+ from the published methodology.
+
+### Initialization
+- **N = 216** by default. Simple-cubic lattice: `n = ceil(N^(1/3))` per axis (= 6 for 216),
+ spacing `pos = L/n`, positions at `((i+0.5)·pos, (j+0.5)·pos, (k+0.5)·pos)`.
+- Velocities: draw each component from a Gaussian, remove center-of-mass velocity, then rescale
+ to the target temperature with `lambda = sqrt(3·(N−1)·Tinit / vSqdSum)` (the `(N−1)` accounts
+ for the removed COM d.o.f.). `Tinit` is in natural units (input Kelvin ÷ `TempFac`).
+- RNG: use `numpy.random.Generator(PCG64(seed))`. Bit-matching C's `rand()` is neither possible
+ nor required.
+
+### Timestep & step count (per gas)
+- Non-helium: `dt = 0.5e-14 / timefac`, `NumTime = 20000`.
+- Helium: `dt = 0.2e-14 / timefac`, `NumTime = 50000`.
+
+---
+
+## 2. Units and per-gas constants
+
+Natural units per noble gas; SI is used only for reporting. **Transcribe these from `MD.cpp`**;
+values below are for cross-checking only.
+
+Global: `NA = 6.022140857e23`, `kBSI = 1.38064852e-23`.
+
+| Gas | VolFac (m³) | PressFac (Pa) | TempFac (K) | timefac (s) |
+|-----|-------------------------|-----------------------|----------------------|-------------------------|
+| He | 1.8399744000000005e-29 | 8152287.336171632 | 10.864459551225972 | 1.7572698825166272e-12 |
+| Ne | 2.0570823999999997e-29 | 27223022.27659913 | 40.560648991243625 | 2.1192341945685407e-12 |
+| Ar | 3.7949992920124995e-29 | 51695201.06691862 | 142.0950000000000 | 2.09618e-12 |
+| Kr | 4.5882712000000004e-29 | 59935428.40275003 | 199.1817584391428 | 8.051563913585078e-13 |
+| Xe | 5.4872e-29 | 70527773.72794868 | 280.30305642163006 | 9.018957925790732e-13 |
+
+Volume from density: input `rho` in **mol/m³**, then `Vol = N/(rho·NA)`, convert to natural
+units `Vol /= VolFac`, box side `L = Vol^(1/3)`. Preserve the reference's guard rails
+(reject `rho ≤ 0`, warn/της on `Vol < N`, i.e. density above ~1 particle per natural-unit volume).
+
+---
+
+## 3. Target API
+
+Ergonomic for both students and an inline LLM assistant. Rich docstrings and type hints
+directly improve autocomplete quality (and reduce wasted AI calls), so treat them as part of
+the deliverable.
+
+```python
+from noblegasmd import run, sweep # package name TBD; check PyPI availability
+
+# Single simulation
+result = run(
+ gas: str = "Ar", # He | Ne | Ar | Kr | Xe
+ T: float = 300.0, # initial temperature, Kelvin
+ rho: float = 40.0, # number density, mol/m^3
+ n_particles: int = 216,
+ n_steps: int | None = None, # default 20000 (50000 for He)
+ n_equil: int = 0, # steps discarded before averaging; 0 = match reference
+ seed: int | None = None,
+ record_trajectory: bool = False,
+)
+# result exposes at least: T_avg (K), P_avg (Pa), Z, gc (PV/nT), V (m^3), N,
+# energy drift diagnostic, and optionally per-step arrays / trajectory.
+
+# Parameter sweep -> tidy DataFrame, one row per (state point, replicate)
+df = sweep(
+ gas="Ar",
+ T=[100, 200, 300, 400],
+ rho=np.linspace(1.0, 5000.0, 25),
+ n_replicates=3, # independent seeds per state point for error bars
+ seed=0,
+)
+# Columns: gas, T_set, rho_set, seed, T_avg, P_avg, Z, gc, V, N, n_steps
+# Figure 4 = df.groupby("T_set") plotted as Z (or gc) vs rho_set (or P_avg).
+```
+
+Keep the heavy loop inside numba; `run`/`sweep` are thin Python wrappers. **Never** place an
+LLM/API call inside a sweep loop — that reintroduces the per-run token cost the whole design
+avoids.
+
+---
+
+## 4. Implementation guidance (numba)
+
+- Put the integrator + force kernel in `@njit(fastmath=True, cache=True)` functions operating on
+ preallocated `(N,3)` float64 arrays. Write the force loop as **explicit nested loops with
+ Newton's third law** — do *not* build N×N numpy temporaries (slower and memory-heavy).
+- Consider `parallel=True` + `prange` on the outer force loop for a free ~2× on Colab's second
+ vCPU; benchmark with and without.
+- **Warm up the JIT** with a throwaway tiny run at import or first call so compile time doesn't
+ land on the user's first timed simulation.
+- No per-step Python callbacks or printing inside the kernel; collect diagnostics into arrays and
+ process after the loop.
+- Ship a small timing harness (`benchmark.py`) that reports wall-clock for a couple of
+ `(N, n_steps)` points, so the team can measure Colab performance directly rather than guess.
+
+---
+
+## 5. Validation protocol
+
+Build the C as a **golden oracle** and compare averages. `MD.cpp` reads its inputs from stdin
+(title, gas, T in K, rho in mol/m³) and writes averages to `_average.txt`; drive it by
+piping inputs, parse `Z`, `P`, `T` from that file.
+
+Required tests (pytest):
+1. **Ideal-gas limit** (oracle-free): at low density, `Z → 1`. Assert `|Z − 1| < ~0.02` for the
+ lowest-density state points. Tune threshold with the domain expert.
+2. **Energy conservation** (oracle-free): elastic walls + conservative forces ⇒ total energy is
+ conserved. Assert relative drift `|E(t) − E(0)| / |E(0)|` stays below a small bound over a full
+ run. (Do **not** test momentum conservation — walls break it by design.)
+3. **Agreement with the C oracle**: for a grid of Ar state points (e.g. T ∈ {100, 200, 300, 400} K
+ across a density range spanning ideal → strongly non-ideal), the Python **ensemble mean** of Z
+ over `n_replicates` seeds must agree with the C value within combined statistical error
+ (start with ~2–3 %, or within 2σ of the ensemble — final tolerance is the domain expert's call).
+ The C is deterministic (unseeded `rand()` ⇒ fixed seed 1), so it yields one value per input;
+ compare it against the Python ensemble.
+4. **Constants round-trip**: unit-test that the per-gas factors and derived quantities match
+ `MD.cpp` to full precision.
+
+Acceptance test / **definition of done**: a `notebooks/reproduce_figure4.ipynb` that calls
+`sweep()` and regenerates the published quasi-isotherms within agreed tolerance, plus all tests
+green and a clean `pip install` in a fresh Colab runtime.
+
+---
+
+## 6. Packaging & publishing
+
+- Pure Python + numba ⇒ a **pure-Python wheel** (`py3-none-any`), no compiled extensions, no
+ `cibuildwheel`/manylinux machinery. Dependencies: `numpy`, `numba`, `pandas` (and `matplotlib`
+ for the example notebook only).
+- `pyproject.toml` with a standard backend (hatchling or setuptools). Target the Python versions
+ Colab currently ships and that numba supports; pin a floor.
+- **License must be GPLv3** (derivative of GPLv3 `MD.cpp`). Include `LICENSE` and proper
+ attribution to the original authors and the *J. Chem. Educ.* article.
+- Verify the package name is free on PyPI before committing to it.
+- **Publishing / credential boundary:** Claude Code can build the sdist/wheel, do a **TestPyPI**
+ dry-run install, and set up a GitHub Actions release workflow using **PyPI Trusted Publishing
+ (OIDC)** so no long-lived token exists. Do **not** hand it a PyPI API token; the human performs
+ the final publish / authorizes the trusted publisher.
+
+### Suggested repo layout
+```
+<pkg>/
+ pyproject.toml
+ README.md
+ LICENSE # GPLv3
+ src/<pkg>/__init__.py
+ src/<pkg>/units.py # per-gas constants (transcribed from MD.cpp)
+ src/<pkg>/engine.py # numba kernels: forces, velocity Verlet, walls
+ src/<pkg>/api.py # run(), sweep(), result object
+ benchmark.py
+ tests/
+ oracle/ # MD.cpp + build script + stdin driver
+ test_units.py
+ test_energy_conservation.py
+ test_ideal_gas_limit.py
+ test_vs_c_oracle.py
+ notebooks/
+ reproduce_figure4.ipynb # the acceptance test
+```
+
+---
+
+## 7. Decisions to confirm with the domain expert
+
+- Final numerical tolerances for oracle agreement and the ideal-gas limit (§5).
+- Whether `n_equil` should default to 0 (faithful to the paper) — recommended yes.
+- Package name.
+- Whether the shipped trajectory format needs to stay VMD-compatible or can be simplified
+ (the browser app already covers live visualization, so the Python package can focus on
+ thermodynamics).
+
diff --git a/PUBLISHING.md b/PUBLISHING.md
new file mode 100644
index 0000000..791c9bf
--- /dev/null
+++ b/PUBLISHING.md
@@ -0,0 +1,48 @@
+# Publishing `noblegasmd`
+
+Claude Code will not touch PyPI credentials or trigger a real publish. This
+document covers what's already prepared and what a human needs to do.
+
+## Already done
+
+- `pyproject.toml` builds a pure-Python wheel (`py3-none-any`) via hatchling.
+- `python -m build` produces `dist/noblegasmd-<version>-py3-none-any.whl` and
+ the sdist; both pass `twine check`.
+- A clean-venv install of the built wheel was verified locally (`pip install
+ dist/*.whl` in a fresh venv, then `import noblegasmd; noblegasmd.run(...)`).
+- `.github/workflows/release.yml` builds the package and publishes via
+ **PyPI Trusted Publishing (OIDC)** — no API token is stored in this repo.
+
+## What a human needs to do
+
+1. **Register the trusted publisher** (one-time, per target):
+ - On [test.pypi.org](https://test.pypi.org) and/or
+ [pypi.org](https://pypi.org), create (or claim) the `noblegasmd` project,
+ then under its "Publishing" settings add a trusted publisher pointing at:
+ - Repository: `jayfoleyiv/MolecularDynamics` (or wherever this ends up)
+ - Workflow: `.github/workflows/release.yml`
+ - Environment: `testpypi` or `pypi` (must match the workflow's
+ `environment:` field for that job)
+ - Create matching GitHub Environments named `testpypi` and `pypi` under
+ the repo's Settings > Environments (optionally with required reviewers
+ for `pypi`, as a manual publish gate).
+2. **TestPyPI dry run**: Actions tab -> "Publish to PyPI" -> "Run workflow" ->
+ target = `testpypi`. Then verify with:
+ ```bash
+ pip install -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ noblegasmd
+ ```
+ (`--extra-index-url` is needed because `numpy`/`numba`/`pandas` aren't on
+ TestPyPI.)
+3. **Real publish**: push a tag `vX.Y.Z` (or run the workflow manually with
+ target = `pypi`) once you're satisfied.
+
+## Local dry run (no GitHub Actions, no upload)
+
+```bash
+python -m pip install --upgrade build twine
+python -m build # writes dist/*.whl and dist/*.tar.gz
+python -m twine check dist/* # metadata sanity check, no network
+python -m venv /tmp/noblegasmd_check
+/tmp/noblegasmd_check/bin/pip install dist/noblegasmd-*-py3-none-any.whl
+/tmp/noblegasmd_check/bin/python -c "from noblegasmd import run; print(run(gas='Ar', T=300, rho=40, n_steps=500, seed=1).Z)"
+```
diff --git a/README.md b/README.md
index f5c7fcf..8c5b16d 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,30 @@
+## Python package (`noblegasmd`)
+
+A `pip install`-able, numba-accelerated Python port of the same NVE Lennard-Jones engine,
+runnable with zero build step (e.g. in Google Colab). `MD.cpp` below remains the authoritative
+physics reference; the port is validated against it (see `tests/`).
+
+```bash
+pip install noblegasmd
+```
+
+```python
+from noblegasmd import run, sweep
+
+result = run(gas="Ar", T=300.0, rho=40.0)
+print(result.Z, result.P_avg, result.T_avg)
+
+df = sweep(gas="Ar", T=[100, 200, 300, 400], rho=[1, 40, 500, 2000], n_replicates=3)
+```
+
+See [`notebooks/reproduce_figure4.ipynb`](notebooks/reproduce_figure4.ipynb) to regenerate the
+published quasi-isotherms, and [`MD_python_port_spec.md`](MD_python_port_spec.md) for the full
+port specification. Development install: `pip install -e ".[dev]"`, then `pytest`.
+
+---
+
+## C++ reference implementation (`MD.cpp`)
+
- Source code for Molecular Dynamics Program - can compile and run on Linux, Windows, and Mac OSX
- More information about this program, including detailed instructions for its use, can be found [here for instructions](https://pubs.acs.org/doi/suppl/10.1021/acs.jchemed.7b00747) and [here for discussion of its use in an undergraduate laboratory setting](https://pubs.acs.org/doi/pdf/10.1021/acs.jchemed.7b00747)
diff --git a/benchmark.py b/benchmark.py
new file mode 100644
index 0000000..932accf
--- /dev/null
+++ b/benchmark.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+"""Timing harness for noblegasmd's numba kernel: reports wall-clock for a few
+(N, n_steps) points, so performance on e.g. Colab can be measured directly
+rather than guessed. Run with: python benchmark.py
+"""
+from __future__ import annotations
+
+import time
+
+from noblegasmd import run
+from noblegasmd.api import _ensure_warm
+
+POINTS = [
+ (108, 2000),
+ (216, 20000), # default production point
+ (216, 50000), # He-equivalent step count
+ (512, 5000),
+]
+
+
+def main() -> None:
+ print("Warming up JIT (not timed)...")
+ t0 = time.perf_counter()
+ _ensure_warm()
+ print(f" warmup: {time.perf_counter() - t0:.2f}s\n")
+
+ print(f"{'N':>6} {'n_steps':>8} {'wall (s)':>10} {'steps/s':>10}")
+ for n_particles, n_steps in POINTS:
+ t0 = time.perf_counter()
+ run(gas="Ar", T=300.0, rho=40.0, n_particles=n_particles,
+ n_steps=n_steps, seed=0)
+ elapsed = time.perf_counter() - t0
+ print(f"{n_particles:>6} {n_steps:>8} {elapsed:>10.3f} {n_steps / elapsed:>10.1f}")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/notebooks/reproduce_figure4.ipynb b/notebooks/reproduce_figure4.ipynb
new file mode 100644
index 0000000..70ae37a
--- /dev/null
+++ b/notebooks/reproduce_figure4.ipynb
@@ -0,0 +1,159 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "57f7fb4e",
+ "metadata": {},
+ "source": [
+ "# Reproducing Figure 4: Ar quasi-isotherms (Z vs. density)\n",
+ "\n",
+ "This notebook is the acceptance test for `noblegasmd`: it calls `sweep()` to regenerate the\n",
+ "published quasi-isotherms of the compressibility factor `Z = PV/(NkT)` vs. number density for\n",
+ "argon at several temperatures, from the *J. Chem. Educ.* article this package ports\n",
+ "(Foley, Sweet, Akinfenwa). Run top-to-bottom in a fresh environment (including Google Colab)\n",
+ "with just `pip install noblegasmd`.\n",
+ "\n",
+ "No LLM/API calls are used anywhere in the simulation loop -- `sweep()` is a thin Python\n",
+ "wrapper around a numba-jitted kernel."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "fabeef7f",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "ModuleNotFoundError",
+ "evalue": "No module named 'noblegasmd'",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
+ "\u001b[31mModuleNotFoundError\u001b[39m Traceback (most recent call last)",
+ "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[1]\u001b[39m\u001b[32m, line 5\u001b[39m\n\u001b[32m 1\u001b[39m \u001b[38;5;66;03m# In Colab: !pip install noblegasmd\u001b[39;00m\n\u001b[32m 2\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m numpy \u001b[38;5;28;01mas\u001b[39;00m np\n\u001b[32m 3\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m matplotlib.pyplot \u001b[38;5;28;01mas\u001b[39;00m plt\n\u001b[32m 4\u001b[39m \n\u001b[32m----> \u001b[39m\u001b[32m5\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m noblegasmd \u001b[38;5;28;01mimport\u001b[39;00m sweep\n",
+ "\u001b[31mModuleNotFoundError\u001b[39m: No module named 'noblegasmd'"
+ ]
+ }
+ ],
+ "source": [
+ "# In Colab: !pip install noblegasmd\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "from noblegasmd import sweep"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d44f564c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Full grid: takes a few minutes on a single Colab CPU (N=216, 20000 steps/run,\n",
+ "# 4 temperatures x 25 densities x 3 replicates = 300 runs).\n",
+ "# Reduce n_replicates or the density grid for a quick smoke test.\n",
+ "df = sweep(\n",
+ " gas=\"Ar\",\n",
+ " T=[100, 200, 300, 400],\n",
+ " rho=np.linspace(1.0, 5000.0, 25),\n",
+ " n_replicates=3,\n",
+ " seed=0,\n",
+ ")\n",
+ "df.head()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9eda5dd7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Ensemble mean +/- standard error across replicates, one isotherm per T.\n",
+ "summary = (\n",
+ " df.groupby([\"T_set\", \"rho_set\"])[\"Z\"]\n",
+ " .agg([\"mean\", \"std\", \"count\"])\n",
+ " .reset_index()\n",
+ ")\n",
+ "summary[\"sem\"] = summary[\"std\"] / np.sqrt(summary[\"count\"])\n",
+ "\n",
+ "fig, ax = plt.subplots(figsize=(7, 5))\n",
+ "for T_set, group in summary.groupby(\"T_set\"):\n",
+ " ax.errorbar(\n",
+ " group[\"rho_set\"], group[\"mean\"], yerr=group[\"sem\"],\n",
+ " marker=\"o\", markersize=3, linewidth=1, capsize=2,\n",
+ " label=f\"T = {T_set:.0f} K\",\n",
+ " )\n",
+ "ax.axhline(1.0, color=\"gray\", linestyle=\"--\", linewidth=1, label=\"ideal gas (Z=1)\")\n",
+ "ax.set_xlabel(\"number density (mol/m$^3$)\")\n",
+ "ax.set_ylabel(\"Z = PV / (N k$_B$ T)\")\n",
+ "ax.set_title(\"Argon quasi-isotherms (reproduction of Figure 4)\")\n",
+ "ax.legend()\n",
+ "fig.tight_layout()\n",
+ "fig.savefig(\"figure4_reproduction.png\", dpi=150)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c1765cfd",
+ "metadata": {},
+ "source": [
+ "## Acceptance check\n",
+ "\n",
+ "Compares against the C oracle (`MD.cpp`) at a handful of state points, matching\n",
+ "`tests/test_vs_c_oracle.py`. This cell requires a C++ compiler and is optional --\n",
+ "skip it if you're only checking the qualitative isotherm shapes above (e.g. on Colab\n",
+ "without a working oracle build)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "94266aac",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sys\n",
+ "from pathlib import Path\n",
+ "\n",
+ "try:\n",
+ " repo_root = Path.cwd().parent\n",
+ " sys.path.insert(0, str(repo_root / \"tests\" / \"oracle\"))\n",
+ " from oracle import run_oracle\n",
+ "\n",
+ " for T_set in [100.0, 200.0, 300.0, 400.0]:\n",
+ " rho_set = 40.0\n",
+ " python_z = df.loc[\n",
+ " (df[\"T_set\"] == T_set) & np.isclose(df[\"rho_set\"], rho_set, atol=210), \"Z\"\n",
+ " ]\n",
+ " oracle_z = run_oracle(\"Ar\", T_set, rho_set, title=f\"nb_check_{T_set:.0f}\").Z\n",
+ " print(f\"T={T_set:6.1f} K oracle Z={oracle_z:.4f} python nearby Z mean={python_z.mean():.4f}\")\n",
+ "except (ImportError, FileNotFoundError, RuntimeError) as exc:\n",
+ " print(f\"Oracle comparison skipped ({exc}); the isotherm plot above is still valid.\")"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.13.13"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..7dab5d3
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,42 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "noblegasmd"
+version = "0.1.0"
+description = "Numba port of the NVE Lennard-Jones noble-gas MD engine from MD.cpp (Foley, Sweet, Akinfenwa)"
+readme = "README.md"
+license = "GPL-3.0-or-later"
+license-files = ["LICENSE"]
+authors = [
+ { name = "Jonathan J. Foley IV" },
+]
+requires-python = ">=3.9"
+classifiers = [
+ "Development Status :: 3 - Alpha",
+ "Intended Audience :: Science/Research",
+ "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
+ "Programming Language :: Python :: 3",
+ "Topic :: Scientific/Engineering :: Chemistry",
+ "Topic :: Scientific/Engineering :: Physics",
+]
+dependencies = [
+ "numpy>=1.23",
+ "numba>=0.58",
+ "pandas>=1.5",
+]
+
+[project.optional-dependencies]
+notebook = ["matplotlib>=3.5"]
+dev = ["pytest>=7.0", "matplotlib>=3.5"]
+
+[project.urls]
+Homepage = "https://github.com/jayfoleyiv/MolecularDynamics"
+Repository = "https://github.com/jayfoleyiv/MolecularDynamics"
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/noblegasmd"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
diff --git a/src/noblegasmd/__init__.py b/src/noblegasmd/__init__.py
new file mode 100644
index 0000000..145f9c4
--- /dev/null
+++ b/src/noblegasmd/__init__.py
@@ -0,0 +1,21 @@
+"""noblegasmd: a numba port of the NVE Lennard-Jones noble-gas MD engine
+originally implemented in MD.cpp by Foley, Sweet, and Akinfenwa (GPLv3).
+
+See https://doi.org/10.1021/acs.jchemed.7b00790 for the original article.
+"""
+from .api import RunResult, run, sweep
+from .units import GAS_CONSTANTS, KB_SI, NA, GasConstants, get_gas_constants
+
+__version__ = "0.1.0"
+
+__all__ = [
+ "run",
+ "sweep",
+ "RunResult",
+ "GAS_CONSTANTS",
+ "GasConstants",
+ "get_gas_constants",
+ "KB_SI",
+ "NA",
+ "__version__",
+]
diff --git a/src/noblegasmd/api.py b/src/noblegasmd/api.py
new file mode 100644
index 0000000..07c4c8f
--- /dev/null
+++ b/src/noblegasmd/api.py
@@ -0,0 +1,185 @@
+"""Public API: run() for a single simulation, sweep() for a parameter grid."""
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from typing import Iterable, Optional
+
+import numpy as np
+import pandas as pd
+
+from . import engine
+from .units import KB_SI, NA, get_gas_constants, steps_for_gas
+
+_warmed_up = False
+
+
+def _ensure_warm() -> None:
+ global _warmed_up
+ if not _warmed_up:
+ engine.warmup()
+ _warmed_up = True
+
+
+@dataclass
+class RunResult:
+ """Result of a single :func:`run` call.
+
+ Averages (``T_avg``, ``P_avg``, ``Z``, ``gc``) follow MD.cpp's convention:
+ with ``n_equil=0`` they are means over all recorded iterations, matching
+ the published methodology. ``n_equil > 0`` discards that many leading
+ iterations before averaging, which is a documented departure from the
+ paper, useful only for pedagogical experiments.
+ """
+
+ gas: str
+ T_set: float
+ rho_set: float
+ N: int
+ n_steps: int
+ n_equil: int
+ seed: Optional[int]
+ T_avg: float # K
+ P_avg: float # Pa
+ Z: float # compressibility factor, PV/(N kB T)
+ gc: float # PV/(nT), J/(mol K)
+ V: float # m^3
+ energy_drift: float # max |E(t)-E(0)| / |E(0)|, natural units
+ instantaneous_T: np.ndarray = field(repr=False) # K, per iteration
+ instantaneous_P: np.ndarray = field(repr=False) # Pa, per iteration
+ kinetic_energy: np.ndarray = field(repr=False) # natural units, per iteration
+ potential_energy: np.ndarray = field(repr=False) # natural units, per iteration
+ trajectory: Optional[np.ndarray] = field(default=None, repr=False) # (n_records, N, 3)
+
+
+def run(
+ gas: str = "Ar",
+ T: float = 300.0,
+ rho: float = 40.0,
+ n_particles: int = 216,
+ n_steps: Optional[int] = None,
+ n_equil: int = 0,
+ seed: Optional[int] = None,
+ record_trajectory: bool = False,
+) -> RunResult:
+ """Run a single NVE Lennard-Jones simulation and return thermodynamic averages.
+
+ Parameters
+ ----------
+ gas : one of "He", "Ne", "Ar", "Kr", "Xe".
+ T : initial temperature, Kelvin. Must be positive.
+ rho : number density, mol/m^3. Must be positive.
+ n_particles : number of particles (simple-cubic lattice init).
+ n_steps : number of Verlet steps; defaults to MD.cpp's convention
+ (50000 for He, 20000 otherwise).
+ n_equil : leading iterations discarded before averaging. Default 0
+ matches the published methodology (no discard); values > 0 are a
+ documented departure, for pedagogical experiments only.
+ seed : seed for the velocity-initialization RNG (numpy PCG64). The
+ reference C uses unseeded (fixed-seed) ``rand()``; this port uses an
+ independent RNG by design, so trajectories will not match the C
+ step-for-step — validate against statistical averages only.
+ record_trajectory : if True, also return the full position trajectory
+ (memory: ~24 bytes * N * (n_steps+1)).
+ """
+ if T <= 0:
+ raise ValueError(f"T must be positive, got {T}")
+ if rho <= 0:
+ raise ValueError(f"rho must be positive, got {rho}")
+
+ gc_const = get_gas_constants(gas)
+
+ Vol = n_particles / (rho * NA)
+ Vol /= gc_const.vol_fac
+ if Vol < n_particles:
+ raise ValueError(
+ f"Density too high: N={n_particles} particles but only {Vol:.4f} "
+ "natural-unit volume available. Simulations with density greater "
+ "than 1 particle/(natural unit of volume) may diverge."
+ )
+ L = Vol ** (1.0 / 3.0)
+
+ dt_numerator, default_n_steps = steps_for_gas(gas)
+ dt = dt_numerator / gc_const.timefac
+ if n_steps is None:
+ n_steps = default_n_steps
+ if n_equil < 0 or n_equil > n_steps:
+ raise ValueError(f"n_equil must be in [0, n_steps], got {n_equil}")
+
+ T_init_nu = T / gc_const.temp_fac
+
+ rng = np.random.Generator(np.random.PCG64(seed))
+ pos0 = engine.initialize_lattice(n_particles, L)
+ vel0 = engine.initialize_velocities(n_particles, T_init_nu, rng)
+
+ _ensure_warm()
+ trajectory = None
+ if record_trajectory:
+ T_nu, P_nu, KE, PE, _pos_f, _vel_f, trajectory = engine.simulate_with_trajectory(
+ pos0, vel0, L, dt, n_steps
+ )
+ else:
+ T_nu, P_nu, KE, PE, _pos_f, _vel_f = engine.simulate(pos0, vel0, L, dt, n_steps)
+
+ T_si = T_nu * gc_const.temp_fac
+ P_si = P_nu * gc_const.press_fac
+
+ # Match MD.cpp exactly at n_equil=0: sum over all NumTime+1 recorded
+ # iterations, divide by NumTime (== n_steps), not by the record count.
+ if n_equil == 0:
+ T_avg = float(np.sum(T_si) / n_steps)
+ P_avg = float(np.sum(P_si) / n_steps)
+ else:
+ T_avg = float(np.mean(T_si[n_equil:]))
+ P_avg = float(np.mean(P_si[n_equil:]))
+
+ V_si = Vol * gc_const.vol_fac
+ Z = P_avg * V_si / (n_particles * KB_SI * T_avg)
+ gc_val = NA * P_avg * V_si / (n_particles * T_avg)
+
+ E_nu = KE + PE
+ energy_drift = float(np.max(np.abs(E_nu - E_nu[0])) / np.abs(E_nu[0]))
+
+ return RunResult(
+ gas=gas, T_set=T, rho_set=rho, N=n_particles, n_steps=n_steps,
+ n_equil=n_equil, seed=seed,
+ T_avg=T_avg, P_avg=P_avg, Z=Z, gc=gc_val, V=V_si,
+ energy_drift=energy_drift,
+ instantaneous_T=T_si, instantaneous_P=P_si,
+ kinetic_energy=KE, potential_energy=PE,
+ trajectory=trajectory,
+ )
+
+
+def sweep(
+ gas: str = "Ar",
+ T: Iterable[float] = (100.0, 200.0, 300.0, 400.0),
+ rho: Iterable[float] = (40.0,),
+ n_particles: int = 216,
+ n_steps: Optional[int] = None,
+ n_equil: int = 0,
+ n_replicates: int = 1,
+ seed: Optional[int] = 0,
+) -> pd.DataFrame:
+ """Run :func:`run` over a grid of (T, rho) state points x replicate seeds.
+
+ Returns a tidy DataFrame with one row per (state point, replicate):
+ gas, T_set, rho_set, seed, T_avg, P_avg, Z, gc, V, N, n_steps.
+ """
+ rows = []
+ base_rng = np.random.default_rng(seed)
+ for T_val in T:
+ for rho_val in rho:
+ for rep in range(n_replicates):
+ run_seed = int(base_rng.integers(0, 2**32 - 1))
+ result = run(
+ gas=gas, T=T_val, rho=rho_val, n_particles=n_particles,
+ n_steps=n_steps, n_equil=n_equil, seed=run_seed,
+ )
+ rows.append({
+ "gas": result.gas, "T_set": T_val, "rho_set": rho_val,
+ "seed": run_seed, "T_avg": result.T_avg, "P_avg": result.P_avg,
+ "Z": result.Z, "gc": result.gc, "V": result.V,
+ "N": result.N, "n_steps": result.n_steps,
+ "energy_drift": result.energy_drift,
+ })
+ return pd.DataFrame(rows)
diff --git a/src/noblegasmd/engine.py b/src/noblegasmd/engine.py
new file mode 100644
index 0000000..1d1ea36
--- /dev/null
+++ b/src/noblegasmd/engine.py
@@ -0,0 +1,237 @@
+"""Numba kernels: Lennard-Jones forces, velocity Verlet, elastic walls.
+
+Physics matches ``MD.cpp`` (Foley, Sweet, Akinfenwa, GPLv3) exactly: natural
+units (sigma = epsilon = m = kB = 1), full O(N^2) pairwise LJ with no cutoff,
+hard elastic walls on [0, L)^3 that reverse velocity without repositioning.
+
+Two step variants are provided:
+
+- ``vv_step_naive`` mirrors the reference's redundant double force-evaluation
+ per step exactly, for equivalence testing only.
+- ``vv_step_cached`` is the production path: one force evaluation per step,
+ reusing the previous step's post-kick acceleration as this step's initial
+ acceleration (numerically identical to the naive version, ~2x faster).
+"""
+from __future__ import annotations
+
+import numpy as np
+from numba import njit
+
+
+@njit(cache=True, fastmath=True)
+def compute_forces(pos: np.ndarray, N: int) -> np.ndarray:
+ """O(N^2) pairwise Lennard-Jones acceleration (m=1 => a=F), no cutoff."""
+ acc = np.zeros((N, 3))
+ for i in range(N - 1):
+ for j in range(i + 1, N):
+ rij0 = pos[i, 0] - pos[j, 0]
+ rij1 = pos[i, 1] - pos[j, 1]
+ rij2 = pos[i, 2] - pos[j, 2]
+ rSqd = rij0 * rij0 + rij1 * rij1 + rij2 * rij2
+ f = 24.0 * (2.0 * rSqd**-7 - rSqd**-4)
+ acc[i, 0] += rij0 * f
+ acc[i, 1] += rij1 * f
+ acc[i, 2] += rij2 * f
+ acc[j, 0] -= rij0 * f
+ acc[j, 1] -= rij1 * f
+ acc[j, 2] -= rij2 * f
+ return acc
+
+
+@njit(cache=True, fastmath=True)
+def compute_potential(pos: np.ndarray, N: int) -> float:
+ """Total LJ potential energy, summed once per distinct pair i<j."""
+ pot = 0.0
+ for i in range(N - 1):
+ for j in range(i + 1, N):
+ rij0 = pos[i, 0] - pos[j, 0]
+ rij1 = pos[i, 1] - pos[j, 1]
+ rij2 = pos[i, 2] - pos[j, 2]
+ rSqd = rij0 * rij0 + rij1 * rij1 + rij2 * rij2
+ rnorm = np.sqrt(rSqd)
+ quot = 1.0 / rnorm
+ term1 = quot**12
+ term2 = quot**6
+ pot += 4.0 * (term1 - term2)
+ return pot
+
+
+@njit(cache=True, fastmath=True)
+def kinetic_energy(vel: np.ndarray, N: int) -> float:
+ kin = 0.0
+ for i in range(N):
+ v2 = vel[i, 0] ** 2 + vel[i, 1] ** 2 + vel[i, 2] ** 2
+ kin += 0.5 * v2
+ return kin
+
+
+@njit(cache=True, fastmath=True)
+def mean_squared_velocity(vel: np.ndarray, N: int) -> float:
+ vx2 = 0.0
+ vy2 = 0.0
+ vz2 = 0.0
+ for i in range(N):
+ vx2 += vel[i, 0] * vel[i, 0]
+ vy2 += vel[i, 1] * vel[i, 1]
+ vz2 += vel[i, 2] * vel[i, 2]
+ return (vx2 + vy2 + vz2) / N
+
+
+@njit(cache=True, fastmath=True)
+def _apply_walls(pos: np.ndarray, vel: np.ndarray, N: int, L: float, dt: float) -> float:
+ """Elastic walls: reverse velocity (not position) on out-of-box coords.
+
+ Returns psum, the summed momentum-transfer contribution to pressure.
+ """
+ psum = 0.0
+ for i in range(N):
+ for k in range(3):
+ if pos[i, k] < 0.0:
+ vel[i, k] *= -1.0
+ psum += 2.0 * abs(vel[i, k]) / dt
+ if pos[i, k] >= L:
+ vel[i, k] *= -1.0
+ psum += 2.0 * abs(vel[i, k]) / dt
+ return psum
+
+
+@njit(cache=True, fastmath=True)
+def vv_step_naive(pos: np.ndarray, vel: np.ndarray, N: int, L: float, dt: float) -> float:
+ """One velocity-Verlet step, recomputing forces twice (matches MD.cpp's
+ ``VelocityVerlet`` exactly, including its redundant leading force call).
+ Mutates pos/vel in place. Returns instantaneous pressure (natural units).
+ """
+ acc = compute_forces(pos, N) # redundant recompute, as in the reference
+ for i in range(N):
+ for k in range(3):
+ pos[i, k] += vel[i, k] * dt + 0.5 * acc[i, k] * dt * dt
+ vel[i, k] += 0.5 * acc[i, k] * dt
+ acc = compute_forces(pos, N)
+ for i in range(N):
+ for k in range(3):
+ vel[i, k] += 0.5 * acc[i, k] * dt
+ psum = _apply_walls(pos, vel, N, L, dt)
+ return psum / (6.0 * L * L)
+
+
+@njit(cache=True, fastmath=True)
+def vv_step_cached(pos: np.ndarray, vel: np.ndarray, acc: np.ndarray, N: int,
+ L: float, dt: float) -> float:
+ """One velocity-Verlet step given the acceleration already evaluated at
+ the current position (from the previous step's second force call, or the
+ pre-loop initial force call). Mutates pos/vel/acc in place; acc is left
+ holding the force at the *new* position for reuse by the next step.
+ Returns instantaneous pressure (natural units).
+ """
+ for i in range(N):
+ for k in range(3):
+ pos[i, k] += vel[i, k] * dt + 0.5 * acc[i, k] * dt * dt
+ vel[i, k] += 0.5 * acc[i, k] * dt
+ new_acc = compute_forces(pos, N)
+ for i in range(N):
+ for k in range(3):
+ vel[i, k] += 0.5 * new_acc[i, k] * dt
+ acc[i, k] = new_acc[i, k]
+ psum = _apply_walls(pos, vel, N, L, dt)
+ return psum / (6.0 * L * L)
+
+
+@njit(cache=True, fastmath=True)
+def simulate(pos0: np.ndarray, vel0: np.ndarray, L: float, dt: float,
+ n_steps: int):
+ """Run n_steps+1 velocity-Verlet iterations (matching MD.cpp's ``i<NumTime+1``
+ loop bound), recording per-iteration instantaneous T, P, KE, PE (all
+ natural units). Returns (T_nu, P_nu, KE, PE, pos_final, vel_final).
+ """
+ N = pos0.shape[0]
+ pos = pos0.copy()
+ vel = vel0.copy()
+ acc = compute_forces(pos, N)
+
+ n_records = n_steps + 1
+ T_nu = np.empty(n_records)
+ P_nu = np.empty(n_records)
+ KE = np.empty(n_records)
+ PE = np.empty(n_records)
+
+ for i in range(n_records):
+ p_inst = vv_step_cached(pos, vel, acc, N, L, dt)
+ mvs = mean_squared_velocity(vel, N)
+ T_nu[i] = mvs / 3.0
+ P_nu[i] = p_inst
+ KE[i] = kinetic_energy(vel, N)
+ PE[i] = compute_potential(pos, N)
+
+ return T_nu, P_nu, KE, PE, pos, vel
+
+
+@njit(cache=True, fastmath=True)
+def simulate_with_trajectory(pos0: np.ndarray, vel0: np.ndarray, L: float, dt: float,
+ n_steps: int):
+ """Same as ``simulate`` but additionally records positions at every
+ iteration (memory-heavy; only used when the caller requests a trajectory).
+ """
+ N = pos0.shape[0]
+ pos = pos0.copy()
+ vel = vel0.copy()
+ acc = compute_forces(pos, N)
+
+ n_records = n_steps + 1
+ T_nu = np.empty(n_records)
+ P_nu = np.empty(n_records)
+ KE = np.empty(n_records)
+ PE = np.empty(n_records)
+ traj = np.empty((n_records, N, 3))
+
+ for i in range(n_records):
+ p_inst = vv_step_cached(pos, vel, acc, N, L, dt)
+ mvs = mean_squared_velocity(vel, N)
+ T_nu[i] = mvs / 3.0
+ P_nu[i] = p_inst
+ KE[i] = kinetic_energy(vel, N)
+ PE[i] = compute_potential(pos, N)
+ traj[i] = pos
+
+ return T_nu, P_nu, KE, PE, pos, vel, traj
+
+
+def initialize_lattice(N: int, L: float) -> np.ndarray:
+ """Simple-cubic lattice positions, matching MD.cpp's ``initialize()``."""
+ n = int(np.ceil(N ** (1.0 / 3.0)))
+ pos = np.empty((N, 3))
+ p = 0
+ spacing = L / n
+ for i in range(n):
+ for j in range(n):
+ for k in range(n):
+ if p < N:
+ pos[p, 0] = (i + 0.5) * spacing
+ pos[p, 1] = (j + 0.5) * spacing
+ pos[p, 2] = (k + 0.5) * spacing
+ p += 1
+ return pos
+
+
+def initialize_velocities(N: int, T_init_nu: float, rng: np.random.Generator) -> np.ndarray:
+ """Gaussian velocities, COM removed, rescaled to T_init_nu (natural units),
+ matching MD.cpp's ``initializeVelocities()`` (including its (N-1) factor).
+ """
+ vel = rng.normal(0.0, 1.0, size=(N, 3))
+ vcm = vel.mean(axis=0)
+ vel -= vcm
+ vsqd_sum = np.sum(vel ** 2)
+ lam = np.sqrt(3.0 * (N - 1) * T_init_nu / vsqd_sum)
+ vel *= lam
+ return vel
+
+
+def warmup() -> None:
+ """Trigger JIT compilation on a throwaway tiny system so it doesn't land
+ on the user's first timed simulation."""
+ pos = initialize_lattice(4, 4.0)
+ vel = np.zeros((4, 3))
+ simulate(pos, vel, 4.0, 0.001, 1)
+ # also warm the naive path used by the equivalence test
+ p2 = pos.copy()
+ v2 = np.zeros((4, 3))
+ vv_step_naive(p2, v2, 4, 4.0, 0.001)
diff --git a/src/noblegasmd/units.py b/src/noblegasmd/units.py
new file mode 100644
index 0000000..1819963
--- /dev/null
+++ b/src/noblegasmd/units.py
@@ -0,0 +1,61 @@
+"""Physical constants and per-gas unit-conversion factors.
+
+Transcribed verbatim from ``MD.cpp`` (Foley, Sweet, Akinfenwa, GPLv3) — the C
+source is the authoritative reference; do not "clean up" these values, they
+must match to full precision (see tests/test_units.py).
+"""
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+# Avogadro's number, mol^-1
+NA = 6.022140857e23
+# Boltzmann constant, m^2 kg / (s^2 K)
+KB_SI = 1.38064852e-23
+
+
+@dataclass(frozen=True)
+class GasConstants:
+ """Natural-unit <-> SI conversion factors for one noble gas."""
+
+ name: str
+ vol_fac: float # m^3 per natural unit of volume
+ press_fac: float # Pa per natural unit of pressure
+ temp_fac: float # K per natural unit of temperature
+ timefac: float # s per natural unit of time
+
+
+GAS_CONSTANTS: dict[str, GasConstants] = {
+ "He": GasConstants("He", 1.8399744000000005e-29, 8152287.336171632,
+ 10.864459551225972, 1.7572698825166272e-12),
+ "Ne": GasConstants("Ne", 2.0570823999999997e-29, 27223022.27659913,
+ 40.560648991243625, 2.1192341945685407e-12),
+ "Ar": GasConstants("Ar", 3.7949992920124995e-29, 51695201.06691862,
+ 142.0950000000000, 2.09618e-12),
+ "Kr": GasConstants("Kr", 4.5882712000000004e-29, 59935428.40275003,
+ 199.1817584391428, 8.051563913585078e-13),
+ "Xe": GasConstants("Xe", 5.4872e-29, 70527773.72794868,
+ 280.30305642163006, 9.018957925790732e-13),
+}
+
+
+def get_gas_constants(gas: str) -> GasConstants:
+ """Look up per-gas unit factors, matching MD.cpp's fallback-to-Ar behavior
+ only for the empty/unset case; unlike the C, an unrecognized gas string is
+ an error here rather than a silent fallback."""
+ try:
+ return GAS_CONSTANTS[gas]
+ except KeyError as exc:
+ valid = ", ".join(GAS_CONSTANTS)
+ raise ValueError(f"Unknown gas {gas!r}; must be one of: {valid}") from exc
+
+
+def steps_for_gas(gas: str) -> tuple[float, int]:
+ """Return (dt_fraction_seconds, NumTime) matching MD.cpp's timestep rule.
+
+ dt is returned as the SI seconds numerator (0.2e-14 or 0.5e-14); divide by
+ timefac to get dt in natural units, exactly as MD.cpp does.
+ """
+ if gas == "He":
+ return 0.2e-14, 50000
+ return 0.5e-14, 20000
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..82a3abc
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,4 @@
+import sys
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parent / "oracle"))
diff --git a/tests/oracle/__init__.py b/tests/oracle/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/oracle/build.sh b/tests/oracle/build.sh
new file mode 100755
index 0000000..1cb1276
--- /dev/null
+++ b/tests/oracle/build.sh
@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+# Compiles the reference MD.cpp (golden oracle) exactly as the repo's own Makefile does.
+set -euo pipefail
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+REPO_ROOT="$(cd "$SCRIPT_DIR/../.." && pwd)"
+OUT="${1:-$SCRIPT_DIR/MD_oracle.exe}"
+g++ -O3 -I"$REPO_ROOT" -L"$REPO_ROOT" "$REPO_ROOT/MD.cpp" -o "$OUT"
+echo "Built oracle binary: $OUT"
diff --git a/tests/oracle/oracle.py b/tests/oracle/oracle.py
new file mode 100644
index 0000000..28b4076
--- /dev/null
+++ b/tests/oracle/oracle.py
@@ -0,0 +1,70 @@
+"""Golden-oracle harness: drives the compiled MD.cpp reference binary via stdin
+and parses its <title>_average.txt output. Used only by tests, never by the package.
+"""
+from __future__ import annotations
+
+import subprocess
+import tempfile
+from dataclasses import dataclass
+from pathlib import Path
+
+REPO_ROOT = Path(__file__).resolve().parents[2]
+ORACLE_SRC = REPO_ROOT / "MD.cpp"
+ORACLE_BIN = Path(__file__).resolve().parent / "MD_oracle.exe"
+
+
+@dataclass(frozen=True)
+class OracleResult:
+ T: float # average temperature, K
+ P: float # average pressure, Pa
+ gc: float # PV/nT, J/(mol K)
+ Z: float # compressibility factor, unitless
+ V: float # volume, m^3
+ N: int
+
+
+def build_oracle(force: bool = False) -> Path:
+ """Compile MD.cpp into ORACLE_BIN if not already built (or force rebuild)."""
+ if ORACLE_BIN.exists() and not force:
+ return ORACLE_BIN
+ subprocess.run(
+ ["g++", "-O3", "-I", str(REPO_ROOT), "-L", str(REPO_ROOT),
+ str(ORACLE_SRC), "-o", str(ORACLE_BIN)],
+ check=True,
+ )
+ return ORACLE_BIN
+
+
+def run_oracle(gas: str, T: float, rho: float, title: str = "oracle_run") -> OracleResult:
+ """Run the C reference MD.cpp for one state point and return its averages.
+
+ MD.cpp prompts on stdin for: title, gas, T (K), rho (mol/m^3) and writes
+ averages to '<title>_average.txt' in the current working directory. The
+ binary's RNG is unseeded `rand()`, which glibc/libstdc++ seed identically
+ to seed 1 on every run, so results are deterministic.
+ """
+ build_oracle()
+ stdin_text = f"{title}\n{gas}\n{T}\n{rho}\n"
+
+ with tempfile.TemporaryDirectory() as tmpdir:
+ subprocess.run(
+ [str(ORACLE_BIN)],
+ input=stdin_text,
+ text=True,
+ cwd=tmpdir,
+ check=True,
+ capture_output=True,
+ timeout=600,
+ )
+ avg_path = Path(tmpdir) / f"{title}_average.txt"
+ return _parse_average_file(avg_path)
+
+
+def _parse_average_file(path: Path) -> OracleResult:
+ lines = path.read_text().splitlines()
+ # Line 0: header, line 1: dashes, line 2: data row.
+ data_line = lines[2]
+ tokens = data_line.split()
+ # time, T, P, gc, Z, V, N
+ _time, T, P, gc, Z, V, N = tokens
+ return OracleResult(T=float(T), P=float(P), gc=float(gc), Z=float(Z), V=float(V), N=int(N))
diff --git a/tests/test_energy_conservation.py b/tests/test_energy_conservation.py
new file mode 100644
index 0000000..cc801c8
--- /dev/null
+++ b/tests/test_energy_conservation.py
@@ -0,0 +1,26 @@
+"""Oracle-free: elastic walls + conservative forces => total energy is conserved.
+
+Do NOT test momentum conservation here: walls inject impulse by design, so
+linear momentum is not conserved and is not a physics bug.
+"""
+from noblegasmd import run
+
+ENERGY_DRIFT_TOLERANCE = 1e-3
+
+
+def test_energy_conserved_over_full_run():
+ result = run(gas="Ar", T=300.0, rho=40.0, seed=1)
+ assert result.energy_drift < ENERGY_DRIFT_TOLERANCE, (
+ f"Relative energy drift {result.energy_drift} exceeds "
+ f"{ENERGY_DRIFT_TOLERANCE}"
+ )
+
+
+def test_energy_conserved_at_higher_density():
+ # Denser system => more frequent close encounters, a harder case for
+ # energy conservation.
+ result = run(gas="Ar", T=300.0, rho=5000.0, seed=2)
+ assert result.energy_drift < ENERGY_DRIFT_TOLERANCE, (
+ f"Relative energy drift {result.energy_drift} exceeds "
+ f"{ENERGY_DRIFT_TOLERANCE}"
+ )
diff --git a/tests/test_ideal_gas_limit.py b/tests/test_ideal_gas_limit.py
new file mode 100644
index 0000000..2ebac75
--- /dev/null
+++ b/tests/test_ideal_gas_limit.py
@@ -0,0 +1,31 @@
+"""Oracle-free: at low density, the compressibility factor Z -> 1 (ideal gas).
+
+The instantaneous pressure estimator (momentum transfer at wall collisions)
+is intrinsically noisy for a single seed -- collisions are rare events, so a
+single trajectory's Z can swing several percent even at low density. We
+therefore average over a handful of replicate seeds, exactly as the published
+methodology and the oracle-agreement test do. rho=40 mol/m^3 is MD.cpp's own
+reference value for "ideal gas at STP" (see the density prompt in MD.cpp).
+"""
+from noblegasmd import sweep
+
+IDEAL_GAS_Z_TOLERANCE = 0.02
+IDEAL_GAS_RHO = 40.0
+N_REPLICATES = 5
+
+
+def test_z_approaches_one_at_low_density():
+ df = sweep(gas="Ar", T=[300.0], rho=[IDEAL_GAS_RHO], n_replicates=N_REPLICATES, seed=42)
+ z_mean = df["Z"].mean()
+ assert abs(z_mean - 1.0) < IDEAL_GAS_Z_TOLERANCE, (
+ f"Expected ensemble-mean Z near 1 at low density, got Z={z_mean}"
+ )
+
+
+def test_z_approaches_one_across_low_density_temperatures():
+ for T in (100.0, 200.0, 300.0, 400.0):
+ df = sweep(gas="Ar", T=[T], rho=[IDEAL_GAS_RHO], n_replicates=N_REPLICATES, seed=int(T))
+ z_mean = df["Z"].mean()
+ assert abs(z_mean - 1.0) < IDEAL_GAS_Z_TOLERANCE, (
+ f"T={T}: expected ensemble-mean Z near 1 at low density, got Z={z_mean}"
+ )
diff --git a/tests/test_units.py b/tests/test_units.py
new file mode 100644
index 0000000..10af0c1
--- /dev/null
+++ b/tests/test_units.py
@@ -0,0 +1,40 @@
+"""Constants round-trip: per-gas factors must match MD.cpp to full precision."""
+from noblegasmd.units import GAS_CONSTANTS, KB_SI, NA, get_gas_constants
+
+# Transcribed directly from the literal values in MD.cpp's main().
+EXPECTED = {
+ "He": dict(vol_fac=1.8399744000000005e-29, press_fac=8152287.336171632,
+ temp_fac=10.864459551225972, timefac=1.7572698825166272e-12),
+ "Ne": dict(vol_fac=2.0570823999999997e-29, press_fac=27223022.27659913,
+ temp_fac=40.560648991243625, timefac=2.1192341945685407e-12),
+ "Ar": dict(vol_fac=3.7949992920124995e-29, press_fac=51695201.06691862,
+ temp_fac=142.0950000000000, timefac=2.09618e-12),
+ "Kr": dict(vol_fac=4.5882712000000004e-29, press_fac=59935428.40275003,
+ temp_fac=199.1817584391428, timefac=8.051563913585078e-13),
+ "Xe": dict(vol_fac=5.4872e-29, press_fac=70527773.72794868,
+ temp_fac=280.30305642163006, timefac=9.018957925790732e-13),
+}
+
+
+def test_global_constants_match_md_cpp():
+ assert NA == 6.022140857e23
+ assert KB_SI == 1.38064852e-23
+
+
+def test_all_five_gases_present():
+ assert set(GAS_CONSTANTS) == set(EXPECTED)
+
+
+def test_per_gas_constants_match_md_cpp_exactly():
+ for gas, expected in EXPECTED.items():
+ gc = get_gas_constants(gas)
+ assert gc.vol_fac == expected["vol_fac"], gas
+ assert gc.press_fac == expected["press_fac"], gas
+ assert gc.temp_fac == expected["temp_fac"], gas
+ assert gc.timefac == expected["timefac"], gas
+
+
+def test_unknown_gas_raises():
+ import pytest
+ with pytest.raises(ValueError):
+ get_gas_constants("Rn")
diff --git a/tests/test_vs_c_oracle.py b/tests/test_vs_c_oracle.py
new file mode 100644
index 0000000..96077c4
--- /dev/null
+++ b/tests/test_vs_c_oracle.py
@@ -0,0 +1,39 @@
+"""Python ensemble-mean Z must agree with the deterministic C oracle (MD.cpp)
+within combined statistical tolerance, across an Ar grid spanning ideal to
+strongly non-ideal densities. MD is chaotic and the Python port uses a
+different RNG, so we validate ensemble-mean thermodynamic averages only,
+never trajectory identity.
+"""
+import numpy as np
+import pytest
+from oracle import run_oracle
+
+from noblegasmd import sweep
+
+Z_AGREEMENT_TOLERANCE = 0.03 # 3% relative, per domain-expert sign-off
+
+TEMPERATURES = [100.0, 200.0, 300.0, 400.0]
+DENSITIES = [500.0, 2000.0] # dense enough that the oracle's single deterministic sample is low-noise;
+# rho=40 (sparse wall collisions) is covered instead by the oracle-free, ensemble-averaged
+# check in test_ideal_gas_limit.py, since at that density even the C oracle's own single
+# run swings several percent from its own long-run mean.
+N_REPLICATES = 3
+
+
+@pytest.mark.parametrize("T_set", TEMPERATURES)
+@pytest.mark.parametrize("rho_set", DENSITIES)
+def test_ensemble_mean_z_agrees_with_c_oracle(T_set, rho_set):
+ df = sweep(
+ gas="Ar", T=[T_set], rho=[rho_set],
+ n_replicates=N_REPLICATES, seed=int(T_set * 100 + rho_set),
+ )
+ python_z_mean = df["Z"].mean()
+
+ oracle_result = run_oracle("Ar", T_set, rho_set, title=f"oracle_{T_set:.0f}_{rho_set:.0f}")
+
+ rel_err = abs(python_z_mean - oracle_result.Z) / abs(oracle_result.Z)
+ assert rel_err < Z_AGREEMENT_TOLERANCE, (
+ f"T={T_set} rho={rho_set}: python Z_mean={python_z_mean:.5f} "
+ f"(std={df['Z'].std():.5f}) vs oracle Z={oracle_result.Z:.5f}, "
+ f"rel_err={rel_err:.4f}"
+ )
diff --git a/tests/test_vv_equivalence.py b/tests/test_vv_equivalence.py
new file mode 100644
index 0000000..14c914f
--- /dev/null
+++ b/tests/test_vv_equivalence.py
@@ -0,0 +1,28 @@
+"""The cached single-force-per-step velocity Verlet must be numerically
+identical to the naive reference implementation that recomputes forces twice
+per step (matching MD.cpp's VelocityVerlet exactly)."""
+import numpy as np
+
+from noblegasmd import engine
+
+
+def test_cached_matches_naive_trajectory():
+ rng = np.random.default_rng(0)
+ N, L, dt, n_steps = 20, 6.0, 1e-3, 50
+
+ pos0 = engine.initialize_lattice(N, L)
+ vel0 = engine.initialize_velocities(N, 1.0, rng)
+
+ pos_naive = pos0.copy()
+ vel_naive = vel0.copy()
+ for _ in range(n_steps):
+ engine.vv_step_naive(pos_naive, vel_naive, N, L, dt)
+
+ pos_cached = pos0.copy()
+ vel_cached = vel0.copy()
+ acc_cached = engine.compute_forces(pos_cached, N)
+ for _ in range(n_steps):
+ engine.vv_step_cached(pos_cached, vel_cached, acc_cached, N, L, dt)
+
+ np.testing.assert_allclose(pos_cached, pos_naive, rtol=1e-9, atol=1e-12)
+ np.testing.assert_allclose(vel_cached, vel_naive, rtol=1e-9, atol=1e-12)