From b8982f7dff77e6bb03698740f084c53efc68e034 Mon Sep 17 00:00:00 2001 From: FoleyLab Date: Sun, 16 Aug 2026 15:52:09 -0400 Subject: [PATCH] python implementation --- .github/workflows/ci.yml | 31 ++ .github/workflows/release.yml | 72 ++++ .gitignore | 12 + LICENSE | 674 ++++++++++++++++++++++++++++++ MD_python_port_spec.md | 240 +++++++++++ PUBLISHING.md | 48 +++ README.md | 27 ++ benchmark.py | 37 ++ notebooks/reproduce_figure4.ipynb | 159 +++++++ pyproject.toml | 42 ++ src/noblegasmd/__init__.py | 21 + src/noblegasmd/api.py | 185 ++++++++ src/noblegasmd/engine.py | 237 +++++++++++ src/noblegasmd/units.py | 61 +++ tests/__init__.py | 0 tests/conftest.py | 4 + tests/oracle/__init__.py | 0 tests/oracle/build.sh | 8 + tests/oracle/oracle.py | 70 ++++ tests/test_energy_conservation.py | 26 ++ tests/test_ideal_gas_limit.py | 31 ++ tests/test_units.py | 40 ++ tests/test_vs_c_oracle.py | 39 ++ tests/test_vv_equivalence.py | 28 ++ 24 files changed, 2092 insertions(+) create mode 100644 .github/workflows/ci.yml create mode 100644 .github/workflows/release.yml create mode 100644 LICENSE create mode 100644 MD_python_port_spec.md create mode 100644 PUBLISHING.md create mode 100644 benchmark.py create mode 100644 notebooks/reproduce_figure4.ipynb create mode 100644 pyproject.toml create mode 100644 src/noblegasmd/__init__.py create mode 100644 src/noblegasmd/api.py create mode 100644 src/noblegasmd/engine.py create mode 100644 src/noblegasmd/units.py create mode 100644 tests/__init__.py create mode 100644 tests/conftest.py create mode 100644 tests/oracle/__init__.py create mode 100755 tests/oracle/build.sh create mode 100644 tests/oracle/oracle.py create mode 100644 tests/test_energy_conservation.py create mode 100644 tests/test_ideal_gas_limit.py create mode 100644 tests/test_units.py create mode 100644 tests/test_vs_c_oracle.py create mode 100644 tests/test_vv_equivalence.py 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)